From 56b9f2d3b655fa3f50058cae49821070193a55a4 Mon Sep 17 00:00:00 2001 From: dkazanc Date: Thu, 4 Oct 2018 13:02:14 +0100 Subject: [PATCH 1/5] reoarganisation and matlab mexing fixed, demos updated --- Wrappers/MATLAB/Model2D_tGeneratorDemo.m | 10 +- Wrappers/MATLAB/Model3DGeneratorDemo.m | 4 +- Wrappers/MATLAB/Model3D_tGeneratorDemo.m | 7 +- Wrappers/MATLAB/Reconstruction_Demo.m | 2 +- Wrappers/MATLAB/SpectralPhantomDemo.m | 4 +- Wrappers/MATLAB/install/compile_mex_linux.m | 15 +- Wrappers/MATLAB/install/compile_mex_windows.m | 22 +- Wrappers/MATLAB/mex_wrappers/TomoP2DModel.c | 304 +++++++++++++ .../MATLAB/mex_wrappers/TomoP2DModelSino.c | 383 ++++++++++++++++ Wrappers/MATLAB/mex_wrappers/TomoP2DObject.c | 96 ++++ .../MATLAB/mex_wrappers/TomoP2DObjectSino.c | 101 +++++ Wrappers/MATLAB/mex_wrappers/TomoP2DSinoNum.c | 70 +++ Wrappers/MATLAB/mex_wrappers/TomoP3DModel.c | 419 ++++++++++++++++++ Wrappers/MATLAB/mex_wrappers/TomoP3DObject.c | 124 ++++++ Wrappers/Python/Demos/DemoModel.py | 2 +- Wrappers/Python/Demos/DemoModel2.py | 12 +- Wrappers/Python/Demos/DemoModel3D.py | 81 ++++ Wrappers/Python/Demos/DemoModel4D.py | 88 ++++ Wrappers/Python/Demos/DemoModel_temporal.py | 2 +- Wrappers/Python/Demos/DemoReconASTRA.py | 2 +- Wrappers/Python/Demos/DemoReconTomoPy.py | 97 ++++ Wrappers/Python/Demos/libraryToDict.py | 2 +- Wrappers/Python/tomophantom/supp/artifacts.py | 16 +- Wrappers/Python/tomophantom/supp/recMod.py | 103 +++++ python/Demos/DemoModel.py | 2 +- python/Demos/DemoModel2.py | 2 +- python/Demos/DemoModel3D.py | 2 +- python/Demos/DemoModel4D.py | 2 +- python/Demos/DemoModel_temporal.py | 2 +- python/Demos/DemoReconASTRA.py | 2 +- python/Demos/DemoReconTomoPy.py | 2 +- 31 files changed, 1929 insertions(+), 51 deletions(-) mode change 100644 => 100755 Wrappers/MATLAB/install/compile_mex_linux.m create mode 100644 Wrappers/MATLAB/mex_wrappers/TomoP2DModel.c create mode 100644 Wrappers/MATLAB/mex_wrappers/TomoP2DModelSino.c create mode 100644 Wrappers/MATLAB/mex_wrappers/TomoP2DObject.c create mode 100644 Wrappers/MATLAB/mex_wrappers/TomoP2DObjectSino.c create mode 100644 Wrappers/MATLAB/mex_wrappers/TomoP2DSinoNum.c create mode 100644 Wrappers/MATLAB/mex_wrappers/TomoP3DModel.c create mode 100644 Wrappers/MATLAB/mex_wrappers/TomoP3DObject.c create mode 100644 Wrappers/Python/Demos/DemoModel3D.py create mode 100644 Wrappers/Python/Demos/DemoModel4D.py create mode 100644 Wrappers/Python/Demos/DemoReconTomoPy.py create mode 100644 Wrappers/Python/tomophantom/supp/recMod.py 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..506cc72 100644 --- a/Wrappers/MATLAB/Reconstruction_Demo.m +++ b/Wrappers/MATLAB/Reconstruction_Demo.m @@ -21,7 +21,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/Wrappers/MATLAB/mex_wrappers/TomoP2DModel.c b/Wrappers/MATLAB/mex_wrappers/TomoP2DModel.c new file mode 100644 index 0000000..a91a447 --- /dev/null +++ b/Wrappers/MATLAB/mex_wrappers/TomoP2DModel.c @@ -0,0 +1,304 @@ +/* + * Copyright 2017 Daniil Kazantsev + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * http://www.apache.org/licenses/LICENSE-2.0 + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + +#include +#include "mex.h" +#include +#include +#include +#include +#include "omp.h" + +#include "TomoP2DModel_core.h" +#include "utils.h" + +#define M_PI 3.14159265358979323846 +#define MAXCHAR 1000 + +/* MATLAB mex-wrapper for TomoP2Dmodel + * + * Input Parameters: + * 1. ModelNo - a model number from Phantom2DLibrary file + * 2. VolumeSize in voxels (N x N) or (N x N x time-frames) + * 3. An absolute path to the file Phantom2DLibrary.dat (see OS-specific syntax-differences) + * + * Output: + * 1. The analytical phantom size of [N x N] or a temporal phantom size of [N xN x Time-frames] + */ + +void mexFunction( + int nlhs, mxArray *plhs[], + int nrhs, const mxArray *prhs[]) + +{ + int ModelSelected, N; + float *A; + + + ModelSelected = (int) mxGetScalar(prhs[0]); /* selected model */ + N = (int) mxGetScalar(prhs[1]); /* choosen dimension (N x N x N) */ + + int Model=0, Components=0, steps = 0, counter=0, ii; + float C0 = 0.0f, x0 = 0.0f, y0 = 0.0f, a = 0.0f, b = 0.0f, psi_gr1; + + char *filename; + FILE * fp; + if(!mxIsChar(prhs[2]) ) { + mexErrMsgTxt("Need filename absolute path string input."); + } + filename = mxArrayToString(prhs[2]); + fp= fopen(filename, "rb"); + mxFree(filename); + if( fp == NULL ) { + mexErrMsgTxt("Cannot open the model library file (Phantom3DLibrary.dat)"); + } + else { + + char str[MAXCHAR]; + char tmpstr1[16]; + char tmpstr2[22]; + char tmpstr3[16]; + char tmpstr4[16]; + char tmpstr5[16]; + char tmpstr6[16]; + char tmpstr7[16]; + char tmpstr8[16]; + + while (fgets(str, MAXCHAR, fp) != NULL) + { + /* work with non-# commented lines */ + if(str[0] != '#') { + sscanf(str, "%15s : %21[^;];", tmpstr1, tmpstr2); + if (strcmp(tmpstr1,"Model")==0) + { + Model = atoi(tmpstr2); + if ((ModelSelected == Model) && (counter == 0)) { + /* check if we have a right model */ + if (fgets(str, MAXCHAR, fp) != NULL) sscanf(str, "%15s : %21[^;];", tmpstr1, tmpstr2); + else { + mexErrMsgTxt("Unexpected the end of the line (Components) in parameters file"); + break; } + if (strcmp(tmpstr1,"Components") == 0) Components = atoi(tmpstr2); + printf("%s %i\n", "Components:", Components); + if (Components <= 0) { + printf("%s %i\n", "Components cannot be negative, the given value is", Components); + mexErrMsgTxt("Components cannot be negative"); + break; } + if (fgets(str, MAXCHAR, fp) != NULL) sscanf(str, "%15s : %21[^;];", tmpstr1, tmpstr2); + else { + mexErrMsgTxt("Unexpected the end of the line (TimeSteps) in parameters file"); + break; } + if (strcmp(tmpstr1,"TimeSteps") == 0) steps = atoi(tmpstr2); + if (steps <= 0) { + printf("%s %i\n", "TimeSteps cannot be negative, the given value is", steps); + mexErrMsgTxt("TimeSteps cannot be negative"); + break; } + printf("%s %i\n", "TimeSteps:", steps); + if (steps == 1) { + /**************************************************/ + printf("\n %s \n", "Stationary model is selected"); + + const mwSize N_dims[2] = {N, N}; /* image dimensions */ + A = (float*)mxGetPr(plhs[0] = mxCreateNumericArray(2, N_dims, mxSINGLE_CLASS, mxREAL)); /*output array*/ + + /* loop over all components */ + for(ii=0; ii 1)) { + printf("%s %f\n", "x0 (object position) must be in [-1,1] range, the given value is", x0); + mexErrMsgTxt("x0 (object position) must be in [-1,1] range"); + break; } + if ((y0 < -1) || (y0 > 1)) { + printf("%s %f\n", "y0 (object position) must be in [-1,1] range, the given value is", y0); + mexErrMsgTxt("y0 (object position) must be in [-1,1] range"); + break; } + if ((a <= 0) || (a > 2)) { + printf("%s %f\n", "a (object size) must be positive in [0,2] range, the given value is", a); + mexErrMsgTxt("a (object size) must be positive in [0,2] range"); + break; } + if ((b <= 0) || (b > 2)) { + printf("%s %f\n", "b (object size) must be positive in [0,2] range, the given value is", b); + mexErrMsgTxt("b (object size) must be positive in [0,2] range"); + break; } + printf("\nObject : %s \nC0 : %f \nx0 : %f \ny0 : %f \na : %f \nb : %f \n", tmpstr2, C0, x0, y0, a, b); + + TomoP2DObject_core(A, N, tmpstr2, C0, x0, y0, b, a, -psi_gr1, 0); /* Matlab */ + } + } + else { + /**************************************************/ + printf("\n %s \n", "Temporal model is selected"); + /* temporal phantom 2D + time (3D) */ + const mwSize N_dims[3] = {N, N, steps}; /* image dimensions */ + A = (float*)mxGetPr(plhs[0] = mxCreateNumericArray(3, N_dims, mxSINGLE_CLASS, mxREAL)); + + float C1 = 0.0f, x1 = 0.0f, y1 = 0.0f, a1 = 0.0f, b1 = 0.0f, psi_gr1_1 = 0.0f; + /* loop over all components */ + for(ii=0; ii 1)) { + printf("%s %f\n", "x0 (object position) must be in [-1,1] range, the given value is", x0); + mexErrMsgTxt("x0 (object position) must be in [-1,1] range"); + break; } + if ((y0 < -1) || (y0 > 1)) { + printf("%s %f\n", "y0 (object position) must be in [-1,1] range, the given value is", y0); + mexErrMsgTxt("y0 (object position) must be in [-1,1] range"); + break; } + if ((a <= 0) || (a > 2)) { + printf("%s %f\n", "a (object size) must be positive in [0,2] range, the given value is", a); + mexErrMsgTxt("a (object size) must be positive in [0,2] range"); + break; } + if ((b <= 0) || (b > 2)) { + printf("%s %f\n", "b (object size) must be positive in [0,2] range, the given value is", b); + mexErrMsgTxt("b (object size) must be positive in [0,2] range"); + break; } + // printf("\nObject : %s \nC0 : %f \nx0 : %f \ny0 : %f \nz0 : %f \na : %f \nb : %f \n", tmpstr2, C0, x0, y0, z0, a, b, c); + + /* check Endvar relatedparameters */ + if (fgets(str, MAXCHAR, fp) != NULL) sscanf(str, "%15s : %15s %15s %15s %15s %15s %15[^;];", tmpstr1, tmpstr3, tmpstr4, tmpstr5, tmpstr6, tmpstr7, tmpstr8); + else { + mexErrMsgTxt("Unexpected the end of the line (Endvar loop) in parameters file"); + break; } + + if (strcmp(tmpstr1,"Endvar") == 0) { + C1 = (float)atof(tmpstr3); /* intensity */ + x1 = (float)atof(tmpstr4); /* x0 position */ + y1 = (float)atof(tmpstr5); /* y0 position */ + a1 = (float)atof(tmpstr6); /* a - size object */ + b1 = (float)atof(tmpstr7); /* b - size object */ + psi_gr1_1 = (float)atof(tmpstr8); /* rotation angle 1*/ + } + else { + printf("%s\n", "Cannot find 'Endvar' string in parameters file"); + break; } + + if (C1 == 0) { + printf("%s %f\n", "Endvar C1 should not be equal to zero, the given value is", C1); + mexErrMsgTxt("Endvar C0 should not be equal to zero"); + break; } + if ((x1 < -1) || (x1 > 1)) { + printf("%s %f\n", "Endvar x1 (object position) must be in [-1,1] range, the given value is", x1); + mexErrMsgTxt("Endvar x0 (object position) must be in [-1,1] range"); + break; } + if ((y1 < -1) || (y1 > 1)) { + printf("%s %f\n", "Endvar y1 (object position) must be in [-1,1] range, the given value is", y1); + mexErrMsgTxt("Endvar y0 (object position) must be in [-1,1] range"); + break; } + if ((a1 <= 0) || (a1 > 2)) { + printf("%s %f\n", "Endvar a1 (object size) must be positive in [0,2] range, the given value is", a1); + mexErrMsgTxt("Endvar a (object size) must be positive in [0,2] range"); + break; } + if ((b1 <= 0) || (b1 > 2)) { + printf("%s %f\n", "Endvar b1 (object size) must be positive in [0,2] range, the given value is", b1); + mexErrMsgTxt("Endvar b (object size) must be positive in [0,2] range"); + break; } + //printf("\nObject : %s \nC0 : %f \nx0 : %f \ny0 : %f \nz0 : %f \na : %f \nb : %f \nc : %f \n", tmpstr2, C0, x0, y0, z0, a1, b1, c1); + + /*now we know the initial parameters of the object and the final ones. We linearly extrapolate to establish steps and coordinates. */ + /* calculating the full distance berween the start and the end points */ + float distance = sqrtf(pow((x1 - x0),2) + pow((y1 - y0),2)); + float d_dist = distance/(steps-1); /*a step over line */ + float C_step = (C1 - C0)/(steps-1); + float a_step = (a1 - a)/(steps-1); + float b_step = (b1 - b)/(steps-1); + float phi_rot_step = (psi_gr1_1 - psi_gr1)/(steps-1); + + int tt; + float x_t, y_t, a_t, b_t, C_t, phi_t, d_step; + /* initialize */ + x_t = x0; y_t = y0; a_t = a; b_t = b; C_t = C0; phi_t = psi_gr1; d_step = d_dist; + /*loop over time frames*/ + for(tt=0; tt < steps; tt++) { + + TomoP2DObject_core(A, N, tmpstr2, C_t, x_t, y_t, b_t, a_t, -phi_t, tt); /* Matlab */ + + /* calculating new coordinates of an object */ + if (distance != 0.0f) { + float t = d_step/distance; + x_t = (1-t)*x0 + t*x1; + y_t = (1-t)*y0 + t*y1; } + else { + x_t = x0; + y_t = y0; } + + d_step += d_dist; + a_t += a_step; + b_t += b_step; + C_t += C_step; + phi_t += phi_rot_step; + } /*time steps*/ + + } /*components loop*/ + } + counter++; + } + } + } + } + } + fclose(fp); + if (counter == 0) { + printf("%s %i %s \n", "Model no. ", ModelSelected, "is not found!"); + mexErrMsgTxt("No object found, check models file"); + } +} diff --git a/Wrappers/MATLAB/mex_wrappers/TomoP2DModelSino.c b/Wrappers/MATLAB/mex_wrappers/TomoP2DModelSino.c new file mode 100644 index 0000000..7908a61 --- /dev/null +++ b/Wrappers/MATLAB/mex_wrappers/TomoP2DModelSino.c @@ -0,0 +1,383 @@ +/* + * Copyright 2017 Daniil Kazantsev + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * http://www.apache.org/licenses/LICENSE-2.0 + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + +#include "mex.h" +#include +#include +#include +#include +#include "omp.h" + +#include "TomoP2DModelSino_core.h" +#include "utils.h" + +#define M_PI 3.14159265358979323846 +#define MAXCHAR 1000 + +/* MATLAB mex-wrapper to generate parallel beam sinograms using TomoP2DModelSino_core.c + * + * Input Parameters: + * 1. Model number (see Phantom2DLibrary.dat) [required] + * 2. ImageSize in pixels (N x N) [required] + * 3. Detector array size P (in pixels) [required] + * 4. Projection angles Theta (in degrees) [required] + * 5. An absolute path to the file Phantom2DLibrary.dat (see OS-specific syntax-differences) [required] + * 6. ImageCentring, choose 'radon' or 'astra' (default) [optional] + * + * Output: + * 1. 2D sinogram size of [length(angles), P] or a temporal sinogram size of [length(angles), P, Time-Frames] + */ + +void mexFunction( + int nlhs, mxArray *plhs[], + int nrhs, const mxArray *prhs[]) + + { + int ModelSelected, N, CenTypeIn, P; + float *A, *Th; + mwSize NStructElems; + char *filename; + FILE * fp; + + /*Handling Matlab input data*/ + if ((nrhs < 5) && (nrhs > 6)) mexErrMsgTxt("Input of 5 or 6 parameters is required: Model, ImageSize, Detector Size, Projection angles, PATH, Centering"); + if (mxGetClassID(prhs[3]) != mxSINGLE_CLASS) {mexErrMsgTxt("The vector of angles must be in a single precision"); } + + ModelSelected = (int) mxGetScalar(prhs[0]); /* selected model */ + N = (int) mxGetScalar(prhs[1]); /* choosen dimension (N x N) */ + P = (int) mxGetScalar(prhs[2]); /* detector size */ + Th = (float*) mxGetData(prhs[3]); /* angles */ + CenTypeIn = 1; /* astra-type centering is the default one */ + + if (nrhs == 6) { + char *CenType; + CenType = mxArrayToString(prhs[5]); /* 'radon' or 'astra' (default) */ + if ((strcmp(CenType, "radon") != 0) && (strcmp(CenType, "astra") != 0)) mexErrMsgTxt("Choose 'radon' or 'astra''"); + if (strcmp(CenType, "radon") == 0) CenTypeIn = 0; /* enable 'radon'-type centaering */ + mxFree(CenType); + } + NStructElems = mxGetNumberOfElements(prhs[3]); + + int Model=0, Components=0, steps = 0, counter=0, ii; + float C0 = 0.0f, x0 = 0.0f, y0 = 0.0f, a = 0.0f, b = 0.0f, psi_gr1; + + if(!mxIsChar(prhs[4]) ) { + mexErrMsgTxt("Need filename absolute path string input."); + } + filename = mxArrayToString(prhs[4]); + fp= fopen(filename, "r"); + mxFree(filename); + if( fp == NULL ) { + mexErrMsgTxt("Cannot open the model library file (Phantom3DLibrary.dat)"); + } + else { + + char str[MAXCHAR]; + char tmpstr1[16]; + char tmpstr2[22]; + char tmpstr3[16]; + char tmpstr4[16]; + char tmpstr5[16]; + char tmpstr6[16]; + char tmpstr7[16]; + char tmpstr8[16]; + + while (fgets(str, MAXCHAR, fp) != NULL) + { + /* work with non-# commented lines */ + if(str[0] != '#') { + sscanf(str, "%15s : %21[^;];", tmpstr1, tmpstr2); + if (strcmp(tmpstr1,"Model")==0) + { + Model = atoi(tmpstr2); + if ((ModelSelected == Model) && (counter == 0)) { + /* check if we have a right model */ + if (fgets(str, MAXCHAR, fp) != NULL) sscanf(str, "%15s : %21[^;];", tmpstr1, tmpstr2); + else { + mexErrMsgTxt("Unexpected the end of the line (Components) in parameters file"); + break; } + if (strcmp(tmpstr1,"Components") == 0) Components = atoi(tmpstr2); + printf("%s %i\n", "Components:", Components); + if (Components <= 0) { + printf("%s %i\n", "Components cannot be negative, the given value is", Components); + mexErrMsgTxt("Components cannot be negative"); + break; } + if (fgets(str, MAXCHAR, fp) != NULL) sscanf(str, "%15s : %21[^;];", tmpstr1, tmpstr2); + else { + mexErrMsgTxt("Unexpected the end of the line (TimeSteps) in parameters file"); + break; } + if (strcmp(tmpstr1,"TimeSteps") == 0) steps = atoi(tmpstr2); + if (steps <= 0) { + printf("%s %i\n", "TimeSteps cannot be negative, the given value is", steps); + mexErrMsgTxt("TimeSteps cannot be negative"); + break; } + printf("%s %i\n", "TimeSteps:", steps); + if (steps == 1) { + /**************************************************/ + printf("\n %s %i %s \n", "Stationary 2D model", ModelSelected, " is selected"); + + const mwSize N_dims[2] = {NStructElems,P}; /*format: X-detectors, Y-angles dim*/ + A = (float*)mxGetPr(plhs[0] = mxCreateNumericArray(2, N_dims, mxSINGLE_CLASS, mxREAL)); + + /* loop over all components */ + for(ii=0; ii 1)) { + printf("%s %f\n", "x0 (object position) must be in [-1,1] range, the given value is", x0); + mexErrMsgTxt("x0 (object position) must be in [-1,1] range"); + break; } + if ((y0 < -1) || (y0 > 1)) { + printf("%s %f\n", "y0 (object position) must be in [-1,1] range, the given value is", y0); + mexErrMsgTxt("y0 (object position) must be in [-1,1] range"); + break; } + if ((a <= 0) || (a > 2)) { + printf("%s %f\n", "a (object size) must be positive in [0,2] range, the given value is", a); + mexErrMsgTxt("a (object size) must be positive in [0,2] range"); + break; } + if ((b <= 0) || (b > 2)) { + printf("%s %f\n", "b (object size) must be positive in [0,2] range, the given value is", b); + mexErrMsgTxt("b (object size) must be positive in [0,2] range"); + break; } + printf("\nObject : %s \nC0 : %f \nx0 : %f \ny0 : %f \na : %f \nb : %f \n", tmpstr2, C0, x0, y0, a, b); + + TomoP2DObjectSino_core(A, N, P, Th, (int)NStructElems, CenTypeIn, tmpstr2, C0, x0, y0, b, a, -psi_gr1, 0); /* Matlab */ + } + } + else { + /**************************************************/ + printf("\n %s %i %s \n", "Temporal 2D+time model", ModelSelected, " is selected"); + /* temporal phantom 2D + time (3D) */ + const mwSize N_dims[3] = {NStructElems, P, steps}; /*format: X-detectors, Y-angles dim, Time-Frames*/ + A = (float*)mxGetPr(plhs[0] = mxCreateNumericArray(3, N_dims, mxSINGLE_CLASS, mxREAL)); + + float C1 = 0.0f, x1 = 0.0f, y1 = 0.0f, a1 = 0.0f, b1 = 0.0f, psi_gr1_1 = 0.0f; + /* loop over all components */ + for(ii=0; ii 1)) { + printf("%s %f\n", "x0 (object position) must be in [-1,1] range, the given value is", x0); + mexErrMsgTxt("x0 (object position) must be in [-1,1] range"); + break; } + if ((y0 < -1) || (y0 > 1)) { + printf("%s %f\n", "y0 (object position) must be in [-1,1] range, the given value is", y0); + mexErrMsgTxt("y0 (object position) must be in [-1,1] range"); + break; } + if ((a <= 0) || (a > 2)) { + printf("%s %f\n", "a (object size) must be positive in [0,2] range, the given value is", a); + mexErrMsgTxt("a (object size) must be positive in [0,2] range"); + break; } + if ((b <= 0) || (b > 2)) { + printf("%s %f\n", "b (object size) must be positive in [0,2] range, the given value is", b); + mexErrMsgTxt("b (object size) must be positive in [0,2] range"); + break; } + // printf("\nObject : %s \nC0 : %f \nx0 : %f \ny0 : %f \nz0 : %f \na : %f \nb : %f \n", tmpstr2, C0, x0, y0, z0, a, b, c); + + /* check Endvar relatedparameters */ + if (fgets(str, MAXCHAR, fp) != NULL) sscanf(str, "%15s : %15s %15s %15s %15s %15s %15[^;];", tmpstr1, tmpstr3, tmpstr4, tmpstr5, tmpstr6, tmpstr7, tmpstr8); + else { + mexErrMsgTxt("Unexpected the end of the line (Endvar loop) in parameters file"); + break; } + + if (strcmp(tmpstr1,"Endvar") == 0) { + C1 = (float)atof(tmpstr3); /* intensity */ + x1 = (float)atof(tmpstr4); /* x0 position */ + y1 = (float)atof(tmpstr5); /* y0 position */ + a1 = (float)atof(tmpstr6); /* a - size object */ + b1 = (float)atof(tmpstr7); /* b - size object */ + psi_gr1_1 = (float)atof(tmpstr8); /* rotation angle 1*/ + } + else { + printf("%s\n", "Cannot find 'Endvar' string in parameters file"); + break; } + + if (C1 == 0) { + printf("%s %f\n", "Endvar C1 should not be equal to zero, the given value is", C1); + mexErrMsgTxt("Endvar C0 should not be equal to zero"); + break; } + if ((x1 < -1) || (x1 > 1)) { + printf("%s %f\n", "Endvar x1 (object position) must be in [-1,1] range, the given value is", x1); + mexErrMsgTxt("Endvar x0 (object position) must be in [-1,1] range"); + break; } + if ((y1 < -1) || (y1 > 1)) { + printf("%s %f\n", "Endvar y1 (object position) must be in [-1,1] range, the given value is", y1); + mexErrMsgTxt("Endvar y0 (object position) must be in [-1,1] range"); + break; } + if ((a1 <= 0) || (a1 > 2)) { + printf("%s %f\n", "Endvar a1 (object size) must be positive in [0,2] range, the given value is", a1); + mexErrMsgTxt("Endvar a (object size) must be positive in [0,2] range"); + break; } + if ((b1 <= 0) || (b1 > 2)) { + printf("%s %f\n", "Endvar b1 (object size) must be positive in [0,2] range, the given value is", b1); + mexErrMsgTxt("Endvar b (object size) must be positive in [0,2] range"); + break; } + //printf("\nObject : %s \nC0 : %f \nx0 : %f \ny0 : %f \nz0 : %f \na : %f \nb : %f \nc : %f \n", tmpstr2, C0, x0, y0, z0, a1, b1, c1); + + /*now we know the initial parameters of the object and the final ones. We linearly extrapolate to establish steps and coordinates. */ + /* calculating the full distance berween the start and the end points */ + float distance = sqrtf(pow((x1 - x0),2) + pow((y1 - y0),2)); + float d_dist = distance/(steps-1); /*a step over line */ + float C_step = (C1 - C0)/(steps-1); + float a_step = (a1 - a)/(steps-1); + float b_step = (b1 - b)/(steps-1); + float phi_rot_step = (psi_gr1_1 - psi_gr1)/(steps-1); + + int tt; + float x_t, y_t, a_t, b_t, C_t, phi_t, d_step; + /* initialize */ + x_t = x0; y_t = y0; a_t = a; b_t = b; C_t = C0; phi_t = psi_gr1; d_step = d_dist; + /*loop over time frames*/ + for(tt=0; tt < steps; tt++) { + + TomoP2DObjectSino_core(A, N, P, Th, (int)NStructElems, CenTypeIn, tmpstr2, C_t, x_t, y_t, b_t, a_t, -phi_t, tt); + + /* calculating new coordinates of an object */ + if (distance != 0.0f) { + float t = d_step/distance; + x_t = (1-t)*x0 + t*x1; + y_t = (1-t)*y0 + t*y1; } + else { + x_t = x0; + y_t = y0; } + + d_step += d_dist; + a_t += a_step; + b_t += b_step; + C_t += C_step; + phi_t += phi_rot_step; + } /*time steps*/ + + } /*components loop*/ + } + counter++; + } + } + } + } + } + fclose(fp); + if (counter == 0) { + printf("%s %i %s \n", "Model no. ", ModelSelected, "is not found!"); + mexErrMsgTxt("No object found, check models file"); + } +} + + + +// /* first to check if the model is stationary or temporal */ +// int *timeFrames = (int *) mxGetData(mxCreateNumericMatrix(1, 1, mxINT16_CLASS, mxREAL)); +// /* extract the number of time-step */ +// extractTimeFrames(timeFrames, ModelSelected, ModelParameters_PATH); +// +// int *params_switch = (int *) mxGetData(mxCreateNumericMatrix(1, 10, mxINT16_CLASS, mxREAL)); +// /* run checks for parameters (only integer switches and the number of time-frames) */ +// checkParams2D(params_switch, ModelSelected, ModelParameters_PATH, timeFrames[0]); +// +// if (params_switch[0] == 0) mexErrMsgTxt("Parameters file (Phantom2DLibrary.dat) cannot be read, check that your PATH to the file is a valid one and the file exists"); +// if (params_switch[1] == 0) { +// printf("%s %i\n", ">>>>> No such model available in the library, the given model is", ModelSelected); +// mexErrMsgTxt("Check the syntax in Phantom2DLibrary.dat"); +// } +// if (params_switch[2] == 0) { +// printf("%s %i\n", ">>>>> Components number cannot be negative for the given model", ModelSelected); +// mexErrMsgTxt("Check the syntax in Phantom2DLibrary.dat"); +// } +// if (params_switch[3] == 0) { +// printf("%s %i\n", ">>>>> TimeSteps value cannot be negative for the given model", ModelSelected); +// mexErrMsgTxt("Check the syntax in Phantom2DLibrary.dat"); +// } +// if (params_switch[4] == 0) { +// printf("%s %i\n", ">>>>> Unknown name of the object for the given model", ModelSelected); +// mexErrMsgTxt("Check the syntax in Phantom2DLibrary.dat"); +// } +// if (params_switch[5] == 0) { +// printf("%s %i\n", ">>>>> C0 should not be equal to zero for the given model", ModelSelected); +// mexErrMsgTxt("Check the syntax in Phantom2DLibrary.dat"); +// } +// if (params_switch[6] == 0) { +// printf("%s %i\n", ">>>>> x0 (object position) must be in [-1,1] range for the given model", ModelSelected); +// mexErrMsgTxt("Check the syntax in Phantom2DLibrary.dat"); +// } +// if (params_switch[7] == 0) { +// printf("%s %i\n", ">>>>> y0 (object position) must be in [-1,1] range for the given model", ModelSelected); +// mexErrMsgTxt("Check the syntax in Phantom2DLibrary.dat"); +// } +// if (params_switch[8] == 0) { +// printf("%s %i\n", ">>>>> a (object size) must be positive in [0,2] range", ModelSelected); +// mexErrMsgTxt("Check the syntax in Phantom2DLibrary.dat"); +// } +// if (params_switch[9] == 0) { +// printf("%s %i\n", ">>>>> b (object size) must be positive in [0,2] range", ModelSelected); +// mexErrMsgTxt("Check the syntax in Phantom2DLibrary.dat"); +// } +// +// /*Handling Matlab output data*/ +// if (params_switch[3] == 1) { +// /* static phantom case */ +// int N_dims[] = {NStructElems,P}; /*format: X-detectors, Y-angles dim*/ +// A = (float*)mxGetPr(plhs[0] = mxCreateNumericArray(2, N_dims, mxSINGLE_CLASS, mxREAL)); } +// else { +// int N_dims[] = {NStructElems, P, params_switch[3]}; /*format: X-detectors, Y-angles dim, Time-Frames*/ +// A = (float*)mxGetPr(plhs[0] = mxCreateNumericArray(3, N_dims, mxSINGLE_CLASS, mxREAL)); } +// TomoP2DModelSino_core(A, ModelSelected, N, P, Th, (int)NStructElems, CenTypeIn, ModelParameters_PATH, 0); +// mxFree(ModelParameters_PATH); +// } diff --git a/Wrappers/MATLAB/mex_wrappers/TomoP2DObject.c b/Wrappers/MATLAB/mex_wrappers/TomoP2DObject.c new file mode 100644 index 0000000..1276b6f --- /dev/null +++ b/Wrappers/MATLAB/mex_wrappers/TomoP2DObject.c @@ -0,0 +1,96 @@ +/* + * Copyright 2017 Daniil Kazantsev + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * http://www.apache.org/licenses/LICENSE-2.0 + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + +#include "mex.h" +#include +#include +#include +#include +#include "omp.h" + +#include "TomoP2DModel_core.h" +#include "utils.h" + +#define M_PI 3.14159265358979323846 + +/* Function to build 2D objects without using the library file + * (MATLAB mex-wrapper) + * + * Input Parameters: + * 1. parameters for 2D object + * 2. Size of the 2D object + * + * VolumeSize in voxels (N x N ) + * + * Output: + * 1. The analytical phantom size of [N x N] + * + */ +#define MAXCHAR 1000 + +void mexFunction( + int nlhs, mxArray *plhs[], + int nrhs, const mxArray *prhs[]) + +{ + int N; + float *A; + char *tmpstr2; + float C0 = 0.0f, x0 = 0.0f, y0 = 0.0f, a = 0.0f, b = 0.0f, psi_gr1 = 0.0f; + + tmpstr2 = mxArrayToString(prhs[0]); /* name of the object */ + C0 = (float) mxGetScalar(prhs[1]); /* intensity */ + x0 = (float) mxGetScalar(prhs[2]); /* x0 position */ + y0 = (float) mxGetScalar(prhs[3]); /* y0 position */ + a = (float) mxGetScalar(prhs[4]); /* a - size object */ + b = (float) mxGetScalar(prhs[5]); /* b - size object */ + psi_gr1 = (float) mxGetScalar(prhs[6]); /* rotation angle 1*/ + N = (int) mxGetScalar(prhs[7]); /* choosen dimension (N x N x N) */ + + /*Handling Matlab input data*/ + if (nrhs != 8) mexErrMsgTxt("Input of 8 parameters is required"); + + const mwSize N_dims[2] = {N, N}; /* image dimensions */ + A = (float*)mxGetPr(plhs[0] = mxCreateNumericArray(2, N_dims, mxSINGLE_CLASS, mxREAL)); /*output array*/ + + if ((strcmp("gaussian",tmpstr2) != 0) && (strcmp("parabola",tmpstr2) != 0) && (strcmp("ellipse",tmpstr2) != 0) && (strcmp("parabola1",tmpstr2) != 0) && (strcmp("cone",tmpstr2) != 0) && (strcmp("rectangle",tmpstr2) != 0) ) { + printf("%s %s\n", "Unknown name of the object, the given name is", tmpstr2); + mexErrMsgTxt("Unknown name of the object"); + } + if (C0 == 0) { + printf("%s %f\n", "C0 should not be equal to zero, the given value is", C0); + mexErrMsgTxt("C0 should not be equal to zero"); + } + if ((x0 < -1) || (x0 > 1)) { + printf("%s %f\n", "x0 (object position) must be in [-1,1] range, the given value is", x0); + mexErrMsgTxt("x0 (object position) must be in [-1,1] range"); + } + if ((y0 < -1) || (y0 > 1)) { + printf("%s %f\n", "y0 (object position) must be in [-1,1] range, the given value is", y0); + mexErrMsgTxt("y0 (object position) must be in [-1,1] range"); + } + if ((a <= 0) || (a > 2)) { + printf("%s %f\n", "a (object size) must be positive in [0,2] range, the given value is", a); + mexErrMsgTxt("a (object size) must be positive in [0,2] range"); + } + if ((b <= 0) || (b > 2)) { + printf("%s %f\n", "b (object size) must be positive in [0,2] range, the given value is", b); + mexErrMsgTxt("b (object size) must be positive in [0,2] range"); + } + printf("\nObject : %s \nC0 : %f \nx0 : %f \ny0 : %f \na : %f \nb : %f \n", tmpstr2, C0, x0, y0, a, b); + + TomoP2DObject_core(A, N, tmpstr2, C0, x0, y0, b, a, -psi_gr1, 0); /* Matlab */ + + mxFree(tmpstr2); +} diff --git a/Wrappers/MATLAB/mex_wrappers/TomoP2DObjectSino.c b/Wrappers/MATLAB/mex_wrappers/TomoP2DObjectSino.c new file mode 100644 index 0000000..d4b8252 --- /dev/null +++ b/Wrappers/MATLAB/mex_wrappers/TomoP2DObjectSino.c @@ -0,0 +1,101 @@ +/* + * Copyright 2017 Daniil Kazantsev + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * http://www.apache.org/licenses/LICENSE-2.0 + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + +#include "mex.h" +#include +#include +#include +#include +#include "omp.h" + +#include "TomoP2DModelSino_core.h" +#include "utils.h" + +/* Function to build 2D sinogram objects without using the library file + * (MATLAB mex-wrapper) + * + * Input Parameters: + * 1. parameters for 2D sinogram + * 3. + * 2. Size of the 2D object + * + * VolumeSize in voxels (N x N ) + * + * Output: + * 1. The analytical phantom size of [N x N] + * + */ +#define MAXCHAR 1000 + +void mexFunction( + int nlhs, mxArray *plhs[], + int nrhs, const mxArray *prhs[]) + +{ + mwSize NStructElems; + int N, P; + float *A, *Th; + char *tmpstr2; + float C0 = 0.0f, x0 = 0.0f, y0 = 0.0f, a = 0.0f, b = 0.0f, psi_gr1 = 0.0f; + + tmpstr2 = mxArrayToString(prhs[0]); /* name of the object */ + C0 = (float) mxGetScalar(prhs[1]); /* intensity */ + x0 = (float) mxGetScalar(prhs[2]); /* x0 position */ + y0 = (float) mxGetScalar(prhs[3]); /* y0 position */ + a = (float) mxGetScalar(prhs[4]); /* a - size object */ + b = (float) mxGetScalar(prhs[5]); /* b - size object */ + psi_gr1 = (float) mxGetScalar(prhs[6]); /* rotation angle 1*/ + N = (int) mxGetScalar(prhs[7]); /* choosen dimension (N x N) */ + P = (int) mxGetScalar(prhs[8]); /* detector size */ + Th = (float*) mxGetData(prhs[9]); /* angles */ + NStructElems = mxGetNumberOfElements(prhs[9]); + + int CenTypeIn = 1; /* astra-type centering is the default one */ + + /*Handling Matlab input data*/ + if (nrhs != 10) mexErrMsgTxt("Input of 10 parameters is required"); + + const mwSize N_dims[2] = {NStructElems,P}; /*format: X-detectors, Y-angles dim*/ + A = (float*)mxGetPr(plhs[0] = mxCreateNumericArray(2, N_dims, mxSINGLE_CLASS, mxREAL)); + + if ((strcmp("gaussian",tmpstr2) != 0) && (strcmp("parabola",tmpstr2) != 0) && (strcmp("ellipse",tmpstr2) != 0) && (strcmp("parabola1",tmpstr2) != 0) && (strcmp("cone",tmpstr2) != 0) && (strcmp("rectangle",tmpstr2) != 0) ) { + printf("%s %s\n", "Unknown name of the object, the given name is", tmpstr2); + mexErrMsgTxt("Unknown name of the object"); + } + if (C0 == 0) { + printf("%s %f\n", "C0 should not be equal to zero, the given value is", C0); + mexErrMsgTxt("C0 should not be equal to zero"); + } + if ((x0 < -1) || (x0 > 1)) { + printf("%s %f\n", "x0 (object position) must be in [-1,1] range, the given value is", x0); + mexErrMsgTxt("x0 (object position) must be in [-1,1] range"); + } + if ((y0 < -1) || (y0 > 1)) { + printf("%s %f\n", "y0 (object position) must be in [-1,1] range, the given value is", y0); + mexErrMsgTxt("y0 (object position) must be in [-1,1] range"); + } + if ((a <= 0) || (a > 2)) { + printf("%s %f\n", "a (object size) must be positive in [0,2] range, the given value is", a); + mexErrMsgTxt("a (object size) must be positive in [0,2] range"); + } + if ((b <= 0) || (b > 2)) { + printf("%s %f\n", "b (object size) must be positive in [0,2] range, the given value is", b); + mexErrMsgTxt("b (object size) must be positive in [0,2] range"); + } + printf("\nObject : %s \nC0 : %f \nx0 : %f \ny0 : %f \na : %f \nb : %f \n", tmpstr2, C0, x0, y0, a, b); + + TomoP2DObjectSino_core(A, N, P, Th, (int)NStructElems, CenTypeIn, tmpstr2, C0, x0, y0, b, a, -psi_gr1, 0); /* Matlab */ + + mxFree(tmpstr2); +} diff --git a/Wrappers/MATLAB/mex_wrappers/TomoP2DSinoNum.c b/Wrappers/MATLAB/mex_wrappers/TomoP2DSinoNum.c new file mode 100644 index 0000000..d170bdb --- /dev/null +++ b/Wrappers/MATLAB/mex_wrappers/TomoP2DSinoNum.c @@ -0,0 +1,70 @@ +/* + * Copyright 2017 Daniil Kazantsev + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * http://www.apache.org/licenses/LICENSE-2.0 + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + +#include "mex.h" +#include +#include +#include +#include +#include +#include +#include "omp.h" + +#include "TomoP2DSinoNum_core.h" +#include "utils.h" + + +/* C OMP implementation of the forward projection (The Radon Transform) + * by rotating a padded image and summing over columns (rotation-based projector + * for parallel beam) + * + * Input Parameters: + * 1. Phantom to calculate projections from [required] + * 2. Detector array size P (in pixels) [required] + * 3. Projection angles Theta (in degrees) [required] + * + * Output: + * Sinogram [No. angles] x [No. detectors] + * + * mex ForwardProjNum.c CFLAGS="\$CFLAGS -fopenmp -Wall -std=c99" LDFLAGS="\$LDFLAGS -fopenmp" + * sinogram = ForwardProjNum(single(G), P, single(angles)); + */ + +void mexFunction( + int nlhs, mxArray *plhs[], + int nrhs, const mxArray *prhs[]) { + int dimX, dimY,ThetaLength, DetSize, sys=1; + const int *dim_array, *dimsTh; + float *Phantom, *Sinogram, *Theta; + + /*Handling Matlab input data*/ + Phantom = (float *) mxGetData(prhs[0]); + DetSize = (int) mxGetScalar(prhs[1]); + Theta = (float *) mxGetData(prhs[2]); + + dim_array = mxGetDimensions(prhs[0]); + dimsTh = mxGetDimensions(prhs[2]); + + dimX = dim_array[0]; dimY = dim_array[1]; /*Image Dimensions (must be squared) */ + if( dimX != dimY) mexErrMsgTxt("The input image must be squared"); + if( DetSize < dimX) mexErrMsgTxt("The detector dimension is smaller than the phantom dimension! "); + + ThetaLength = dimsTh[1]; /* angles dimension */ + + /*Handling Matlab output data*/ + const mwSize N_dims[2] = {ThetaLength, DetSize}; /*format: Y-angles x X-detectors dim*/ + Sinogram = (float*)mxGetPr(plhs[0] = mxCreateNumericArray(2, N_dims, mxSINGLE_CLASS, mxREAL)); + + TomoP2DSinoNum_core(Sinogram, Phantom, dimX, DetSize, Theta, ThetaLength, sys); +} diff --git a/Wrappers/MATLAB/mex_wrappers/TomoP3DModel.c b/Wrappers/MATLAB/mex_wrappers/TomoP3DModel.c new file mode 100644 index 0000000..30393b0 --- /dev/null +++ b/Wrappers/MATLAB/mex_wrappers/TomoP3DModel.c @@ -0,0 +1,419 @@ +/* + * Copyright 2017 Daniil Kazantsev + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * http://www.apache.org/licenses/LICENSE-2.0 + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + +#include "mex.h" +#include +#include +#include +#include +#include "omp.h" + +#include "TomoP3DModel_core.h" +#include "utils.h" + +/* Function to read from the file Phantom3DLibrary.dat the required parameters to build 3D or 3D + time analytical models + * (MATLAB mex-wrapper) + * + * Input Parameters: + * 1. ModelNo - a model number from Phantom3DLibrary file + * 2. DIM volume dimensions [N1,N2,N3] in voxels (N1 x N2 x N3) or (N1 x N2 x N3 x time-frames) + * 3. An absolute path to the file Phantom3DLibrary.dat (ssee OS-specific syntax-differences) + * + * Output: + * 1. The analytical phantom size of [N1 x N2 x N3] or temporal 4D phantom (N1 x N2 x N3 x time-frames) + * + */ +#define MAXCHAR 1000 + +void mexFunction( + int nlhs, mxArray *plhs[], + int nrhs, const mxArray *prhs[]) + +{ + int ModelSelected; + float *A; + double *DimAr; + mwSize N1, N2, N3; + const mwSize *dim_array; + + /*Handling Matlab input data*/ + if (nrhs != 3) mexErrMsgTxt("Input of 3 parameters is required: model, DIMensions, PATH"); + + dim_array = mxGetDimensions(prhs[1]); /* get dimensions of DIM */ + + //printf("%i \n", dim_array[1]); + + N1 = 0; N2 = 0; N3 = 0; + if (dim_array[1] == 1) { + N1 = (long) mxGetScalar(prhs[1]); + N2 = N1; + N3 = N1; + } + else if (dim_array[1] == 3) { + DimAr = (double*) mxGetData(prhs[1]); + N1 = (long)(DimAr[0]); + N2 = (long)(DimAr[1]); + N3 = (long)(DimAr[2]); + } + else {mexErrMsgTxt("DIM must be scalar [N] or a vector with 3 elements [N1, N2, N3]");} + + ModelSelected = (int) mxGetScalar(prhs[0]); /* selected model from Phantom3DLibrary */ + + int Model=0, Components=0, steps = 0, counter=0, ii; + float C0 = 0.0f, x0 = 0.0f, y0 = 0.0f, z0 = 0.0f, a = 0.0f, b = 0.0f, c = 0.0f, psi_gr1 = 0.0f, psi_gr2 = 0.0f, psi_gr3 = 0.0f; + + char *filename; + FILE * fp; + if(!mxIsChar(prhs[2]) ) { + mexErrMsgTxt("Need filename absolute path string input."); + } + filename = mxArrayToString(prhs[2]); + fp= fopen(filename, "rb"); + mxFree(filename); + if( fp == NULL ) { + mexErrMsgTxt("Cannot open the model library file (Phantom3DLibrary.dat)"); + } + else { + + char str[MAXCHAR]; + char tmpstr1[16]; + char tmpstr2[22]; + char tmpstr3[16]; + char tmpstr4[16]; + char tmpstr5[16]; + char tmpstr6[16]; + char tmpstr7[16]; + char tmpstr8[16]; + char tmpstr9[16]; + char tmpstr10[16]; + char tmpstr11[16]; + char tmpstr12[16]; + + while (fgets(str, MAXCHAR, fp) != NULL) + { + /* work with non-# commented lines */ + if(str[0] != '#') { + sscanf(str, "%15s : %21[^;];", tmpstr1, tmpstr2); + if (strcmp(tmpstr1,"Model")==0) + { + Model = atoi(tmpstr2); + if ((ModelSelected == Model) && (counter == 0)) { + /* check if we have a right model */ + if (fgets(str, MAXCHAR, fp) != NULL) sscanf(str, "%15s : %21[^;];", tmpstr1, tmpstr2); + else { + mexErrMsgTxt("Unexpected the end of the line (Components) in parameters file"); + break; } + if (strcmp(tmpstr1,"Components") == 0) Components = atoi(tmpstr2); + printf("%s %i\n", "Components:", Components); + if (Components <= 0) { + printf("%s %i\n", "Components cannot be negative, the given value is", Components); + mexErrMsgTxt("Components cannot be negative"); + break; } + if (fgets(str, MAXCHAR, fp) != NULL) sscanf(str, "%15s : %21[^;];", tmpstr1, tmpstr2); + else { + mexErrMsgTxt("Unexpected the end of the line (TimeSteps) in parameters file"); + break; } + if (strcmp(tmpstr1,"TimeSteps") == 0) steps = atoi(tmpstr2); + if (steps <= 0) { + printf("%s %i\n", "TimeSteps cannot be negative, the given value is", steps); + mexErrMsgTxt("TimeSteps cannot be negative"); + break; } + printf("%s %i\n", "TimeSteps:", steps); + if (steps == 1) { + /**************************************************/ + printf("\n %s %i %s \n", "Stationary 3D model", ModelSelected, " is selected"); + + const mwSize N_dims[3] = {N1, N2, N3}; /* image dimensions */ + A = (float*)mxGetPr(plhs[0] = mxCreateNumericArray(3, N_dims, mxSINGLE_CLASS, mxREAL)); /*output array*/ + + /* loop over all components */ + for(ii=0; ii 1)) { + printf("%s %f\n", "x0 (object position) must be in [-1,1] range, the given value is", x0); + mexErrMsgTxt("x0 (object position) must be in [-1,1] range"); + break; } + if ((y0 < -1) || (y0 > 1)) { + printf("%s %f\n", "y0 (object position) must be in [-1,1] range, the given value is", y0); + mexErrMsgTxt("y0 (object position) must be in [-1,1] range"); + break; } + if ((z0 < -1) || (z0 > 1)) { + printf("%s %f\n", "z0 (object position) must be in [-1,1] range, the given value is", z0); + mexErrMsgTxt("z0 (object position) must be in [-1,1] range"); + break; } + if ((a <= 0) || (a > 2)) { + printf("%s %f\n", "a (object size) must be positive in [0,2] range, the given value is", a); + mexErrMsgTxt("a (object size) must be positive in [0,2] range"); + break; } + if ((b <= 0) || (b > 2)) { + printf("%s %f\n", "b (object size) must be positive in [0,2] range, the given value is", b); + mexErrMsgTxt("b (object size) must be positive in [0,2] range"); + break; } + if ((c <= 0) || (c > 2)) { + printf("%s %f\n", "c (object size) must be positive in [0,2] range, the given value is", c); + mexErrMsgTxt("c (object size) must be positive in [0,2] range"); + break; } + printf("\nObject : %s \nC0 : %f \nx0 : %f \ny0 : %f \nz0 : %f \na : %f \nb : %f \nc : %f \n", tmpstr2, C0, x0, y0, z0, a, b, c); + + TomoP3DObject_core(A, N1, N2, N3, 0l, N3, tmpstr2, C0, x0, y0, z0, b, a, c, -psi_gr1, psi_gr2, psi_gr3, 0l); /* Matlab */ + } + } + else { + /**************************************************/ + printf("\n %s %i %s \n", "Temporal 3D+time model", ModelSelected, " is selected"); + /* temporal phantom 3D + time (4D) */ + const mwSize N_dims[4] = {N1, N2, N3, steps}; /* image dimensions */ + A = (float*)mxGetPr(plhs[0] = mxCreateNumericArray(4, N_dims, mxSINGLE_CLASS, mxREAL)); + + float C1 = 0.0f, x1 = 0.0f, y1 = 0.0f, z1 = 0.0f, a1 = 0.0f, b1 = 0.0f, c1 = 0.0f, psi_gr1_1 = 0.0f, psi_gr2_1 = 0.0f, psi_gr3_1 = 0.0f; + /* loop over all components */ + for(ii=0; ii 1)) { + printf("%s %f\n", "x0 (object position) must be in [-1,1] range, the given value is", x0); + mexErrMsgTxt("x0 (object position) must be in [-1,1] range"); + break; } + if ((y0 < -1) || (y0 > 1)) { + printf("%s %f\n", "y0 (object position) must be in [-1,1] range, the given value is", y0); + mexErrMsgTxt("y0 (object position) must be in [-1,1] range"); + break; } + if ((z0 < -1) || (z0 > 1)) { + printf("%s %f\n", "z0 (object position) must be in [-1,1] range, the given value is", z0); + mexErrMsgTxt("z0 (object position) must be in [-1,1] range"); + break; } + if ((a <= 0) || (a > 2)) { + printf("%s %f\n", "a (object size) must be positive in [0,2] range, the given value is", a); + mexErrMsgTxt("a (object size) must be positive in [0,2] range"); + break; } + if ((b <= 0) || (b > 2)) { + printf("%s %f\n", "b (object size) must be positive in [0,2] range, the given value is", b); + mexErrMsgTxt("b (object size) must be positive in [0,2] range"); + break; } + if ((c <= 0) || (c > 2)) { + printf("%s %f\n", "c (object size) must be positive in [0,2] range, the given value is", c); + mexErrMsgTxt("c (object size) must be positive in [0,2] range"); + break; } + // printf("\nObject : %s \nC0 : %f \nx0 : %f \ny0 : %f \nz0 : %f \na : %f \nb : %f \n", tmpstr2, C0, x0, y0, z0, a, b, c); + + /* check Endvar relatedparameters */ + if (fgets(str, MAXCHAR, fp) != NULL) sscanf(str, "%15s : %15s %15s %15s %15s %15s %15s %15s %15s %15s %15[^;];", tmpstr1, tmpstr3, tmpstr4, tmpstr5, tmpstr6, tmpstr7, tmpstr8, tmpstr9, tmpstr10, tmpstr11, tmpstr12); + else { + mexErrMsgTxt("Unexpected the end of the line (Endvar loop) in parameters file"); + break; } + + if (strcmp(tmpstr1,"Endvar") == 0) { + C1 = (float)atof(tmpstr3); /* intensity */ + x1 = (float)atof(tmpstr4); /* x0 position */ + y1 = (float)atof(tmpstr5); /* y0 position */ + z1 = (float)atof(tmpstr6); /* z0 position */ + a1 = (float)atof(tmpstr7); /* a - size object */ + b1 = (float)atof(tmpstr8); /* b - size object */ + c1 = (float)atof(tmpstr9); /* c - size object */ + psi_gr1_1 = (float)atof(tmpstr10); /* rotation angle 1*/ + psi_gr2_1 = (float)atof(tmpstr11); /* rotation angle 2*/ + psi_gr3_1 = (float)atof(tmpstr12); /* rotation angle 3*/ + } + else { + printf("%s\n", "Cannot find 'Endvar' string in parameters file"); + break; } + + if (C1 == 0) { + printf("%s %f\n", "Endvar C1 should not be equal to zero, the given value is", C1); + mexErrMsgTxt("Endvar C0 should not be equal to zero"); + break; } + if ((x1 < -1) || (x1 > 1)) { + printf("%s %f\n", "Endvar x1 (object position) must be in [-1,1] range, the given value is", x1); + mexErrMsgTxt("Endvar x0 (object position) must be in [-1,1] range"); + break; } + if ((y1 < -1) || (y1 > 1)) { + printf("%s %f\n", "Endvar y1 (object position) must be in [-1,1] range, the given value is", y1); + mexErrMsgTxt("Endvar y0 (object position) must be in [-1,1] range"); + break; } + if ((z1 < -1) || (z1 > 1)) { + printf("%s %f\n", "Endvar z1 (object position) must be in [-1,1] range, the given value is", z1); + mexErrMsgTxt("Endvar z0 (object position) must be in [-1,1] range"); + break; } + if ((a1 <= 0) || (a1 > 2)) { + printf("%s %f\n", "Endvar a1 (object size) must be positive in [0,2] range, the given value is", a1); + mexErrMsgTxt("Endvar a (object size) must be positive in [0,2] range"); + break; } + if ((b1 <= 0) || (b1 > 2)) { + printf("%s %f\n", "Endvar b1 (object size) must be positive in [0,2] range, the given value is", b1); + mexErrMsgTxt("Endvar b (object size) must be positive in [0,2] range"); + break; } + if ((c1 <= 0) || (c1 > 2)) { + printf("%s %f\n", "Endvar c1 (object size) must be positive in [0,2] range, the given value is", c1); + mexErrMsgTxt("Endvar c (object size) must be positive in [0,2] range"); + break; } + //printf("\nObject : %s \nC0 : %f \nx0 : %f \ny0 : %f \nz0 : %f \na : %f \nb : %f \nc : %f \n", tmpstr2, C0, x0, y0, z0, a1, b1, c1); + + /*now we know the initial parameters of the object and the final ones. We linearly extrapolate to establish steps and coordinates. */ + /* calculating the full distance berween the start and the end points */ + float distance = sqrtf(pow((x1 - x0),2) + pow((y1 - y0),2) + pow((z1 - z0),2)); + float d_dist = distance/(steps-1); /*a step over line */ + float C_step = (C1 - C0)/(steps-1); + float a_step = (a1 - a)/(steps-1); + float b_step = (b1 - b)/(steps-1); + float c_step = (c1 - c)/(steps-1); + float phi_rot_step1 = (psi_gr1_1 - psi_gr1)/(steps-1); + float phi_rot_step2 = (psi_gr2_1 - psi_gr2)/(steps-1); + float phi_rot_step3 = (psi_gr3_1 - psi_gr3)/(steps-1); + + int tt; + float x_t, y_t, z_t, a_t, b_t, c_t, C_t, phi1_t, phi2_t, phi3_t, d_step; + /* initialize */ + x_t = x0; y_t = y0; z_t = z0; a_t = a; b_t = b; c_t = c; C_t = C0; phi1_t = psi_gr1; phi2_t = psi_gr2; phi3_t = psi_gr3; d_step = d_dist; + + /*loop over time frames*/ + for(tt=0; tt < steps; tt++) { + + TomoP3DObject_core(A, N1, N2, N3, 0l, N3, tmpstr2, C_t, x_t, y_t, z_t, b_t, a_t, c_t, -phi1_t, phi2_t, phi3_t, tt); /* Matlab */ + + /* calculating new coordinates of an object */ + if (distance != 0.0f) { + float t = d_step/distance; + x_t = (1-t)*x0 + t*x1; + y_t = (1-t)*y0 + t*y1; + z_t = (1-t)*z0 + t*z1;} + else { + x_t = x0; + y_t = y0; + z_t = z0; } + + d_step += d_dist; + a_t += a_step; + b_t += b_step; + c_t += c_step; + C_t += C_step; + phi1_t += phi_rot_step1; + phi2_t += phi_rot_step2; + phi3_t += phi_rot_step3; + } /*time steps*/ + + } /*components loop*/ + } + counter++; + } + } + } + } + } + fclose(fp); + if (counter == 0) { + printf("%s %i %s \n", "Model no. ", ModelSelected, "is not found!"); + mexErrMsgTxt("No object found, check models file"); + } +} +// if (params_switch[0] == 0) mexErrMsgTxt("Parameters file (Phantom3DLibrary.dat) cannot be read OR some major syntax errors in it; Check your PATH to the file is a valid one and the file exists"); +// if (params_switch[1] == 0) { +// printf("%s %i\n", ">>>>> No such model available in the library, the given model is", ModelSelected); +// mexErrMsgTxt("Check the syntax in Phantom3DLibrary.dat"); +// } +// if (params_switch[2] == 0) { +// printf("%s %i\n", ">>>>> Components number cannot be negative for the given model", ModelSelected); +// mexErrMsgTxt("Check the syntax in Phantom3DLibrary.dat"); +// } +// if (params_switch[3] == 0) { +// printf("%s %i\n", ">>>>> TimeSteps value cannot be negative for the given model", ModelSelected); +// mexErrMsgTxt("Check the syntax in Phantom3DLibrary.dat"); +// } +// if (params_switch[4] == 0) { +// printf("%s %i\n", ">>>>> Unknown name of the object for the given model", ModelSelected); +// mexErrMsgTxt("Check the syntax in Phantom3DLibrary.dat"); +// } +// if (params_switch[5] == 0) { +// printf("%s %i\n", ">>>>> C0 should not be equal to zero for the given model", ModelSelected); +// mexErrMsgTxt("Check the syntax in Phantom3DLibrary.dat"); +// } +// if (params_switch[6] == 0) { +// printf("%s %i\n", ">>>>> x0 (object position) must be in [-1,1] range for the given model", ModelSelected); +// mexErrMsgTxt("Check the syntax in Phantom3DLibrary.dat"); +// } +// if (params_switch[7] == 0) { +// printf("%s %i\n", ">>>>> y0 (object position) must be in [-1,1] range for the given model", ModelSelected); +// mexErrMsgTxt("Check the syntax in Phantom3DLibrary.dat"); +// } +// if (params_switch[8] == 0) { +// printf("%s %i\n", ">>>>> z0 (object position) must be in [-1,1] range for the given model", ModelSelected); +// mexErrMsgTxt("Check the syntax in Phantom3DLibrary.dat"); +// } +// if (params_switch[9] == 0) { +// printf("%s %i\n", ">>>>> a (object size) must be positive in [0,2] range for the given model", ModelSelected); +// mexErrMsgTxt("Check the syntax in Phantom3DLibrary.dat"); +// } +// if (params_switch[10] == 0) { +// printf("%s %i\n", ">>>>> b (object size) must be positive in [0,2] range for the given model", ModelSelected); +// mexErrMsgTxt("Check the syntax in Phantom3DLibrary.dat"); +// } +// if (params_switch[11] == 0) { +// printf("%s %i\n", ">>>>> c (object size) must be positive in [0,2] range for the given model", ModelSelected); +// mexErrMsgTxt("Check the syntax in Phantom3DLibrary.dat"); +// } diff --git a/Wrappers/MATLAB/mex_wrappers/TomoP3DObject.c b/Wrappers/MATLAB/mex_wrappers/TomoP3DObject.c new file mode 100644 index 0000000..9f0f9b0 --- /dev/null +++ b/Wrappers/MATLAB/mex_wrappers/TomoP3DObject.c @@ -0,0 +1,124 @@ +/* + * Copyright 2017 Daniil Kazantsev + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * http://www.apache.org/licenses/LICENSE-2.0 + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + +#include "mex.h" +#include +#include +#include +#include +#include "omp.h" + +#include "TomoP3DModel_core.h" +#include "utils.h" + +#define M_PI 3.14159265358979323846 + +/* Function to build 3D objects without using the library file + * (MATLAB mex-wrapper) + * + * Input Parameters: + * 1. parameters for 3D object + * 2. DIM volume dimensions [N1,N2,N3] in voxels (N1 x N2 x N3) + * + * VolumeSize in voxels (N x N x N) + * + * Output: + * 1. The analytical phantom size of [N1 x N2 x N3] + */ +#define MAXCHAR 1000 + +void mexFunction( + int nlhs, mxArray *plhs[], + int nrhs, const mxArray *prhs[]) + +{ + float *A; + char *tmpstr2; + float C0 = 0.0f, x0 = 0.0f, y0 = 0.0f, z0 = 0.0f, a = 0.0f, b = 0.0f, c = 0.0f, psi_gr1 = 0.0f, psi_gr2 = 0.0f, psi_gr3 = 0.0f; + double *DimAr; + mwSize N1, N2, N3; + const mwSize *dim_array; + + tmpstr2 = mxArrayToString(prhs[0]); /* name of the object */ + C0 = (float) mxGetScalar(prhs[1]); /* intensity */ + x0 = (float) mxGetScalar(prhs[2]); /* x0 position */ + y0 = (float) mxGetScalar(prhs[3]); /* y0 position */ + z0 = (float) mxGetScalar(prhs[4]); /* z0 position */ + a = (float) mxGetScalar(prhs[5]); /* a - size object */ + b = (float) mxGetScalar(prhs[6]); /* b - size object */ + c = (float) mxGetScalar(prhs[7]); /* c - size object */ + psi_gr1 = (float) mxGetScalar(prhs[8]); /* rotation angle 1*/ + psi_gr2 = (float) mxGetScalar(prhs[9]); /* rotation angle 2*/ + psi_gr3 = (float) mxGetScalar(prhs[10]); /* rotation angle 3*/ + + dim_array = mxGetDimensions(prhs[11]); /* get dimensions of DIM */ + + N1 = 0; N2 = 0; N3 = 0; + if (dim_array[1] == 1) { + N1 = (long) mxGetScalar(prhs[11]); + N2 = N1; + N3 = N1; + } + else if (dim_array[1] == 3) { + DimAr = (double*) mxGetData(prhs[11]); + N1 = (long)(DimAr[0]); + N2 = (long)(DimAr[1]); + N3 = (long)(DimAr[2]); + } + else {mexErrMsgTxt("DIM must be scalar [N] or a vector with 3 elements [N1, N2, N3]");} + + /*Handling Matlab input data*/ + if (nrhs != 12) mexErrMsgTxt("Input of 12 parameters is required"); + + const mwSize N_dims[3] = {N1, N2, N3}; /* image dimensions */ + A = (float*)mxGetPr(plhs[0] = mxCreateNumericArray(3, N_dims, mxSINGLE_CLASS, mxREAL)); /*output array*/ + + + if ((strcmp("gaussian",tmpstr2) != 0) && (strcmp("paraboloid",tmpstr2) != 0) && (strcmp("ellipsoid",tmpstr2) != 0) && (strcmp("cone",tmpstr2) != 0) && (strcmp("cuboid",tmpstr2) != 0) && (strcmp("ellipticalcylinder",tmpstr2) != 0) ) { + printf("%s %s\n", "Unknown name of the object, the given name is", tmpstr2); + mexErrMsgTxt("Unknown name of the object"); + } + if (C0 == 0) { + printf("%s %f\n", "C0 should not be equal to zero, the given value is", C0); + mexErrMsgTxt("C0 should not be equal to zero"); + } + if ((x0 < -1) || (x0 > 1)) { + printf("%s %f\n", "x0 (object position) must be in [-1,1] range, the given value is", x0); + mexErrMsgTxt("x0 (object position) must be in [-1,1] range"); + } + if ((y0 < -1) || (y0 > 1)) { + printf("%s %f\n", "y0 (object position) must be in [-1,1] range, the given value is", y0); + mexErrMsgTxt("y0 (object position) must be in [-1,1] range"); + } + if ((z0 < -1) || (z0 > 1)) { + printf("%s %f\n", "z0 (object position) must be in [-1,1] range, the given value is", z0); + mexErrMsgTxt("z0 (object position) must be in [-1,1] range"); + } + if ((a <= 0) || (a > 2)) { + printf("%s %f\n", "a (object size) must be positive in [0,2] range, the given value is", a); + mexErrMsgTxt("a (object size) must be positive in [0,2] range"); + } + if ((b <= 0) || (b > 2)) { + printf("%s %f\n", "b (object size) must be positive in [0,2] range, the given value is", b); + mexErrMsgTxt("b (object size) must be positive in [0,2] range"); + } + if ((c <= 0) || (c > 2)) { + printf("%s %f\n", "c (object size) must be positive in [0,2] range, the given value is", c); + mexErrMsgTxt("c (object size) must be positive in [0,2] range"); + } + printf("\nObject : %s \nC0 : %f \nx0 : %f \ny0 : %f \nz0 : %f \na : %f \nb : %f \nc : %f \n", tmpstr2, C0, x0, y0, z0, a, b, c); + + TomoP3DObject_core(A, N1, N2, N3, 0l, N3, tmpstr2, C0, x0, y0, z0, b, a, c, -psi_gr1, psi_gr2, psi_gr3, 0); /* Matlab */ + mxFree(tmpstr2); +} 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/DemoReconTomoPy.py b/Wrappers/Python/Demos/DemoReconTomoPy.py new file mode 100644 index 0000000..66c2268 --- /dev/null +++ b/Wrappers/Python/Demos/DemoReconTomoPy.py @@ -0,0 +1,97 @@ +#!/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 analytical phantoms and their sinograms with added noise and artifacts +Sinograms then reconstructed using TomoPy package + +>>>> Prerequisites: TomoPy package <<<<< + +This demo demonstrates frequent inaccuracies which are accosiated with X-ray imaging: +zingers, rings and noise + +@author: Daniil Kazantsev +""" +import numpy as np +import matplotlib.pyplot as plt +from tomophantom import TomoP2D + +model = 4 # select a model +N_size = 512 +#specify a full path to the parameters file +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) + +plt.close('all') +plt.figure(1) +plt.rcParams.update({'font.size': 21}) +plt.imshow(phantom_2D, vmin=0, vmax=1, cmap="gray") +plt.colorbar(ticks=[0, 0.5, 1], orientation='vertical') +plt.title('{}''{}'.format('2D Phantom using model no.',model)) + +# 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="gray") +plt.colorbar(ticks=[0, 150, 250], orientation='vertical') +plt.title('{}''{}'.format('Analytical sinogram of model no.',model)) + +#%% +# Adding artifacts and noise +from tomophantom.supp.artifacts import ArtifactsClass + +# adding noise +artifacts_add =ArtifactsClass(sino_an) +#noisy_sino = artifacts_add.noise(sigma=0.1,noisetype='Gaussian') +noisy_sino = artifacts_add.noise(sigma=10000,noisetype='Poisson') + +# adding zingers +artifacts_add =ArtifactsClass(noisy_sino) +noisy_zing = artifacts_add.zingers(percentage=0.25, modulus = 10) + +#adding stripes +artifacts_add =ArtifactsClass(noisy_zing) +noisy_zing_stripe = artifacts_add.stripes(percentage=1, maxthickness = 1) +noisy_zing_stripe[noisy_zing_stripe < 0] = 0 + +plt.figure() +plt.rcParams.update({'font.size': 21}) +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 gridrec (TomoPy)...") +print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") +import tomopy + +sinoTP = np.zeros((angles_num,1,P),dtype='float32') +sinoTP[:,0,:] = sino_an +rot_center = tomopy.find_center(sinoTP, angles_rad, init=290, ind=0, tol=0.5) +reconTomoPy_ideal = tomopy.recon(sinoTP, angles_rad, center=rot_center, algorithm='gridrec') +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") +plt.colorbar(ticks=[0, 0.5, 1], orientation='vertical') +plt.title('Ideal reconstruction (TomoPy)') +plt.subplot(122) +plt.imshow(reconTomoPy_noisy[0,:,:], vmin=0, vmax=1, cmap="gray") +plt.colorbar(ticks=[0, 0.5, 1], orientation='vertical') +plt.title('Erroneous Reconstruction (TomoPy)') +plt.show() +#%% \ No newline at end of file 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/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/python/Demos/DemoModel.py b/python/Demos/DemoModel.py index e8c2cb6..f2110cc 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) diff --git a/python/Demos/DemoModel2.py b/python/Demos/DemoModel2.py index e99ec16..29e387e 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..0ad8cfc 100644 --- a/python/Demos/DemoModel3D.py +++ b/python/Demos/DemoModel3D.py @@ -21,7 +21,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 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() diff --git a/python/Demos/DemoModel4D.py b/python/Demos/DemoModel4D.py index 5d3c996..5066b60 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() 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..2de570e 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..66c2268 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) From 7522e3b55ec68a26d608288a29096b49e16009bd Mon Sep 17 00:00:00 2001 From: dkazanc Date: Thu, 4 Oct 2018 15:11:20 +0100 Subject: [PATCH 2/5] 3D utiles fixed, demos and a bash script to use cmake --- Core/TomoP3DModel_core.c | 8 +- Core/TomoP3DModel_core.h | 2 - Core/utils.c | 96 ++++-------- Core/utils.h | 6 +- Wrappers/Python/Demos/Reconstruction_Demo.py | 125 --------------- .../Demos/__pycache__/astraOP.cpython-35.pyc | Bin 2978 -> 0 bytes .../__pycache__/libraryToDict.cpython-35.pyc | Bin 1395 -> 0 bytes Wrappers/Python/Demos/model4.txt | 41 ----- Wrappers/Python/Demos/test_ModelObject.py | 147 ------------------ matlab/Model2DGeneratorDemo.m | 4 +- python/Demos/DemoModel3D.py | 4 +- python/Demos/DemoModel4D.py | 2 +- run.sh | 11 ++ 13 files changed, 51 insertions(+), 395 deletions(-) delete mode 100644 Wrappers/Python/Demos/Reconstruction_Demo.py delete mode 100644 Wrappers/Python/Demos/__pycache__/astraOP.cpython-35.pyc delete mode 100644 Wrappers/Python/Demos/__pycache__/libraryToDict.cpython-35.pyc delete mode 100644 Wrappers/Python/Demos/model4.txt delete mode 100644 Wrappers/Python/Demos/test_ModelObject.py create mode 100644 run.sh 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/Python/Demos/Reconstruction_Demo.py b/Wrappers/Python/Demos/Reconstruction_Demo.py deleted file mode 100644 index d2d802c..0000000 --- a/Wrappers/Python/Demos/Reconstruction_Demo.py +++ /dev/null @@ -1,125 +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 analytical phantoms and their sinograms with added noise and artifacts -Sinograms then reconstructed using ASTRA TOOLBOX and TomoPy packages - ->>>> Prerequisites: ASTRA toolbox and TomoPy packages <<<<< - -This demo demonstrates frequent inaccuracies which are accosiated with X-ray imaging: -zingers, rings and noise - -@author: Daniil Kazantsev -""" -import numpy as np -import matplotlib.pyplot as plt -from tomophantom import TomoP2D - -model = 4 # select a model -N_size = 512 -#specify a full path to the parameters file -pathTP = '../../functions/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) - -plt.close('all') -plt.figure(1) -plt.rcParams.update({'font.size': 21}) -plt.imshow(phantom_2D, vmin=0, vmax=1, cmap="gray") -plt.colorbar(ticks=[0, 0.5, 1], orientation='vertical') -plt.title('{}''{}'.format('2D Phantom using model no.',model)) - -# 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="gray") -plt.colorbar(ticks=[0, 150, 250], orientation='vertical') -plt.title('{}''{}'.format('Analytical sinogram of model no.',model)) - -#%% -# Adding artifacts and noise -from tomophantom.supp.artifacts import ArtifactsClass - -# adding noise -artifacts_add =ArtifactsClass(sino_an) -#noisy_sino = artifacts_add.noise(sigma=0.1,noisetype='Gaussian') -noisy_sino = artifacts_add.noise(sigma=10000,noisetype='Poisson') - -# adding zingers -artifacts_add =ArtifactsClass(noisy_sino) -noisy_zing = artifacts_add.zingers(percentage=0.25, modulus = 10) - -#adding stripes -artifacts_add =ArtifactsClass(noisy_zing) -noisy_zing_stripe = artifacts_add.stripes(percentage=1, maxthickness = 1) -noisy_zing_stripe[noisy_zing_stripe < 0] = 0 - -plt.figure() -plt.rcParams.update({'font.size': 21}) -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 ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") -print ("Reconstructing analytical sinogram using gridrec (TomoPy)...") -print ("%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%") -import tomopy - -sinoTP = np.zeros((angles_num,1,P),dtype='float32') -sinoTP[:,0,:] = sino_an -rot_center = tomopy.find_center(sinoTP, angles_rad, init=290, ind=0, tol=0.5) -reconTomoPy_ideal = tomopy.recon(sinoTP, angles_rad, center=rot_center, algorithm='gridrec') -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") -plt.colorbar(ticks=[0, 0.5, 1], orientation='vertical') -plt.title('Ideal reconstruction (TomoPy)') -plt.subplot(122) -plt.imshow(reconTomoPy_noisy[0,:,:], vmin=0, vmax=1, cmap="gray") -plt.colorbar(ticks=[0, 0.5, 1], orientation='vertical') -plt.title('Erroneous Reconstruction (TomoPy)') -plt.show() -#%% \ No newline at end of file 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 ffac956246adea0bce3fb300cb01a7c2e05817c3..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 2978 zcmbuB%Wm676o!YSNJ*4ZJ8^D}+myXJ0_qg-CQV^9Q7j}4P$Y)!bU_FpC~~aO;YFRH zl3H*Uh5G_s_W_DNL7$*6!EIOBb=_6}GZb~PQ*=?vJUTOThT=K%%|E_&;evhdmtWnh z?=bc!TeuAP_i@V(OvsoEtFcgLvBq4DJqisL*4R(zt+7*`xjK#Qcg)GpcGn4FUrJ}- zOCCB|>WpQS4xR4)!Cu$Nvn(EDKj7OLadDW5#2Ja~hzEI;rC$vE7xR|mr=er}vCJHq zCEO7_$Wob$v5Iy4c$kSOf1XHtbJGdNLhv+?kDaZaGxCKW$2@ih+)wN`$bdZw$X1sQ z=RR(^2g8vaQiDYj$PSsoF(^5=M)`%N@=X?6EUYUZPJ@L_<=0tgv#_Onbhc533x(aq zoB9XHR!(-ij*OB~%+I~sa_+%m@@HzeRK38&Ytu#FbnoDnf54n*?1V87j;? z@2gpsn)Q1%3y(tl`B{x>md)}O74uYM8kS+dV?QYx)fx_w6|TST6-F=`7qvJ_`5O&A zt}&!`Fb;h|?a8fu9)lw+dj5*@)Y(A%hLZ_r1A!Y3SaU-DG75O7U9@{VM{Fi!FG`AL zHyy@Ye#e8N{&aA(A5C~+`m})VY9P3ubC0m}hCD;?%yn0$gD14{evvYcAxuco+bU!asE6@`#qjyvahmwx_f(cTr@l{ zN~7HKF5oxhbr_~K)n{G1qP4Y~xRuo>Yjc`!SF-}zM0C|~lUyhWAHd16+yG3Wcm&LJ z*4No<{26LUaY*=Kxhb@4!zsG}?0!d^tanGF_?UzP3-`(w;lf_&n4*Q>aa+^+Ltioh z1!JhKpZi;?#xYj7Qxg}d+l&%_$cq*Q4uzk1ScXQ`P&F>i8U;7#2Whd^(UP<*OkdYl zwKZ+24pr3utivqtc`26aaBWtHfxjpaRSm98Yp^`3sv)l6)x}kqX&uD7@CyqIi-@qQh}M}1wCCtz2HMDNigMY;qoI!IIU(r^yRal!o#wGZd z^i^P|w%zK2{HVzR=?z`9J7?VK|H@q{JXjs|Zwi4n-p*ZhYJn1*a_6|iFlfUV zafc_t19x-U6kBtaKpT!(;G>tyR|2Q-d0k%sPC>1-@$G;qh!@7^z|{dKsI!!w9Qa|CISJ)oI+DmQC(%SE| iEqrdMk&Yogrznb09dCO^(bz6u>1z^GNmN^UTmJ_WxR~Jp 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 02d88dbb6dbac9d8aa7f011a9417592f7bfcec8c..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 1395 zcmZuwO>g5w7=Fi&oy7Uwo?&z?Cuqy;ZqL>9k1ThX~ld!xCb=CJMI*vB!%k^V*i1i%(l*%g6nfXq<*GZVIG zwU!L&8<1IW`4emn*qYElc`HO1K(v-$L6{(NutIgY*G?GosDXl3wn~dR8|aL=?mVqq zPy7EoEf9HCOl=7MaI=n()!m#zH*A#0Sh*P5<>F&Fy>mN9Nj6;m38DbooXR%(`gh1~ zgSUCCy(((u=ul9ngD3if0Sz|~SrIY^vJzxv$SROkA*;dV+aMgw@)l%u*w$fHf^AW6 zL!XQ<`a~D6-SDW4~Yl@m*)9f%y zkXLKu;%_+J+IVB?FFBqk?3CHmUofAB%|Y|K&Y;uFo;#a|uZAp2SWIJ);{1oM|D^T& z@#f&wB#7vVpvk@2ow_bN6OAz*30F%N_u_LZtuzS(!BtIZed~p54Db@n z=}T)7GEcPjWML76H1WhzT9HyF$wHDqI`cQ0rbU55QYB>SFn#gEAQh6Vq~+6e&Zjoc zrh@Z47UBhEipn_3l$0r>nHr4!4E$Y$|4fGi_p6lh)V+OM?ugCTk=jCc%$a^%?r=B` z+|L5>`TESAu!xP9UMyJTCRg|_V|PR&mbzhZ#yx&DVWVI!+R2sN(=wE*-bgUTRPjiU z(SY&5w8>33NrQZZ-O;^D2&odsu#J6kK;9=6QX>|wl<)O-2!D(@WMMXoyda*G{DI1~k9t^gqpPho_D1@= v-js^Hh|XrxR!1}$^M@Fv_wb&1PqX?2i?2iasj8Ez4Y!ed-6(%hMpgd-nTSa( 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/matlab/Model2DGeneratorDemo.m b/matlab/Model2DGeneratorDemo.m index c49553f..44dcbb9 100644 --- a/matlab/Model2DGeneratorDemo.m +++ b/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 = 1; % 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/python/Demos/DemoModel3D.py b/python/Demos/DemoModel3D.py index 0ad8cfc..cff302b 100644 --- a/python/Demos/DemoModel3D.py +++ b/python/Demos/DemoModel3D.py @@ -17,7 +17,7 @@ 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 @@ -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 5066b60..8dedabe 100644 --- a/python/Demos/DemoModel4D.py +++ b/python/Demos/DemoModel4D.py @@ -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/run.sh b/run.sh new file mode 100644 index 0000000..3624419 --- /dev/null +++ b/run.sh @@ -0,0 +1,11 @@ +#!/bin/bash +echo "Building Tomophantom using CMake" +#rm -r build +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 From 9a9b9e375b13c33bd22bb5b0cb642bb335304f2d Mon Sep 17 00:00:00 2001 From: dkazanc Date: Thu, 4 Oct 2018 15:13:42 +0100 Subject: [PATCH 3/5] matlab DIR deleted --- matlab/Model2DGeneratorDemo.m | 77 ----- matlab/Model2D_tGeneratorDemo.m | 96 ------ matlab/Model3DGeneratorDemo.m | 87 ----- matlab/Model3D_tGeneratorDemo.m | 55 ---- matlab/Object2DGeneratorDemo.m | 62 ---- matlab/Object3DGeneratorDemo.m | 39 --- matlab/Reconstruction_Demo.m | 60 ---- matlab/SpectralPhantomDemo.m | 57 ---- matlab/install/compile_mex_linux.m | 37 --- matlab/install/compile_mex_windows.m | 65 ---- matlab/mex_wrappers/CMakeLists.txt | 102 ------ matlab/mex_wrappers/TomoP2DModel.c | 304 ----------------- matlab/mex_wrappers/TomoP2DModelSino.c | 383 ---------------------- matlab/mex_wrappers/TomoP2DObject.c | 96 ------ matlab/mex_wrappers/TomoP2DObjectSino.c | 101 ------ matlab/mex_wrappers/TomoP2DSinoNum.c | 70 ---- matlab/mex_wrappers/TomoP3DModel.c | 419 ------------------------ matlab/mex_wrappers/TomoP3DObject.c | 124 ------- matlab/supplem/add_noise.m | 20 -- matlab/supplem/add_stripes.m | 30 -- matlab/supplem/add_zingers.m | 31 -- matlab/supplem/rec2Dastra.m | 40 --- matlab/supplem/sino2Dastra.m | 16 - 23 files changed, 2371 deletions(-) delete mode 100644 matlab/Model2DGeneratorDemo.m delete mode 100644 matlab/Model2D_tGeneratorDemo.m delete mode 100644 matlab/Model3DGeneratorDemo.m delete mode 100644 matlab/Model3D_tGeneratorDemo.m delete mode 100644 matlab/Object2DGeneratorDemo.m delete mode 100644 matlab/Object3DGeneratorDemo.m delete mode 100644 matlab/Reconstruction_Demo.m delete mode 100644 matlab/SpectralPhantomDemo.m delete mode 100755 matlab/install/compile_mex_linux.m delete mode 100644 matlab/install/compile_mex_windows.m delete mode 100644 matlab/mex_wrappers/CMakeLists.txt delete mode 100644 matlab/mex_wrappers/TomoP2DModel.c delete mode 100644 matlab/mex_wrappers/TomoP2DModelSino.c delete mode 100644 matlab/mex_wrappers/TomoP2DObject.c delete mode 100644 matlab/mex_wrappers/TomoP2DObjectSino.c delete mode 100644 matlab/mex_wrappers/TomoP2DSinoNum.c delete mode 100644 matlab/mex_wrappers/TomoP3DModel.c delete mode 100644 matlab/mex_wrappers/TomoP3DObject.c delete mode 100644 matlab/supplem/add_noise.m delete mode 100644 matlab/supplem/add_stripes.m delete mode 100644 matlab/supplem/add_zingers.m delete mode 100644 matlab/supplem/rec2Dastra.m delete mode 100644 matlab/supplem/sino2Dastra.m diff --git a/matlab/Model2DGeneratorDemo.m b/matlab/Model2DGeneratorDemo.m deleted file mode 100644 index 44dcbb9..0000000 --- a/matlab/Model2DGeneratorDemo.m +++ /dev/null @@ -1,77 +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 = '/'; -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 '..' 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; -%% -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/mex_wrappers/TomoP2DModel.c b/matlab/mex_wrappers/TomoP2DModel.c deleted file mode 100644 index a91a447..0000000 --- a/matlab/mex_wrappers/TomoP2DModel.c +++ /dev/null @@ -1,304 +0,0 @@ -/* - * Copyright 2017 Daniil Kazantsev - * - * Licensed under the Apache License, Version 2.0 (the "License"); - * you may not use this file except in compliance with the License. - * You may obtain a copy of the License at - * http://www.apache.org/licenses/LICENSE-2.0 - * Unless required by applicable law or agreed to in writing, software - * distributed under the License is distributed on an "AS IS" BASIS, - * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. - * See the License for the specific language governing permissions and - * limitations under the License. - */ - -#include -#include "mex.h" -#include -#include -#include -#include -#include "omp.h" - -#include "TomoP2DModel_core.h" -#include "utils.h" - -#define M_PI 3.14159265358979323846 -#define MAXCHAR 1000 - -/* MATLAB mex-wrapper for TomoP2Dmodel - * - * Input Parameters: - * 1. ModelNo - a model number from Phantom2DLibrary file - * 2. VolumeSize in voxels (N x N) or (N x N x time-frames) - * 3. An absolute path to the file Phantom2DLibrary.dat (see OS-specific syntax-differences) - * - * Output: - * 1. The analytical phantom size of [N x N] or a temporal phantom size of [N xN x Time-frames] - */ - -void mexFunction( - int nlhs, mxArray *plhs[], - int nrhs, const mxArray *prhs[]) - -{ - int ModelSelected, N; - float *A; - - - ModelSelected = (int) mxGetScalar(prhs[0]); /* selected model */ - N = (int) mxGetScalar(prhs[1]); /* choosen dimension (N x N x N) */ - - int Model=0, Components=0, steps = 0, counter=0, ii; - float C0 = 0.0f, x0 = 0.0f, y0 = 0.0f, a = 0.0f, b = 0.0f, psi_gr1; - - char *filename; - FILE * fp; - if(!mxIsChar(prhs[2]) ) { - mexErrMsgTxt("Need filename absolute path string input."); - } - filename = mxArrayToString(prhs[2]); - fp= fopen(filename, "rb"); - mxFree(filename); - if( fp == NULL ) { - mexErrMsgTxt("Cannot open the model library file (Phantom3DLibrary.dat)"); - } - else { - - char str[MAXCHAR]; - char tmpstr1[16]; - char tmpstr2[22]; - char tmpstr3[16]; - char tmpstr4[16]; - char tmpstr5[16]; - char tmpstr6[16]; - char tmpstr7[16]; - char tmpstr8[16]; - - while (fgets(str, MAXCHAR, fp) != NULL) - { - /* work with non-# commented lines */ - if(str[0] != '#') { - sscanf(str, "%15s : %21[^;];", tmpstr1, tmpstr2); - if (strcmp(tmpstr1,"Model")==0) - { - Model = atoi(tmpstr2); - if ((ModelSelected == Model) && (counter == 0)) { - /* check if we have a right model */ - if (fgets(str, MAXCHAR, fp) != NULL) sscanf(str, "%15s : %21[^;];", tmpstr1, tmpstr2); - else { - mexErrMsgTxt("Unexpected the end of the line (Components) in parameters file"); - break; } - if (strcmp(tmpstr1,"Components") == 0) Components = atoi(tmpstr2); - printf("%s %i\n", "Components:", Components); - if (Components <= 0) { - printf("%s %i\n", "Components cannot be negative, the given value is", Components); - mexErrMsgTxt("Components cannot be negative"); - break; } - if (fgets(str, MAXCHAR, fp) != NULL) sscanf(str, "%15s : %21[^;];", tmpstr1, tmpstr2); - else { - mexErrMsgTxt("Unexpected the end of the line (TimeSteps) in parameters file"); - break; } - if (strcmp(tmpstr1,"TimeSteps") == 0) steps = atoi(tmpstr2); - if (steps <= 0) { - printf("%s %i\n", "TimeSteps cannot be negative, the given value is", steps); - mexErrMsgTxt("TimeSteps cannot be negative"); - break; } - printf("%s %i\n", "TimeSteps:", steps); - if (steps == 1) { - /**************************************************/ - printf("\n %s \n", "Stationary model is selected"); - - const mwSize N_dims[2] = {N, N}; /* image dimensions */ - A = (float*)mxGetPr(plhs[0] = mxCreateNumericArray(2, N_dims, mxSINGLE_CLASS, mxREAL)); /*output array*/ - - /* loop over all components */ - for(ii=0; ii 1)) { - printf("%s %f\n", "x0 (object position) must be in [-1,1] range, the given value is", x0); - mexErrMsgTxt("x0 (object position) must be in [-1,1] range"); - break; } - if ((y0 < -1) || (y0 > 1)) { - printf("%s %f\n", "y0 (object position) must be in [-1,1] range, the given value is", y0); - mexErrMsgTxt("y0 (object position) must be in [-1,1] range"); - break; } - if ((a <= 0) || (a > 2)) { - printf("%s %f\n", "a (object size) must be positive in [0,2] range, the given value is", a); - mexErrMsgTxt("a (object size) must be positive in [0,2] range"); - break; } - if ((b <= 0) || (b > 2)) { - printf("%s %f\n", "b (object size) must be positive in [0,2] range, the given value is", b); - mexErrMsgTxt("b (object size) must be positive in [0,2] range"); - break; } - printf("\nObject : %s \nC0 : %f \nx0 : %f \ny0 : %f \na : %f \nb : %f \n", tmpstr2, C0, x0, y0, a, b); - - TomoP2DObject_core(A, N, tmpstr2, C0, x0, y0, b, a, -psi_gr1, 0); /* Matlab */ - } - } - else { - /**************************************************/ - printf("\n %s \n", "Temporal model is selected"); - /* temporal phantom 2D + time (3D) */ - const mwSize N_dims[3] = {N, N, steps}; /* image dimensions */ - A = (float*)mxGetPr(plhs[0] = mxCreateNumericArray(3, N_dims, mxSINGLE_CLASS, mxREAL)); - - float C1 = 0.0f, x1 = 0.0f, y1 = 0.0f, a1 = 0.0f, b1 = 0.0f, psi_gr1_1 = 0.0f; - /* loop over all components */ - for(ii=0; ii 1)) { - printf("%s %f\n", "x0 (object position) must be in [-1,1] range, the given value is", x0); - mexErrMsgTxt("x0 (object position) must be in [-1,1] range"); - break; } - if ((y0 < -1) || (y0 > 1)) { - printf("%s %f\n", "y0 (object position) must be in [-1,1] range, the given value is", y0); - mexErrMsgTxt("y0 (object position) must be in [-1,1] range"); - break; } - if ((a <= 0) || (a > 2)) { - printf("%s %f\n", "a (object size) must be positive in [0,2] range, the given value is", a); - mexErrMsgTxt("a (object size) must be positive in [0,2] range"); - break; } - if ((b <= 0) || (b > 2)) { - printf("%s %f\n", "b (object size) must be positive in [0,2] range, the given value is", b); - mexErrMsgTxt("b (object size) must be positive in [0,2] range"); - break; } - // printf("\nObject : %s \nC0 : %f \nx0 : %f \ny0 : %f \nz0 : %f \na : %f \nb : %f \n", tmpstr2, C0, x0, y0, z0, a, b, c); - - /* check Endvar relatedparameters */ - if (fgets(str, MAXCHAR, fp) != NULL) sscanf(str, "%15s : %15s %15s %15s %15s %15s %15[^;];", tmpstr1, tmpstr3, tmpstr4, tmpstr5, tmpstr6, tmpstr7, tmpstr8); - else { - mexErrMsgTxt("Unexpected the end of the line (Endvar loop) in parameters file"); - break; } - - if (strcmp(tmpstr1,"Endvar") == 0) { - C1 = (float)atof(tmpstr3); /* intensity */ - x1 = (float)atof(tmpstr4); /* x0 position */ - y1 = (float)atof(tmpstr5); /* y0 position */ - a1 = (float)atof(tmpstr6); /* a - size object */ - b1 = (float)atof(tmpstr7); /* b - size object */ - psi_gr1_1 = (float)atof(tmpstr8); /* rotation angle 1*/ - } - else { - printf("%s\n", "Cannot find 'Endvar' string in parameters file"); - break; } - - if (C1 == 0) { - printf("%s %f\n", "Endvar C1 should not be equal to zero, the given value is", C1); - mexErrMsgTxt("Endvar C0 should not be equal to zero"); - break; } - if ((x1 < -1) || (x1 > 1)) { - printf("%s %f\n", "Endvar x1 (object position) must be in [-1,1] range, the given value is", x1); - mexErrMsgTxt("Endvar x0 (object position) must be in [-1,1] range"); - break; } - if ((y1 < -1) || (y1 > 1)) { - printf("%s %f\n", "Endvar y1 (object position) must be in [-1,1] range, the given value is", y1); - mexErrMsgTxt("Endvar y0 (object position) must be in [-1,1] range"); - break; } - if ((a1 <= 0) || (a1 > 2)) { - printf("%s %f\n", "Endvar a1 (object size) must be positive in [0,2] range, the given value is", a1); - mexErrMsgTxt("Endvar a (object size) must be positive in [0,2] range"); - break; } - if ((b1 <= 0) || (b1 > 2)) { - printf("%s %f\n", "Endvar b1 (object size) must be positive in [0,2] range, the given value is", b1); - mexErrMsgTxt("Endvar b (object size) must be positive in [0,2] range"); - break; } - //printf("\nObject : %s \nC0 : %f \nx0 : %f \ny0 : %f \nz0 : %f \na : %f \nb : %f \nc : %f \n", tmpstr2, C0, x0, y0, z0, a1, b1, c1); - - /*now we know the initial parameters of the object and the final ones. We linearly extrapolate to establish steps and coordinates. */ - /* calculating the full distance berween the start and the end points */ - float distance = sqrtf(pow((x1 - x0),2) + pow((y1 - y0),2)); - float d_dist = distance/(steps-1); /*a step over line */ - float C_step = (C1 - C0)/(steps-1); - float a_step = (a1 - a)/(steps-1); - float b_step = (b1 - b)/(steps-1); - float phi_rot_step = (psi_gr1_1 - psi_gr1)/(steps-1); - - int tt; - float x_t, y_t, a_t, b_t, C_t, phi_t, d_step; - /* initialize */ - x_t = x0; y_t = y0; a_t = a; b_t = b; C_t = C0; phi_t = psi_gr1; d_step = d_dist; - /*loop over time frames*/ - for(tt=0; tt < steps; tt++) { - - TomoP2DObject_core(A, N, tmpstr2, C_t, x_t, y_t, b_t, a_t, -phi_t, tt); /* Matlab */ - - /* calculating new coordinates of an object */ - if (distance != 0.0f) { - float t = d_step/distance; - x_t = (1-t)*x0 + t*x1; - y_t = (1-t)*y0 + t*y1; } - else { - x_t = x0; - y_t = y0; } - - d_step += d_dist; - a_t += a_step; - b_t += b_step; - C_t += C_step; - phi_t += phi_rot_step; - } /*time steps*/ - - } /*components loop*/ - } - counter++; - } - } - } - } - } - fclose(fp); - if (counter == 0) { - printf("%s %i %s \n", "Model no. ", ModelSelected, "is not found!"); - mexErrMsgTxt("No object found, check models file"); - } -} diff --git a/matlab/mex_wrappers/TomoP2DModelSino.c b/matlab/mex_wrappers/TomoP2DModelSino.c deleted file mode 100644 index 7908a61..0000000 --- a/matlab/mex_wrappers/TomoP2DModelSino.c +++ /dev/null @@ -1,383 +0,0 @@ -/* - * Copyright 2017 Daniil Kazantsev - * - * Licensed under the Apache License, Version 2.0 (the "License"); - * you may not use this file except in compliance with the License. - * You may obtain a copy of the License at - * http://www.apache.org/licenses/LICENSE-2.0 - * Unless required by applicable law or agreed to in writing, software - * distributed under the License is distributed on an "AS IS" BASIS, - * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. - * See the License for the specific language governing permissions and - * limitations under the License. - */ - -#include "mex.h" -#include -#include -#include -#include -#include "omp.h" - -#include "TomoP2DModelSino_core.h" -#include "utils.h" - -#define M_PI 3.14159265358979323846 -#define MAXCHAR 1000 - -/* MATLAB mex-wrapper to generate parallel beam sinograms using TomoP2DModelSino_core.c - * - * Input Parameters: - * 1. Model number (see Phantom2DLibrary.dat) [required] - * 2. ImageSize in pixels (N x N) [required] - * 3. Detector array size P (in pixels) [required] - * 4. Projection angles Theta (in degrees) [required] - * 5. An absolute path to the file Phantom2DLibrary.dat (see OS-specific syntax-differences) [required] - * 6. ImageCentring, choose 'radon' or 'astra' (default) [optional] - * - * Output: - * 1. 2D sinogram size of [length(angles), P] or a temporal sinogram size of [length(angles), P, Time-Frames] - */ - -void mexFunction( - int nlhs, mxArray *plhs[], - int nrhs, const mxArray *prhs[]) - - { - int ModelSelected, N, CenTypeIn, P; - float *A, *Th; - mwSize NStructElems; - char *filename; - FILE * fp; - - /*Handling Matlab input data*/ - if ((nrhs < 5) && (nrhs > 6)) mexErrMsgTxt("Input of 5 or 6 parameters is required: Model, ImageSize, Detector Size, Projection angles, PATH, Centering"); - if (mxGetClassID(prhs[3]) != mxSINGLE_CLASS) {mexErrMsgTxt("The vector of angles must be in a single precision"); } - - ModelSelected = (int) mxGetScalar(prhs[0]); /* selected model */ - N = (int) mxGetScalar(prhs[1]); /* choosen dimension (N x N) */ - P = (int) mxGetScalar(prhs[2]); /* detector size */ - Th = (float*) mxGetData(prhs[3]); /* angles */ - CenTypeIn = 1; /* astra-type centering is the default one */ - - if (nrhs == 6) { - char *CenType; - CenType = mxArrayToString(prhs[5]); /* 'radon' or 'astra' (default) */ - if ((strcmp(CenType, "radon") != 0) && (strcmp(CenType, "astra") != 0)) mexErrMsgTxt("Choose 'radon' or 'astra''"); - if (strcmp(CenType, "radon") == 0) CenTypeIn = 0; /* enable 'radon'-type centaering */ - mxFree(CenType); - } - NStructElems = mxGetNumberOfElements(prhs[3]); - - int Model=0, Components=0, steps = 0, counter=0, ii; - float C0 = 0.0f, x0 = 0.0f, y0 = 0.0f, a = 0.0f, b = 0.0f, psi_gr1; - - if(!mxIsChar(prhs[4]) ) { - mexErrMsgTxt("Need filename absolute path string input."); - } - filename = mxArrayToString(prhs[4]); - fp= fopen(filename, "r"); - mxFree(filename); - if( fp == NULL ) { - mexErrMsgTxt("Cannot open the model library file (Phantom3DLibrary.dat)"); - } - else { - - char str[MAXCHAR]; - char tmpstr1[16]; - char tmpstr2[22]; - char tmpstr3[16]; - char tmpstr4[16]; - char tmpstr5[16]; - char tmpstr6[16]; - char tmpstr7[16]; - char tmpstr8[16]; - - while (fgets(str, MAXCHAR, fp) != NULL) - { - /* work with non-# commented lines */ - if(str[0] != '#') { - sscanf(str, "%15s : %21[^;];", tmpstr1, tmpstr2); - if (strcmp(tmpstr1,"Model")==0) - { - Model = atoi(tmpstr2); - if ((ModelSelected == Model) && (counter == 0)) { - /* check if we have a right model */ - if (fgets(str, MAXCHAR, fp) != NULL) sscanf(str, "%15s : %21[^;];", tmpstr1, tmpstr2); - else { - mexErrMsgTxt("Unexpected the end of the line (Components) in parameters file"); - break; } - if (strcmp(tmpstr1,"Components") == 0) Components = atoi(tmpstr2); - printf("%s %i\n", "Components:", Components); - if (Components <= 0) { - printf("%s %i\n", "Components cannot be negative, the given value is", Components); - mexErrMsgTxt("Components cannot be negative"); - break; } - if (fgets(str, MAXCHAR, fp) != NULL) sscanf(str, "%15s : %21[^;];", tmpstr1, tmpstr2); - else { - mexErrMsgTxt("Unexpected the end of the line (TimeSteps) in parameters file"); - break; } - if (strcmp(tmpstr1,"TimeSteps") == 0) steps = atoi(tmpstr2); - if (steps <= 0) { - printf("%s %i\n", "TimeSteps cannot be negative, the given value is", steps); - mexErrMsgTxt("TimeSteps cannot be negative"); - break; } - printf("%s %i\n", "TimeSteps:", steps); - if (steps == 1) { - /**************************************************/ - printf("\n %s %i %s \n", "Stationary 2D model", ModelSelected, " is selected"); - - const mwSize N_dims[2] = {NStructElems,P}; /*format: X-detectors, Y-angles dim*/ - A = (float*)mxGetPr(plhs[0] = mxCreateNumericArray(2, N_dims, mxSINGLE_CLASS, mxREAL)); - - /* loop over all components */ - for(ii=0; ii 1)) { - printf("%s %f\n", "x0 (object position) must be in [-1,1] range, the given value is", x0); - mexErrMsgTxt("x0 (object position) must be in [-1,1] range"); - break; } - if ((y0 < -1) || (y0 > 1)) { - printf("%s %f\n", "y0 (object position) must be in [-1,1] range, the given value is", y0); - mexErrMsgTxt("y0 (object position) must be in [-1,1] range"); - break; } - if ((a <= 0) || (a > 2)) { - printf("%s %f\n", "a (object size) must be positive in [0,2] range, the given value is", a); - mexErrMsgTxt("a (object size) must be positive in [0,2] range"); - break; } - if ((b <= 0) || (b > 2)) { - printf("%s %f\n", "b (object size) must be positive in [0,2] range, the given value is", b); - mexErrMsgTxt("b (object size) must be positive in [0,2] range"); - break; } - printf("\nObject : %s \nC0 : %f \nx0 : %f \ny0 : %f \na : %f \nb : %f \n", tmpstr2, C0, x0, y0, a, b); - - TomoP2DObjectSino_core(A, N, P, Th, (int)NStructElems, CenTypeIn, tmpstr2, C0, x0, y0, b, a, -psi_gr1, 0); /* Matlab */ - } - } - else { - /**************************************************/ - printf("\n %s %i %s \n", "Temporal 2D+time model", ModelSelected, " is selected"); - /* temporal phantom 2D + time (3D) */ - const mwSize N_dims[3] = {NStructElems, P, steps}; /*format: X-detectors, Y-angles dim, Time-Frames*/ - A = (float*)mxGetPr(plhs[0] = mxCreateNumericArray(3, N_dims, mxSINGLE_CLASS, mxREAL)); - - float C1 = 0.0f, x1 = 0.0f, y1 = 0.0f, a1 = 0.0f, b1 = 0.0f, psi_gr1_1 = 0.0f; - /* loop over all components */ - for(ii=0; ii 1)) { - printf("%s %f\n", "x0 (object position) must be in [-1,1] range, the given value is", x0); - mexErrMsgTxt("x0 (object position) must be in [-1,1] range"); - break; } - if ((y0 < -1) || (y0 > 1)) { - printf("%s %f\n", "y0 (object position) must be in [-1,1] range, the given value is", y0); - mexErrMsgTxt("y0 (object position) must be in [-1,1] range"); - break; } - if ((a <= 0) || (a > 2)) { - printf("%s %f\n", "a (object size) must be positive in [0,2] range, the given value is", a); - mexErrMsgTxt("a (object size) must be positive in [0,2] range"); - break; } - if ((b <= 0) || (b > 2)) { - printf("%s %f\n", "b (object size) must be positive in [0,2] range, the given value is", b); - mexErrMsgTxt("b (object size) must be positive in [0,2] range"); - break; } - // printf("\nObject : %s \nC0 : %f \nx0 : %f \ny0 : %f \nz0 : %f \na : %f \nb : %f \n", tmpstr2, C0, x0, y0, z0, a, b, c); - - /* check Endvar relatedparameters */ - if (fgets(str, MAXCHAR, fp) != NULL) sscanf(str, "%15s : %15s %15s %15s %15s %15s %15[^;];", tmpstr1, tmpstr3, tmpstr4, tmpstr5, tmpstr6, tmpstr7, tmpstr8); - else { - mexErrMsgTxt("Unexpected the end of the line (Endvar loop) in parameters file"); - break; } - - if (strcmp(tmpstr1,"Endvar") == 0) { - C1 = (float)atof(tmpstr3); /* intensity */ - x1 = (float)atof(tmpstr4); /* x0 position */ - y1 = (float)atof(tmpstr5); /* y0 position */ - a1 = (float)atof(tmpstr6); /* a - size object */ - b1 = (float)atof(tmpstr7); /* b - size object */ - psi_gr1_1 = (float)atof(tmpstr8); /* rotation angle 1*/ - } - else { - printf("%s\n", "Cannot find 'Endvar' string in parameters file"); - break; } - - if (C1 == 0) { - printf("%s %f\n", "Endvar C1 should not be equal to zero, the given value is", C1); - mexErrMsgTxt("Endvar C0 should not be equal to zero"); - break; } - if ((x1 < -1) || (x1 > 1)) { - printf("%s %f\n", "Endvar x1 (object position) must be in [-1,1] range, the given value is", x1); - mexErrMsgTxt("Endvar x0 (object position) must be in [-1,1] range"); - break; } - if ((y1 < -1) || (y1 > 1)) { - printf("%s %f\n", "Endvar y1 (object position) must be in [-1,1] range, the given value is", y1); - mexErrMsgTxt("Endvar y0 (object position) must be in [-1,1] range"); - break; } - if ((a1 <= 0) || (a1 > 2)) { - printf("%s %f\n", "Endvar a1 (object size) must be positive in [0,2] range, the given value is", a1); - mexErrMsgTxt("Endvar a (object size) must be positive in [0,2] range"); - break; } - if ((b1 <= 0) || (b1 > 2)) { - printf("%s %f\n", "Endvar b1 (object size) must be positive in [0,2] range, the given value is", b1); - mexErrMsgTxt("Endvar b (object size) must be positive in [0,2] range"); - break; } - //printf("\nObject : %s \nC0 : %f \nx0 : %f \ny0 : %f \nz0 : %f \na : %f \nb : %f \nc : %f \n", tmpstr2, C0, x0, y0, z0, a1, b1, c1); - - /*now we know the initial parameters of the object and the final ones. We linearly extrapolate to establish steps and coordinates. */ - /* calculating the full distance berween the start and the end points */ - float distance = sqrtf(pow((x1 - x0),2) + pow((y1 - y0),2)); - float d_dist = distance/(steps-1); /*a step over line */ - float C_step = (C1 - C0)/(steps-1); - float a_step = (a1 - a)/(steps-1); - float b_step = (b1 - b)/(steps-1); - float phi_rot_step = (psi_gr1_1 - psi_gr1)/(steps-1); - - int tt; - float x_t, y_t, a_t, b_t, C_t, phi_t, d_step; - /* initialize */ - x_t = x0; y_t = y0; a_t = a; b_t = b; C_t = C0; phi_t = psi_gr1; d_step = d_dist; - /*loop over time frames*/ - for(tt=0; tt < steps; tt++) { - - TomoP2DObjectSino_core(A, N, P, Th, (int)NStructElems, CenTypeIn, tmpstr2, C_t, x_t, y_t, b_t, a_t, -phi_t, tt); - - /* calculating new coordinates of an object */ - if (distance != 0.0f) { - float t = d_step/distance; - x_t = (1-t)*x0 + t*x1; - y_t = (1-t)*y0 + t*y1; } - else { - x_t = x0; - y_t = y0; } - - d_step += d_dist; - a_t += a_step; - b_t += b_step; - C_t += C_step; - phi_t += phi_rot_step; - } /*time steps*/ - - } /*components loop*/ - } - counter++; - } - } - } - } - } - fclose(fp); - if (counter == 0) { - printf("%s %i %s \n", "Model no. ", ModelSelected, "is not found!"); - mexErrMsgTxt("No object found, check models file"); - } -} - - - -// /* first to check if the model is stationary or temporal */ -// int *timeFrames = (int *) mxGetData(mxCreateNumericMatrix(1, 1, mxINT16_CLASS, mxREAL)); -// /* extract the number of time-step */ -// extractTimeFrames(timeFrames, ModelSelected, ModelParameters_PATH); -// -// int *params_switch = (int *) mxGetData(mxCreateNumericMatrix(1, 10, mxINT16_CLASS, mxREAL)); -// /* run checks for parameters (only integer switches and the number of time-frames) */ -// checkParams2D(params_switch, ModelSelected, ModelParameters_PATH, timeFrames[0]); -// -// if (params_switch[0] == 0) mexErrMsgTxt("Parameters file (Phantom2DLibrary.dat) cannot be read, check that your PATH to the file is a valid one and the file exists"); -// if (params_switch[1] == 0) { -// printf("%s %i\n", ">>>>> No such model available in the library, the given model is", ModelSelected); -// mexErrMsgTxt("Check the syntax in Phantom2DLibrary.dat"); -// } -// if (params_switch[2] == 0) { -// printf("%s %i\n", ">>>>> Components number cannot be negative for the given model", ModelSelected); -// mexErrMsgTxt("Check the syntax in Phantom2DLibrary.dat"); -// } -// if (params_switch[3] == 0) { -// printf("%s %i\n", ">>>>> TimeSteps value cannot be negative for the given model", ModelSelected); -// mexErrMsgTxt("Check the syntax in Phantom2DLibrary.dat"); -// } -// if (params_switch[4] == 0) { -// printf("%s %i\n", ">>>>> Unknown name of the object for the given model", ModelSelected); -// mexErrMsgTxt("Check the syntax in Phantom2DLibrary.dat"); -// } -// if (params_switch[5] == 0) { -// printf("%s %i\n", ">>>>> C0 should not be equal to zero for the given model", ModelSelected); -// mexErrMsgTxt("Check the syntax in Phantom2DLibrary.dat"); -// } -// if (params_switch[6] == 0) { -// printf("%s %i\n", ">>>>> x0 (object position) must be in [-1,1] range for the given model", ModelSelected); -// mexErrMsgTxt("Check the syntax in Phantom2DLibrary.dat"); -// } -// if (params_switch[7] == 0) { -// printf("%s %i\n", ">>>>> y0 (object position) must be in [-1,1] range for the given model", ModelSelected); -// mexErrMsgTxt("Check the syntax in Phantom2DLibrary.dat"); -// } -// if (params_switch[8] == 0) { -// printf("%s %i\n", ">>>>> a (object size) must be positive in [0,2] range", ModelSelected); -// mexErrMsgTxt("Check the syntax in Phantom2DLibrary.dat"); -// } -// if (params_switch[9] == 0) { -// printf("%s %i\n", ">>>>> b (object size) must be positive in [0,2] range", ModelSelected); -// mexErrMsgTxt("Check the syntax in Phantom2DLibrary.dat"); -// } -// -// /*Handling Matlab output data*/ -// if (params_switch[3] == 1) { -// /* static phantom case */ -// int N_dims[] = {NStructElems,P}; /*format: X-detectors, Y-angles dim*/ -// A = (float*)mxGetPr(plhs[0] = mxCreateNumericArray(2, N_dims, mxSINGLE_CLASS, mxREAL)); } -// else { -// int N_dims[] = {NStructElems, P, params_switch[3]}; /*format: X-detectors, Y-angles dim, Time-Frames*/ -// A = (float*)mxGetPr(plhs[0] = mxCreateNumericArray(3, N_dims, mxSINGLE_CLASS, mxREAL)); } -// TomoP2DModelSino_core(A, ModelSelected, N, P, Th, (int)NStructElems, CenTypeIn, ModelParameters_PATH, 0); -// mxFree(ModelParameters_PATH); -// } diff --git a/matlab/mex_wrappers/TomoP2DObject.c b/matlab/mex_wrappers/TomoP2DObject.c deleted file mode 100644 index 1276b6f..0000000 --- a/matlab/mex_wrappers/TomoP2DObject.c +++ /dev/null @@ -1,96 +0,0 @@ -/* - * Copyright 2017 Daniil Kazantsev - * - * Licensed under the Apache License, Version 2.0 (the "License"); - * you may not use this file except in compliance with the License. - * You may obtain a copy of the License at - * http://www.apache.org/licenses/LICENSE-2.0 - * Unless required by applicable law or agreed to in writing, software - * distributed under the License is distributed on an "AS IS" BASIS, - * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. - * See the License for the specific language governing permissions and - * limitations under the License. - */ - -#include "mex.h" -#include -#include -#include -#include -#include "omp.h" - -#include "TomoP2DModel_core.h" -#include "utils.h" - -#define M_PI 3.14159265358979323846 - -/* Function to build 2D objects without using the library file - * (MATLAB mex-wrapper) - * - * Input Parameters: - * 1. parameters for 2D object - * 2. Size of the 2D object - * - * VolumeSize in voxels (N x N ) - * - * Output: - * 1. The analytical phantom size of [N x N] - * - */ -#define MAXCHAR 1000 - -void mexFunction( - int nlhs, mxArray *plhs[], - int nrhs, const mxArray *prhs[]) - -{ - int N; - float *A; - char *tmpstr2; - float C0 = 0.0f, x0 = 0.0f, y0 = 0.0f, a = 0.0f, b = 0.0f, psi_gr1 = 0.0f; - - tmpstr2 = mxArrayToString(prhs[0]); /* name of the object */ - C0 = (float) mxGetScalar(prhs[1]); /* intensity */ - x0 = (float) mxGetScalar(prhs[2]); /* x0 position */ - y0 = (float) mxGetScalar(prhs[3]); /* y0 position */ - a = (float) mxGetScalar(prhs[4]); /* a - size object */ - b = (float) mxGetScalar(prhs[5]); /* b - size object */ - psi_gr1 = (float) mxGetScalar(prhs[6]); /* rotation angle 1*/ - N = (int) mxGetScalar(prhs[7]); /* choosen dimension (N x N x N) */ - - /*Handling Matlab input data*/ - if (nrhs != 8) mexErrMsgTxt("Input of 8 parameters is required"); - - const mwSize N_dims[2] = {N, N}; /* image dimensions */ - A = (float*)mxGetPr(plhs[0] = mxCreateNumericArray(2, N_dims, mxSINGLE_CLASS, mxREAL)); /*output array*/ - - if ((strcmp("gaussian",tmpstr2) != 0) && (strcmp("parabola",tmpstr2) != 0) && (strcmp("ellipse",tmpstr2) != 0) && (strcmp("parabola1",tmpstr2) != 0) && (strcmp("cone",tmpstr2) != 0) && (strcmp("rectangle",tmpstr2) != 0) ) { - printf("%s %s\n", "Unknown name of the object, the given name is", tmpstr2); - mexErrMsgTxt("Unknown name of the object"); - } - if (C0 == 0) { - printf("%s %f\n", "C0 should not be equal to zero, the given value is", C0); - mexErrMsgTxt("C0 should not be equal to zero"); - } - if ((x0 < -1) || (x0 > 1)) { - printf("%s %f\n", "x0 (object position) must be in [-1,1] range, the given value is", x0); - mexErrMsgTxt("x0 (object position) must be in [-1,1] range"); - } - if ((y0 < -1) || (y0 > 1)) { - printf("%s %f\n", "y0 (object position) must be in [-1,1] range, the given value is", y0); - mexErrMsgTxt("y0 (object position) must be in [-1,1] range"); - } - if ((a <= 0) || (a > 2)) { - printf("%s %f\n", "a (object size) must be positive in [0,2] range, the given value is", a); - mexErrMsgTxt("a (object size) must be positive in [0,2] range"); - } - if ((b <= 0) || (b > 2)) { - printf("%s %f\n", "b (object size) must be positive in [0,2] range, the given value is", b); - mexErrMsgTxt("b (object size) must be positive in [0,2] range"); - } - printf("\nObject : %s \nC0 : %f \nx0 : %f \ny0 : %f \na : %f \nb : %f \n", tmpstr2, C0, x0, y0, a, b); - - TomoP2DObject_core(A, N, tmpstr2, C0, x0, y0, b, a, -psi_gr1, 0); /* Matlab */ - - mxFree(tmpstr2); -} diff --git a/matlab/mex_wrappers/TomoP2DObjectSino.c b/matlab/mex_wrappers/TomoP2DObjectSino.c deleted file mode 100644 index d4b8252..0000000 --- a/matlab/mex_wrappers/TomoP2DObjectSino.c +++ /dev/null @@ -1,101 +0,0 @@ -/* - * Copyright 2017 Daniil Kazantsev - * - * Licensed under the Apache License, Version 2.0 (the "License"); - * you may not use this file except in compliance with the License. - * You may obtain a copy of the License at - * http://www.apache.org/licenses/LICENSE-2.0 - * Unless required by applicable law or agreed to in writing, software - * distributed under the License is distributed on an "AS IS" BASIS, - * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. - * See the License for the specific language governing permissions and - * limitations under the License. - */ - -#include "mex.h" -#include -#include -#include -#include -#include "omp.h" - -#include "TomoP2DModelSino_core.h" -#include "utils.h" - -/* Function to build 2D sinogram objects without using the library file - * (MATLAB mex-wrapper) - * - * Input Parameters: - * 1. parameters for 2D sinogram - * 3. - * 2. Size of the 2D object - * - * VolumeSize in voxels (N x N ) - * - * Output: - * 1. The analytical phantom size of [N x N] - * - */ -#define MAXCHAR 1000 - -void mexFunction( - int nlhs, mxArray *plhs[], - int nrhs, const mxArray *prhs[]) - -{ - mwSize NStructElems; - int N, P; - float *A, *Th; - char *tmpstr2; - float C0 = 0.0f, x0 = 0.0f, y0 = 0.0f, a = 0.0f, b = 0.0f, psi_gr1 = 0.0f; - - tmpstr2 = mxArrayToString(prhs[0]); /* name of the object */ - C0 = (float) mxGetScalar(prhs[1]); /* intensity */ - x0 = (float) mxGetScalar(prhs[2]); /* x0 position */ - y0 = (float) mxGetScalar(prhs[3]); /* y0 position */ - a = (float) mxGetScalar(prhs[4]); /* a - size object */ - b = (float) mxGetScalar(prhs[5]); /* b - size object */ - psi_gr1 = (float) mxGetScalar(prhs[6]); /* rotation angle 1*/ - N = (int) mxGetScalar(prhs[7]); /* choosen dimension (N x N) */ - P = (int) mxGetScalar(prhs[8]); /* detector size */ - Th = (float*) mxGetData(prhs[9]); /* angles */ - NStructElems = mxGetNumberOfElements(prhs[9]); - - int CenTypeIn = 1; /* astra-type centering is the default one */ - - /*Handling Matlab input data*/ - if (nrhs != 10) mexErrMsgTxt("Input of 10 parameters is required"); - - const mwSize N_dims[2] = {NStructElems,P}; /*format: X-detectors, Y-angles dim*/ - A = (float*)mxGetPr(plhs[0] = mxCreateNumericArray(2, N_dims, mxSINGLE_CLASS, mxREAL)); - - if ((strcmp("gaussian",tmpstr2) != 0) && (strcmp("parabola",tmpstr2) != 0) && (strcmp("ellipse",tmpstr2) != 0) && (strcmp("parabola1",tmpstr2) != 0) && (strcmp("cone",tmpstr2) != 0) && (strcmp("rectangle",tmpstr2) != 0) ) { - printf("%s %s\n", "Unknown name of the object, the given name is", tmpstr2); - mexErrMsgTxt("Unknown name of the object"); - } - if (C0 == 0) { - printf("%s %f\n", "C0 should not be equal to zero, the given value is", C0); - mexErrMsgTxt("C0 should not be equal to zero"); - } - if ((x0 < -1) || (x0 > 1)) { - printf("%s %f\n", "x0 (object position) must be in [-1,1] range, the given value is", x0); - mexErrMsgTxt("x0 (object position) must be in [-1,1] range"); - } - if ((y0 < -1) || (y0 > 1)) { - printf("%s %f\n", "y0 (object position) must be in [-1,1] range, the given value is", y0); - mexErrMsgTxt("y0 (object position) must be in [-1,1] range"); - } - if ((a <= 0) || (a > 2)) { - printf("%s %f\n", "a (object size) must be positive in [0,2] range, the given value is", a); - mexErrMsgTxt("a (object size) must be positive in [0,2] range"); - } - if ((b <= 0) || (b > 2)) { - printf("%s %f\n", "b (object size) must be positive in [0,2] range, the given value is", b); - mexErrMsgTxt("b (object size) must be positive in [0,2] range"); - } - printf("\nObject : %s \nC0 : %f \nx0 : %f \ny0 : %f \na : %f \nb : %f \n", tmpstr2, C0, x0, y0, a, b); - - TomoP2DObjectSino_core(A, N, P, Th, (int)NStructElems, CenTypeIn, tmpstr2, C0, x0, y0, b, a, -psi_gr1, 0); /* Matlab */ - - mxFree(tmpstr2); -} diff --git a/matlab/mex_wrappers/TomoP2DSinoNum.c b/matlab/mex_wrappers/TomoP2DSinoNum.c deleted file mode 100644 index d170bdb..0000000 --- a/matlab/mex_wrappers/TomoP2DSinoNum.c +++ /dev/null @@ -1,70 +0,0 @@ -/* - * Copyright 2017 Daniil Kazantsev - * - * Licensed under the Apache License, Version 2.0 (the "License"); - * you may not use this file except in compliance with the License. - * You may obtain a copy of the License at - * http://www.apache.org/licenses/LICENSE-2.0 - * Unless required by applicable law or agreed to in writing, software - * distributed under the License is distributed on an "AS IS" BASIS, - * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. - * See the License for the specific language governing permissions and - * limitations under the License. - */ - -#include "mex.h" -#include -#include -#include -#include -#include -#include -#include "omp.h" - -#include "TomoP2DSinoNum_core.h" -#include "utils.h" - - -/* C OMP implementation of the forward projection (The Radon Transform) - * by rotating a padded image and summing over columns (rotation-based projector - * for parallel beam) - * - * Input Parameters: - * 1. Phantom to calculate projections from [required] - * 2. Detector array size P (in pixels) [required] - * 3. Projection angles Theta (in degrees) [required] - * - * Output: - * Sinogram [No. angles] x [No. detectors] - * - * mex ForwardProjNum.c CFLAGS="\$CFLAGS -fopenmp -Wall -std=c99" LDFLAGS="\$LDFLAGS -fopenmp" - * sinogram = ForwardProjNum(single(G), P, single(angles)); - */ - -void mexFunction( - int nlhs, mxArray *plhs[], - int nrhs, const mxArray *prhs[]) { - int dimX, dimY,ThetaLength, DetSize, sys=1; - const int *dim_array, *dimsTh; - float *Phantom, *Sinogram, *Theta; - - /*Handling Matlab input data*/ - Phantom = (float *) mxGetData(prhs[0]); - DetSize = (int) mxGetScalar(prhs[1]); - Theta = (float *) mxGetData(prhs[2]); - - dim_array = mxGetDimensions(prhs[0]); - dimsTh = mxGetDimensions(prhs[2]); - - dimX = dim_array[0]; dimY = dim_array[1]; /*Image Dimensions (must be squared) */ - if( dimX != dimY) mexErrMsgTxt("The input image must be squared"); - if( DetSize < dimX) mexErrMsgTxt("The detector dimension is smaller than the phantom dimension! "); - - ThetaLength = dimsTh[1]; /* angles dimension */ - - /*Handling Matlab output data*/ - const mwSize N_dims[2] = {ThetaLength, DetSize}; /*format: Y-angles x X-detectors dim*/ - Sinogram = (float*)mxGetPr(plhs[0] = mxCreateNumericArray(2, N_dims, mxSINGLE_CLASS, mxREAL)); - - TomoP2DSinoNum_core(Sinogram, Phantom, dimX, DetSize, Theta, ThetaLength, sys); -} diff --git a/matlab/mex_wrappers/TomoP3DModel.c b/matlab/mex_wrappers/TomoP3DModel.c deleted file mode 100644 index 30393b0..0000000 --- a/matlab/mex_wrappers/TomoP3DModel.c +++ /dev/null @@ -1,419 +0,0 @@ -/* - * Copyright 2017 Daniil Kazantsev - * - * Licensed under the Apache License, Version 2.0 (the "License"); - * you may not use this file except in compliance with the License. - * You may obtain a copy of the License at - * http://www.apache.org/licenses/LICENSE-2.0 - * Unless required by applicable law or agreed to in writing, software - * distributed under the License is distributed on an "AS IS" BASIS, - * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. - * See the License for the specific language governing permissions and - * limitations under the License. - */ - -#include "mex.h" -#include -#include -#include -#include -#include "omp.h" - -#include "TomoP3DModel_core.h" -#include "utils.h" - -/* Function to read from the file Phantom3DLibrary.dat the required parameters to build 3D or 3D + time analytical models - * (MATLAB mex-wrapper) - * - * Input Parameters: - * 1. ModelNo - a model number from Phantom3DLibrary file - * 2. DIM volume dimensions [N1,N2,N3] in voxels (N1 x N2 x N3) or (N1 x N2 x N3 x time-frames) - * 3. An absolute path to the file Phantom3DLibrary.dat (ssee OS-specific syntax-differences) - * - * Output: - * 1. The analytical phantom size of [N1 x N2 x N3] or temporal 4D phantom (N1 x N2 x N3 x time-frames) - * - */ -#define MAXCHAR 1000 - -void mexFunction( - int nlhs, mxArray *plhs[], - int nrhs, const mxArray *prhs[]) - -{ - int ModelSelected; - float *A; - double *DimAr; - mwSize N1, N2, N3; - const mwSize *dim_array; - - /*Handling Matlab input data*/ - if (nrhs != 3) mexErrMsgTxt("Input of 3 parameters is required: model, DIMensions, PATH"); - - dim_array = mxGetDimensions(prhs[1]); /* get dimensions of DIM */ - - //printf("%i \n", dim_array[1]); - - N1 = 0; N2 = 0; N3 = 0; - if (dim_array[1] == 1) { - N1 = (long) mxGetScalar(prhs[1]); - N2 = N1; - N3 = N1; - } - else if (dim_array[1] == 3) { - DimAr = (double*) mxGetData(prhs[1]); - N1 = (long)(DimAr[0]); - N2 = (long)(DimAr[1]); - N3 = (long)(DimAr[2]); - } - else {mexErrMsgTxt("DIM must be scalar [N] or a vector with 3 elements [N1, N2, N3]");} - - ModelSelected = (int) mxGetScalar(prhs[0]); /* selected model from Phantom3DLibrary */ - - int Model=0, Components=0, steps = 0, counter=0, ii; - float C0 = 0.0f, x0 = 0.0f, y0 = 0.0f, z0 = 0.0f, a = 0.0f, b = 0.0f, c = 0.0f, psi_gr1 = 0.0f, psi_gr2 = 0.0f, psi_gr3 = 0.0f; - - char *filename; - FILE * fp; - if(!mxIsChar(prhs[2]) ) { - mexErrMsgTxt("Need filename absolute path string input."); - } - filename = mxArrayToString(prhs[2]); - fp= fopen(filename, "rb"); - mxFree(filename); - if( fp == NULL ) { - mexErrMsgTxt("Cannot open the model library file (Phantom3DLibrary.dat)"); - } - else { - - char str[MAXCHAR]; - char tmpstr1[16]; - char tmpstr2[22]; - char tmpstr3[16]; - char tmpstr4[16]; - char tmpstr5[16]; - char tmpstr6[16]; - char tmpstr7[16]; - char tmpstr8[16]; - char tmpstr9[16]; - char tmpstr10[16]; - char tmpstr11[16]; - char tmpstr12[16]; - - while (fgets(str, MAXCHAR, fp) != NULL) - { - /* work with non-# commented lines */ - if(str[0] != '#') { - sscanf(str, "%15s : %21[^;];", tmpstr1, tmpstr2); - if (strcmp(tmpstr1,"Model")==0) - { - Model = atoi(tmpstr2); - if ((ModelSelected == Model) && (counter == 0)) { - /* check if we have a right model */ - if (fgets(str, MAXCHAR, fp) != NULL) sscanf(str, "%15s : %21[^;];", tmpstr1, tmpstr2); - else { - mexErrMsgTxt("Unexpected the end of the line (Components) in parameters file"); - break; } - if (strcmp(tmpstr1,"Components") == 0) Components = atoi(tmpstr2); - printf("%s %i\n", "Components:", Components); - if (Components <= 0) { - printf("%s %i\n", "Components cannot be negative, the given value is", Components); - mexErrMsgTxt("Components cannot be negative"); - break; } - if (fgets(str, MAXCHAR, fp) != NULL) sscanf(str, "%15s : %21[^;];", tmpstr1, tmpstr2); - else { - mexErrMsgTxt("Unexpected the end of the line (TimeSteps) in parameters file"); - break; } - if (strcmp(tmpstr1,"TimeSteps") == 0) steps = atoi(tmpstr2); - if (steps <= 0) { - printf("%s %i\n", "TimeSteps cannot be negative, the given value is", steps); - mexErrMsgTxt("TimeSteps cannot be negative"); - break; } - printf("%s %i\n", "TimeSteps:", steps); - if (steps == 1) { - /**************************************************/ - printf("\n %s %i %s \n", "Stationary 3D model", ModelSelected, " is selected"); - - const mwSize N_dims[3] = {N1, N2, N3}; /* image dimensions */ - A = (float*)mxGetPr(plhs[0] = mxCreateNumericArray(3, N_dims, mxSINGLE_CLASS, mxREAL)); /*output array*/ - - /* loop over all components */ - for(ii=0; ii 1)) { - printf("%s %f\n", "x0 (object position) must be in [-1,1] range, the given value is", x0); - mexErrMsgTxt("x0 (object position) must be in [-1,1] range"); - break; } - if ((y0 < -1) || (y0 > 1)) { - printf("%s %f\n", "y0 (object position) must be in [-1,1] range, the given value is", y0); - mexErrMsgTxt("y0 (object position) must be in [-1,1] range"); - break; } - if ((z0 < -1) || (z0 > 1)) { - printf("%s %f\n", "z0 (object position) must be in [-1,1] range, the given value is", z0); - mexErrMsgTxt("z0 (object position) must be in [-1,1] range"); - break; } - if ((a <= 0) || (a > 2)) { - printf("%s %f\n", "a (object size) must be positive in [0,2] range, the given value is", a); - mexErrMsgTxt("a (object size) must be positive in [0,2] range"); - break; } - if ((b <= 0) || (b > 2)) { - printf("%s %f\n", "b (object size) must be positive in [0,2] range, the given value is", b); - mexErrMsgTxt("b (object size) must be positive in [0,2] range"); - break; } - if ((c <= 0) || (c > 2)) { - printf("%s %f\n", "c (object size) must be positive in [0,2] range, the given value is", c); - mexErrMsgTxt("c (object size) must be positive in [0,2] range"); - break; } - printf("\nObject : %s \nC0 : %f \nx0 : %f \ny0 : %f \nz0 : %f \na : %f \nb : %f \nc : %f \n", tmpstr2, C0, x0, y0, z0, a, b, c); - - TomoP3DObject_core(A, N1, N2, N3, 0l, N3, tmpstr2, C0, x0, y0, z0, b, a, c, -psi_gr1, psi_gr2, psi_gr3, 0l); /* Matlab */ - } - } - else { - /**************************************************/ - printf("\n %s %i %s \n", "Temporal 3D+time model", ModelSelected, " is selected"); - /* temporal phantom 3D + time (4D) */ - const mwSize N_dims[4] = {N1, N2, N3, steps}; /* image dimensions */ - A = (float*)mxGetPr(plhs[0] = mxCreateNumericArray(4, N_dims, mxSINGLE_CLASS, mxREAL)); - - float C1 = 0.0f, x1 = 0.0f, y1 = 0.0f, z1 = 0.0f, a1 = 0.0f, b1 = 0.0f, c1 = 0.0f, psi_gr1_1 = 0.0f, psi_gr2_1 = 0.0f, psi_gr3_1 = 0.0f; - /* loop over all components */ - for(ii=0; ii 1)) { - printf("%s %f\n", "x0 (object position) must be in [-1,1] range, the given value is", x0); - mexErrMsgTxt("x0 (object position) must be in [-1,1] range"); - break; } - if ((y0 < -1) || (y0 > 1)) { - printf("%s %f\n", "y0 (object position) must be in [-1,1] range, the given value is", y0); - mexErrMsgTxt("y0 (object position) must be in [-1,1] range"); - break; } - if ((z0 < -1) || (z0 > 1)) { - printf("%s %f\n", "z0 (object position) must be in [-1,1] range, the given value is", z0); - mexErrMsgTxt("z0 (object position) must be in [-1,1] range"); - break; } - if ((a <= 0) || (a > 2)) { - printf("%s %f\n", "a (object size) must be positive in [0,2] range, the given value is", a); - mexErrMsgTxt("a (object size) must be positive in [0,2] range"); - break; } - if ((b <= 0) || (b > 2)) { - printf("%s %f\n", "b (object size) must be positive in [0,2] range, the given value is", b); - mexErrMsgTxt("b (object size) must be positive in [0,2] range"); - break; } - if ((c <= 0) || (c > 2)) { - printf("%s %f\n", "c (object size) must be positive in [0,2] range, the given value is", c); - mexErrMsgTxt("c (object size) must be positive in [0,2] range"); - break; } - // printf("\nObject : %s \nC0 : %f \nx0 : %f \ny0 : %f \nz0 : %f \na : %f \nb : %f \n", tmpstr2, C0, x0, y0, z0, a, b, c); - - /* check Endvar relatedparameters */ - if (fgets(str, MAXCHAR, fp) != NULL) sscanf(str, "%15s : %15s %15s %15s %15s %15s %15s %15s %15s %15s %15[^;];", tmpstr1, tmpstr3, tmpstr4, tmpstr5, tmpstr6, tmpstr7, tmpstr8, tmpstr9, tmpstr10, tmpstr11, tmpstr12); - else { - mexErrMsgTxt("Unexpected the end of the line (Endvar loop) in parameters file"); - break; } - - if (strcmp(tmpstr1,"Endvar") == 0) { - C1 = (float)atof(tmpstr3); /* intensity */ - x1 = (float)atof(tmpstr4); /* x0 position */ - y1 = (float)atof(tmpstr5); /* y0 position */ - z1 = (float)atof(tmpstr6); /* z0 position */ - a1 = (float)atof(tmpstr7); /* a - size object */ - b1 = (float)atof(tmpstr8); /* b - size object */ - c1 = (float)atof(tmpstr9); /* c - size object */ - psi_gr1_1 = (float)atof(tmpstr10); /* rotation angle 1*/ - psi_gr2_1 = (float)atof(tmpstr11); /* rotation angle 2*/ - psi_gr3_1 = (float)atof(tmpstr12); /* rotation angle 3*/ - } - else { - printf("%s\n", "Cannot find 'Endvar' string in parameters file"); - break; } - - if (C1 == 0) { - printf("%s %f\n", "Endvar C1 should not be equal to zero, the given value is", C1); - mexErrMsgTxt("Endvar C0 should not be equal to zero"); - break; } - if ((x1 < -1) || (x1 > 1)) { - printf("%s %f\n", "Endvar x1 (object position) must be in [-1,1] range, the given value is", x1); - mexErrMsgTxt("Endvar x0 (object position) must be in [-1,1] range"); - break; } - if ((y1 < -1) || (y1 > 1)) { - printf("%s %f\n", "Endvar y1 (object position) must be in [-1,1] range, the given value is", y1); - mexErrMsgTxt("Endvar y0 (object position) must be in [-1,1] range"); - break; } - if ((z1 < -1) || (z1 > 1)) { - printf("%s %f\n", "Endvar z1 (object position) must be in [-1,1] range, the given value is", z1); - mexErrMsgTxt("Endvar z0 (object position) must be in [-1,1] range"); - break; } - if ((a1 <= 0) || (a1 > 2)) { - printf("%s %f\n", "Endvar a1 (object size) must be positive in [0,2] range, the given value is", a1); - mexErrMsgTxt("Endvar a (object size) must be positive in [0,2] range"); - break; } - if ((b1 <= 0) || (b1 > 2)) { - printf("%s %f\n", "Endvar b1 (object size) must be positive in [0,2] range, the given value is", b1); - mexErrMsgTxt("Endvar b (object size) must be positive in [0,2] range"); - break; } - if ((c1 <= 0) || (c1 > 2)) { - printf("%s %f\n", "Endvar c1 (object size) must be positive in [0,2] range, the given value is", c1); - mexErrMsgTxt("Endvar c (object size) must be positive in [0,2] range"); - break; } - //printf("\nObject : %s \nC0 : %f \nx0 : %f \ny0 : %f \nz0 : %f \na : %f \nb : %f \nc : %f \n", tmpstr2, C0, x0, y0, z0, a1, b1, c1); - - /*now we know the initial parameters of the object and the final ones. We linearly extrapolate to establish steps and coordinates. */ - /* calculating the full distance berween the start and the end points */ - float distance = sqrtf(pow((x1 - x0),2) + pow((y1 - y0),2) + pow((z1 - z0),2)); - float d_dist = distance/(steps-1); /*a step over line */ - float C_step = (C1 - C0)/(steps-1); - float a_step = (a1 - a)/(steps-1); - float b_step = (b1 - b)/(steps-1); - float c_step = (c1 - c)/(steps-1); - float phi_rot_step1 = (psi_gr1_1 - psi_gr1)/(steps-1); - float phi_rot_step2 = (psi_gr2_1 - psi_gr2)/(steps-1); - float phi_rot_step3 = (psi_gr3_1 - psi_gr3)/(steps-1); - - int tt; - float x_t, y_t, z_t, a_t, b_t, c_t, C_t, phi1_t, phi2_t, phi3_t, d_step; - /* initialize */ - x_t = x0; y_t = y0; z_t = z0; a_t = a; b_t = b; c_t = c; C_t = C0; phi1_t = psi_gr1; phi2_t = psi_gr2; phi3_t = psi_gr3; d_step = d_dist; - - /*loop over time frames*/ - for(tt=0; tt < steps; tt++) { - - TomoP3DObject_core(A, N1, N2, N3, 0l, N3, tmpstr2, C_t, x_t, y_t, z_t, b_t, a_t, c_t, -phi1_t, phi2_t, phi3_t, tt); /* Matlab */ - - /* calculating new coordinates of an object */ - if (distance != 0.0f) { - float t = d_step/distance; - x_t = (1-t)*x0 + t*x1; - y_t = (1-t)*y0 + t*y1; - z_t = (1-t)*z0 + t*z1;} - else { - x_t = x0; - y_t = y0; - z_t = z0; } - - d_step += d_dist; - a_t += a_step; - b_t += b_step; - c_t += c_step; - C_t += C_step; - phi1_t += phi_rot_step1; - phi2_t += phi_rot_step2; - phi3_t += phi_rot_step3; - } /*time steps*/ - - } /*components loop*/ - } - counter++; - } - } - } - } - } - fclose(fp); - if (counter == 0) { - printf("%s %i %s \n", "Model no. ", ModelSelected, "is not found!"); - mexErrMsgTxt("No object found, check models file"); - } -} -// if (params_switch[0] == 0) mexErrMsgTxt("Parameters file (Phantom3DLibrary.dat) cannot be read OR some major syntax errors in it; Check your PATH to the file is a valid one and the file exists"); -// if (params_switch[1] == 0) { -// printf("%s %i\n", ">>>>> No such model available in the library, the given model is", ModelSelected); -// mexErrMsgTxt("Check the syntax in Phantom3DLibrary.dat"); -// } -// if (params_switch[2] == 0) { -// printf("%s %i\n", ">>>>> Components number cannot be negative for the given model", ModelSelected); -// mexErrMsgTxt("Check the syntax in Phantom3DLibrary.dat"); -// } -// if (params_switch[3] == 0) { -// printf("%s %i\n", ">>>>> TimeSteps value cannot be negative for the given model", ModelSelected); -// mexErrMsgTxt("Check the syntax in Phantom3DLibrary.dat"); -// } -// if (params_switch[4] == 0) { -// printf("%s %i\n", ">>>>> Unknown name of the object for the given model", ModelSelected); -// mexErrMsgTxt("Check the syntax in Phantom3DLibrary.dat"); -// } -// if (params_switch[5] == 0) { -// printf("%s %i\n", ">>>>> C0 should not be equal to zero for the given model", ModelSelected); -// mexErrMsgTxt("Check the syntax in Phantom3DLibrary.dat"); -// } -// if (params_switch[6] == 0) { -// printf("%s %i\n", ">>>>> x0 (object position) must be in [-1,1] range for the given model", ModelSelected); -// mexErrMsgTxt("Check the syntax in Phantom3DLibrary.dat"); -// } -// if (params_switch[7] == 0) { -// printf("%s %i\n", ">>>>> y0 (object position) must be in [-1,1] range for the given model", ModelSelected); -// mexErrMsgTxt("Check the syntax in Phantom3DLibrary.dat"); -// } -// if (params_switch[8] == 0) { -// printf("%s %i\n", ">>>>> z0 (object position) must be in [-1,1] range for the given model", ModelSelected); -// mexErrMsgTxt("Check the syntax in Phantom3DLibrary.dat"); -// } -// if (params_switch[9] == 0) { -// printf("%s %i\n", ">>>>> a (object size) must be positive in [0,2] range for the given model", ModelSelected); -// mexErrMsgTxt("Check the syntax in Phantom3DLibrary.dat"); -// } -// if (params_switch[10] == 0) { -// printf("%s %i\n", ">>>>> b (object size) must be positive in [0,2] range for the given model", ModelSelected); -// mexErrMsgTxt("Check the syntax in Phantom3DLibrary.dat"); -// } -// if (params_switch[11] == 0) { -// printf("%s %i\n", ">>>>> c (object size) must be positive in [0,2] range for the given model", ModelSelected); -// mexErrMsgTxt("Check the syntax in Phantom3DLibrary.dat"); -// } diff --git a/matlab/mex_wrappers/TomoP3DObject.c b/matlab/mex_wrappers/TomoP3DObject.c deleted file mode 100644 index 9f0f9b0..0000000 --- a/matlab/mex_wrappers/TomoP3DObject.c +++ /dev/null @@ -1,124 +0,0 @@ -/* - * Copyright 2017 Daniil Kazantsev - * - * Licensed under the Apache License, Version 2.0 (the "License"); - * you may not use this file except in compliance with the License. - * You may obtain a copy of the License at - * http://www.apache.org/licenses/LICENSE-2.0 - * Unless required by applicable law or agreed to in writing, software - * distributed under the License is distributed on an "AS IS" BASIS, - * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. - * See the License for the specific language governing permissions and - * limitations under the License. - */ - -#include "mex.h" -#include -#include -#include -#include -#include "omp.h" - -#include "TomoP3DModel_core.h" -#include "utils.h" - -#define M_PI 3.14159265358979323846 - -/* Function to build 3D objects without using the library file - * (MATLAB mex-wrapper) - * - * Input Parameters: - * 1. parameters for 3D object - * 2. DIM volume dimensions [N1,N2,N3] in voxels (N1 x N2 x N3) - * - * VolumeSize in voxels (N x N x N) - * - * Output: - * 1. The analytical phantom size of [N1 x N2 x N3] - */ -#define MAXCHAR 1000 - -void mexFunction( - int nlhs, mxArray *plhs[], - int nrhs, const mxArray *prhs[]) - -{ - float *A; - char *tmpstr2; - float C0 = 0.0f, x0 = 0.0f, y0 = 0.0f, z0 = 0.0f, a = 0.0f, b = 0.0f, c = 0.0f, psi_gr1 = 0.0f, psi_gr2 = 0.0f, psi_gr3 = 0.0f; - double *DimAr; - mwSize N1, N2, N3; - const mwSize *dim_array; - - tmpstr2 = mxArrayToString(prhs[0]); /* name of the object */ - C0 = (float) mxGetScalar(prhs[1]); /* intensity */ - x0 = (float) mxGetScalar(prhs[2]); /* x0 position */ - y0 = (float) mxGetScalar(prhs[3]); /* y0 position */ - z0 = (float) mxGetScalar(prhs[4]); /* z0 position */ - a = (float) mxGetScalar(prhs[5]); /* a - size object */ - b = (float) mxGetScalar(prhs[6]); /* b - size object */ - c = (float) mxGetScalar(prhs[7]); /* c - size object */ - psi_gr1 = (float) mxGetScalar(prhs[8]); /* rotation angle 1*/ - psi_gr2 = (float) mxGetScalar(prhs[9]); /* rotation angle 2*/ - psi_gr3 = (float) mxGetScalar(prhs[10]); /* rotation angle 3*/ - - dim_array = mxGetDimensions(prhs[11]); /* get dimensions of DIM */ - - N1 = 0; N2 = 0; N3 = 0; - if (dim_array[1] == 1) { - N1 = (long) mxGetScalar(prhs[11]); - N2 = N1; - N3 = N1; - } - else if (dim_array[1] == 3) { - DimAr = (double*) mxGetData(prhs[11]); - N1 = (long)(DimAr[0]); - N2 = (long)(DimAr[1]); - N3 = (long)(DimAr[2]); - } - else {mexErrMsgTxt("DIM must be scalar [N] or a vector with 3 elements [N1, N2, N3]");} - - /*Handling Matlab input data*/ - if (nrhs != 12) mexErrMsgTxt("Input of 12 parameters is required"); - - const mwSize N_dims[3] = {N1, N2, N3}; /* image dimensions */ - A = (float*)mxGetPr(plhs[0] = mxCreateNumericArray(3, N_dims, mxSINGLE_CLASS, mxREAL)); /*output array*/ - - - if ((strcmp("gaussian",tmpstr2) != 0) && (strcmp("paraboloid",tmpstr2) != 0) && (strcmp("ellipsoid",tmpstr2) != 0) && (strcmp("cone",tmpstr2) != 0) && (strcmp("cuboid",tmpstr2) != 0) && (strcmp("ellipticalcylinder",tmpstr2) != 0) ) { - printf("%s %s\n", "Unknown name of the object, the given name is", tmpstr2); - mexErrMsgTxt("Unknown name of the object"); - } - if (C0 == 0) { - printf("%s %f\n", "C0 should not be equal to zero, the given value is", C0); - mexErrMsgTxt("C0 should not be equal to zero"); - } - if ((x0 < -1) || (x0 > 1)) { - printf("%s %f\n", "x0 (object position) must be in [-1,1] range, the given value is", x0); - mexErrMsgTxt("x0 (object position) must be in [-1,1] range"); - } - if ((y0 < -1) || (y0 > 1)) { - printf("%s %f\n", "y0 (object position) must be in [-1,1] range, the given value is", y0); - mexErrMsgTxt("y0 (object position) must be in [-1,1] range"); - } - if ((z0 < -1) || (z0 > 1)) { - printf("%s %f\n", "z0 (object position) must be in [-1,1] range, the given value is", z0); - mexErrMsgTxt("z0 (object position) must be in [-1,1] range"); - } - if ((a <= 0) || (a > 2)) { - printf("%s %f\n", "a (object size) must be positive in [0,2] range, the given value is", a); - mexErrMsgTxt("a (object size) must be positive in [0,2] range"); - } - if ((b <= 0) || (b > 2)) { - printf("%s %f\n", "b (object size) must be positive in [0,2] range, the given value is", b); - mexErrMsgTxt("b (object size) must be positive in [0,2] range"); - } - if ((c <= 0) || (c > 2)) { - printf("%s %f\n", "c (object size) must be positive in [0,2] range, the given value is", c); - mexErrMsgTxt("c (object size) must be positive in [0,2] range"); - } - printf("\nObject : %s \nC0 : %f \nx0 : %f \ny0 : %f \nz0 : %f \na : %f \nb : %f \nc : %f \n", tmpstr2, C0, x0, y0, z0, a, b, c); - - TomoP3DObject_core(A, N1, N2, N3, 0l, N3, tmpstr2, C0, x0, y0, z0, b, a, c, -psi_gr1, psi_gr2, psi_gr3, 0); /* Matlab */ - mxFree(tmpstr2); -} 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); - From 969940fd75bb7869fcff1948368021431cf63be8 Mon Sep 17 00:00:00 2001 From: dkazanc Date: Thu, 4 Oct 2018 15:57:25 +0100 Subject: [PATCH 4/5] demos fixed --- Wrappers/Python/setup.py.in | 6 +++--- python/Demos/DemoModel.py | 6 ++++-- python/Demos/DemoModel2.py | 2 +- python/Demos/DemoModel3D.py | 4 ++-- python/Demos/DemoModel4D.py | 4 ++-- python/Demos/DemoReconASTRA.py | 2 +- python/Demos/DemoReconTomoPy.py | 2 +- 7 files changed, 14 insertions(+), 12 deletions(-) diff --git a/Wrappers/Python/setup.py.in b/Wrappers/Python/setup.py.in index 78c6015..93f2e1d 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 = [] @@ -55,7 +55,7 @@ setup( extra_link_args = extra_link_args)]), zip_safe = False, include_package_data=True, - package_data={'tomophantom':['*.dat','../functions/models/*.dat']}, + package_data={'tomophantom':['*.dat','../PhantomLibrary/models/*.dat']}, packages = {'tomophantom'} ) @@ -83,6 +83,6 @@ setup( gdb_debug=True), zip_safe = False, include_package_data=True, - package_data={'tomophantom':['*.dat','../functions/models/*.dat']}, + package_data={'tomophantom':['*.dat','../PhantomLibrary/models/*.dat']}, packages = {'tomophantom'} ) diff --git a/python/Demos/DemoModel.py b/python/Demos/DemoModel.py index f2110cc..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 = '../../../PhantomLibrary/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 29e387e..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 = '../../../PhantomLibrary/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 cff302b..07c3d4d 100644 --- a/python/Demos/DemoModel3D.py +++ b/python/Demos/DemoModel3D.py @@ -21,7 +21,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 = '../../../PhantomLibrary/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 = '../../../PhantomLibrary/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 8dedabe..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 = '../../../PhantomLibrary/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 = '../../../PhantomLibrary/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/DemoReconASTRA.py b/python/Demos/DemoReconASTRA.py index 2de570e..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 = '../../../PhantomLibrary/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 66c2268..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 = '../../../PhantomLibrary/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) From 9c136aed567d76c261d25d8bf3b413012c85c0ac Mon Sep 17 00:00:00 2001 From: dkazanc Date: Thu, 4 Oct 2018 16:40:41 +0100 Subject: [PATCH 5/5] minor matlab demo bits --- Wrappers/MATLAB/Model2DGeneratorDemo.m | 4 +--- Wrappers/MATLAB/Reconstruction_Demo.m | 2 -- run.sh | 3 ++- 3 files changed, 3 insertions(+), 6 deletions(-) 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/Reconstruction_Demo.m b/Wrappers/MATLAB/Reconstruction_Demo.m index 506cc72..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 diff --git a/run.sh b/run.sh index 3624419..1105a16 100644 --- a/run.sh +++ b/run.sh @@ -1,8 +1,9 @@ #!/bin/bash echo "Building Tomophantom using CMake" #rm -r build +# pip install cython mkdir build -cd build +cd build/ #make clean cmake ../ -DBUILD_PYTHON_WRAPPER=ON -DCMAKE_BUILD_TYPE=Release -DCMAKE_INSTALL_PREFIX=./install make install