diff --git a/Core/TomoP3DModel_core.c b/Core/TomoP3DModel_core.c index fd93ce3..9e93433 100644 --- a/Core/TomoP3DModel_core.c +++ b/Core/TomoP3DModel_core.c @@ -99,10 +99,10 @@ float TomoP3DObject_core(float *A, long N1, long N2, long N3, long Z1, long Z2, psi2 = psi_gr2*((float)M_PI / 180.0f); psi3 = psi_gr3*((float)M_PI / 180.0f); - float bs[3][3] = { - { 0.0f,0.0f,0.0f }, - { 0.0f,0.0f,0.0f }, - { 0.0f,0.0f,0.0f } }; + float bs[3][3] = { + {0.0f,0.0f,0.0f}, + {0.0f,0.0f,0.0f}, + {0.0f,0.0f,0.0f} }; float xh[3] = { 0.0f, 0.0f, 0.0f }; float xh3[3] = { 0.0f, 0.0f, 0.0f }; diff --git a/Core/TomoP3DModel_core.h b/Core/TomoP3DModel_core.h index 5e1ae51..8883198 100644 --- a/Core/TomoP3DModel_core.h +++ b/Core/TomoP3DModel_core.h @@ -21,8 +21,6 @@ limitations under the License. #include "CCPiDefines.h" #include "utils.h" - -//CCPI_EXPORT float TomoP3DModel_core(float *A, int ModelSelected, int N, char *ModelParametersFilename); CCPI_EXPORT float TomoP3DModel_core(float *A, int ModelSelected, long N1, long N2, long N3, long Z1, long Z2, char *ModelParametersFilename); CCPI_EXPORT float TomoP3DObject_core(float *A, long N1, long N2, long N3, long Z1, long Z2, char *Object, float C0, /* intensity */ diff --git a/Core/utils.c b/Core/utils.c index 747d8cb..fcb8124 100644 --- a/Core/utils.c +++ b/Core/utils.c @@ -502,82 +502,44 @@ float checkParams3D(int *params_switch, int ModelSelected, char *ModelParameters } /***********************************************************************************************/ -/* rotation matrix routine */ -float matrot3(float *Ad, float psi1, float psi2, float psi3) +/* rotation matrix */ +float matrot3(float Ad[3][3], float psi1, float psi2, float psi3) { - Ad[0] = cosf(psi1)*cosf(psi2)*cosf(psi3)-sinf(psi1)*sinf(psi3); - Ad[1] = sinf(psi1)*cosf(psi2)*cosf(psi3)+cosf(psi1)*sinf(psi3); - Ad[2] = -sinf(psi2)*cosf(psi3); - Ad[3] = -cosf(psi1)*cosf(psi2)*sinf(psi3)-sinf(psi1)*cosf(psi3); - Ad[4] = -sinf(psi1)*cosf(psi2)*sinf(psi3)+cosf(psi1)*cosf(psi3); - Ad[5] = sinf(psi2)*sinf(psi3); - Ad[6] = cosf(psi1)*sinf(psi2); - Ad[7] = sinf(psi1)*sinf(psi2); - Ad[8] = cosf(psi2); - return *Ad; + Ad[0][0]=cosf(psi1)*cosf(psi2)*cosf(psi3)-sinf(psi1)*sinf(psi3); + Ad[0][1]=sinf(psi1)*cosf(psi2)*cosf(psi3)+cosf(psi1)*sinf(psi3); + Ad[0][2]=-sinf(psi2)*cosf(psi3); + Ad[1][0]=-cosf(psi1)*cosf(psi2)*sinf(psi3)-sinf(psi1)*cosf(psi3); + Ad[1][1]=-sinf(psi1)*cosf(psi2)*sinf(psi3)+cosf(psi1)*cosf(psi3); + Ad[1][2]=sinf(psi2)*sinf(psi3); + Ad[2][0]=cosf(psi1)*sinf(psi2); + Ad[2][1]=sinf(psi1)*sinf(psi2); + Ad[2][2]=cosf(psi2); + return 1; } /*matrix-vector multiplication*/ -float matvet3(float *Ad, float *V1, float *V2) +float matvet3(float Ad[3][3], float V1[3], float V2[3]) { - int l, m, counterT; - - counterT = 0; + int l, m; for(l=0; l<3; l++) { V2[l] = 0.0f; for(m=0; m<3; m++) { - V2[l] += Ad[counterT]*V1[m]; - counterT++; + V2[l] += Ad[l][m]*V1[m]; }} - return *V2; + return 1; } /*matrix-matrix multiplication*/ -// float matmat3(float *A, float *B, float *C) -// { -// float Am[3][3]; -// float Bm[3][3]; -// float Cm[3][3]; -// -// Am[0][0]=A[0]; -// Am[0][1]=A[1]; -// Am[0][2]=A[2]; -// Am[1][0]=A[3]; -// Am[1][1]=A[4]; -// Am[1][2]=A[5]; -// Am[2][0]=A[6]; -// Am[2][1]=A[7]; -// Am[2][2]=A[8]; -// -// Bm[0][0]=B[0]; -// Bm[0][1]=B[1]; -// Bm[0][2]=B[2]; -// Bm[1][0]=B[3]; -// Bm[1][1]=B[4]; -// Bm[1][2]=B[5]; -// Bm[2][0]=B[6]; -// Bm[2][1]=B[7]; -// Bm[2][2]=B[8]; -// -// int i, j, k; -// for(i=0; i<3; i++) { -// for(j=0; j<3; j++) { -// Cm[i][j]=0.0f; -// for(k=0; k<3; k++) { -// Cm[i][j] += Am[i][k]*Bm[k][j]; -// } -// }} -// -// C[0] = Cm[0][0]; -// C[1] = Cm[0][1]; -// C[2] = Cm[0][2]; -// C[3] = Cm[1][0]; -// C[4] = Cm[1][1]; -// C[5] = Cm[1][2]; -// C[6] = Cm[2][0]; -// C[7] = Cm[2][1]; -// C[8] = Cm[2][2]; -// -// -// return *C; -// } +/*float matmat3(float Am[3][3], float Bm[3][3], float Cm[3][3]) +{ + int i, j, k; + for(i=0; i<3; i++) { + for(j=0; j<3; j++) { + Cm[i][j]=0.0f; + for(k=0; k<3; k++) { + Cm[i][j] += Am[i][k]*Bm[k][j]; + } + }} + return 1; +} +*/ diff --git a/Core/utils.h b/Core/utils.h index d5baa62..771e1a0 100644 --- a/Core/utils.h +++ b/Core/utils.h @@ -23,9 +23,9 @@ extern "C" { #endif CCPI_EXPORT float checkParams2D(int *params_switch, int ModelSelected, char *ModelParametersFilename); CCPI_EXPORT float checkParams3D(int *params_switch, int ModelSelected, char *ModelParametersFilename); -CCPI_EXPORT float matrot3(float *Ad, float psi1, float psi2, float psi3); -CCPI_EXPORT float matvet3(float *Ad, float *V1, float *V2); -// float matmat3(float *A, float *B, float *C); +CCPI_EXPORT float matrot3(float Ad[3][3], float psi1, float psi2, float psi3); +CCPI_EXPORT float matvet3(float Ad[3][3], float V1[3], float V2[3]); +// CCPI_EXPORT float matmat3(float Am[3][3], float Bm[3][3], float Cm[3][3]); #ifdef __cplusplus } #endif diff --git a/Wrappers/MATLAB/Model2DGeneratorDemo.m b/Wrappers/MATLAB/Model2DGeneratorDemo.m index 43c83aa..c1a85e7 100644 --- a/Wrappers/MATLAB/Model2DGeneratorDemo.m +++ b/Wrappers/MATLAB/Model2DGeneratorDemo.m @@ -10,8 +10,6 @@ close all;clc;clear; % adding paths fsep = '/'; -pathtoModels = sprintf(['..' fsep 'functions' fsep 'models' fsep], 1i); -addpath(pathtoModels); addpath('compiled'); addpath('supplem'); ModelNo = 4; % Select a model from Phantom2DLibrary.dat @@ -21,7 +19,7 @@ % Generate 2D phantom: curDir = pwd; mainDir = fileparts(curDir); -pathtoLibrary = sprintf([fsep 'functions' fsep 'models' fsep 'Phantom2DLibrary.dat'], 1i); +pathtoLibrary = sprintf([fsep '..' fsep 'PhantomLibrary' fsep 'models' fsep 'Phantom2DLibrary.dat'], 1i); pathTP = strcat(mainDir, pathtoLibrary); % path to TomoPhantom parameters file [G] = TomoP2DModel(ModelNo,N,pathTP); figure; imagesc(G, [0 1]); daspect([1 1 1]); colormap hot; diff --git a/Wrappers/MATLAB/Model2D_tGeneratorDemo.m b/Wrappers/MATLAB/Model2D_tGeneratorDemo.m index 66cd559..36e70b3 100644 --- a/Wrappers/MATLAB/Model2D_tGeneratorDemo.m +++ b/Wrappers/MATLAB/Model2D_tGeneratorDemo.m @@ -10,8 +10,6 @@ close all;clc;clear; % adding paths fsep = '/'; -pathtoModels = sprintf(['..' fsep 'functions' fsep 'models' fsep], 1i); -addpath(pathtoModels); addpath('compiled'); addpath('supplem'); ModelNo = 100; % Select a model from Phantom2DLibrary.dat @@ -21,10 +19,10 @@ % Generate 2D+t phantom: curDir = pwd; mainDir = fileparts(curDir); -pathtoLibrary = sprintf([fsep 'functions' fsep 'models' fsep 'Phantom2DLibrary.dat'], 1i); +pathtoLibrary = sprintf([fsep '..' fsep 'PhantomLibrary' fsep 'models' fsep 'Phantom2DLibrary.dat'], 1i); pathTP = strcat(mainDir, pathtoLibrary); % path to TomoPhantom parameters file [G] = TomoP2DModel(ModelNo,N,pathTP); -figure(1); imagesc(G, [0 1]); daspect([1 1 1]); title('2D+t model, t=3 here'); colormap hot; +figure(1); imagesc(G(:,:,3), [0 1]); daspect([1 1 1]); title('2D+t model, t=3 here'); colormap hot; %% % Lets look at more finely discretized temporal model and related sinograms ModelNo = 101; % Select a model from Phantom2DLibrary.dat @@ -34,7 +32,7 @@ % Generate 2D phantom: curDir = pwd; mainDir = fileparts(curDir); -pathtoLibrary = sprintf([fsep 'functions' fsep 'models' fsep 'Phantom2DLibrary.dat'], 1i); +pathtoLibrary = sprintf([fsep '..' fsep 'PhantomLibrary' fsep 'models' fsep 'Phantom2DLibrary.dat'], 1i); pathTP = strcat(mainDir, pathtoLibrary); % path to TomoPhantom parameters file [G] = TomoP2DModel(ModelNo,N,pathTP); @@ -58,7 +56,7 @@ timeFrames = 25; %must be the same as in model curDir = pwd; mainDir = fileparts(curDir); -pathtoLibrary = sprintf([fsep 'functions' fsep 'models' fsep 'Phantom2DLibrary.dat'], 1i); +pathtoLibrary = sprintf([fsep '..' fsep 'PhantomLibrary' fsep 'models' fsep 'Phantom2DLibrary.dat'], 1i); pathTP = strcat(mainDir, pathtoLibrary); % path to TomoPhantom parameters file [G] = TomoP2DModel(ModelNo,N,pathTP); diff --git a/Wrappers/MATLAB/Model3DGeneratorDemo.m b/Wrappers/MATLAB/Model3DGeneratorDemo.m index 5726402..96730c6 100644 --- a/Wrappers/MATLAB/Model3DGeneratorDemo.m +++ b/Wrappers/MATLAB/Model3DGeneratorDemo.m @@ -8,8 +8,6 @@ close all;clc;clear; % adding paths fsep = '/'; -pathtoModels = sprintf(['..' fsep 'functions' fsep 'models' fsep], 1i); -addpath(pathtoModels); addpath('compiled'); addpath('supplem'); ModelNo = 10; % Select a model @@ -19,7 +17,7 @@ % generate 3D phantom (modify your PATH bellow): curDir = pwd; mainDir = fileparts(curDir); -pathtoLibrary = sprintf([fsep 'functions' fsep 'models' fsep 'Phantom3DLibrary.dat'], 1i); +pathtoLibrary = sprintf([fsep '..' fsep 'PhantomLibrary' fsep 'models' fsep 'Phantom3DLibrary.dat'], 1i); pathTP = strcat(mainDir, pathtoLibrary); % path to TomoPhantom parameters file tic; [G] = TomoP3DModel(ModelNo,N,pathTP); toc; diff --git a/Wrappers/MATLAB/Model3D_tGeneratorDemo.m b/Wrappers/MATLAB/Model3D_tGeneratorDemo.m index 782e4e8..b46a1c1 100644 --- a/Wrappers/MATLAB/Model3D_tGeneratorDemo.m +++ b/Wrappers/MATLAB/Model3D_tGeneratorDemo.m @@ -8,8 +8,6 @@ close all;clc;clear; % adding paths fsep = '/'; -pathtoModels = sprintf(['..' fsep 'functions' fsep 'models' fsep], 1i); -addpath(pathtoModels); addpath('compiled'); addpath('supplem'); ModelNo = 100; % Select a model @@ -19,7 +17,7 @@ % generate 4D phantom (modify your PATH bellow): curDir = pwd; mainDir = fileparts(curDir); -pathtoLibrary = sprintf([fsep 'functions' fsep 'models' fsep 'Phantom3DLibrary.dat'], 1i); +pathtoLibrary = sprintf([fsep '..' fsep 'PhantomLibrary' fsep 'models' fsep 'Phantom3DLibrary.dat'], 1i); pathTP = strcat(mainDir, pathtoLibrary); % path to TomoPhantom parameters file [G] = TomoP3DModel(ModelNo,N,pathTP); @@ -38,11 +36,10 @@ % generate 4D phantom (modify your PATH bellow): curDir = pwd; mainDir = fileparts(curDir); -pathtoLibrary = sprintf([fsep 'functions' fsep 'models' fsep 'Phantom3DLibrary.dat'], 1i); +pathtoLibrary = sprintf([fsep '..' fsep 'PhantomLibrary' fsep 'models' fsep 'Phantom3DLibrary.dat'], 1i); pathTP = strcat(mainDir, pathtoLibrary); % path to TomoPhantom parameters file [G] = TomoP3DModel(ModelNo,N,pathTP); -%% sliceM = round(0.5*N); figure(2); ll = 10; diff --git a/Wrappers/MATLAB/Reconstruction_Demo.m b/Wrappers/MATLAB/Reconstruction_Demo.m index 991b1c5..a85d075 100644 --- a/Wrappers/MATLAB/Reconstruction_Demo.m +++ b/Wrappers/MATLAB/Reconstruction_Demo.m @@ -10,8 +10,6 @@ close all;clc;clear; % adding paths fsep = '/'; -pathtoModels = sprintf(['..' fsep 'functions' fsep 'models' fsep], 1i); -addpath(pathtoModels); addpath('compiled'); addpath('supplem'); ModelNo = 4; % Select a model from Phantom2DLibrary.dat @@ -21,7 +19,7 @@ % Generate 2D phantom: curDir = pwd; mainDir = fileparts(curDir); -pathtoLibrary = sprintf([fsep 'functions' fsep 'models' fsep 'Phantom2DLibrary.dat'], 1i); +pathtoLibrary = sprintf([fsep '..' fsep 'PhantomLibrary' fsep 'models' fsep 'Phantom2DLibrary.dat'], 1i); pathTP = strcat(mainDir, pathtoLibrary); % path to TomoPhantom parameters file [G] = TomoP2DModel(ModelNo,N,pathTP); figure; imagesc(G, [0 1]); daspect([1 1 1]); colormap hot; diff --git a/Wrappers/MATLAB/SpectralPhantomDemo.m b/Wrappers/MATLAB/SpectralPhantomDemo.m index 6e9ee4d..1d3c5e5 100644 --- a/Wrappers/MATLAB/SpectralPhantomDemo.m +++ b/Wrappers/MATLAB/SpectralPhantomDemo.m @@ -10,8 +10,6 @@ close all;clc;clear; % adding paths fsep = '/'; -pathtoModels = sprintf(['..' fsep 'functions' fsep 'models' fsep], 1i); -addpath(pathtoModels); addpath('compiled'); addpath('supplem'); @@ -22,7 +20,7 @@ % generate the 2D phantom: curDir = pwd; mainDir = fileparts(curDir); -pathtoLibrary = sprintf([fsep 'functions' fsep 'models' fsep 'Phantom2DLibrary.dat'], 1i); +pathtoLibrary = sprintf([fsep '..' fsep 'PhantomLibrary' fsep 'models' fsep 'Phantom2DLibrary.dat'], 1i); pathTP = strcat(mainDir, pathtoLibrary); % path to TomoPhantom parameters file [G] = TomoP2DModel(ModelNo,N,pathTP); diff --git a/Wrappers/MATLAB/install/compile_mex_linux.m b/Wrappers/MATLAB/install/compile_mex_linux.m old mode 100644 new mode 100755 index be5c5ba..a9a669b --- a/Wrappers/MATLAB/install/compile_mex_linux.m +++ b/Wrappers/MATLAB/install/compile_mex_linux.m @@ -6,10 +6,13 @@ mkdir compiled -PathFunc = sprintf(['..' fsep 'functions' fsep], 1i); -cd(PathFunc); +PathFunc = sprintf(['..' fsep '..' fsep 'Core' fsep], 1i); -Pathmove = sprintf(['..' fsep 'matlab' fsep 'compiled' fsep], 1i); +copyfile(PathFunc, 'mex_wrappers'); + +cd('mex_wrappers'); + +Pathmove = sprintf(['..' fsep 'compiled' fsep], 1i); fprintf('%s \n', 'Building functions...'); mex TomoP2DModel.c TomoP2DModel_core.c utils.c CFLAGS="\$CFLAGS -fopenmp -Wall -std=c99" LDFLAGS="\$LDFLAGS -fopenmp" @@ -24,7 +27,11 @@ movefile('TomoP3DModel.mex*',Pathmove); mex TomoP3DObject.c TomoP3DModel_core.c utils.c CFLAGS="\$CFLAGS -fopenmp -Wall -std=c99" LDFLAGS="\$LDFLAGS -fopenmp" movefile('TomoP3DObject.mex*',Pathmove); +mex TomoP2DSinoNum.c TomoP2DSinoNum_core.c utils.c CFLAGS="\$CFLAGS -fopenmp -Wall -std=c99" LDFLAGS="\$LDFLAGS -fopenmp" +movefile('TomoP2DSinoNum.mex*',Pathmove); + +delete TomoP2DModel_core* TomoP2DModelSino_core* TomoP3DModel_core* TomoP2DSinoNum_core* CCPiDefines.h utils* CMakeLists.txt; + fprintf('%s \n', 'All compiled!'); cd(UpPath); -cd matlab \ No newline at end of file diff --git a/Wrappers/MATLAB/install/compile_mex_windows.m b/Wrappers/MATLAB/install/compile_mex_windows.m index 574372a..2a3ce5c 100644 --- a/Wrappers/MATLAB/install/compile_mex_windows.m +++ b/Wrappers/MATLAB/install/compile_mex_windows.m @@ -14,10 +14,13 @@ mkdir compiled -PathFunc = sprintf(['..' fsep 'functions' fsep], 1i); -cd(PathFunc); -Pathmove = sprintf(['..' fsep 'matlab' fsep 'compiled' fsep], 1i); +PathFunc = sprintf(['..' fsep '..' fsep 'Core' fsep], 1i); +copyfile(PathFunc, 'mex_wrappers'); + +cd('mex_wrappers'); + +Pathmove = sprintf(['..' fsep 'compiled' fsep], 1i); % One way to compile using Matlab-installed MinGW: fprintf('%s \n', 'Building functions...'); mex TomoP2DModel.c TomoP2DModel_core.c utils.c COMPFLAGS="\$COMPFLAGS -fopenmp -Wall -std=c99" @@ -32,6 +35,11 @@ movefile('TomoP3DModel.mex*',Pathmove); mex TomoP3DObject.c TomoP3DModel_core.c utils.c COMPFLAGS="\$COMPFLAGS -fopenmp -Wall -std=c99" movefile('TomoP3DObject.mex*',Pathmove); +mex TomoP2DSinoNum.c TomoP2DSinoNum_core.c utils.c COMPFLAGS="\$COMPFLAGS -fopenmp -Wall -std=c99" +movefile('TomoP2DSinoNum.mex*',Pathmove); + + +delete TomoP2DModel_core* TomoP2DModelSino_core* TomoP3DModel_core* TomoP2DSinoNum_core* CCPiDefines.h utils* CMakeLists.txt; %%% The second approach to compile using TDM-GCC which follows this %%% discussion: @@ -55,7 +63,9 @@ % movefile('TomoP3DModel.mex*',Pathmove); % mex C:\TDMGCC\lib\gcc\x86_64-w64-mingw32\5.1.0\libgomp.a CXXFLAGS="$CXXFLAGS -std=c++11 -fopenmp" TomoP3DObject.c TomoP3DModel_core.c utils.c % movefile('TomoP3DObject.mex*',Pathmove); -% fprintf('%s \n', 'All compiled!'); +% mex C:\TDMGCC\lib\gcc\x86_64-w64-mingw32\5.1.0\libgomp.a CXXFLAGS="$CXXFLAGS -std=c++11 -fopenmp" TomoP2DSinoNum.c TomoP2DSinoNum_core.c utils.c +% movefile('TomoP2DSinoNum.mex*',Pathmove); +% delete TomoP2DModel_core* TomoP2DModelSino_core* TomoP3DModel_core* TomoP2DSinoNum_core* CCPiDefines.h utils* CMakeLists.txt; -cd(UpPath); -cd matlab +fprintf('%s \n', 'All compiled!'); +cd(UpPath); \ No newline at end of file diff --git a/matlab/mex_wrappers/TomoP2DModel.c b/Wrappers/MATLAB/mex_wrappers/TomoP2DModel.c similarity index 100% rename from matlab/mex_wrappers/TomoP2DModel.c rename to Wrappers/MATLAB/mex_wrappers/TomoP2DModel.c diff --git a/matlab/mex_wrappers/TomoP2DModelSino.c b/Wrappers/MATLAB/mex_wrappers/TomoP2DModelSino.c similarity index 100% rename from matlab/mex_wrappers/TomoP2DModelSino.c rename to Wrappers/MATLAB/mex_wrappers/TomoP2DModelSino.c diff --git a/matlab/mex_wrappers/TomoP2DObject.c b/Wrappers/MATLAB/mex_wrappers/TomoP2DObject.c similarity index 100% rename from matlab/mex_wrappers/TomoP2DObject.c rename to Wrappers/MATLAB/mex_wrappers/TomoP2DObject.c diff --git a/matlab/mex_wrappers/TomoP2DObjectSino.c b/Wrappers/MATLAB/mex_wrappers/TomoP2DObjectSino.c similarity index 100% rename from matlab/mex_wrappers/TomoP2DObjectSino.c rename to Wrappers/MATLAB/mex_wrappers/TomoP2DObjectSino.c diff --git a/matlab/mex_wrappers/TomoP2DSinoNum.c b/Wrappers/MATLAB/mex_wrappers/TomoP2DSinoNum.c similarity index 100% rename from matlab/mex_wrappers/TomoP2DSinoNum.c rename to Wrappers/MATLAB/mex_wrappers/TomoP2DSinoNum.c diff --git a/matlab/mex_wrappers/TomoP3DModel.c b/Wrappers/MATLAB/mex_wrappers/TomoP3DModel.c similarity index 100% rename from matlab/mex_wrappers/TomoP3DModel.c rename to Wrappers/MATLAB/mex_wrappers/TomoP3DModel.c diff --git a/matlab/mex_wrappers/TomoP3DObject.c b/Wrappers/MATLAB/mex_wrappers/TomoP3DObject.c similarity index 100% rename from matlab/mex_wrappers/TomoP3DObject.c rename to Wrappers/MATLAB/mex_wrappers/TomoP3DObject.c diff --git a/Wrappers/Python/Demos/DemoModel.py b/Wrappers/Python/Demos/DemoModel.py index e8c2cb6..f2110cc 100644 --- a/Wrappers/Python/Demos/DemoModel.py +++ b/Wrappers/Python/Demos/DemoModel.py @@ -22,7 +22,7 @@ model = 1 # selecting a model N_size = 512 #specify a full path to the parameters file -pathTP = '../../functions/models/Phantom2DLibrary.dat' +pathTP = '../../../PhantomLibrary/models/Phantom2DLibrary.dat' #objlist = modelfile2Dtolist(pathTP, model) # one can extract parameters #This will generate a N_size x N_size phantom (2D) phantom_2D = TomoP2D.Model(model, N_size, pathTP) diff --git a/Wrappers/Python/Demos/DemoModel2.py b/Wrappers/Python/Demos/DemoModel2.py index 392f667..29e387e 100644 --- a/Wrappers/Python/Demos/DemoModel2.py +++ b/Wrappers/Python/Demos/DemoModel2.py @@ -10,9 +10,11 @@ >>>> Prerequisites: ASTRA toolbox, if one needs to do reconstruction <<<<< Main difference from DemoModel.py is that we extract all parameters from the -library file using Python and then pass it to the Object function instead of model +library file using Python and then pass it to the Object function instead of model. +This can be helpful if one would like to avoid using library files and can +pass parameters directly into object function -Run demo from the folder "Demos" +!Run script from "Demos" folder in order to ensure a correct path to *dat file! @author: Daniil Kazantsev """ @@ -24,7 +26,7 @@ model = 11 N_size = 512 #specify a full path to the parameters file -pathTP = '../../functions/models/Phantom2DLibrary.dat' +pathTP = '../../../PhantomLibrary/models/Phantom2DLibrary.dat' objlist = modelfile2Dtolist(pathTP, model) # extract parameters using Python #This will generate a N_size x N_size phantom (2D) phantom_2D = TomoP2D.Object(N_size, objlist) @@ -38,8 +40,8 @@ # create sinogram analytically angles_num = int(0.5*np.pi*N_size); # angles number -angles = np.linspace(0,180,angles_num,dtype='float32') -angles_rad = angles*(np.pi/180) +angles = np.linspace(0.0,179.9,angles_num,dtype='float32') +angles_rad = angles*(np.pi/180.0) P = int(np.sqrt(2)*N_size) #detectors sino_an = TomoP2D.ObjectSino(N_size, P, angles, objlist) diff --git a/Wrappers/Python/Demos/DemoModel3D.py b/Wrappers/Python/Demos/DemoModel3D.py new file mode 100644 index 0000000..0ad8cfc --- /dev/null +++ b/Wrappers/Python/Demos/DemoModel3D.py @@ -0,0 +1,81 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +GPLv3 license (ASTRA toolbox) +Note that the TomoPhantom package is released under Apache License, Version 2.0 + +Script to generate 3D analytical phantoms (wip: generation of 3D projection data ) +If one needs to modify/add phantoms, please edit Phantom3DLibrary.dat + +!Run script from "Demos" folder in order to ensure a correct path to *dat file! + +@author: Daniil Kazantsev +""" +import timeit +from tomophantom import TomoP3D +import matplotlib.pyplot as plt + +print ("Building 3D phantom using TomoPhantom software") +tic=timeit.default_timer() +model = 2 +# Define phantom dimensions using a scalar (cubic) or a tuple [N1, N2, N3] +N_size = 256 # or as a tuple of a custom size (256,256,256) +#specify a full path to the parameters file +pathTP3 = '../../../PhantomLibrary/models/Phantom3DLibrary.dat' +#This will generate a N_size x N_size x N_size phantom (3D) or non-cubic phantom +phantom_tm = TomoP3D.Model(model, N_size, pathTP3) +toc=timeit.default_timer() +Run_time = toc - tic +print("Phantom has been built in {} seconds".format(Run_time)) +sliceSel = int(0.5*N_size) +#plt.gray() +plt.figure(1) +plt.subplot(131) +plt.imshow(phantom_tm[sliceSel,:,:],vmin=0, vmax=1) +plt.title('3D Phantom, axial view') + +plt.subplot(132) +plt.imshow(phantom_tm[:,sliceSel,:],vmin=0, vmax=1) +plt.title('3D Phantom, coronal view') + +plt.subplot(133) +plt.imshow(phantom_tm[:,:,sliceSel],vmin=0, vmax=1) +plt.title('3D Phantom, sagittal view') +plt.show() + +#%% +# The capability of building a subset of vertical slices out of 3D phantom (faster) +import timeit +from tomophantom import TomoP3D +import matplotlib.pyplot as plt + +print ("Building a subset of 3D phantom using TomoPhantom software") +tic=timeit.default_timer() +model = 3 +# Define phantom dimensions using a scalar (cubic) or a tuple [Z, Y, X] +DIM = (256,256,256) # full dimension of required phantom (z, y, x) +DIM_z = (94, 158) # selected vertical subset (a slab) of the phantom +#specify a full path to the parameters file +pathTP3 = '../../functions/models/Phantom3DLibrary.dat' +#This will generate a N1 x N2 x N_slab phantom (3D) +phantom_tm = TomoP3D.ModelSub(model, DIM, DIM_z, pathTP3) +#phantom_tm = TomoP3D.Model(model, DIM, pathTP3) +toc=timeit.default_timer() +Run_time = toc - tic +print("Phantom has been built in {} seconds".format(Run_time)) +sliceSel = 35 +#plt.gray() +plt.figure(2) +plt.subplot(131) +plt.imshow(phantom_tm[sliceSel,:,:],vmin=0, vmax=1) +plt.title('3D Phantom, axial view') + +plt.subplot(132) +plt.imshow(phantom_tm[:,128,:],vmin=0, vmax=1) +plt.title('3D Phantom, coronal view') + +plt.subplot(133) +plt.imshow(phantom_tm[:,:,128],vmin=0, vmax=1) +plt.title('3D Phantom, sagittal view') +plt.show() +#%% \ No newline at end of file diff --git a/Wrappers/Python/Demos/DemoModel4D.py b/Wrappers/Python/Demos/DemoModel4D.py new file mode 100644 index 0000000..5066b60 --- /dev/null +++ b/Wrappers/Python/Demos/DemoModel4D.py @@ -0,0 +1,88 @@ +#!/usr/bin/env python3 +# -*- coding: utf-8 -*- +""" +GPLv3 license (ASTRA toolbox) +Note that the TomoPhantom package is released under Apache License, Version 2.0 + +Script to generate 4D analytical phantoms (wip: generation of 4D projection data ) +If one needs to modify/add phantoms, please edit Phantom3DLibrary.dat + +!Run script from "Demos" folder in order to ensure a correct path to *dat file! + +@author: Daniil Kazantsev +""" +import timeit +from tomophantom import TomoP3D +import matplotlib.pyplot as plt +import numpy as np + +print ("Building 4D phantom using TomoPhantom software") +tic=timeit.default_timer() +model = 100 +# Define phantom dimensions using a scalar (cubic) or a tuple [N1, N2, N3] +N_size = 256 # or as a tuple of a custom size (256,256,256) +#specify a full path to the parameters file +pathTP3 = '../../../PhantomLibrary/models/Phantom3DLibrary.dat' +#This will generate a Time x N_size x N_size x N_size phantom (4D) +phantom_tm = TomoP3D.ModelTemporal(model, N_size, pathTP3) +toc=timeit.default_timer() +Run_time = toc - tic +print("Phantom has been built in {} seconds".format(Run_time)) + +for i in range(0,np.size(phantom_tm,0)): + sliceSel = int(0.5*N_size) + #plt.gray() + plt.figure(1) + plt.subplot(131) + plt.imshow(phantom_tm[i,sliceSel,:,:],vmin=0, vmax=1) + plt.title('3D Phantom, axial view') + + plt.subplot(132) + plt.imshow(phantom_tm[i,:,sliceSel,:],vmin=0, vmax=1) + plt.title('3D Phantom, coronal view') + + plt.subplot(133) + plt.imshow(phantom_tm[i,:,:,sliceSel],vmin=0, vmax=1) + plt.title('3D Phantom, sagittal view') + plt.show() + plt.pause(0.3) + +#%% +# The capability of building a subset of vertical slices out of 4D phantom (faster) +import timeit +from tomophantom import TomoP3D +import matplotlib.pyplot as plt + +print ("Building a subset of 3D phantom using TomoPhantom software") +tic=timeit.default_timer() +model = 101 +# Define phantom dimensions using a scalar (cubic) or a tuple [Z, Y, X] +DIM = (256,256,256) # full dimension of required phantom (z, y, x) +DIM_z = (94, 158) # selected vertical subset (a slab) of the phantom +#specify a full path to the parameters file +pathTP3 = '../../functions/models/Phantom3DLibrary.dat' +#This will generate a Time x N1 x N2 x N_slab phantom (4D) +phantom_tm = TomoP3D.ModelTemporalSub(model, DIM, DIM_z, pathTP3) +#phantom_tm = TomoP3D.Model(model, DIM, pathTP3) +toc=timeit.default_timer() +Run_time = toc - tic +print("Phantom has been built in {} seconds".format(Run_time)) + +for i in range(0,np.size(phantom_tm,0)): + sliceSel = 32 + #plt.gray() + plt.figure(1) + plt.subplot(131) + plt.imshow(phantom_tm[i,sliceSel,:,:],vmin=0, vmax=1) + plt.title('3D Phantom, axial view') + + plt.subplot(132) + plt.imshow(phantom_tm[i,:,70,:],vmin=0, vmax=1) + plt.title('3D Phantom, coronal view') + + plt.subplot(133) + plt.imshow(phantom_tm[i,:,:,70],vmin=0, vmax=1) + plt.title('3D Phantom, sagittal view') + plt.show() + plt.pause(0.5) +#%% \ No newline at end of file diff --git a/Wrappers/Python/Demos/DemoModel_temporal.py b/Wrappers/Python/Demos/DemoModel_temporal.py index 5119955..20f32b4 100644 --- a/Wrappers/Python/Demos/DemoModel_temporal.py +++ b/Wrappers/Python/Demos/DemoModel_temporal.py @@ -20,7 +20,7 @@ N_size = 512 timeframes = 25 #specify a full path to the parameters file -pathTP = '../../functions/models/Phantom2DLibrary.dat' +pathTP = '../../../PhantomLibrary/models/Phantom2DLibrary.dat' #This will generate a N_size x N_size x Time frames phantom (2D + time) phantom_2Dt = TomoP2D.ModelTemporal(model, N_size, pathTP) diff --git a/Wrappers/Python/Demos/DemoReconASTRA.py b/Wrappers/Python/Demos/DemoReconASTRA.py index 277edbc..2de570e 100644 --- a/Wrappers/Python/Demos/DemoReconASTRA.py +++ b/Wrappers/Python/Demos/DemoReconASTRA.py @@ -22,7 +22,7 @@ model = 4 # select a model N_size = 512 #specify a full path to the parameters file -pathTP = '../../functions/models/Phantom2DLibrary.dat' +pathTP = '../../../PhantomLibrary/models/Phantom2DLibrary.dat' #objlist = modelfile2Dtolist(pathTP, model) # one can extract parameters #This will generate a N_size x N_size phantom (2D) phantom_2D = TomoP2D.Model(model, N_size, pathTP) diff --git a/Wrappers/Python/Demos/Reconstruction_Demo.py b/Wrappers/Python/Demos/DemoReconTomoPy.py similarity index 71% rename from Wrappers/Python/Demos/Reconstruction_Demo.py rename to Wrappers/Python/Demos/DemoReconTomoPy.py index d2d802c..66c2268 100644 --- a/Wrappers/Python/Demos/Reconstruction_Demo.py +++ b/Wrappers/Python/Demos/DemoReconTomoPy.py @@ -5,9 +5,9 @@ Note that the TomoPhantom package is released under Apache License, Version 2.0 Script to generate 2D analytical phantoms and their sinograms with added noise and artifacts -Sinograms then reconstructed using ASTRA TOOLBOX and TomoPy packages +Sinograms then reconstructed using TomoPy package ->>>> Prerequisites: ASTRA toolbox and TomoPy packages <<<<< +>>>> Prerequisites: TomoPy package <<<<< This demo demonstrates frequent inaccuracies which are accosiated with X-ray imaging: zingers, rings and noise @@ -21,7 +21,7 @@ model = 4 # select a model N_size = 512 #specify a full path to the parameters file -pathTP = '../../functions/models/Phantom2DLibrary.dat' +pathTP = '../../../PhantomLibrary/models/Phantom2DLibrary.dat' #objlist = modelfile2Dtolist(pathTP, model) # one can extract parameters #This will generate a N_size x N_size phantom (2D) phantom_2D = TomoP2D.Model(model, N_size, pathTP) @@ -70,33 +70,6 @@ plt.imshow(noisy_zing_stripe,cmap="gray") plt.colorbar(ticks=[0, 150, 250], orientation='vertical') plt.title('{}''{}'.format('Analytical noisy sinogram with artifacts.',model)) -#%% -print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") -print ("Reconstructing analytical sinogram using FBP (ASTRA-TOOLBOX)...") -print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") -from tomophantom.supp.astraOP import AstraTools - -Atools = AstraTools(P, angles_rad, N_size, 'cpu') # initiate a class object - -FBPrec_ideal = Atools.fbp2D(sino_an) # ideal reconstruction -FBPrec_error = Atools.fbp2D(noisy_zing_stripe) # error reconstruction - -plt.figure() -plt.subplot(121) -plt.imshow(FBPrec_ideal, vmin=0, vmax=1, cmap="gray") -plt.colorbar(ticks=[0, 0.5, 1], orientation='vertical') -plt.title('Ideal reconstruction (ASTRA)') -plt.subplot(122) -plt.imshow(FBPrec_error, vmin=0, vmax=1, cmap="gray") -plt.colorbar(ticks=[0, 0.5, 1], orientation='vertical') -plt.title('Erroneous Reconstruction (ASTRA)') -plt.show() - -plt.figure() -plt.imshow(abs(FBPrec_ideal-FBPrec_error), vmin=0, vmax=1, cmap="gray") -plt.colorbar(ticks=[0, 0.5, 1], orientation='vertical') -plt.title('FBP reconsrtuction differences') -rmse2 = np.linalg.norm(FBPrec_ideal-FBPrec_error)/np.linalg.norm(FBPrec_error) #%% print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") @@ -111,7 +84,6 @@ sinoTP[:,0,:] = noisy_zing_stripe reconTomoPy_noisy = tomopy.recon(sinoTP, angles_rad, center=rot_center, algorithm='gridrec') - plt.figure() plt.subplot(121) plt.imshow(reconTomoPy_ideal[0,:,:], vmin=0, vmax=1, cmap="gray") diff --git a/Wrappers/Python/Demos/__pycache__/astraOP.cpython-35.pyc b/Wrappers/Python/Demos/__pycache__/astraOP.cpython-35.pyc deleted file mode 100644 index ffac956..0000000 Binary files a/Wrappers/Python/Demos/__pycache__/astraOP.cpython-35.pyc and /dev/null differ diff --git a/Wrappers/Python/Demos/__pycache__/libraryToDict.cpython-35.pyc b/Wrappers/Python/Demos/__pycache__/libraryToDict.cpython-35.pyc deleted file mode 100644 index 02d88db..0000000 Binary files a/Wrappers/Python/Demos/__pycache__/libraryToDict.cpython-35.pyc and /dev/null differ diff --git a/Wrappers/Python/Demos/libraryToDict.py b/Wrappers/Python/Demos/libraryToDict.py index 72d7350..774b910 100644 --- a/Wrappers/Python/Demos/libraryToDict.py +++ b/Wrappers/Python/Demos/libraryToDict.py @@ -26,7 +26,7 @@ def modelfile2Dtolist(filepath, model): # print (f[i]) break i += 1 - print (i) + #print (i) objectlist = [] diff --git a/Wrappers/Python/Demos/model4.txt b/Wrappers/Python/Demos/model4.txt deleted file mode 100644 index 2601737..0000000 --- a/Wrappers/Python/Demos/model4.txt +++ /dev/null @@ -1,41 +0,0 @@ -Object : ellipse .3e00 0e00 0.0e00 .9500 .95e00 00.e00; -Object : ellipse 1.0 0.0 0.45 0.1 0.1 00; -Object : ellipse 1.0 0.3 0.45 0.1 0.1 00; -Object : ellipse 1.0 -0.3 0.45 0.1 0.1 00; -Object : ellipse 1.0 -0.6 0.45 0.1 0.1 00; -Object : ellipse 1.0 0.6 0.45 0.1 0.1 00; -Object : ellipse 0.98 -0.6 0.1 0.06 0.06 00; -Object : ellipse 0.98 -0.8 0.1 0.06 0.06 00; -Object : ellipse 0.98 -0.4 0.1 0.06 0.06 00; -Object : ellipse 0.98 -0.2 0.1 0.06 0.06 00; -Object : ellipse 0.98 -0.0 0.1 0.06 0.06 00; -Object : ellipse 0.98 0.2 0.1 0.06 0.06 00; -Object : ellipse 0.98 0.4 0.1 0.06 0.06 00; -Object : ellipse 0.98 0.6 0.1 0.06 0.06 00; -Object : ellipse 0.98 0.8 0.1 0.06 0.06 00; -Object : ellipse 0.97 0.7 -0.1 0.05 0.05 00; -Object : ellipse 0.97 0.5 -0.1 0.05 0.05 00; -Object : ellipse 0.97 0.3 -0.1 0.05 0.05 00; -Object : ellipse 0.97 0.1 -0.1 0.05 0.05 00; -Object : ellipse 0.97 -0.1 -0.1 0.05 0.05 00; -Object : ellipse 0.97 -0.3 -0.1 0.05 0.05 00; -Object : ellipse 0.97 -0.5 -0.1 0.05 0.05 00; -Object : ellipse 0.97 -0.7 -0.1 0.05 0.05 00; -Object : ellipse 0.96 -0.6 -0.3 0.04 0.04 00; -Object : ellipse 0.96 -0.4 -0.3 0.04 0.04 00; -Object : ellipse 0.96 -0.2 -0.3 0.04 0.04 00; -Object : ellipse 0.96 0.0 -0.3 0.04 0.04 00; -Object : ellipse 0.96 0.2 -0.3 0.04 0.04 00; -Object : ellipse 0.96 0.4 -0.3 0.04 0.04 00; -Object : ellipse 0.96 0.6 -0.3 0.04 0.04 00; -Object : ellipse 0.95 0.5 -0.5 0.03 0.03 00; -Object : ellipse 0.95 0.3 -0.5 0.03 0.03 00; -Object : ellipse 0.95 0.1 -0.5 0.03 0.03 00; -Object : ellipse 0.95 -0.1 -0.5 0.03 0.03 00; -Object : ellipse 0.95 -0.3 -0.5 0.03 0.03 00; -Object : ellipse 0.95 -0.5 -0.5 0.03 0.03 00; -Object : ellipse 0.94 -0.4 -0.7 0.015 0.015 00; -Object : ellipse 0.94 -0.2 -0.7 0.015 0.015 00; -Object : ellipse 0.94 0.0 -0.7 0.015 0.015 00; -Object : ellipse 0.94 0.2 -0.7 0.015 0.015 00; -Object : ellipse 0.94 0.4 -0.7 0.015 0.015 00; \ No newline at end of file diff --git a/Wrappers/Python/Demos/test_ModelObject.py b/Wrappers/Python/Demos/test_ModelObject.py deleted file mode 100644 index b403caf..0000000 --- a/Wrappers/Python/Demos/test_ModelObject.py +++ /dev/null @@ -1,147 +0,0 @@ -#!/usr/bin/env python3 -# -*- coding: utf-8 -*- -""" -GPLv3 license (ASTRA toolbox) -Note that the TomoPhantom package is released under Apache License, Version 2.0 - -Script to generate 2D/3D analytical phantoms and their sinograms -If one needs to modify/add phantoms, please edit Phantom2DLibrary.dat or -Phantom3DLibrary.dat ->>>> Prerequisites: ASTRA toolbox, if one needs to do reconstruction <<<<< - -Run demo from the folder "Demos" - -@author: Daniil Kazantsev -""" -import numpy as np -import matplotlib.pyplot as plt -from tomophantom import TomoP2D -from astraOP import AstraTools -from libraryToDict import modelfile2Dtolist - - -#%% -model = 10 -N_size = 512 -#specify a full path to the parameters file -pathTP = '../../functions/models/Phantom2DLibrary.dat' -#This will generate a N_size x N_size phantom (2D) -phantom_2D = TomoP2D.Model(model, N_size, pathTP) -objlist = modelfile2Dtolist(pathTP, model) -print (objlist) -phantom_2Da = TomoP2D.Object(N_size, objlist) - -rows = 1 -cols = 2 -current = 1 -plt.close('all') -fig = plt.figure(1) -a=fig.add_subplot(rows,cols,current) -plt.rcParams.update({'font.size': 21}) -plt.imshow(phantom_2D, vmin=0, vmax=1, cmap="BuPu") -plt.colorbar(ticks=[0, 0.5, 1], orientation='vertical') -plt.title('{}''{}'.format('using dat',model)) - -a=fig.add_subplot(rows,cols,current+1) -plt.rcParams.update({'font.size': 21}) -plt.imshow(phantom_2Da, vmin=0, vmax=1, cmap="BuPu") -plt.colorbar(ticks=[0, 0.5, 1], orientation='vertical') -plt.title('{}''{}'.format('using dicts',model)) - -plt.show() - -""" -#%% -# create sinogram analytically -angles_num = int(0.5*np.pi*N_size); # angles number -angles = np.linspace(0,180,angles_num,dtype='float32') -angles_rad = angles*(np.pi/180) -P = int(np.sqrt(2)*N_size) #detectors - -sino_an = TomoP2D.ModelSino(model, N_size, P, angles, pathTP) - -plt.figure(2) -plt.rcParams.update({'font.size': 21}) -plt.imshow(sino_an, cmap="BuPu") -plt.colorbar(ticks=[0, 150, 250], orientation='vertical') -plt.title('{}''{}'.format('Analytical sinogram of model no.',model)) -#%% -Atools = AstraTools(P, angles_rad - 0.5*np.pi, N_size, 'cpu') # initiate a class object -sino_num_ASTRA = Atools.forwproj(phantom_2D) # generate numerical sino (Ax) -#x = Atools.backproj(sino_an) # generate backprojection (A'b) - -plt.figure(2) -plt.subplot(121) -plt.imshow(sino_an,cmap="BuPu") -plt.title('Analytical sinogram') -plt.subplot(122) -plt.imshow(sino_num_ASTRA,cmap="BuPu") -plt.title('Numerical sinogram') -plt.show() -#%% -print ("Reconstructing analytical sinogram using FBP (astra TB)...") -FBPrec1 = Atools.fbp2D(sino_an) - -plt.figure(3) -plt.imshow(FBPrec1, vmin=0, vmax=1, cmap="BuPu") -plt.colorbar(ticks=[0, 0.5, 1], orientation='vertical') -plt.title('FBP Reconstructed Phantom (analyt)') -#%% -print ("Reconstructing numerical sinogram using FBP (astra TB)...") -FBPrec2 = Atools.fbp2D(sino_num_ASTRA) - -plt.figure(4) -plt.imshow(FBPrec2, vmin=0, vmax=1, cmap="BuPu") -plt.colorbar(ticks=[0, 0.5, 1], orientation='vertical') -plt.title('FBP Reconstructed Phantom (numeric)') -#%% -plt.figure(5) -plt.imshow(abs(FBPrec1-FBPrec2), vmin=0, vmax=0.05, cmap="BuPu") -plt.colorbar(ticks=[0, 0.02, 0.05], orientation='vertical') -plt.title('FBP rec differences') -#%% - -""" - -""" -print ("Reconstructing using SIRT...") -SIRTrec = Atools.sirt2D(sino_an, 100) - -plt.figure(4) -plt.imshow(SIRTrec, vmin=0, vmax=1,cmap="BuPu") -plt.title('SIRT Reconstructed Phantom') -""" -#%% -""" -import timeit -from tomophantom import TomoP3D -import matplotlib.pyplot as plt - -print ("Building 3D phantom using TomoPhantom software") -tic=timeit.default_timer() -model = 1 -N_size = 256 -#specify a full path to the parameters file -pathTP3 = '../../functions/models/Phantom3DLibrary.dat' -#This will generate a N_size x N_size x N_size phantom (3D) -phantom_tm = TomoP3D.Model(model, N_size, pathTP3) -toc=timeit.default_timer() -Run_time = toc - tic -print("Phantom has been built in {} seconds".format(Run_time)) - -sliceSel = int(0.5*N_size) -#plt.gray() -plt.figure(5) -plt.subplot(131) -plt.imshow(phantom_tm[sliceSel,:,:],vmin=0, vmax=1) -plt.title('3D Phantom, axial view') - -plt.subplot(132) -plt.imshow(phantom_tm[:,sliceSel,:],vmin=0, vmax=1) -plt.title('3D Phantom, coronal view') - -plt.subplot(133) -plt.imshow(phantom_tm[:,:,sliceSel],vmin=0, vmax=1) -plt.title('3D Phantom, sagittal view') -plt.show() -""" diff --git a/Wrappers/Python/setup.py.in b/Wrappers/Python/setup.py.in index 70a5bb7..e169e52 100644 --- a/Wrappers/Python/setup.py.in +++ b/Wrappers/Python/setup.py.in @@ -24,7 +24,7 @@ import numpy import platform #import sys -version = '1.0' +version = '1.1' extra_include_dirs = [numpy.get_include(), '@CMAKE_SOURCE_DIR@/Core'] extra_library_dirs = ['@LIBRARY_LIB@'] extra_compile_args = [] @@ -45,8 +45,6 @@ setup( cmdclass = {'build_ext': build_ext}, ext_modules = cythonize([ Extension("tomophantom.TomoP3D", sources = [ "src/TomoP3D.pyx", - #"@CMAKE_SOURCE_DIR@/Core/TomoP3DModel_core.c", - #"@CMAKE_SOURCE_DIR@/Core/utils.c" ], include_dirs = extra_include_dirs, library_dirs = extra_library_dirs, @@ -55,9 +53,6 @@ setup( extra_link_args = extra_link_args), Extension("tomophantom.TomoP2D", sources = [ "src/TomoP2D.pyx", - #"@CMAKE_SOURCE_DIR@/Core/TomoP2DModel_core.c", - #"@CMAKE_SOURCE_DIR@/Core/TomoP2DModelSino_core.c", - #"@CMAKE_SOURCE_DIR@/Core/utils.c" ], include_dirs = extra_include_dirs, library_dirs = extra_library_dirs, diff --git a/Wrappers/Python/tomophantom/supp/artifacts.py b/Wrappers/Python/tomophantom/supp/artifacts.py index 26a32cb..7b73af1 100644 --- a/Wrappers/Python/tomophantom/supp/artifacts.py +++ b/Wrappers/Python/tomophantom/supp/artifacts.py @@ -26,14 +26,16 @@ def noise(self, sigma, noisetype): if noisetype == 'Gaussian': # add normal Gaussian noise sino_noisy += np.random.normal(loc = 0 ,scale = sigma * sino_noisy, size = np.shape(sino_noisy)) + sino_noisy[sino_noisy<0] = 0 elif noisetype == 'Poisson': # add Poisson noise - maxSino = np.max(sino_noisy) + maxSino = np.max(self.sinogram) if maxSino > 0: - sino_noisy = sino_noisy/maxSino + sino_noisy = (self.sinogram)/maxSino dataExp = sigma*np.exp(-sino_noisy) # noiseless raw data sino_noisy = np.random.poisson(dataExp) #adding Poisson noise - sino_noisy = -np.log(sino_noisy/np.max(sino_noisy))*maxSino #log corrected data -> sinogram + div_res = np.float32(sino_noisy)/np.max(sino_noisy) + sino_noisy = -np.log(div_res)*maxSino #log corrected data -> sinogram sino_noisy[sino_noisy<0] = 0 else: print ("Select 'Gaussian' or 'Poisson' for noise type") @@ -53,7 +55,7 @@ def zingers(self, percentage, modulus): raise ("Modulus integer must be positive") sino_zingers = self.sinogram length_sino = np.size(sino_zingers) - num_values = int((length_sino)*(percentage/100)) + num_values = int((length_sino)*(np.float32(percentage)/100.0)) sino_zingers_fl = sino_zingers.flatten() for x in range(num_values): randind = random.randint(0,length_sino) # generate random index @@ -77,13 +79,13 @@ def stripes(self, percentage, maxthickness): pass else: raise ("percentage must be larger than zero but smaller than 100") - if 0 <= maxthickness <= 5: + if 0 <= maxthickness <= 10: pass else: - raise ("maximum thickness must be in [0,5] range") + raise ("maximum thickness must be in [0,10] range") sino_stripes = self.sinogram max_intensity = np.max(sino_stripes) - range_detect = int((self.DetectorsDim)*(percentage/100)) + range_detect = int((np.float32(self.DetectorsDim))*(np.float32(percentage)/100.0)) for x in range(range_detect): randind = random.randint(0,self.DetectorsDim) # generate random index randthickness = random.randint(0,maxthickness) #generate random thickness diff --git a/Wrappers/Python/tomophantom/supp/recMod.py b/Wrappers/Python/tomophantom/supp/recMod.py new file mode 100644 index 0000000..c701463 --- /dev/null +++ b/Wrappers/Python/tomophantom/supp/recMod.py @@ -0,0 +1,103 @@ +#!/usr/bin/env python2 +# -*- coding: utf-8 -*- +""" +Created on Thu Sep 27 12:28:29 2018 +A reconstruction class for TomoPhantom, currently includes +-- Fourier Slice Theorem reconstruction (adopted from Tim Day's code) +-- Filtered Back Projection + +@author: Daniil Kazantsev +""" +import numpy as np + +class RecTools: + """ Class for reconstruction using TomoPhantom """ + def __init__(self, DetectorsDim, AnglesVec, ObjSize): + self.DetectorsDim = DetectorsDim # detector dimension + self.AnglesVec = AnglesVec # angles arrray in radians + if ObjSize is tuple: + raise (" Reconstruction is currently for square or cubic objects only ") + else: + self.ObjSize = ObjSize # size of the object + + def fourier(self, sinogram, method): + """ + 2D Reconstruction using Fourier slice theorem (scipy required) + for griddata interpolation module choose nearest, linear or cubic + """ + if sinogram.ndim == 3: + raise ("Fourier method is currently for 2D data only, use FBP if 3D needed ") + else: + pass + if ((method == 'linear') or (method == 'nearest') or (method == 'cubic')): + pass + else: + raise ("For griddata interpolation module choose nearest, linear or cubic ") + import scipy.interpolate + import scipy.fftpack + import scipy.misc + import scipy.ndimage.interpolation + + # Fourier transform the rows of the sinogram, move the DC component to the row's centre + sinogram_fft_rows=scipy.fftpack.fftshift(scipy.fftpack.fft(scipy.fftpack.ifftshift(sinogram,axes=1)),axes=1) + + """ + V = 100 + plt.figure() + plt.subplot(121) + plt.title("Sinogram rows FFT (real)") + plt.imshow(np.real(sinogram_fft_rows),vmin=-V,vmax=V) + plt.subplot(122) + plt.title("Sinogram rows FFT (imag)") + plt.imshow(np.imag(sinogram_fft_rows),vmin=-V,vmax=V) + """ + # Coordinates of sinogram FFT-ed rows' samples in 2D FFT space + a = -self.AnglesVec + r=np.arange(self.DetectorsDim) - self.DetectorsDim/2 + r,a=np.meshgrid(r,a) + r=r.flatten() + a=a.flatten() + srcx=(self.DetectorsDim /2)+r*np.cos(a) + srcy=(self.DetectorsDim /2)+r*np.sin(a) + + # Coordinates of regular grid in 2D FFT space + dstx,dsty=np.meshgrid(np.arange(self.DetectorsDim),np.arange(self.DetectorsDim)) + dstx=dstx.flatten() + dsty=dsty.flatten() + + """ + V = 100 + plt.figure() + plt.title("Sinogram samples in 2D FFT (abs)") + plt.scatter(srcx, srcy,c=np.absolute(sinogram_fft_rows.flatten()), marker='.', edgecolor='none', vmin=-V, vmax=V) + """ + # Interpolate the 2D Fourier space grid from the transformed sinogram rows + fft2=scipy.interpolate.griddata((srcy,srcx), sinogram_fft_rows.flatten(), (dsty,dstx), method, fill_value=0.0).reshape((self.DetectorsDim,self.DetectorsDim)) + """ + plt.figure() + plt.suptitle("FFT2 space") + plt.subplot(221) + plt.title("Recon (real)") + plt.imshow(np.real(fft2),vmin=-V,vmax=V) + plt.subplot(222) + plt.title("Recon (imag)") + plt.imshow(np.imag(fft2),vmin=-V,vmax=V) + """ + + """ + # Show 2D FFT of target, just for comparison + expected_fft2=scipy.fftpack.fftshift(scipy.fftpack.fft2(scipy.fftpack.ifftshift(phantom_2D))) + + plt.subplot(223) + plt.title("Expected (real)") + plt.imshow(np.real(expected_fft2),vmin=-V,vmax=V) + plt.subplot(224) + plt.title("Expected (imag)") + plt.imshow(np.imag(expected_fft2),vmin=-V,vmax=V) + """ + # Transform from 2D Fourier space back to a reconstruction of the target + recon=np.real(scipy.fftpack.fftshift(scipy.fftpack.ifft2(scipy.fftpack.ifftshift(fft2)))) + + # Cropping reconstruction to size of the original image + image = recon[int(((self.DetectorsDim-self.ObjSize)/2)+1):self.DetectorsDim-int(((self.DetectorsDim-self.ObjSize)/2)-1),int(((self.DetectorsDim-self.ObjSize)/2)):self.DetectorsDim-int(((self.DetectorsDim-self.ObjSize)/2))] + return image diff --git a/matlab/Model2DGeneratorDemo.m b/matlab/Model2DGeneratorDemo.m deleted file mode 100644 index c49553f..0000000 --- a/matlab/Model2DGeneratorDemo.m +++ /dev/null @@ -1,79 +0,0 @@ -% GPLv3 license (ASTRA toolbox) -% Note that the TomoPhantom package is released under Apache License, Version 2.0 - -% Script to generate 2D analytical phantoms and sinograms (parallel beam) -% that can be used to test reconstruction algorithms without the "Inverse -% Crime" -% If one needs to modify/add phantoms, please edit Phantom2DLibrary.dat -% >>>> Prerequisites: ASTRA toolbox, if one needs to do reconstruction <<<<< - -close all;clc;clear; -% adding paths -fsep = '/'; -pathtoModels = sprintf(['..' fsep 'functions' fsep 'models' fsep], 1i); -addpath(pathtoModels); -addpath('compiled'); addpath('supplem'); - -ModelNo = 1; % Select a model from Phantom2DLibrary.dat -% Define phantom dimensions -N = 512; % x-y size (squared image) - -% Generate 2D phantom: -curDir = pwd; -mainDir = fileparts(curDir); -pathtoLibrary = sprintf([fsep 'functions' fsep 'models' fsep 'Phantom2DLibrary.dat'], 1i); -pathTP = strcat(mainDir, pathtoLibrary); % path to TomoPhantom parameters file -[G] = TomoP2DModel(ModelNo,N,pathTP); -figure; imagesc(G, [0 1]); daspect([1 1 1]); colormap hot; -%% -fprintf('%s \n', 'Generating sinogram analytically and numerically with Matlab (radon)...'); -% generate angles -angles = linspace(0,179.9,round(1.3*N)); % projection angles - -% lets use Matlab's radon function -[F_d,xp] = radon(G,angles); % discrete sinogram -P = size(F_d,1); %detectors dimension -F_d = F_d'; - -% [F_num] = TomoP2DSinoNum(G, P, single(angles)); % numerical sinogram - -% generate the 2D analytical parallel beam sinogram -[F_a] = TomoP2DModelSino(ModelNo, N, P, single(angles), pathTP, 'radon'); - -figure; -subplot(1,2,1); imshow(F_a, []); title('Analytical Sinogram'); -subplot(1,2,2); imshow(F_d, []); title('Numerical Sinogram'); - -% calculate residiual norm (the error is expected since projection models not the same) -err_diff = norm(F_a(:) - F_d(:))./norm(F_d(:)); -fprintf('%s %.4f\n', 'NMSE for sino residuals:', err_diff); - -% reconstructing with FBP (iradon) -FBP_F_a = iradon(F_a',angles,N); -FBP_F_d = iradon(F_d',angles,N); -figure; -subplot(1,2,1); imagesc(FBP_F_a, [0 1]); title('Analytical Sinogram Reconstruction'); daspect([1 1 1]); colormap hot; -subplot(1,2,2); imagesc(FBP_F_d, [0 1]); title('Numerical Sinogram (radon) Reconstruction'); daspect([1 1 1]); colormap hot; -%% -fprintf('%s \n', 'Use the ASTRA-toolbox to generate numerical sinogram...'); -% generate 2D analytical parallel beam sinogram (note the 'astra' opton) -[F_a] = TomoP2DModelSino(ModelNo, N, P, single(angles), pathTP, 'astra'); -[F_num_astra] = sino2Dastra(G, (angles*pi/180), P, N, 'cpu'); - -% calculate residiual norm (the error is expected since projection models not the same) -err_diff = norm(F_a(:) - F_num_astra(:))./norm(F_num_astra(:)); -fprintf('%s %.4f\n', 'NMSE for sino residuals:', err_diff); - -figure; -subplot(1,2,1); imshow(F_a, []); title('Analytical Sinogram'); -subplot(1,2,2); imshow(F_num_astra, []); title('Numerical Sinogram (ASTRA)'); -%% -fprintf('%s \n', 'Reconstruction using the ASTRA-toolbox (FBP)...'); - -rec_an = rec2Dastra(F_a, (angles*pi/180), P, N, 'cpu'); -rec_num = rec2Dastra(F_num_astra, (angles*pi/180), P, N, 'cpu'); - -figure; -subplot(1,2,1); imagesc(rec_an, [0 1]); daspect([1 1 1]); colormap hot; title('Analytical Sinogram Reconstruction [ASTRA]'); -subplot(1,2,2); imagesc(rec_num, [0 1]); daspect([1 1 1]); colormap hot; title('Numerical Sinogram Reconstruction [ASTRA]'); -%% diff --git a/matlab/Model2D_tGeneratorDemo.m b/matlab/Model2D_tGeneratorDemo.m deleted file mode 100644 index 66cd559..0000000 --- a/matlab/Model2D_tGeneratorDemo.m +++ /dev/null @@ -1,96 +0,0 @@ -% GPLv3 license (ASTRA toolbox) -% Note that the TomoPhantom package is released under Apache License, Version 2.0 - -% Script to generate 2D analytical temporal (2D + time) phantoms and sinograms (parallel beam) -% that can be used to test reconstruction algorithms without the "Inverse -% Crime" -% If one needs to modify/add phantoms, please edit Phantom2DLibrary.dat -% >>>> Prerequisites: ASTRA toolbox, if one needs to do reconstruction <<<<< - -close all;clc;clear; -% adding paths -fsep = '/'; -pathtoModels = sprintf(['..' fsep 'functions' fsep 'models' fsep], 1i); -addpath(pathtoModels); -addpath('compiled'); addpath('supplem'); - -ModelNo = 100; % Select a model from Phantom2DLibrary.dat -% Define phantom dimensions -N = 512; % x-y size (squared image) - -% Generate 2D+t phantom: -curDir = pwd; -mainDir = fileparts(curDir); -pathtoLibrary = sprintf([fsep 'functions' fsep 'models' fsep 'Phantom2DLibrary.dat'], 1i); -pathTP = strcat(mainDir, pathtoLibrary); % path to TomoPhantom parameters file -[G] = TomoP2DModel(ModelNo,N,pathTP); -figure(1); imagesc(G, [0 1]); daspect([1 1 1]); title('2D+t model, t=3 here'); colormap hot; -%% -% Lets look at more finely discretized temporal model and related sinograms -ModelNo = 101; % Select a model from Phantom2DLibrary.dat -% Define phantom dimensions -N = 512; % x-y size (squared image) - -% Generate 2D phantom: -curDir = pwd; -mainDir = fileparts(curDir); -pathtoLibrary = sprintf([fsep 'functions' fsep 'models' fsep 'Phantom2DLibrary.dat'], 1i); -pathTP = strcat(mainDir, pathtoLibrary); % path to TomoPhantom parameters file -[G] = TomoP2DModel(ModelNo,N,pathTP); - -angles = linspace(0,180,N); % projection angles -P = round(sqrt(2)*N); -F = TomoP2DModelSino(ModelNo, N, P, single(angles), pathTP, 'radon'); - -figure(2); -for i = 1:350 - subplot(1,2,1); imagesc(G(:,:,i), [0 1]); daspect([1 1 1]); title('2D+t phantom'); colormap hot; - subplot(1,2,2); imagesc(F(:,:,i), [0 165]); daspect([1 1 1]); title('Corresponding sinogram'); colormap hot; - pause(0.01); -end -%% -% another temporal phantom with stationary and dynamic features -ModelNo = 102; % Select a model from Phantom2DLibrary.dat -% Define phantom dimensions -N = 512; % x-y size (squared image) - -% Generate 2D+t phantom: -timeFrames = 25; %must be the same as in model -curDir = pwd; -mainDir = fileparts(curDir); -pathtoLibrary = sprintf([fsep 'functions' fsep 'models' fsep 'Phantom2DLibrary.dat'], 1i); -pathTP = strcat(mainDir, pathtoLibrary); % path to TomoPhantom parameters file -[G] = TomoP2DModel(ModelNo,N,pathTP); - -angles = linspace(0,180,N); % projection angles -P = round(sqrt(2)*N); -F = TomoP2DModelSino(ModelNo, N, P, single(angles), pathTP, 'radon'); - -% reconstruct -FBP_F_a = zeros(N,N,timeFrames); -for i = 1:timeFrames - FBP_F_a(:,:,i) = iradon(F(:,:,i)',angles,N); -end - -fig_num = 2; -for i = 1:timeFrames - figure(fig_num); - subplot(1,3,1); imagesc(G(:,:,i), [0 1]); daspect([1 1 1]); title('2D + time phantom'); colormap hot; - subplot(1,3,2); imagesc(F(:,:,i), [0 175]); daspect([1 1 1]); title('Corresponding sinogram'); colormap hot; - subplot(1,3,3); imagesc(FBP_F_a(:,:,i), [0 1]); daspect([1 1 1]); title('FBP reconstruction'); colormap hot; - pause(0.15); - - % if one needs animated gif here: -% filename = 'animat.gif'; -% del = 0.1; % time between animation frames -% drawnow -% frame = getframe(fig_num); -% im = frame2im(frame); -% [imind,cm] = rgb2ind(im,256); -% if (i == 1) -% imwrite(imind,cm,filename,'gif','Loopcount',inf,'DelayTime',del); -% else -% imwrite(imind,cm,filename,'gif','WriteMode','append','DelayTime',del); -% end -end -%% diff --git a/matlab/Model3DGeneratorDemo.m b/matlab/Model3DGeneratorDemo.m deleted file mode 100644 index 49433a3..0000000 --- a/matlab/Model3DGeneratorDemo.m +++ /dev/null @@ -1,87 +0,0 @@ -% GPLv3 license (ASTRA toolbox) -% Note that the TomoPhantom package is released under Apache License, Version 2.0 - -% Script to generate 3D analytical phantoms -% If one needs to modify/add phantoms, please edit Phantom3DLibrary.dat -% >>>> Prerequisites: ASTRA toolbox, if one needs to do reconstruction <<<<< - -close all;clc;clear; -% adding paths -fsep = '/'; -pathtoModels = sprintf(['..' fsep 'functions' fsep 'models' fsep], 1i); -addpath(pathtoModels); -addpath('compiled'); addpath('supplem'); - -ModelNo = 2; % Select a model -% Define phantom dimensions using a scalar (cubic) or a tuple [N1, N2, N3] -N = 256; % x-y-z size (cubic phantom), one can also pass DIM = [N N N] instead - -% generate 3D phantom (modify your PATH bellow): -curDir = pwd; -mainDir = fileparts(curDir); -pathtoLibrary = sprintf([fsep 'functions' fsep 'models' fsep 'Phantom3DLibrary.dat'], 1i); -pathTP = strcat(mainDir, pathtoLibrary); % path to TomoPhantom parameters file -tic; [G] = TomoP3DModel(ModelNo,N,pathTP); toc; - -% check 3 projections -figure; -slice = round(0.5*N); -subplot(1,3,1); imagesc(G(:,:,slice), [0 1]); daspect([1 1 1]); colormap hot; title('Axial Slice'); -subplot(1,3,2); imagesc(squeeze(G(:,slice,:)), [0 1]); daspect([1 1 1]); colormap hot; title('Y-Slice'); -subplot(1,3,3); imagesc(squeeze(G(slice,:,:)), [0 1]); daspect([1 1 1]); colormap hot; title('X-Slice'); - -% visualise/save the whole 3D Phantom -% % figure(2); -% filename = strcat('ModelNo',num2str(ModelNo)); -% counter = 1; -% for i = 1:N -% imshow(G(:,:,i), [0 1]); -% pause(0.01); -% % -% % % % write png images -% IM = im2uint16(G(:,:,i)); -% setStrNo = num2str(counter); -% if (counter < 10) -% filename_save = strcat(filename,'_','000',setStrNo, '.png'); -% elseif ((counter >= 10) && (counter < 100)) -% filename_save = strcat(filename,'_','00',setStrNo,'.png'); -% elseif ((counter >= 100) && (counter < 1000)) -% filename_save = strcat(filename,'_','0',setStrNo, '.png'); -% else -% filename_save = strcat(filename,'_',setStrNo, '.png'); -% end -% imwrite(IM,filename_save,'png'); -% counter = counter + 1; -% end -% close (figure(2)); -%% -fprintf('%s \n', 'Calculating 3D parallel-beam sinogram of a phantom using ASTRA-toolbox...'); -angles = linspace(0,pi,N); % projection angles -det = round(sqrt(2)*N); -sino_astra3D = zeros(length(angles),det,N,'single'); - -tic; -for i = 1:N -sino_astra3D(:,:,i) = sino2Dastra(G(:,:,i), angles, det, N, 'gpu'); -end -toc; - -% calculate residiual norm (the error is expected since projection models not the same) -% err_diff = norm(sino_tomophan3D(:) - sino_astra3D(:))./norm(sino_astra3D(:)); -% fprintf('%s %.4f\n', 'NRMSE for sino residuals:', err_diff); -% figure; -% subplot(1,2,1); imagesc(sino_tomophan3D(:,:,slice)', [0 70]); colormap hot; colorbar; daspect([1 1 1]); title('Exact sinogram'); -% subplot(1,2,2); imagesc(sino_astra3D(:,:,slice)', [0 70]); colormap hot; colorbar; daspect([1 1 1]); title('Discrete sinogram'); -%% -fprintf('%s \n', 'Reconstruction using ASTRA-toolbox (FBP)...'); -FBP3D = zeros(N,N,N,'single'); - -% choose which sinogram to substitute -for i = 1:N -FBP3D(:,:,i) = rec2Dastra(sino_astra3D(:,:,i), angles, det, N); -end -figure(3); -subplot(1,3,1); imagesc(FBP3D(:,:,slice), [0 1]); daspect([1 1 1]); colormap hot; title('Axial Slice'); -subplot(1,3,2); imagesc(squeeze(FBP3D(:,slice,:)), [0 1]); daspect([1 1 1]); colormap hot; title('Y-Slice'); -subplot(1,3,3); imagesc(squeeze(FBP3D(slice,:,:)), [0 1]); daspect([1 1 1]); colormap hot; title('X-Slice'); -%% diff --git a/matlab/Model3D_tGeneratorDemo.m b/matlab/Model3D_tGeneratorDemo.m deleted file mode 100644 index 76ee6d7..0000000 --- a/matlab/Model3D_tGeneratorDemo.m +++ /dev/null @@ -1,55 +0,0 @@ -% GPLv3 license (ASTRA toolbox) -% Note that the TomoPhantom package is released under Apache License, Version 2.0 - -% Script to generate 4D analytical phantoms (3D + time) -% If one needs to modify/add phantoms, please edit Phantom3DLibrary.dat -% >>>> Prerequisites: ASTRA toolbox, if one needs to do reconstruction <<<<< - -close all;clc;clear; -% adding paths -fsep = '/'; -pathtoModels = sprintf(['..' fsep 'functions' fsep 'models' fsep], 1i); -addpath(pathtoModels); -addpath('compiled'); addpath('supplem'); - -ModelNo = 100; % Select a model -% Define phantom dimensions using a scalar (cubic) or a tuple [N1, N2, N3] -N = 256; % x-y-z size (cubic phantom), one can also pass DIM = [N N N] instead - -% generate 4D phantom (modify your PATH bellow): -curDir = pwd; -mainDir = fileparts(curDir); -pathtoLibrary = sprintf([fsep 'functions' fsep 'models' fsep 'Phantom3DLibrary.dat'], 1i); -pathTP = strcat(mainDir, pathtoLibrary); % path to TomoPhantom parameters file -[G] = TomoP3DModel(ModelNo,N,pathTP); - -sliceM = round(0.5*N); -figure(1); -for i = 1:5 - imagesc(G(:,:,sliceM, i), [0 1]); daspect([1 1 1]); title('3D+t phantom'); colormap hot; - pause(0.1); -end -%% -% another 3D + time model -ModelNo = 101; % Select a model -% Define phantom dimensions -N = 256; % x-y-z size (cubic image) - -% generate 4D phantom (modify your PATH bellow): -curDir = pwd; -mainDir = fileparts(curDir); -pathtoLibrary = sprintf([fsep 'functions' fsep 'models' fsep 'Phantom3DLibrary.dat'], 1i); -pathTP = strcat(mainDir, pathtoLibrary); % path to TomoPhantom parameters file -[G] = TomoP3DModel(ModelNo,N,pathTP); - -%% -sliceM = round(0.5*N); -figure(2); -ll = 10; -counter = 1; -for ll = 1:5 - imagesc(G(:,:,sliceM, ll), [0 1]); daspect([1 1 1]); title('3D+t phantom'); colormap hot; - pause(0.1); -end -%% - diff --git a/matlab/Object2DGeneratorDemo.m b/matlab/Object2DGeneratorDemo.m deleted file mode 100644 index f648432..0000000 --- a/matlab/Object2DGeneratorDemo.m +++ /dev/null @@ -1,62 +0,0 @@ -% GPLv3 license (ASTRA toolbox) -% Note that the TomoPhantom package is released under Apache License, Version 2.0 - -% Script to generate 2D analytical object and the corresponding sinogram -% without using Phantom2DLibrary.dat - -close all;clc;clear; -% adding paths -fsep = '/'; -pathtoModels = sprintf(['..' fsep 'functions' fsep 'models' fsep], 1i); -addpath(pathtoModels); -addpath('compiled'); addpath('supplem'); - -% Define object dimensions -N = 256; % x-y-z size (cubic image) - -% define parameters -paramsObject.Ob = 'gaussian'; -paramsObject.C0 = 1; -paramsObject.x0 = -0.25; -paramsObject.y0 = 0.1; -paramsObject.a = 0.2; -paramsObject.b = 0.35; -paramsObject.phi1 = 30; - -% generate 2D phantom [N x N ]: -[G1] = TomoP2DObject(paramsObject.Ob,paramsObject.C0, paramsObject.x0, paramsObject.y0, paramsObject.a, paramsObject.b, paramsObject.phi1, N); -figure; imagesc(G1, [0 1]); daspect([1 1 1]); colormap hot; - -% generate corresponding 2D sinogram -angles = single(linspace(0,180,N)); % projection angles -P = round(sqrt(2)*N); - -[F1] = TomoP2DObjectSino(paramsObject.Ob,paramsObject.C0, paramsObject.x0, paramsObject.y0, paramsObject.a, paramsObject.b, paramsObject.phi1, N, P, angles); -figure; imagesc(F1, [0 50]); colormap hot; - -% generate another object -paramsObject.Ob = 'rectangle'; -paramsObject.C0 = 0.75; -paramsObject.x0 = 0.25; -paramsObject.y0 = -0.1; -paramsObject.a = 0.2; -paramsObject.b = 0.4; -paramsObject.phi1 = -45; - -[G2] = TomoP2DObject(paramsObject.Ob,paramsObject.C0, paramsObject.x0, paramsObject.y0, paramsObject.a, paramsObject.b, paramsObject.phi1, N); -figure; imagesc(G2, [0 1]); daspect([1 1 1]); colormap hot; -% generate corresponding 2D sinogram -[F2] = TomoP2DObjectSino(paramsObject.Ob,paramsObject.C0, paramsObject.x0, paramsObject.y0, paramsObject.a, paramsObject.b, paramsObject.phi1, N, P, angles); -figure; imagesc(F2, [0 50]); colormap hot; - -% build a new model -G = G1 + G2; -F = F1 + F2; - -figure; -subplot(1,2,1); imagesc(G, [0 1]); daspect([1 1 1]); colormap hot; title('New model'); -subplot(1,2,2); imagesc(F, [0 70]); daspect([1 1 1]); colormap hot; title('Sinogram'); - -% reconstruct the model -REC = rec2Dastra(F, double(angles*pi/180), P, N, 'cpu'); -figure; imagesc(REC, [0 1]); daspect([1 1 1]); colormap hot; \ No newline at end of file diff --git a/matlab/Object3DGeneratorDemo.m b/matlab/Object3DGeneratorDemo.m deleted file mode 100644 index c940d86..0000000 --- a/matlab/Object3DGeneratorDemo.m +++ /dev/null @@ -1,39 +0,0 @@ -% GPLv3 license (ASTRA toolbox) -% Note that the TomoPhantom package is released under Apache License, Version 2.0 - -% Script to generate 3D analytical object without using Phantom3DLibrary.dat - - -close all;clc;clear; -% adding paths -fsep = '/'; -pathtoModels = sprintf(['..' fsep 'functions' fsep 'models' fsep], 1i); -addpath(pathtoModels); -addpath('compiled'); addpath('supplem'); - - -% Define phantom dimensions using a scalar (cubic) or a tuple [N1, N2, N3] -N = 256; % x-y-z size (cubic phantom), one can also pass DIM = [N N N] instead - -% define parameters -paramsObject.Ob = 'gaussian'; -paramsObject.C0 = 1; -paramsObject.x0 = -0.25; -paramsObject.y0 = 0.1; -paramsObject.z0 = 0.0; -paramsObject.a = 0.2; -paramsObject.b = 0.35; -paramsObject.c = 0.7; -paramsObject.phi1 = 30; -paramsObject.phi2 = 60; -paramsObject.phi3 = -25; - -% Generate 3D phantom -tic; [G] = TomoP3DObject(paramsObject.Ob,paramsObject.C0, paramsObject.x0, paramsObject.y0, paramsObject.z0, paramsObject.a, paramsObject.b, paramsObject.c, paramsObject.phi1, paramsObject.phi2, paramsObject.phi3, N); toc; - -% check 3 projections -figure; -slice = round(0.5*N); -subplot(1,3,1); imagesc(G(:,:,slice), [0 1]); daspect([1 1 1]); colormap hot; title('Axial Slice'); -subplot(1,3,2); imagesc(squeeze(G(:,slice,:)), [0 1]); daspect([1 1 1]); colormap hot; title('Y-Slice'); -subplot(1,3,3); imagesc(squeeze(G(slice,:,:)), [0 1]); daspect([1 1 1]); colormap hot; title('X-Slice'); \ No newline at end of file diff --git a/matlab/Reconstruction_Demo.m b/matlab/Reconstruction_Demo.m deleted file mode 100644 index 991b1c5..0000000 --- a/matlab/Reconstruction_Demo.m +++ /dev/null @@ -1,60 +0,0 @@ -% GPLv3 license (ASTRA toolbox) -% Note that the TomoPhantom package is released under Apache License, Version 2.0 - -% Script to generate 2D analytical phantoms and their sinograms with added noise and artifacts -% Sinograms then reconstructed using ASTRA TOOLBOX - -% If one needs to modify/add phantoms, please edit Phantom2DLibrary.dat -% >>>> Prerequisites: ASTRA toolbox, if one needs to do reconstruction <<<<< - -close all;clc;clear; -% adding paths -fsep = '/'; -pathtoModels = sprintf(['..' fsep 'functions' fsep 'models' fsep], 1i); -addpath(pathtoModels); -addpath('compiled'); addpath('supplem'); - -ModelNo = 4; % Select a model from Phantom2DLibrary.dat -% Define phantom dimensions -N = 512; % x-y size (squared image) - -% Generate 2D phantom: -curDir = pwd; -mainDir = fileparts(curDir); -pathtoLibrary = sprintf([fsep 'functions' fsep 'models' fsep 'Phantom2DLibrary.dat'], 1i); -pathTP = strcat(mainDir, pathtoLibrary); % path to TomoPhantom parameters file -[G] = TomoP2DModel(ModelNo,N,pathTP); -figure; imagesc(G, [0 1]); daspect([1 1 1]); colormap hot; -%% -fprintf('%s \n', 'Generating sinogram analytically and numerically using ASTRA-toolbox...'); -angles = linspace(0,180,N); % projection angles -P = round(sqrt(2)*N); %detectors dimension -% generate 2D analytical parallel beam sinogram (note the 'astra' opton) -[F_a] = TomoP2DModelSino(ModelNo, N, P, single(angles), pathTP, 'astra'); -[F_num_astra] = sino2Dastra(G, (angles*pi/180), P, N, 'cpu'); - -% calculate residiual norm (the error is expected since projection models not the same) -err_diff = norm(F_a(:) - F_num_astra(:))./norm(F_num_astra(:)); -fprintf('%s %.4f\n', 'NMSE for sino residuals:', err_diff); - -figure; -subplot(1,2,1); imshow(F_a, []); title('Analytical Sinogram'); -subplot(1,2,2); imshow(F_num_astra, []); title('Numerical Sinogram (ASTRA)'); -%% -fprintf('%s \n', 'Adding noise and artifacts to analytical sinogram...'); -dose = 1e4; % photon flux (controls noise level) -[sino_noise] = add_noise(F_a, dose, 'Poisson'); % adding Poisson noise -[sino_noise_zingers] = add_zingers(sino_noise, 0.5, 10); % adding zingers -[sino_noise_zingers_stripes] = add_stripes(sino_noise_zingers, 1, 1); % adding stripes - -figure; imshow(sino_noise_zingers_stripes, []); title('Analytical Sinogram degraded with noise and artifacts'); -%% -fprintf('%s \n', 'Reconstruction using ASTRA-toolbox (FBP)...'); - -rec_an = rec2Dastra(sino_noise_zingers_stripes, (angles*pi/180), P, N, 'cpu'); -rec_num = rec2Dastra(F_num_astra, (angles*pi/180), P, N, 'cpu'); - -figure; -subplot(1,2,1); imagesc(rec_an, [0 1]); daspect([1 1 1]); colormap hot; title('Degraded analytical sinogram reconstruction'); -subplot(1,2,2); imagesc(rec_num, [0 1]); daspect([1 1 1]); colormap hot; title('Numerical sinogram reconstruction'); -%% \ No newline at end of file diff --git a/matlab/SpectralPhantomDemo.m b/matlab/SpectralPhantomDemo.m deleted file mode 100644 index 6e9ee4d..0000000 --- a/matlab/SpectralPhantomDemo.m +++ /dev/null @@ -1,57 +0,0 @@ -% Note that the TomoPhantom package is released under Apache License, Version 2.0 - -% Script to generate material-specific phantom - -% This phantom has been used in paper: -% Joint image reconstruction method with correlative multi-channel prior -% for X-ray spectral computed tomography by D. Kazantsev et al. -% Inverse Problems 2018. - -close all;clc;clear; -% adding paths -fsep = '/'; -pathtoModels = sprintf(['..' fsep 'functions' fsep 'models' fsep], 1i); -addpath(pathtoModels); -addpath('compiled'); addpath('supplem'); - - -ModelNo = 11; -% Define phantom dimension -N = 512; % x-y size (squared image) - -% generate the 2D phantom: -curDir = pwd; -mainDir = fileparts(curDir); -pathtoLibrary = sprintf([fsep 'functions' fsep 'models' fsep 'Phantom2DLibrary.dat'], 1i); -pathTP = strcat(mainDir, pathtoLibrary); % path to TomoPhantom parameters file -[G] = TomoP2DModel(ModelNo,N,pathTP); - -% create 4 phantoms with dedicated materials -G(G >= 1.29) = 0.8; -G1 = zeros(N,N); -G2 = zeros(N,N); -G3 = zeros(N,N); -G4 = zeros(N,N); - -G1((G > 0.2) & (G < 0.4)) = 1; -G2((G > 0.7) & (G < 0.9)) = 2; -G3((G > 1.0) & (G < 1.15)) = 3; -G4(G > 1.17) = 4; -G1(G > 1) = 0; - -figure; -subplot(2,2,1) % add first plot in 2 x 2 grid -imagesc(G1, [0 1]); daspect([1 1 1]); colormap hot; -title('Material 1') - -subplot(2,2,2) % add second plot in 2 x 2 grid -imagesc(G2, [0 2]); daspect([1 1 1]); colormap hot; -title('Material 2') - -subplot(2,2,3) % add third plot in 2 x 2 grid -imagesc(G3, [0 3]); daspect([1 1 1]); colormap hot; -title('Material 3') - -subplot(2,2,4) % add fourth plot in 2 x 2 grid -imagesc(G4, [0 4]); daspect([1 1 1]); colormap hot; -title('Material 4') diff --git a/matlab/install/compile_mex_linux.m b/matlab/install/compile_mex_linux.m deleted file mode 100755 index fef7223..0000000 --- a/matlab/install/compile_mex_linux.m +++ /dev/null @@ -1,37 +0,0 @@ -% compile mex files once from Matlab running this script - -fsep = '/'; -UpPath = sprintf(['..' fsep], 1i); -cd(UpPath); - -mkdir compiled - -PathFunc = sprintf(['..' fsep 'Core' fsep], 1i); - -copyfile(PathFunc, 'mex_wrappers'); - -cd('mex_wrappers'); - -Pathmove = sprintf(['..' fsep 'compiled' fsep], 1i); - -fprintf('%s \n', 'Building functions...'); -mex TomoP2DModel.c TomoP2DModel_core.c utils.c CFLAGS="\$CFLAGS -fopenmp -Wall -std=c99" LDFLAGS="\$LDFLAGS -fopenmp" -movefile('TomoP2DModel.mex*',Pathmove); -mex TomoP2DObject.c TomoP2DModel_core.c utils.c CFLAGS="\$CFLAGS -fopenmp -Wall -std=c99" LDFLAGS="\$LDFLAGS -fopenmp" -movefile('TomoP2DObject.mex*',Pathmove); -mex TomoP2DModelSino.c TomoP2DModelSino_core.c utils.c CFLAGS="\$CFLAGS -fopenmp -Wall -std=c99" LDFLAGS="\$LDFLAGS -fopenmp" -movefile('TomoP2DModelSino.mex*',Pathmove); -mex TomoP2DObjectSino.c TomoP2DModelSino_core.c utils.c CFLAGS="\$CFLAGS -fopenmp -Wall -std=c99" LDFLAGS="\$LDFLAGS -fopenmp" -movefile('TomoP2DObjectSino.mex*',Pathmove); -mex TomoP3DModel.c TomoP3DModel_core.c utils.c CFLAGS="\$CFLAGS -fopenmp -Wall -std=c99" LDFLAGS="\$LDFLAGS -fopenmp" -movefile('TomoP3DModel.mex*',Pathmove); -mex TomoP3DObject.c TomoP3DModel_core.c utils.c CFLAGS="\$CFLAGS -fopenmp -Wall -std=c99" LDFLAGS="\$LDFLAGS -fopenmp" -movefile('TomoP3DObject.mex*',Pathmove); -mex TomoP2DSinoNum.c TomoP2DSinoNum_core.c utils.c CFLAGS="\$CFLAGS -fopenmp -Wall -std=c99" LDFLAGS="\$LDFLAGS -fopenmp" -movefile('TomoP2DSinoNum.mex*',Pathmove); - -delete TomoP2DModel_core* TomoP2DModelSino_core* TomoP3DModel_core* TomoP2DSinoNum_core* CCPiDefines.h utils*; - -fprintf('%s \n', 'All compiled!'); - -cd(UpPath); \ No newline at end of file diff --git a/matlab/install/compile_mex_windows.m b/matlab/install/compile_mex_windows.m deleted file mode 100644 index 13c1f59..0000000 --- a/matlab/install/compile_mex_windows.m +++ /dev/null @@ -1,65 +0,0 @@ -% compile mex files once from Matlab running this script - -% >>>>>>>>>>>>>>>>>>>>>>>>>>>>> -% I've been able to compile on Windows 7 with MinGW and Matlab 2016b, however, -% not sure if openmp is enabled after the compilation. - -% Here I present two ways how software can be compiled, if you have some -% other suggestions please contact me at dkazanc@hotmail.com -% >>>>>>>>>>>>>>>>>>>>>>>>>>>>> - -fsep = '/'; -UpPath = sprintf(['..' fsep], 1i); -cd(UpPath); - -mkdir compiled - -PathFunc = sprintf(['..' fsep 'functions' fsep], 1i); -cd(PathFunc); -Pathmove = sprintf(['..' fsep 'matlab' fsep 'compiled' fsep], 1i); - -% One way to compile using Matlab-installed MinGW: -fprintf('%s \n', 'Building functions...'); -mex TomoP2DModel.c TomoP2DModel_core.c utils.c COMPFLAGS="\$COMPFLAGS -fopenmp -Wall -std=c99" -movefile('TomoP2DModel.mex*',Pathmove); -mex TomoP2DObject.c TomoP2DModel_core.c utils.c COMPFLAGS="\$COMPFLAGS -fopenmp -Wall -std=c99" -movefile('TomoP2DObject.mex*',Pathmove); -mex TomoP2DModelSino.c TomoP2DModelSino_core.c utils.c COMPFLAGS="\$COMPFLAGS -fopenmp -Wall -std=c99" -movefile('TomoP2DModelSino.mex*',Pathmove); -mex TomoP2DObjectSino.c TomoP2DModelSino_core.c utils.c COMPFLAGS="\$COMPFLAGS -fopenmp -Wall -std=c99" -movefile('TomoP2DObjectSino.mex*',Pathmove); -mex TomoP3DModel.c TomoP3DModel_core.c utils.c COMPFLAGS="\$COMPFLAGS -fopenmp -Wall -std=c99" -movefile('TomoP3DModel.mex*',Pathmove); -mex TomoP3DObject.c TomoP3DModel_core.c utils.c COMPFLAGS="\$COMPFLAGS -fopenmp -Wall -std=c99" -movefile('TomoP3DObject.mex*',Pathmove); -mex TomoP2DSinoNum.c TomoP2DSinoNum_core.c utils.c COMPFLAGS="\$COMPFLAGS -fopenmp -Wall -std=c99" -movefile('TomoP2DSinoNum.mex*',Pathmove); - -%%% The second approach to compile using TDM-GCC which follows this -%%% discussion: -%%% https://uk.mathworks.com/matlabcentral/answers/279171-using-mingw-compiler-and-open-mp#comment_359122 -%%% 1. Install TDM-GCC independently from http://tdm-gcc.tdragon.net/ (I installed 5.1.0) -%%% Install openmp version: http://sourceforge.net/projects/tdm-gcc/files/TDM-GCC%205%20series/5.1.0-tdm64-1/gcc-5.1.0-tdm64-1-openmp.zip/download -%%% 2. Link til libgomp.a in that installation when compilling your mex file. - -%%% assuming you unzipped TDM GCC (OpenMp) in folder TDMGCC on C drive, uncomment -%%% bellow -% fprintf('%s \n', 'Building functions...'); -% mex C:\TDMGCC\lib\gcc\x86_64-w64-mingw32\5.1.0\libgomp.a CXXFLAGS="$CXXFLAGS -std=c++11 -fopenmp" TomoP2DModel.c TomoP2DModel_core.c utils.c -% movefile('TomoP2DModel.mex*',Pathmove); -% mex C:\TDMGCC\lib\gcc\x86_64-w64-mingw32\5.1.0\libgomp.a CXXFLAGS="$CXXFLAGS -std=c++11 -fopenmp" TomoP2DObject.c TomoP2DModel_core.c utils.c -% movefile('TomoP2DObject.mex*',Pathmove); -% mex C:\TDMGCC\lib\gcc\x86_64-w64-mingw32\5.1.0\libgomp.a CXXFLAGS="$CXXFLAGS -std=c++11 -fopenmp" TomoP2DModelSino.c TomoP2DModelSino_core.c utils.c -% movefile('TomoP2DModelSino.mex*',Pathmove); -% mex C:\TDMGCC\lib\gcc\x86_64-w64-mingw32\5.1.0\libgomp.a CXXFLAGS="$CXXFLAGS -std=c++11 -fopenmp" TomoP2DObjectSino.c TomoP2DModelSino_core.c utils.c -% movefile('TomoP2DObjectSino.mex*',Pathmove); -% mex C:\TDMGCC\lib\gcc\x86_64-w64-mingw32\5.1.0\libgomp.a CXXFLAGS="$CXXFLAGS -std=c++11 -fopenmp" TomoP3DModel.c TomoP3DModel_core.c utils.c -% movefile('TomoP3DModel.mex*',Pathmove); -% mex C:\TDMGCC\lib\gcc\x86_64-w64-mingw32\5.1.0\libgomp.a CXXFLAGS="$CXXFLAGS -std=c++11 -fopenmp" TomoP3DObject.c TomoP3DModel_core.c utils.c -% movefile('TomoP3DObject.mex*',Pathmove); -% mex C:\TDMGCC\lib\gcc\x86_64-w64-mingw32\5.1.0\libgomp.a CXXFLAGS="$CXXFLAGS -std=c++11 -fopenmp" TomoP2DSinoNum.c TomoP2DSinoNum_core.c utils.c -% movefile('TomoP2DSinoNum.mex*',Pathmove); -% fprintf('%s \n', 'All compiled!'); - -cd(UpPath); -cd matlab \ No newline at end of file diff --git a/matlab/mex_wrappers/CMakeLists.txt b/matlab/mex_wrappers/CMakeLists.txt deleted file mode 100644 index e071c25..0000000 --- a/matlab/mex_wrappers/CMakeLists.txt +++ /dev/null @@ -1,102 +0,0 @@ -# Copyright 2018 Edoardo Pasca -#cmake_minimum_required (VERSION 3.0) - -project(RGL_core) -#https://stackoverflow.com/questions/13298504/using-cmake-with-setup-py - -# The version number. - -#set (CIL_VERSION $ENV{CIL_VERSION} CACHE INTERNAL "Core Imaging Library version" FORCE) - -# conda orchestrated build -message("CIL_VERSION ${CIL_VERSION}") -#include (GenerateExportHeader) - - -find_package(OpenMP REQUIRED) -if (OPENMP_FOUND) - set (CMAKE_C_FLAGS "${CMAKE_C_FLAGS} ${OpenMP_C_FLAGS}") - set (CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${OpenMP_CXX_FLAGS}") - set (CMAKE_EXE_LINKER_FLAGS "${CMAKE_EXE_LINKER_FLAGS} ${OpenMP_EXE_LINKER_FLAGS} ${OpenMP_CXX_FLAGS}") - set (CMAKE_SHARED_LINKER_FLAGS "${CMAKE_EXE_LINKER_FLAGS} ${OpenMP_SHARED_LINKER_FLAGS} ${OpenMP_CXX_FLAGS}") - set (CMAKE_STATIC_LINKER_FLAGS "${CMAKE_EXE_LINKER_FLAGS} ${OpenMP_STATIC_LINKER_FLAGS} ${OpenMP_CXX_FLAGS}") - -endif() - -## Build the regularisers package as a library -message("Creating Regularisers as a shared library") - -message("CMAKE_CXX_FLAGS ${CMAKE_CXX_FLAGS}") -message("CMAKE_C_FLAGS ${CMAKE_C_FLAGS}") -message("CMAKE_EXE_LINKER_FLAGS ${CMAKE_EXE_LINKER_FLAGS}") -message("CMAKE_SHARED_LINKER_FLAGS ${CMAKE_SHARED_LINKER_FLAGS}") -message("CMAKE_STATIC_LINKER_FLAGS ${CMAKE_STATIC_LINKER_FLAGS}") - -set(CMAKE_BUILD_TYPE "Release") - -if(WIN32) - set (FLAGS "/DWIN32 /EHsc /DCCPi_EXPORTS") - set (CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${FLAGS}") - set (CMAKE_C_FLAGS "${CMAKE_CXX_FLAGS} ${FLAGS}") - set (CMAKE_EXE_LINKER_FLAGS "${CMAKE_EXE_LINKER_FLAGS} /NODEFAULTLIB:MSVCRT.lib") - - set (EXTRA_LIBRARIES) - - message("library lib: ${LIBRARY_LIB}") - -elseif(UNIX) - set (FLAGS "-O2 -funsigned-char -Wall -Wl,--no-undefined -DCCPiReconstructionIterative_EXPORTS ") - set (CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${FLAGS}") - set (CMAKE_C_FLAGS "${CMAKE_CXX_FLAGS} ${FLAGS}") - if (OPENMP_FOUND) - set (EXTRA_LIBRARIES - "gomp" - "m" - ) - else() - set (EXTRA_LIBRARIES - "m" - ) - endif() -endif() -message("CMAKE_CXX_FLAGS ${CMAKE_CXX_FLAGS}") - -## Build the regularisers package as a library -message("Adding regularisers as a shared library") - -#set(CMAKE_C_COMPILER /apps/pgi/linux86-64/17.4/bin/pgcc) -#set(CMAKE_C_FLAGS "-acc -Minfo -ta=tesla:cc20 -openmp") -#set(CMAKE_C_FLAGS "-acc -Minfo -ta=multicore -openmp -fPIC") -add_library(ctomophantom "SHARED" - ${CMAKE_CURRENT_SOURCE_DIR}/TomoP2DModel_core.c - ${CMAKE_CURRENT_SOURCE_DIR}/TomoP2DModelSino_core.c - ${CMAKE_CURRENT_SOURCE_DIR}/TomoP3DModel_core.c - ${CMAKE_CURRENT_SOURCE_DIR}/utils.c - ${CMAKE_CURRENT_SOURCE_DIR}/TomoP2DSinoNum_core.c - ) -#add_library(ctomophantom_static "STATIC" -# ${CMAKE_CURRENT_SOURCE_DIR}/TomoP2DModel_core.c -# ${CMAKE_CURRENT_SOURCE_DIR}/TomoP2DModelSino_core.c -# ${CMAKE_CURRENT_SOURCE_DIR}/TomoP3DModel_core.c -# ${CMAKE_CURRENT_SOURCE_DIR}/utils.c -# ${CMAKE_CURRENT_SOURCE_DIR}/TomoP2DSinoNum_core.c -# ) -target_link_libraries(ctomophantom ${EXTRA_LIBRARIES} ) -include_directories(ctomophantom PUBLIC - ${LIBRARY_INC}/include - ${CMAKE_CURRENT_SOURCE_DIR} - ) - -## Install - -set(INSTALL_BIN_DIR ${CMAKE_INSTALL_PREFIX}/bin) -set(INSTALL_LIB_DIR ${CMAKE_INSTALL_PREFIX}/lib) -set(INSTALL_INCLUDE_DIR ${CMAKE_INSTALL_PREFIX}/include) - -install(TARGETS ctomophantom - LIBRARY DESTINATION "${INSTALL_LIB_DIR}" COMPONENT lib - PUBLIC_HEADER DESTINATION "${INSTALL_INCLUDE_DIR}" COMPONENT dev - RUNTIME DESTINATION "${INSTALL_BIN_DIR}" COMPONENT bin - ARCHIVE DESTINATION "${INSTALL_LIB_DIR}" COMPONENT lib - CONFIGURATIONS ${CMAKE_BUILD_TYPE} - ) diff --git a/matlab/supplem/add_noise.m b/matlab/supplem/add_noise.m deleted file mode 100644 index 0969ac4..0000000 --- a/matlab/supplem/add_noise.m +++ /dev/null @@ -1,20 +0,0 @@ -function [noisy_sino,N_noise] = add_noise(b, sigma, noisetype) -% function which adds Gaussian or Poisson noise to data -% return data after negative log and also raw data - -if (strcmp(noisetype, 'Gaussian') == 1) - E = randn(size(b)); - noisy_sino = b + sigma*norm(b,'fro')*E/norm(E,'fro'); % adding normal noise to the sinogram - noisy_sino(noisy_sino<0) = 0; -elseif (strcmp(noisetype, 'Poisson') == 1) - maxVal = max(b(:)); - b = b./maxVal; - dataExp = sigma.*exp(-b); % noiseless raw data - N_noise = poissrnd(dataExp); % adding Poisson noise to the sinogram - noisy_sino = -log(N_noise./max(N_noise(:)))*maxVal; %log corrected data -> sinogram - noisy_sino(noisy_sino<0) = 0; -else - error('Please select Gaussian or Poisson noise'); -end - -return; diff --git a/matlab/supplem/add_stripes.m b/matlab/supplem/add_stripes.m deleted file mode 100644 index ff130bf..0000000 --- a/matlab/supplem/add_stripes.m +++ /dev/null @@ -1,30 +0,0 @@ -function [stripes_sino] = add_stripes(b, percentage, maxthickness) -% function to add stripes (constant offsets) to sinogram which results in rings in the reconstructed image -% - percentage defines the density of stripes -% - maxthickness defines maximal thickness of a stripe - -if ((percentage <= 0) || (percentage > 100)) -error('Percentage must be larger than zero but smaller than 100'); -end -if ((maxthickness < 0) || (maxthickness > 5)) -error('Maximum thickness must be in [0,5] range'); -end - -DetectorsDim = size(b,2); -range_detect = round((DetectorsDim)*(percentage/100)); -stripes_sino = b; -max_intensity = max(b(:)); - -for jj = 1:range_detect - randind = randi([1 DetectorsDim],1); % generate random index of a faulty detector - randthickness = randi([0 maxthickness],1); % generate random thickness - randintens = -1 + (0.5 + 1).*rand(1,1); % generate random multiplier - intensity = max_intensity*randintens; - for ll = -randthickness:randthickness - if (((randind+ll) > 0) && ((randind+ll) < DetectorsDim)) - stripes_sino(:,randind+ll) = stripes_sino(:,randind+ll) + intensity; - end - end -end -stripes_sino(stripes_sino < 0) = 0; -return; diff --git a/matlab/supplem/add_zingers.m b/matlab/supplem/add_zingers.m deleted file mode 100644 index e79c1cc..0000000 --- a/matlab/supplem/add_zingers.m +++ /dev/null @@ -1,31 +0,0 @@ -function [zingers_sino] = add_zingers(b, percentage, modulus) -% adding zingers (zero pixels or small 4 pixels clusters) to data -% - percentage - the amount of zingers to be added -% - modulus controls the amount of 4 pixel clusters to be added - -if ((percentage <= 0) || (percentage > 100)) -error('Percentage must be larger than zero but smaller than 100'); -end -if (modulus < 0) -error('Modulus integer must be positive'); -end - -length_sino = length(b(:)); -num_values = round((length_sino)*(percentage/100)); -zingers_sino = b; - -for jj = 1:num_values - randind = randi([1 length_sino],1); % generate random indeces - zingers_sino(randind) = 0; - if (mod(jj,modulus) == 0) - %make a cluster of 4 pixels - if ((randind > size(b,1)) && (randind < length_sino-size(b,1))) - zingers_sino(randind+1) = 0; - zingers_sino(randind-1) = 0; - zingers_sino(randind+size(b,1)) = 0; - zingers_sino(randind-size(b,1)) = 0; - end - end -end - -return; diff --git a/matlab/supplem/rec2Dastra.m b/matlab/supplem/rec2Dastra.m deleted file mode 100644 index 7cb7732..0000000 --- a/matlab/supplem/rec2Dastra.m +++ /dev/null @@ -1,40 +0,0 @@ -function [reconstructon] = rec2Dastra(sino, angles, detsize, N, device) - - -proj_geom = astra_create_proj_geom('parallel', 1, detsize, angles); -vol_geom = astra_create_vol_geom(N,N); - -if (strcmp(device, 'cpu') == 1) - cfg = astra_struct('FBP'); - proj_id = astra_create_projector('strip', proj_geom, vol_geom); - cfg.ProjectorId = proj_id; - fprintf('%s \n', 'CPU FBP...'); -elseif (strcmp(device, 'gpu') == 1) - cfg = astra_struct('FBP_CUDA'); - fprintf('%s \n', 'GPU FBP...'); -else - error('Device is not selected, please choose between cpu and gpu') -end - -rec_id = astra_mex_data2d('create', '-vol', vol_geom); -sinogram_id = astra_mex_data2d('create', '-sino', proj_geom, sino); - -cfg.FilterType = 'Ram-Lak'; -cfg.ReconstructionDataId = rec_id; -cfg.ProjectionDataId = sinogram_id; - -alg_id = astra_mex_algorithm('create', cfg); -astra_mex_algorithm('run', alg_id); - -% Get the result -reconstructon = astra_mex_data2d('get', rec_id); - -% % Clean up. Note that GPU memory is tied up in the algorithm object, -% % and main RAM in the data objects. -astra_mex_algorithm('delete', alg_id); -astra_mex_data2d('delete', sinogram_id); -astra_mex_data2d('delete', rec_id); - -if (strcmp(device, 'cpu') == 1) -astra_mex_data2d('delete', proj_id); -end \ No newline at end of file diff --git a/matlab/supplem/sino2Dastra.m b/matlab/supplem/sino2Dastra.m deleted file mode 100644 index da10d14..0000000 --- a/matlab/supplem/sino2Dastra.m +++ /dev/null @@ -1,16 +0,0 @@ -function [sino] = sino2Dastra(phantom, angles, detsize, N, device) - -proj_geom = astra_create_proj_geom('parallel', 1, detsize, angles); -vol_geom = astra_create_vol_geom(N,N); - -if (strcmp(device, 'cpu') == 1) - proj_id = astra_create_projector('linear', proj_geom, vol_geom); - [sinogram_id, sino] = astra_create_sino(phantom, proj_id); - astra_mex_data2d('delete', proj_id); -elseif (strcmp(device, 'gpu') == 1) - [sinogram_id, sino] = astra_create_sino_cuda(phantom, proj_geom, vol_geom); -else - error('Device is not selected, please choose between cpu and gpu') -end -astra_mex_data2d('delete', sinogram_id); - diff --git a/python/Demos/DemoModel.py b/python/Demos/DemoModel.py index e8c2cb6..32a7936 100644 --- a/python/Demos/DemoModel.py +++ b/python/Demos/DemoModel.py @@ -22,7 +22,7 @@ model = 1 # selecting a model N_size = 512 #specify a full path to the parameters file -pathTP = '../../functions/models/Phantom2DLibrary.dat' +pathTP = '../../PhantomLibrary/models/Phantom2DLibrary.dat' #objlist = modelfile2Dtolist(pathTP, model) # one can extract parameters #This will generate a N_size x N_size phantom (2D) phantom_2D = TomoP2D.Model(model, N_size, pathTP) @@ -56,8 +56,9 @@ #%% ################################################################### # get numerical sinogram +""" tic=timeit.default_timer() -sino_num = TomoP2D.SinoNum (phantom_2D, P, angles) +sino_num = TomoP2D.SinoNum(phantom_2D, P, angles) toc=timeit.default_timer() Run_time = toc - tic print("Numerical sinogram has been generated in {} seconds".format(Run_time)) @@ -67,6 +68,7 @@ plt.imshow(sino_num, cmap="BuPu") plt.colorbar(ticks=[0, 150, 250], orientation='vertical') plt.title('{}''{}'.format('Numerical sinogram of model no.',model)) +""" #%% ################################################################### # get numerical sinogram (ASTRA-toolbox) diff --git a/python/Demos/DemoModel2.py b/python/Demos/DemoModel2.py index e99ec16..9db7938 100644 --- a/python/Demos/DemoModel2.py +++ b/python/Demos/DemoModel2.py @@ -26,7 +26,7 @@ model = 11 N_size = 512 #specify a full path to the parameters file -pathTP = '../../functions/models/Phantom2DLibrary.dat' +pathTP = '../../PhantomLibrary/models/Phantom2DLibrary.dat' objlist = modelfile2Dtolist(pathTP, model) # extract parameters using Python #This will generate a N_size x N_size phantom (2D) phantom_2D = TomoP2D.Object(N_size, objlist) diff --git a/python/Demos/DemoModel3D.py b/python/Demos/DemoModel3D.py index c854fb3..07c3d4d 100644 --- a/python/Demos/DemoModel3D.py +++ b/python/Demos/DemoModel3D.py @@ -17,11 +17,11 @@ print ("Building 3D phantom using TomoPhantom software") tic=timeit.default_timer() -model = 2 +model = 10 # Define phantom dimensions using a scalar (cubic) or a tuple [N1, N2, N3] N_size = 256 # or as a tuple of a custom size (256,256,256) #specify a full path to the parameters file -pathTP3 = '../../functions/models/Phantom3DLibrary.dat' +pathTP3 = '../../PhantomLibrary/models/Phantom3DLibrary.dat' #This will generate a N_size x N_size x N_size phantom (3D) or non-cubic phantom phantom_tm = TomoP3D.Model(model, N_size, pathTP3) toc=timeit.default_timer() @@ -56,7 +56,7 @@ DIM = (256,256,256) # full dimension of required phantom (z, y, x) DIM_z = (94, 158) # selected vertical subset (a slab) of the phantom #specify a full path to the parameters file -pathTP3 = '../../functions/models/Phantom3DLibrary.dat' +pathTP3 = '../../PhantomLibrary/models/Phantom3DLibrary.dat' #This will generate a N1 x N2 x N_slab phantom (3D) phantom_tm = TomoP3D.ModelSub(model, DIM, DIM_z, pathTP3) #phantom_tm = TomoP3D.Model(model, DIM, pathTP3) diff --git a/python/Demos/DemoModel4D.py b/python/Demos/DemoModel4D.py index 5d3c996..23a9788 100644 --- a/python/Demos/DemoModel4D.py +++ b/python/Demos/DemoModel4D.py @@ -22,7 +22,7 @@ # Define phantom dimensions using a scalar (cubic) or a tuple [N1, N2, N3] N_size = 256 # or as a tuple of a custom size (256,256,256) #specify a full path to the parameters file -pathTP3 = '../../functions/models/Phantom3DLibrary.dat' +pathTP3 = '../../PhantomLibrary/models/Phantom3DLibrary.dat' #This will generate a Time x N_size x N_size x N_size phantom (4D) phantom_tm = TomoP3D.ModelTemporal(model, N_size, pathTP3) toc=timeit.default_timer() @@ -60,7 +60,7 @@ DIM = (256,256,256) # full dimension of required phantom (z, y, x) DIM_z = (94, 158) # selected vertical subset (a slab) of the phantom #specify a full path to the parameters file -pathTP3 = '../../functions/models/Phantom3DLibrary.dat' +pathTP3 = '../../PhantomLibrary/models/Phantom3DLibrary.dat' #This will generate a Time x N1 x N2 x N_slab phantom (4D) phantom_tm = TomoP3D.ModelTemporalSub(model, DIM, DIM_z, pathTP3) #phantom_tm = TomoP3D.Model(model, DIM, pathTP3) diff --git a/python/Demos/DemoModel_temporal.py b/python/Demos/DemoModel_temporal.py index 5119955..20f32b4 100644 --- a/python/Demos/DemoModel_temporal.py +++ b/python/Demos/DemoModel_temporal.py @@ -20,7 +20,7 @@ N_size = 512 timeframes = 25 #specify a full path to the parameters file -pathTP = '../../functions/models/Phantom2DLibrary.dat' +pathTP = '../../../PhantomLibrary/models/Phantom2DLibrary.dat' #This will generate a N_size x N_size x Time frames phantom (2D + time) phantom_2Dt = TomoP2D.ModelTemporal(model, N_size, pathTP) diff --git a/python/Demos/DemoReconASTRA.py b/python/Demos/DemoReconASTRA.py index 277edbc..c8fc948 100644 --- a/python/Demos/DemoReconASTRA.py +++ b/python/Demos/DemoReconASTRA.py @@ -22,7 +22,7 @@ model = 4 # select a model N_size = 512 #specify a full path to the parameters file -pathTP = '../../functions/models/Phantom2DLibrary.dat' +pathTP = '../../PhantomLibrary/models/Phantom2DLibrary.dat' #objlist = modelfile2Dtolist(pathTP, model) # one can extract parameters #This will generate a N_size x N_size phantom (2D) phantom_2D = TomoP2D.Model(model, N_size, pathTP) diff --git a/python/Demos/DemoReconTomoPy.py b/python/Demos/DemoReconTomoPy.py index 0a73aa9..4446e47 100644 --- a/python/Demos/DemoReconTomoPy.py +++ b/python/Demos/DemoReconTomoPy.py @@ -21,7 +21,7 @@ model = 4 # select a model N_size = 512 #specify a full path to the parameters file -pathTP = '../../functions/models/Phantom2DLibrary.dat' +pathTP = '../../PhantomLibrary/models/Phantom2DLibrary.dat' #objlist = modelfile2Dtolist(pathTP, model) # one can extract parameters #This will generate a N_size x N_size phantom (2D) phantom_2D = TomoP2D.Model(model, N_size, pathTP) diff --git a/run.sh b/run.sh new file mode 100644 index 0000000..1105a16 --- /dev/null +++ b/run.sh @@ -0,0 +1,12 @@ +#!/bin/bash +echo "Building Tomophantom using CMake" +#rm -r build +# pip install cython +mkdir build +cd build/ +#make clean +cmake ../ -DBUILD_PYTHON_WRAPPER=ON -DCMAKE_BUILD_TYPE=Release -DCMAKE_INSTALL_PREFIX=./install +make install +cd install/python/ +export LD_LIBRARY_PATH=${LD_LIBRARY_PATH}:. +#spyder