Major reorganization, updated routines

3 files changed, 400 insertions, 0 deletions

diff --git a/demos/Demo1.m b/demos/Demo1.m
new file mode 100644
index 0000000..486b97c
--- /dev/null
+++ b/demos/Demo1.m
@@ -0,0 +1,174 @@
+% 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('../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_tight(1,2,1, [0.05 0.05]); imshow(X_FISTA,[0 0.6]); title('FISTA-PWLS 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-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_tight(1,2,1, [0.05 0.05]); imshow(X_FISTA_TV,[0 0.6]); title('FISTA-PWLS-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.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_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.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_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_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/DemoRD1.m b/demos/DemoRD1.m
new file mode 100644
index 0000000..c25bb3e
--- /dev/null
+++ b/demos/DemoRD1.m
@@ -0,0 +1,100 @@
+% 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');
+%%
+% set projection/reconstruction geometry here
+Z_slices = 1;
+det_row_count = Z_slices;
+proj_geom = astra_create_proj_geom('parallel3d', 1, 1, det_row_count, size_det, angles);
+vol_geom = astra_create_vol_geom(recon_size,recon_size,Z_slices);
+%%
+fprintf('%s\n', 'Reconstruction using FISTA-LS without regularization...');
+clear params
+params.proj_geom = proj_geom; % pass geometry to the function
+params.vol_geom = vol_geom;
+params.sino = sino_basalt;
+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.proj_geom = proj_geom; % pass geometry to the function
+params.vol_geom = vol_geom;
+params.sino = sino_basalt;
+params.iterFISTA  = 60;
+params.Regul_LambdaTV = 0.0003; % TV regularization parameter
+params.show = 0;
+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.proj_geom = proj_geom; % pass geometry to the function
+params.vol_geom = vol_geom;
+params.sino = sino_basalt;
+params.iterFISTA  = 60;
+params.Regul_LambdaTV = 0.0003; % TV regularization parameter
+params.Ring_LambdaR_L1 = 0.001; % Soft-Thresh L1 ring variable parameter
+params.Ring_Alpha = 20; % acceleration for ring variable
+params.show = 0;
+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.proj_geom = proj_geom; % pass geometry to the function
+params.vol_geom = vol_geom;
+params.sino = sino_basalt;
+params.iterFISTA  = 50;
+params.L_const = 3500; % Lipshitz constant
+params.Regul_LambdaTV = 0.0003; % TV regularization parameter
+params.fidelity = 'student'; % choosing Student t penalty
+params.show = 1;
+params.initilize = 1; % warm start with SIRT
+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/demos/DemoRD2.m b/demos/DemoRD2.m
new file mode 100644
index 0000000..ab4da96
--- /dev/null
+++ b/demos/DemoRD2.m
@@ -0,0 +1,126 @@
+% Demonstration of tomographic 3D reconstruction from X-ray synchrotron
+% dataset (dendrites) using various data fidelities
+% warning: can take up to 15-20 minutes to run for the whole 3D data
+clear all
+close all
+%%
+% % adding paths
+addpath('../data/');
+addpath('../main_func/');
+addpath('../supp/');
+
+load('DendrRawData.mat') % load raw data of 3D dendritic set
+angles_rad = angles*(pi/180); % conversion to radians
+size_det = size(data_raw3D,1); % detectors dim
+angSize = size(data_raw3D, 2); % angles dim
+slices_tot = size(data_raw3D, 3); % no of slices
+recon_size = 950; % reconstruction size
+
+Sino3D = zeros(size_det, angSize, slices_tot, 'single'); % log-corrected sino
+% normalizing the data
+for  jj = 1:slices_tot
+    sino = data_raw3D(:,:,jj);
+    for ii = 1:angSize
+        Sino3D(:,ii,jj) = log((flats_ar(:,jj)-darks_ar(:,jj))./(single(sino(:,ii)) - darks_ar(:,jj)));
+    end
+end
+
+Sino3D = Sino3D.*1000;
+Weights3D = single(data_raw3D); % weights for PW model
+clear data_raw3D
+%%
+% set projection/reconstruction geometry here
+Z_slices = 20;
+det_row_count = Z_slices;
+proj_geom = astra_create_proj_geom('parallel3d', 1, 1, det_row_count, size_det, angles_rad);
+vol_geom = astra_create_vol_geom(recon_size,recon_size,Z_slices);
+%%
+fprintf('%s\n', 'Reconstruction using FBP...');
+FBP = iradon(Sino3D(:,:,10), angles,recon_size);
+figure; imshow(FBP , [0, 3]); title ('FBP reconstruction');
+%%
+fprintf('%s\n', 'Reconstruction using FISTA-PWLS without regularization...');
+clear params
+params.proj_geom = proj_geom; % pass geometry to the function
+params.vol_geom = vol_geom;
+params.sino = Sino3D;
+params.L_const = 7.6789e+08; % found quickly for one slice first
+params.iterFISTA  = 30;
+params.weights = Weights3D;
+params.show = 1;
+params.maxvalplot = 2.5; params.slice = 4;
+
+tic; [X_fista, output] = FISTA_REC(params); toc;
+figure; imshow(X_fista(:,:,1) , [0, 2.5]); title ('FISTA-PWLS reconstruction');
+%%
+fprintf('%s\n', 'Reconstruction using FISTA-PWLS-TV...');
+clear params
+params.proj_geom = proj_geom; % pass geometry to the function
+params.vol_geom = vol_geom;
+params.sino = Sino3D;
+params.iterFISTA  = 40;
+params.L_const = 7.6789e+08; 
+params.Regul_LambdaTV = 0.005; % TV regularization parameter for FISTA-TV
+params.weights = Weights3D;
+params.show = 1;
+params.maxvalplot = 2.5; params.slice = 10;
+
+tic; [X_fista_TV] = FISTA_REC(params); toc;
+figure; imshow(X_fista_TV(:,:,1) , [0, 2.5]); title ('FISTA-PWLS-TV reconstruction');
+%%
+%%
+fprintf('%s\n', 'Reconstruction using FISTA-GH-TV...');
+clear params
+params.proj_geom = proj_geom; % pass geometry to the function
+params.vol_geom = vol_geom;
+params.sino = Sino3D;
+params.iterFISTA  = 40;
+params.Regul_LambdaTV = 0.005; % TV regularization parameter for FISTA-TV
+params.Ring_LambdaR_L1 = 0.002; % Soft-Thresh L1 ring variable parameter
+params.Ring_Alpha = 21; % to boost ring removal procedure
+params.weights = Weights3D;
+params.show = 1;
+params.maxvalplot = 2.5; params.slice = 10;
+
+tic; [X_fista_GH_TV] = FISTA_REC(params); toc;
+figure; imshow(X_fista_GH_TV(:,:,1) , [0, 2.5]); title ('FISTA-GH-TV reconstruction');
+%%
+%%
+fprintf('%s\n', 'Reconstruction using FISTA-GH-TV-LLT...');
+clear params
+params.proj_geom = proj_geom; % pass geometry to the function
+params.vol_geom = vol_geom;
+params.sino = Sino3D;
+params.iterFISTA  = 40;
+params.Regul_LambdaTV = 0.005; % TV regularization parameter for FISTA-TV
+params.Regul_LambdaHO = 200;  % regularization parameter for LLT problem
+params.Regul_tauHO = 0.0005; % time-step parameter for the explicit scheme
+params.Regul_iterHO = 250; % the max number of TV iterations
+params.Ring_LambdaR_L1 = 0.002; % Soft-Thresh L1 ring variable parameter
+params.Ring_Alpha = 21; % to boost ring removal procedure
+params.weights = Weights3D;
+params.show = 1;
+params.maxvalplot = 2.5; params.slice = 10;
+
+tic; [X_fista_GH_TVLLT] = FISTA_REC(params); toc;
+figure; imshow(X_fista_GH_TVLLT(:,:,1) , [0, 2.5]); title ('FISTA-GH-TV-LLT reconstruction');
+%%
+%%
+% fprintf('%s\n', 'Reconstruction using FISTA-Student-TV...');
+% clear params
+% params.sino = Sino3D(:,:,10);
+% params.N  = recon_size;
+% params.angles = angles_rad;
+% params.iterFISTA  = 100;
+% params.L_const = 0.01; % Lipshitz constant
+% params.lambdaTV = 0.006; % TV regularization parameter for FISTA-TV
+% params.tol = 1.0e-04;
+% params.iterTV = 20;
+% params.fidelity = 'student'; % choosing Student t penalty
+% params.weights = Weights3D(:,:,10);
+% params.show = 0;
+% params.maxvalplot = 2.5; params.slice = 1;
+% 
+% tic; [X_fistaStudentTV] = FISTA_REC(params); toc;
+% figure; imshow(X_fistaStudentTV(:,:,1), [0, 2.5]); title ('FISTA-Student-TV reconstruction');
+%%