Initial revision

4 files changed, 545 insertions, 0 deletions

diff --git a/demo/Demo1.m b/demo/Demo1.m
new file mode 100644
index 0000000..08d46e1
--- /dev/null
+++ b/demo/Demo1.m
@@ -0,0 +1,160 @@
+% 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
+clear all
+close all 
+
+% adding paths
+addpath('data/'); 
+addpath('main_func/'); 
+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);
+
+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-LS without regularization...');
+clear params
+% define parameters
+params.sino = sino_zing_rings;
+params.N = N; % image size 
+params.angles = theta_rad; % angles in radians 
+params.iterFISTA = 180; %max number of outer iterations
+params.X_ideal = phantom; % ideal phantom
+params.ROI = ROI; % phantom region-of-interest
+params.show = 0; % visualize reconstruction on each iteration
+params.slice = 1; params.maxvalplot = 0.6; 
+params.weights = Dweights; % statistical weighting 
+tic; [X_FISTA, error_FISTA, obj_FISTA, sinoFISTA] = FISTA_REC(params); toc; 
+
+fprintf('%s %.4f\n', 'Min RMSE for FISTA-LS reconstruction is:', min(error_FISTA(:)));
+
+figure(2); clf
+%set(gcf, 'Position', get(0,'Screensize'));
+subplot_tight(1,2,1, [0.05 0.05]); imshow(X_FISTA,[0 0.6]); title('FISTA-LS reconstruction'); colorbar;
+subplot_tight(1,2,2, [0.05 0.05]); imshow((phantom - X_FISTA).^2,[0 0.1]);  title('residual'); colorbar;
+colormap(cmapnew); 
+figure(3); clf
+subplot_tight(1,2,1, [0.05 0.05]); plot(error_FISTA);  title('RMSE plot'); colorbar;
+subplot_tight(1,2,2, [0.05 0.05]); plot(obj_FISTA);  title('Objective plot'); colorbar;
+colormap(cmapnew); 
+%%
+fprintf('%s\n', 'Reconstruction using FISTA-LS-TV...');
+clear params
+% define parameters
+params.sino = sino_zing_rings;
+params.N = N; % image size 
+params.angles = theta_rad; % angles in radians 
+params.iterFISTA = 200; % max number of outer iterations
+params.lambdaTV = 5.39e-05; % regularization parameter for TV problem
+params.tol = 1.0e-04; % tolerance to terminate TV iterations
+params.iterTV = 20; % the max number of TV iterations
+params.X_ideal = phantom; % ideal phantom
+params.ROI = ROI; % phantom region-of-interest
+params.weights = Dweights; % statistical weighting 
+params.show = 0; % visualize reconstruction on each iteration
+params.slice = 1; params.maxvalplot = 0.6; 
+tic; [X_FISTA_TV, error_FISTA_TV, obj_FISTA_TV, sinoFISTA_TV] = FISTA_REC(params); toc; 
+
+fprintf('%s %.4f\n', 'Min RMSE for FISTA-LS-TV reconstruction is:', min(error_FISTA_TV(:)));
+
+figure(4); clf
+subplot_tight(1,2,1, [0.05 0.05]); imshow(X_FISTA_TV,[0 0.6]); title('FISTA-LS-TV reconstruction'); colorbar;
+subplot_tight(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_tight(1,2,1, [0.05 0.05]); plot(error_FISTA_TV);  title('RMSE plot'); colorbar;
+subplot_tight(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.sino = sino_zing_rings;
+params.N = N; % image size 
+params.angles = theta_rad; % angles in radians 
+params.iterFISTA = 60; % max number of outer iterations
+params.lambdaTV = 0.002526;  % regularization parameter for TV problem
+params.tol = 1.0e-04; % tolerance to terminate TV iterations
+params.iterTV = 20; % the max number of TV iterations
+params.X_ideal = phantom; % ideal phantom
+params.ROI = ROI; % phantom region-of-interest
+params.weights = Dweights; % statistical weighting 
+params.lambdaR_L1 = 0.002; % parameter to sparsify the "rings vector"
+params.show = 0; % visualize reconstruction on each iteration
+params.slice = 1; params.maxvalplot = 0.6; 
+tic; [X_FISTA_GH_TV, error_FISTA_GH_TV, obj_FISTA_GH_TV, sinoFISTA_GH_TV] = FISTA_REC(params); toc; 
+
+fprintf('%s %.4f\n', 'Min RMSE for FISTA-GH-TV reconstruction is:', min(error_FISTA_GH_TV(:)));
+
+figure(6); clf
+subplot_tight(1,2,1, [0.05 0.05]); imshow(X_FISTA_GH_TV,[0 0.6]); title('FISTA-GH-TV reconstruction'); colorbar;
+subplot_tight(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_tight(1,2,1, [0.05 0.05]);  plot(error_FISTA_GH_TV);  title('RMSE plot'); colorbar;
+subplot_tight(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.sino = sino_zing_rings;
+params.N = N; % image size 
+params.angles = theta_rad; % angles in radians 
+params.iterFISTA = 67; % max number of outer iterations
+%params.L_const = 80000; % Lipshitz constant (can be chosen manually to accelerate convergence)
+params.lambdaTV = 0.00152;  % regularization parameter for TV problem
+params.tol = 1.0e-04; % tolerance to terminate TV iterations
+params.iterTV = 20; % the max number of TV iterations
+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 = 0; % visualize reconstruction on each iteration
+params.slice = 1; params.maxvalplot = 0.6; 
+tic; [X_FISTA_student_TV, error_FISTA_student_TV, obj_FISTA_student_TV, sinoFISTA_student_TV] = FISTA_REC(params); toc; 
+
+fprintf('%s %.4f\n', 'Min RMSE for FISTA-Student-TV reconstruction is:', min(error_FISTA_student_TV(:)));
+
+figure(8); 
+set(gcf, 'Position', get(0,'Screensize'));
+subplot_tight(1,2,1, [0.05 0.05]); imshow(X_FISTA_student_TV,[0 0.6]); title('FISTA-Student-TV reconstruction'); colorbar;
+subplot_tight(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_tight(1,2,1, [0.05 0.05]); plot(error_FISTA_student_TV);  title('RMSE plot'); colorbar;
+subplot_tight(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_2(:), phantom(:)));
+fprintf('%s %.4f\n', 'Min RMSE for FISTA-LS reconstruction:', min(error_FISTA(:)));
+fprintf('%s %.4f\n', 'Min RMSE for FISTA-LS-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/demo/Demo2.m b/demo/Demo2.m
new file mode 100644
index 0000000..3c1592c
--- /dev/null
+++ b/demo/Demo2.m
@@ -0,0 +1,156 @@
+% Demonstration of tomographic reconstruction from noisy and corrupted by 
+% artifacts undersampled projection data using Students t penalty 
+% This is the missing wedge demo, run it after DemoFISTA_StudT
+
+% see ReadMe file for instructions
+% clear all
+% close all 
+
+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.0);
+
+add_wedges % apply the missing wedge mask
+
+%% 
+fprintf('%s\n', 'Direct reconstruction using FBP...');
+FBP_1 = iradon(MW_sino_artifacts', theta, N); 
+
+fprintf('%s %.4f\n', 'RMSE for FBP reconstruction:', RMSE(FBP_1(:), phantom(:)));
+
+figure(1); 
+% set(gcf, 'Position', get(0,'Screensize'));
+subplot_tight(1,2,1, [0.05 0.05]);  imshow(FBP_1,[-2 0.8]);  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-LS without regularization...');
+clear params
+% define parameters
+params.sino = MW_sino_artifacts;
+params.N = N; % image size 
+params.angles = theta_rad; % angles in radians 
+params.iterFISTA = 132; %max number of outer iterations
+params.X_ideal = phantom; % ideal phantom
+params.ROI = ROI; % phantom region-of-interest
+params.show = 0; % visualize reconstruction on each iteration
+params.slice = 1; params.maxvalplot = 0.6; 
+params.weights = Dweights; % statistical weighting 
+tic; [X_FISTA, error_FISTA, obj_FISTA, sinoFISTA] = FISTA_REC(params); toc; 
+
+fprintf('%s %.4f\n', 'Min RMSE for FISTA-LS reconstruction:', min(error_FISTA(:)));
+
+figure(2); clf
+%set(gcf, 'Position', get(0,'Screensize'));
+subplot_tight(1,2,1, [0.05 0.05]); imshow(X_FISTA,[0 0.6]); title('FISTA-LS reconstruction'); colorbar;
+subplot_tight(1,2,2, [0.05 0.05]); imshow((phantom - X_FISTA).^2,[0 0.1]);  title('residual'); colorbar;
+colormap(cmapnew); 
+figure(3); clf
+subplot_tight(1,2,1, [0.05 0.05]); plot(error_FISTA);  title('RMSE plot'); colorbar;
+subplot_tight(1,2,2, [0.05 0.05]); plot(obj_FISTA);  title('Objective plot'); colorbar;
+colormap(cmapnew); 
+%%
+fprintf('%s\n', 'Reconstruction using FISTA-LS-TV...');
+clear params
+% define parameters
+params.sino = MW_sino_artifacts;
+params.N = N; % image size 
+params.angles = theta_rad; % angles in radians 
+params.iterFISTA = 200; % max number of outer iterations
+params.lambdaTV = 5.39e-05; % regularization parameter for TV problem
+params.tol = 1.0e-04; % tolerance to terminate TV iterations
+params.iterTV = 20; % the max number of TV iterations
+params.X_ideal = phantom; % ideal phantom
+params.ROI = ROI; % phantom region-of-interest
+params.weights = Dweights; % statistical weighting 
+params.show = 0; % visualize reconstruction on each iteration
+params.slice = 1; params.maxvalplot = 0.6; 
+tic; [X_FISTA_TV, error_FISTA_TV, obj_FISTA_TV, sinoFISTA_TV] = FISTA_REC(params); toc; 
+
+fprintf('%s %.4f\n', 'Min RMSE for FISTA-LS-TV reconstruction:', min(error_FISTA_TV(:)));
+
+figure(4); clf
+subplot_tight(1,2,1, [0.05 0.05]); imshow(X_FISTA_TV,[0 0.6]); title('FISTA-LS-TV reconstruction'); colorbar;
+subplot_tight(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_tight(1,2,1, [0.05 0.05]); plot(error_FISTA_TV);  title('RMSE plot'); colorbar;
+subplot_tight(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.sino = MW_sino_artifacts;
+params.N = N; % image size 
+params.angles = theta_rad; % angles in radians 
+params.iterFISTA = 250; % max number of outer iterations
+params.lambdaTV = 0.0019;  % regularization parameter for TV problem
+params.tol = 1.0e-04; % tolerance to terminate TV iterations
+params.iterTV = 20; % the max number of TV iterations
+params.X_ideal = phantom; % ideal phantom
+params.ROI = ROI; % phantom region-of-interest
+params.weights = Dweights; % statistical weighting 
+params.lambdaR_L1 = 0.002; % parameter to sparsify the "rings vector"
+params.show = 0; % visualize reconstruction on each iteration
+params.slice = 1; params.maxvalplot = 0.6; 
+tic; [X_FISTA_GH_TV, error_FISTA_GH_TV, obj_FISTA_GH_TV, sinoFISTA_GH_TV] = FISTA_REC(params); toc; 
+
+fprintf('%s %.4f\n', 'Min RMSE for FISTA-GH-TV reconstruction:', min(error_FISTA_GH_TV(:)));
+
+figure(6); clf
+subplot_tight(1,2,1, [0.05 0.05]); imshow(X_FISTA_GH_TV,[0 0.6]); title('FISTA-GH-TV reconstruction'); colorbar;
+subplot_tight(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_tight(1,2,1, [0.05 0.05]);  plot(error_FISTA_GH_TV);  title('RMSE plot'); colorbar;
+subplot_tight(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.sino = MW_sino_artifacts;
+params.N = N; % image size 
+params.angles = theta_rad; % angles in radians 
+params.iterFISTA = 80; % max number of outer iterations
+% params.L_const = 80000; % Lipshitz constant (can be chosen manually to accelerate convergence)
+params.lambdaTV = 0.0016;  % regularization parameter for TV problem
+params.tol = 1.0e-04; % tolerance to terminate TV iterations
+params.iterTV = 20; % the max number of TV iterations
+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 = 0; % visualize reconstruction on each iteration
+params.slice = 1; params.maxvalplot = 0.6; 
+tic; [X_FISTA_student_TV, error_FISTA_student_TV, obj_FISTA_student_TV, sinoFISTA_student_TV] = FISTA_REC(params); toc; 
+
+fprintf('%s %.4f\n', 'Min RMSE for FISTA-Student-TV reconstruction:', min(error_FISTA_student_TV(:)));
+
+figure(8); 
+set(gcf, 'Position', get(0,'Screensize'));
+subplot_tight(1,2,1, [0.05 0.05]); imshow(X_FISTA_student_TV,[0 0.6]); title('FISTA-Student-TV reconstruction'); colorbar;
+subplot_tight(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_tight(1,2,1, [0.05 0.05]); plot(error_FISTA_student_TV);  title('RMSE plot'); colorbar;
+subplot_tight(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_2(:), phantom(:)));
+fprintf('%s %.4f\n', 'Min RMSE for FISTA-LS reconstruction:', min(error_FISTA(:)));
+fprintf('%s %.4f\n', 'Min RMSE for FISTA-LS-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/demo/DemoRD1.m b/demo/DemoRD1.m
new file mode 100644
index 0000000..9a43cb5
--- /dev/null
+++ b/demo/DemoRD1.m
@@ -0,0 +1,99 @@
+% Demonstration of tomographic reconstruction from neutron tomography
+% dataset (basalt sample) using Student t data fidelity 
+clear all
+close all 
+
+% adding paths
+addpath('data/'); 
+addpath('main_func/'); 
+addpath('supp/'); 
+
+load('sino_basalt.mat') % load real neutron data
+
+size_det = size(sino_basalt, 1); % detector size
+angSize = size(sino_basalt,2); % angles dim
+recon_size = 650; % reconstruction size
+
+FBP = iradon(sino_basalt, rad2deg(angles),recon_size);
+figure; imshow(FBP , [0, 0.45]); title ('FBP reconstruction');
+
+%%
+fprintf('%s\n', 'Reconstruction using FISTA-LS without regularization...');
+clear params
+params.sino = sino_basalt';
+params.N  = recon_size;
+params.angles = angles;
+params.iterFISTA  = 50;
+params.show = 0;
+params.maxvalplot = 0.6; params.slice = 1;
+
+tic; [X_fista] = FISTA_REC(params); toc;
+figure; imshow(X_fista , [0, 0.45]); title ('FISTA-LS reconstruction');
+%%
+fprintf('%s\n', 'Reconstruction using FISTA-LS-TV...');
+clear params
+params.sino = sino_basalt';
+params.N  = recon_size;
+params.angles = angles;
+params.iterFISTA  = 150;
+params.lambdaTV = 0.0003; % TV regularization parameter
+params.tol = 1.0e-04;
+params.iterTV = 20;
+params.show = 1;
+params.maxvalplot = 0.6; params.slice = 1;
+
+tic; [X_fista_TV] = FISTA_REC(params); toc;
+figure; imshow(X_fista_TV , [0, 0.45]); title ('FISTA-LS-TV reconstruction');
+%%
+%%
+fprintf('%s\n', 'Reconstruction using FISTA-GH-TV...');
+clear params
+params.sino = sino_basalt';
+params.N  = recon_size;
+params.angles = angles;
+params.iterFISTA  = 350;
+params.lambdaTV = 0.0003; % TV regularization parameter
+params.tol = 1.0e-04;
+params.iterTV = 20;
+params.lambdaR_L1 = 0.001; % Soft-Thresh L1 ring variable parameter
+params.show = 1;
+params.maxvalplot = 0.6; params.slice = 1;
+
+tic; [X_fista_GH_TV] = FISTA_REC(params); toc;
+figure; imshow(X_fista_GH_TV , [0, 0.45]); title ('FISTA-GH-TV reconstruction');
+%%
+%%
+fprintf('%s\n', 'Reconstruction using FISTA-Student-TV...');
+clear params
+params.sino = sino_basalt';
+params.N  = recon_size;
+params.angles = angles;
+params.iterFISTA  = 350;
+params.L_const = 7000; % Lipshitz constant
+params.lambdaTV = 0.0003; % TV regularization parameter
+params.tol = 1.0e-04;
+params.iterTV = 20;
+params.fidelity = 'student'; % choosing Student t penalty
+params.show = 1;
+params.maxvalplot = 0.6; params.slice = 1;
+
+tic; [X_fistaStudentTV] = FISTA_REC(params); toc;
+figure; imshow(X_fistaStudentTV , [0, 0.45]); title ('FISTA-Student-TV reconstruction');
+%%
+
+fprintf('%s\n', 'Segmentation using OTSU method ...');
+level = graythresh(X_fista);
+Segm_FISTA = im2bw(X_fista,level);
+figure; imshow(Segm_FISTA, []); title ('Segmented FISTA-LS reconstruction');
+
+level = graythresh(X_fista_TV);
+Segm_FISTA_TV = im2bw(X_fista_TV,level);
+figure; imshow(Segm_FISTA_TV, []); title ('Segmented FISTA-LS-TV reconstruction');
+
+level = graythresh(X_fista_GH_TV);
+BW_FISTA_GH_TV = im2bw(X_fista_GH_TV,level);
+figure; imshow(BW_FISTA_GH_TV, []); title ('Segmented FISTA-GH-TV reconstruction');
+
+level = graythresh(X_fistaStudentTV);
+BW_FISTA_Student_TV = im2bw(X_fistaStudentTV,level);
+figure; imshow(BW_FISTA_Student_TV, []); title ('Segmented FISTA-Student-LS reconstruction');
\ No newline at end of file
diff --git a/demo/DemoRD2.m b/demo/DemoRD2.m
new file mode 100644
index 0000000..a8ac2ca
--- /dev/null
+++ b/demo/DemoRD2.m
@@ -0,0 +1,130 @@
+% Demonstration of tomographic 3D reconstruction from X-ray synchrotron 
+% dataset (dendrites) using various data fidelities 
+% clear all
+% close all 
+% 
+% % adding paths
+ addpath('data/'); 
+ addpath('main_func/'); 
+ addpath('supp/'); 
+
+load('sino3D_dendrites.mat') % load 3D normalized sinogram
+angles_rad = angles*(pi/180); % conversion to radians
+
+angSize = size(Sino3D,1); % angles dim
+size_det = size(Sino3D, 2); % detector size
+recon_size = 850; % reconstruction size
+
+FBP = iradon(Sino3D(:,:,10)', angles,recon_size);
+figure; imshow(FBP , [0, 3]); title ('FBP reconstruction');
+
+%%
+fprintf('%s\n', 'Reconstruction using FISTA-LS without regularization...');
+clear params
+params.sino = Sino3D;
+params.N  = recon_size;
+params.angles = angles_rad;
+params.iterFISTA  = 80;
+params.precondition  = 1; % switch on preconditioning
+params.show = 0;
+params.maxvalplot = 2.5; params.slice = 10;
+
+tic; [X_fista] = FISTA_REC(params); toc;
+figure; imshow(X_fista(:,:,10) , [0, 2.5]); title ('FISTA-LS reconstruction');
+%%
+fprintf('%s\n', 'Reconstruction using FISTA-LS-TV...');
+clear params
+params.sino = Sino3D;
+params.N  = recon_size;
+params.angles = angles_rad;
+params.iterFISTA  = 100;
+params.lambdaTV = 0.001; % TV regularization parameter for FISTA-TV
+params.tol = 1.0e-04;
+params.iterTV = 20;
+params.precondition  = 1; % switch on preconditioning
+params.show = 0;
+params.maxvalplot = 2.5; params.slice = 10;
+
+tic; [X_fista_TV] = FISTA_REC(params); toc;
+figure; imshow(X_fista_TV(:,:,10) , [0, 2.5]); title ('FISTA-LS-TV reconstruction');
+%%
+%%
+fprintf('%s\n', 'Reconstruction using FISTA-GH-TV...');
+clear params
+params.sino = Sino3D;
+params.N  = recon_size;
+params.angles = angles_rad;
+params.iterFISTA  = 100;
+params.lambdaTV = 0.001; % TV regularization parameter for FISTA-TV
+params.tol = 1.0e-04;
+params.iterTV = 20;
+params.lambdaR_L1 = 0.001; % Soft-Thresh L1 ring variable parameter
+params.alpha_ring = 20; % to boost ring removal procedure
+params.precondition  = 1; % switch on preconditioning
+params.show = 0;
+params.maxvalplot = 2.5; params.slice = 10;
+
+tic; [X_fista_GH_TV] = FISTA_REC(params); toc;
+figure; imshow(X_fista_GH_TV(:,:,10) , [0, 2.5]); title ('FISTA-GH-TV reconstruction');
+%%
+%%
+fprintf('%s\n', 'Reconstruction using FISTA-GH-TV-LLT...');
+clear params
+params.sino = Sino3D;
+params.N  = recon_size;
+params.angles = angles_rad;
+params.iterFISTA  = 100;
+params.lambdaTV = 0.001; % TV regularization parameter for FISTA-TV
+params.tol = 1.0e-04;
+params.iterTV = 20;
+params.lambdaHO = 35;  % regularization parameter for LLT problem
+params.tauHO = 0.00011; % time-step parameter for explicit scheme
+params.iterHO = 70; % the max number of TV iterations   
+params.lambdaR_L1 = 0.001; % Soft-Thresh L1 ring variable parameter
+params.alpha_ring = 20; % to boost ring removal procedure
+params.precondition  = 1; % switch on preconditioning
+params.show = 0;
+params.maxvalplot = 2.5; params.slice = 10;
+
+tic; [X_fista_GH_TVLLT] = FISTA_REC(params); toc;
+figure; imshow(X_fista_GH_TVLLT(:,:,10) , [0, 2.5]); title ('FISTA-GH-TV-LLT reconstruction');
+%%
+%%
+% fprintf('%s\n', 'Reconstruction using FISTA-Student-TV...');
+% %%%%<<<< Not stable with this dataset! Requires more work >>>> %%%%%
+% clear params
+% params.sino = Sino3D(:,:,15);
+% params.N  = 950;
+% params.angles = angles_rad;
+% params.iterFISTA  = 150;
+% params.L_const = 30; % Lipshitz constant
+% params.lambdaTV = 0.009; % TV regularization parameter for FISTA-TV
+% params.tol = 1.0e-04;
+% params.iterTV = 20;
+% params.fidelity = 'student'; % choosing Student t penalty
+% % params.precondition  = 1; % switch on preconditioning
+% params.show = 1;
+% params.maxvalplot = 2.5; params.slice = 1;
+% 
+% tic; [X_fistaStudentTV] = FISTA_REC(params); toc;
+% figure; imshow(X_fistaStudentTV , [0, 2.5]); title ('FISTA-Student-TV reconstruction');
+%%
+slice = 10; % if 3D reconstruction
+
+fprintf('%s\n', 'Segmentation using OTSU method ...');
+level = graythresh(X_fista(:,:,slice));
+Segm_FISTA = im2bw(X_fista(:,:,slice),level);
+figure; imshow(Segm_FISTA, []); title ('Segmented FISTA-LS reconstruction');
+
+level = graythresh(X_fista_TV(:,:,slice));
+Segm_FISTA_TV = im2bw(X_fista_TV(:,:,slice),level);
+figure; imshow(Segm_FISTA_TV, []); title ('Segmented FISTA-LS-TV reconstruction');
+
+level = graythresh(X_fista_GH_TV(:,:,slice));
+BW_FISTA_GH_TV = im2bw(X_fista_GH_TV(:,:,slice),level);
+figure; imshow(BW_FISTA_GH_TV, []); title ('Segmented FISTA-GH-TV reconstruction');
+
+level = graythresh(X_fista_GH_TVLLT(:,:,slice));
+BW_FISTA_GH_TVLLT = im2bw(X_fista_GH_TVLLT(:,:,slice),level);
+figure; imshow(BW_FISTA_GH_TVLLT, []); title ('Segmented FISTA-GH-TV-LLT reconstruction');
+%%
\ No newline at end of file