Merge branch 'master' of https://github.com/vais-ral/CCPi-FISTA_Reconstruction into pythonize

3 files changed, 115 insertions, 174 deletions

diff --git a/demos/Demo1.m b/demos/Demo1.m
deleted file mode 100644
index 15e2e5b..0000000
--- a/demos/Demo1.m
+++ /dev/null
@@ -1,174 +0,0 @@
-% Demonstration of tomographic reconstruction from noisy and corrupted by 

-% artifacts undersampled projection data using Students't penalty 

-% Optimisation problem is solved using FISTA algorithm (see Beck & Teboulle)

-

-% see Readme file for instructions

-%%

-% compile MEX-files ones

-% cd ..

-% cd main_func

-% compile_mex

-% cd .. 

-% cd demos

-%%

-

-close all;clc;clear all;

-% adding paths

-addpath('../data/'); 

-addpath('../main_func/'); addpath('../main_func/regularizers_CPU/');

-addpath('../supp/'); 

-

-load phantom_bone512.mat % load the phantom

-load my_red_yellowMAP.mat % load the colormap

-% load sino1.mat; % load noisy sinogram

-

-N = 512; %  the size of the tomographic image NxN

-theta = 1:1:180; % acquisition angles (in parallel beam from 0 to Pi)

-theta_rad = theta*(pi/180); % conversion to radians

-P = 2*ceil(N/sqrt(2))+1; % the size of the detector array

-ROI = find(phantom > 0);

-

-% using ASTRA to set the projection geometry

-% potentially parallel geometry can be replaced with a divergent one

-Z_slices = 1;

-det_row_count = Z_slices;

-proj_geom = astra_create_proj_geom('parallel3d', 1, 1, det_row_count, P, theta_rad);

-vol_geom = astra_create_vol_geom(N,N,Z_slices);

-

-zing_rings_add;    % generating data, adding zingers and stripes

-%% 

-fprintf('%s\n', 'Direct reconstruction using FBP...');

-FBP_1 = iradon(sino_zing_rings', theta, N); 

-

-fprintf('%s %.4f\n', 'RMSE for FBP reconstruction:', RMSE(FBP_1(:), phantom(:)));

-

-figure(1); 

-subplot_tight(1,2,1, [0.05 0.05]); imshow(FBP_1,[0 0.6]);  title('FBP reconstruction of noisy and corrupted by artifacts sinogram'); colorbar;

-subplot_tight(1,2,2, [0.05 0.05]); imshow((phantom - FBP_1).^2,[0 0.1]);  title('residual: (ideal phantom - FBP)^2'); colorbar;

-colormap(cmapnew); 

-

-%%

-fprintf('%s\n', 'Reconstruction using FISTA-PWLS without regularization...');

-clear params

-% define parameters

-params.proj_geom = proj_geom; % pass geometry to the function

-params.vol_geom = vol_geom;

-params.sino = sino_zing_rings; % sinogram

-params.iterFISTA = 45; %max number of outer iterations

-params.X_ideal = phantom; % ideal phantom

-params.ROI = ROI; % phantom region-of-interest

-params.show = 1; % visualize reconstruction on each iteration

-params.slice = 1; params.maxvalplot = 0.6; 

-params.weights = Dweights; % statistical weighting 

-tic; [X_FISTA, output] = FISTA_REC(params); toc; 

-

-fprintf('%s %.4f\n', 'Min RMSE for FISTA-PWLS reconstruction is:', min(error_FISTA(:)));

-error_FISTA = output.Resid_error; obj_FISTA = output.objective;

-

-figure(2); clf

-%set(gcf, 'Position', get(0,'Screensize'));

-subplot(1,2,1, [0.05 0.05]); imshow(X_FISTA,[0 0.6]); title('FISTA-PWLS reconstruction'); colorbar;

-subplot(1,2,2, [0.05 0.05]); imshow((phantom - X_FISTA).^2,[0 0.1]);  title('residual'); colorbar;

-colormap(cmapnew); 

-figure(3); clf

-subplot(1,2,1, [0.05 0.05]); plot(error_FISTA);  title('RMSE plot'); colorbar;

-subplot(1,2,2, [0.05 0.05]); plot(obj_FISTA);  title('Objective plot'); colorbar;

-colormap(cmapnew); 

-%%

-fprintf('%s\n', 'Reconstruction using FISTA-PWLS-TV...');

-clear params

-% define parameters

-params.proj_geom = proj_geom; % pass geometry to the function

-params.vol_geom = vol_geom;

-params.sino = sino_zing_rings;

-params.iterFISTA = 45; % max number of outer iterations

-params.Regul_LambdaTV = 0.0015; % regularization parameter for TV problem

-params.X_ideal = phantom; % ideal phantom

-params.ROI = ROI; % phantom region-of-interest

-params.weights = Dweights; % statistical weighting 

-params.show = 1; % visualize reconstruction on each iteration

-params.slice = 1; params.maxvalplot = 0.6; 

-tic; [X_FISTA_TV, output] = FISTA_REC(params); toc; 

-

-fprintf('%s %.4f\n', 'Min RMSE for FISTA-PWLS-TV reconstruction is:', min(error_FISTA_TV(:)));

-error_FISTA_TV = output.Resid_error; obj_FISTA_TV = output.objective;

-

-figure(4); clf

-subplot(1,2,1, [0.05 0.05]); imshow(X_FISTA_TV,[0 0.6]); title('FISTA-PWLS-TV reconstruction'); colorbar;

-subplot(1,2,2, [0.05 0.05]); imshow((phantom - X_FISTA_TV).^2,[0 0.1]);  title('residual'); colorbar;

-colormap(cmapnew); 

-figure(5); clf

-subplot(1,2,1, [0.05 0.05]); plot(error_FISTA_TV);  title('RMSE plot'); colorbar;

-subplot(1,2,2, [0.05 0.05]); plot(obj_FISTA_TV);  title('Objective plot'); colorbar;

-colormap(cmapnew); 

-%%

-fprintf('%s\n', 'Reconstruction using FISTA-GH-TV...');

-clear params

-% define parameters

-params.proj_geom = proj_geom; % pass geometry to the function

-params.vol_geom = vol_geom;

-params.sino = sino_zing_rings;

-params.iterFISTA = 50; % max number of outer iterations

-params.Regul_LambdaTV = 0.0015;  % regularization parameter for TV problem

-params.X_ideal = phantom; % ideal phantom

-params.ROI = ROI; % phantom region-of-interest

-params.weights = Dweights; % statistical weighting 

-params.Ring_LambdaR_L1 = 0.002; % parameter to sparsify the "rings vector"

-params.Ring_Alpha = 20; % to accelerate ring-removal procedure

-params.show = 0; % visualize reconstruction on each iteration

-params.slice = 1; params.maxvalplot = 0.6; 

-tic; [X_FISTA_GH_TV, output] = FISTA_REC(params); toc; 

-

-fprintf('%s %.4f\n', 'Min RMSE for FISTA-GH-TV reconstruction is:', min(error_FISTA_GH_TV(:)));

-error_FISTA_GH_TV = output.Resid_error; obj_FISTA_GH_TV = output.objective;

-

-figure(6); clf

-subplot(1,2,1, [0.05 0.05]); imshow(X_FISTA_GH_TV,[0 0.6]); title('FISTA-GH-TV reconstruction'); colorbar;

-subplot(1,2,2, [0.05 0.05]);imshow((phantom - X_FISTA_GH_TV).^2,[0 0.1]);  title('residual'); colorbar;

-colormap(cmapnew); 

-

-figure(7); clf

-subplot(1,2,1, [0.05 0.05]);  plot(error_FISTA_GH_TV);  title('RMSE plot'); colorbar;

-subplot(1,2,2, [0.05 0.05]);  plot(obj_FISTA_GH_TV);  title('Objective plot'); colorbar;

-colormap(cmapnew); 

-%%

-fprintf('%s\n', 'Reconstruction using FISTA-Student-TV...');

-clear params

-% define parameters

-params.proj_geom = proj_geom; % pass geometry to the function

-params.vol_geom = vol_geom;

-params.sino = sino_zing_rings;

-params.iterFISTA = 55; % max number of outer iterations

-params.L_const = 0.1; % Lipshitz constant (can be chosen manually to accelerate convergence)

-params.Regul_LambdaTV = 0.00152;  % regularization parameter for TV problem

-params.X_ideal = phantom; % ideal phantom

-params.ROI = ROI; % phantom region-of-interest

-params.weights = Dweights; % statistical weighting 

-params.fidelity = 'student'; % selecting students t fidelity

-params.show = 1; % visualize reconstruction on each iteration

-params.slice = 1; params.maxvalplot = 0.6; 

-params.initilize = 1; % warm start with SIRT

-tic; [X_FISTA_student_TV, output] = FISTA_REC(params); toc; 

-

-fprintf('%s %.4f\n', 'Min RMSE for FISTA-Student-TV reconstruction is:', min(error_FISTA_student_TV(:)));

-error_FISTA_student_TV = output.Resid_error; obj_FISTA_student_TV = output.objective;

-

-figure(8); 

-set(gcf, 'Position', get(0,'Screensize'));

-subplot(1,2,1, [0.05 0.05]); imshow(X_FISTA_student_TV,[0 0.6]); title('FISTA-Student-TV reconstruction'); colorbar;

-subplot(1,2,2, [0.05 0.05]); imshow((phantom - X_FISTA_student_TV).^2,[0 0.1]);  title('residual'); colorbar;

-colormap(cmapnew); 

-

-figure(9); 

-subplot(1,2,1, [0.05 0.05]); plot(error_FISTA_student_TV);  title('RMSE plot'); colorbar;

-subplot(1,2,2, [0.05 0.05]); plot(obj_FISTA_student_TV);  title('Objective plot'); colorbar;

-colormap(cmapnew); 

-%%

-% print all RMSE's

-fprintf('%s\n', '--------------------------------------------');

-fprintf('%s %.4f\n', 'RMSE for FBP reconstruction:', RMSE(FBP_1(:), phantom(:)));

-fprintf('%s %.4f\n', 'Min RMSE for FISTA-PWLS reconstruction:', min(error_FISTA(:)));

-fprintf('%s %.4f\n', 'Min RMSE for FISTA-PWLS-TV reconstruction:', min(error_FISTA_TV(:)));

-fprintf('%s %.4f\n', 'Min RMSE for FISTA-GH-TV reconstruction:', min(error_FISTA_GH_TV(:)));

-fprintf('%s %.4f\n', 'Min RMSE for FISTA-Student-TV reconstruction:', min(error_FISTA_student_TV(:)));

-%
\ No newline at end of file
diff --git a/demos/Demo_Phantom3D_Cone.m b/demos/Demo_Phantom3D_Cone.m
new file mode 100644
index 0000000..6419386
--- /dev/null
+++ b/demos/Demo_Phantom3D_Cone.m
@@ -0,0 +1,66 @@
+% A demo script to reconstruct 3D synthetic data using FISTA method for
+% CONE BEAM geometry
+% requirements: ASTRA-toolbox and TomoPhantom toolbox
+
+close all;clc;clear all;
+% adding paths
+addpath('../data/');
+addpath('../main_func/'); addpath('../main_func/regularizers_CPU/'); addpath('../main_func/regularizers_GPU/NL_Regul/'); addpath('../main_func/regularizers_GPU/Diffus_HO/');
+addpath('../supp/');
+
+
+%%
+% build 3D phantom using TomoPhantom 
+modelNo = 3; % see Phantom3DLibrary.dat file in TomoPhantom
+N = 256; % x-y-z size (cubic image)
+angles = 0:1.5:360; % angles vector in degrees
+angles_rad = angles*(pi/180); % conversion to radians
+det_size = round(sqrt(2)*N); % detector size
+% in order to run functions you have to go to the directory:
+cd /home/algol/Documents/MATLAB/TomoPhantom/functions/
+TomoPhantom = buildPhantom3D(modelNo,N); % generate 3D phantom
+%%
+% using ASTRA-toolbox to set the projection geometry (cone beam)
+% eg: astra.create_proj_geom('cone', 1.0 (resol), 1.0 (resol), detectorRowCount, detectorColCount, angles, originToSource, originToDetector)
+vol_geom = astra_create_vol_geom(N,N,N);
+proj_geom = astra_create_proj_geom('cone', 1.0, 1.0, N, det_size, angles_rad, 2000, 2160);
+%% 
+% do forward projection using ASTRA
+% inverse crime data generation
+[sino_id, SinoCone3D] = astra_create_sino3d_cuda(TomoPhantom, proj_geom, vol_geom);
+astra_mex_data3d('delete', sino_id);
+%%
+fprintf('%s\n', 'Reconstructing with CGLS using ASTRA-toolbox ...');
+vol_id = astra_mex_data3d('create', '-vol', vol_geom, 0);
+proj_id = astra_mex_data3d('create', '-proj3d', proj_geom, SinoCone3D); 
+cfg = astra_struct('CGLS3D_CUDA');
+cfg.ProjectionDataId = proj_id;
+cfg.ReconstructionDataId = vol_id;
+cfg.option.MinConstraint = 0;
+alg_id = astra_mex_algorithm('create', cfg);
+astra_mex_algorithm('iterate', alg_id, 15);
+reconASTRA_3D = astra_mex_data3d('get', vol_id);
+%%
+fprintf('%s\n', 'Reconstruction using FISTA-LS without regularization...');
+clear params
+% define parameters
+params.proj_geom = proj_geom; % pass geometry to the function
+params.vol_geom = vol_geom;
+params.sino = single(SinoCone3D); % sinogram
+params.iterFISTA = 30; %max number of outer iterations
+params.X_ideal = TomoPhantom; % ideal phantom
+params.show = 1; % visualize reconstruction on each iteration
+params.slice = round(N/2); params.maxvalplot = 1; 
+tic; [X_FISTA, output] = FISTA_REC(params); toc; 
+
+error_FISTA = output.Resid_error; obj_FISTA = output.objective;
+fprintf('%s %.4f\n', 'Min RMSE for FISTA-LS reconstruction is:', min(error_FISTA(:)));
+
+Resid3D = (TomoPhantom - X_FISTA).^2;
+figure(2);
+subplot(1,2,1); imshow(X_FISTA(:,:,params.slice),[0 params.maxvalplot]); title('FISTA-LS reconstruction'); colorbar;
+subplot(1,2,2); imshow(Resid3D(:,:,params.slice),[0 0.1]);  title('residual'); colorbar;
+figure(3);
+subplot(1,2,1); plot(error_FISTA);  title('RMSE plot'); colorbar;
+subplot(1,2,2); plot(obj_FISTA);  title('Objective plot'); colorbar;
+%%
\ No newline at end of file
diff --git a/demos/Demo_Phantom3D_Parallel.m b/demos/Demo_Phantom3D_Parallel.m
new file mode 100644
index 0000000..fd8096a
--- /dev/null
+++ b/demos/Demo_Phantom3D_Parallel.m
@@ -0,0 +1,49 @@
+% A demo script to reconstruct 3D synthetic data using FISTA method for

+% PARALLEL BEAM geometry

+% requirements: ASTRA-toolbox and TomoPhantom toolbox

+

+close all;clc;clear all;

+% adding paths

+addpath('../data/');

+addpath('../main_func/'); addpath('../main_func/regularizers_CPU/'); addpath('../main_func/regularizers_GPU/NL_Regul/'); addpath('../main_func/regularizers_GPU/Diffus_HO/');

+addpath('../supp/');

+

+%%

+% build 3D phantom using TomoPhantom and generate projection data

+modelNo = 3; % see Phantom3DLibrary.dat file in TomoPhantom

+N = 256; % x-y-z size (cubic image)

+angles = 1:0.5:180; % angles vector in degrees

+angles_rad = angles*(pi/180); % conversion to radians

+det_size = round(sqrt(2)*N); % detector size

+% in order to run functions you have to go to the directory:

+cd /home/algol/Documents/MATLAB/TomoPhantom/functions/

+TomoPhantom = buildPhantom3D(modelNo,N); % generate 3D phantom

+sino_tomophan3D = buildSino3D(modelNo, N, det_size, single(angles)); % generate data

+%%

+% using ASTRA-toolbox to set the projection geometry (parallel beam)

+proj_geom = astra_create_proj_geom('parallel', 1, det_size, angles_rad);

+vol_geom = astra_create_vol_geom(N,N);

+%% 

+fprintf('%s\n', 'Reconstruction using FISTA-LS without regularization...');

+clear params

+% define parameters

+params.proj_geom = proj_geom; % pass geometry to the function

+params.vol_geom = vol_geom;

+params.sino = single(sino_tomophan3D); % sinogram

+params.iterFISTA = 5; %max number of outer iterations

+params.X_ideal = TomoPhantom; % ideal phantom

+params.show = 1; % visualize reconstruction on each iteration

+params.slice = round(N/2); params.maxvalplot = 1; 

+tic; [X_FISTA, output] = FISTA_REC(params); toc; 

+

+error_FISTA = output.Resid_error; obj_FISTA = output.objective;

+fprintf('%s %.4f\n', 'Min RMSE for FISTA-PWLS reconstruction is:', min(error_FISTA(:)));

+

+Resid3D = (TomoPhantom - X_FISTA).^2;

+figure(2);

+subplot(1,2,1); imshow(X_FISTA(:,:,params.slice),[0 params.maxvalplot]); title('FISTA-LS reconstruction'); colorbar;

+subplot(1,2,2); imshow(Resid3D(:,:,params.slice),[0 0.1]);  title('residual'); colorbar;

+figure(3);

+subplot(1,2,1); plot(error_FISTA);  title('RMSE plot'); colorbar;

+subplot(1,2,2); plot(obj_FISTA);  title('Objective plot'); colorbar;

+%%
\ No newline at end of file