all Matlab related stuff have been moved to wrappers

4 files changed, 409 insertions, 0 deletions

diff --git a/Wrappers/Matlab/demos/Demo_Phantom3D_Cone.m b/Wrappers/Matlab/demos/Demo_Phantom3D_Cone.m
new file mode 100644
index 0000000..a8f2c92
--- /dev/null
+++ b/Wrappers/Matlab/demos/Demo_Phantom3D_Cone.m
@@ -0,0 +1,67 @@
+% 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
+
+%---------TomoPhantom routines---------%
+pathTP = '/home/algol/Documents/MATLAB/TomoPhantom/functions/models/Phantom3DLibrary.dat'; % path to TomoPhantom parameters file
+TomoPhantom = buildPhantom3D(modelNo,N,pathTP); % 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/Wrappers/Matlab/demos/Demo_Phantom3D_Parallel.m b/Wrappers/Matlab/demos/Demo_Phantom3D_Parallel.m
new file mode 100644
index 0000000..4219bd1
--- /dev/null
+++ b/Wrappers/Matlab/demos/Demo_Phantom3D_Parallel.m
@@ -0,0 +1,121 @@
+% 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;

+% 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/');

+

+%%

+% Main reconstruction/data generation parameters 

+modelNo = 2; % 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

+

+%---------TomoPhantom routines---------%

+pathTP = '/home/algol/Documents/MATLAB/TomoPhantom/functions/models/Phantom3DLibrary.dat'; % path to TomoPhantom parameters file

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

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

+%--------------------------------------%

+% Adding noise and distortions if required

+sino_tomophan3D = sino_add_artifacts(sino_tomophan3D,'rings');

+% adding Poisson noise

+dose =  3e9; % photon flux (controls noise level)

+multifactor = max(sino_tomophan3D(:));

+dataExp = dose.*exp(-sino_tomophan3D/multifactor); % noiseless raw data

+dataRaw = astra_add_noise_to_sino(dataExp, dose); % pre-log noisy raw data (weights)

+sino3D_log = log(dose./max(dataRaw,1))*multifactor; %log corrected data -> sinogram

+clear dataExp sino_tomophan3D

+%

+%%

+%-------------Astra toolbox------------%

+% one can generate data using ASTRA toolbox

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

+vol_geom = astra_create_vol_geom(N,N);

+sino_ASTRA3D = zeros(det_size, length(angles), N, 'single');

+for i = 1:N

+[sino_id, sinoT] = astra_create_sino_cuda(TomoPhantom(:,:,i), proj_geom, vol_geom);

+sino_ASTRA3D(:,:,i) = sinoT';

+astra_mex_data2d('delete', sino_id);

+end

+%--------------------------------------%

+%%

+% 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', 'Reconstructing with FBP using ASTRA-toolbox ...');

+reconASTRA_3D = zeros(size(TomoPhantom),'single');

+for k = 1:N

+vol_id = astra_mex_data2d('create', '-vol', vol_geom, 0);

+proj_id = astra_mex_data2d('create', '-sino', proj_geom, sino3D_log(:,:,k)'); 

+cfg = astra_struct('FBP_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, 1);

+rec = astra_mex_data2d('get', vol_id);

+reconASTRA_3D(:,:,k) = single(rec);

+end

+figure; imshow(reconASTRA_3D(:,:,128), [0 1.3]);

+%% 

+%% 

+fprintf('%s\n', 'Reconstruction using OS-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 = single(sino3D_log); % sinogram

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

+params.X_ideal = TomoPhantom; % ideal phantom

+params.weights = dataRaw./max(dataRaw(:)); % statistical weight for PWLS

+params.subsets = 12; % the number of subsets

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

+params.slice = 128; params.maxvalplot = 1.3; 

+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'); 

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

+%%

+%% 

+fprintf('%s\n', 'Reconstruction using OS-FISTA-GH with FGP-TV regularization...');

+clear params

+% define parameters

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

+params.vol_geom = vol_geom;

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

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

+params.X_ideal = TomoPhantom; % ideal phantom

+params.weights = dataRaw./max(dataRaw(:)); % statistical weights for PWLS

+params.subsets = 12; % the number of subsets

+params.Regul_Lambda_FGPTV = 100; % TV regularization parameter for FGP-TV

+params.Ring_LambdaR_L1 = 0.02; % Soft-Thresh L1 ring variable parameter

+params.Ring_Alpha = 21; % to boost ring removal procedure

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

+params.slice = 128; params.maxvalplot = 1.3; 

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

+

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

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

+

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

+figure(2);

+subplot(1,2,1); imshow(X_FISTA_GH_TV(:,:,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_GH_TV);  title('RMSE plot'); 

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

+%%
\ No newline at end of file
diff --git a/Wrappers/Matlab/demos/Demo_RealData3D_Parallel.m b/Wrappers/Matlab/demos/Demo_RealData3D_Parallel.m
new file mode 100644
index 0000000..f82e0b0
--- /dev/null
+++ b/Wrappers/Matlab/demos/Demo_RealData3D_Parallel.m
@@ -0,0 +1,186 @@
+% Demonstration of tomographic 3D reconstruction from X-ray synchrotron
+% dataset (dendrites) using various data fidelities
+% ! It is advisable not to run the whole script, it will take lots of time to reconstruct the whole 3D data using many algorithms !
+clear
+close 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/');
+
+load('DendrRawData.mat') % load raw data of 3D dendritic set
+angles_rad = angles*(pi/180); % conversion to radians
+det_size = 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(det_size, 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
+proj_geom = astra_create_proj_geom('parallel', 1, det_size, angles_rad);
+vol_geom = astra_create_vol_geom(recon_size,recon_size);
+%%
+fprintf('%s\n', 'Reconstruction using FBP...');
+FBP = iradon(Sino3D(:,:,10), angles,recon_size);
+figure; imshow(FBP , [0, 3]); title ('FBP reconstruction');
+
+%--------FISTA_REC modular reconstruction alogrithms---------
+%%
+fprintf('%s\n', 'Reconstruction using FISTA-OS-PWLS without regularization...');
+clear params
+params.proj_geom = proj_geom; % pass geometry to the function
+params.vol_geom = vol_geom;
+params.sino = Sino3D;
+params.iterFISTA  = 18;
+params.weights = Weights3D;
+params.subsets = 8; % the number of ordered subsets 
+params.show = 1;
+params.maxvalplot = 2.5; params.slice = 1;
+
+tic; [X_fista, outputFISTA] = FISTA_REC(params); toc;
+figure; imshow(X_fista(:,:,params.slice) , [0, 2.5]); title ('FISTA-OS-PWLS reconstruction');
+%%
+fprintf('%s\n', 'Reconstruction using FISTA-OS-PWLS-TV...');
+clear params
+params.proj_geom = proj_geom; % pass geometry to the function
+params.vol_geom = vol_geom;
+params.sino = Sino3D;
+params.iterFISTA  = 18;
+params.Regul_Lambda_FGPTV = 5.0000e+6; % TV regularization parameter for FGP-TV
+params.weights = Weights3D;
+params.subsets = 8; % the number of ordered subsets 
+params.show = 1;
+params.maxvalplot = 2.5; params.slice = 10;
+
+tic; [X_fista_TV, outputTV] = FISTA_REC(params); toc;
+figure; imshow(X_fista_TV(:,:,params.slice) , [0, 2.5]); title ('FISTA-OS-PWLS-TV reconstruction');
+%%
+fprintf('%s\n', 'Reconstruction using FISTA-OS-GH-TV...');
+clear params
+params.proj_geom = proj_geom; % pass geometry to the function
+params.vol_geom = vol_geom;
+params.sino = Sino3D(:,:,10);
+params.iterFISTA  = 18;
+params.Regul_Lambda_FGPTV = 5.0000e+6; % TV regularization parameter for FGP-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(:,:,10);
+params.subsets = 8; % the number of ordered subsets 
+params.show = 1;
+params.maxvalplot = 2.5; params.slice = 1;
+
+tic; [X_fista_GH_TV, outputGHTV] = FISTA_REC(params); toc;
+figure; imshow(X_fista_GH_TV(:,:,params.slice) , [0, 2.5]); title ('FISTA-OS-GH-TV reconstruction');
+%%
+fprintf('%s\n', 'Reconstruction using FISTA-OS-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  = 12;
+params.Regul_Lambda_FGPTV = 5.0000e+6; % TV regularization parameter for FGP-TV
+params.Regul_LambdaLLT = 100;  % regularization parameter for LLT problem
+params.Regul_tauLLT = 0.0005; % time-step parameter for the explicit scheme
+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.subsets = 16; % the number of ordered subsets 
+params.show = 1;
+params.maxvalplot = 2.5; params.slice = 2;
+
+tic; [X_fista_GH_TVLLT, outputGH_TVLLT] = FISTA_REC(params); toc;
+figure; imshow(X_fista_GH_TVLLT(:,:,params.slice) , [0, 2.5]); title ('FISTA-OS-GH-TV-LLT reconstruction');
+
+%%
+fprintf('%s\n', 'Reconstruction using FISTA-OS-GH-HigherOrderDiffusion...');
+% !GPU version!
+clear params
+params.proj_geom = proj_geom; % pass geometry to the function
+params.vol_geom = vol_geom;
+params.sino = Sino3D(:,:,1:5);
+params.iterFISTA  = 25;
+params.Regul_LambdaDiffHO = 2; % DiffHO regularization parameter 
+params.Regul_DiffHO_EdgePar = 0.05; % threshold parameter
+params.Regul_Iterations = 150;
+params.Ring_LambdaR_L1 = 0.002; % Soft-Thresh L1 ring variable parameter
+params.Ring_Alpha = 21; % to boost ring removal procedure
+params.weights = Weights3D(:,:,1:5);
+params.subsets = 16; % the number of ordered subsets 
+params.show = 1;
+params.maxvalplot = 2.5; params.slice = 1;
+ 
+tic; [X_fista_GH_HO, outputHO] = FISTA_REC(params); toc;
+figure; imshow(X_fista_GH_HO(:,:,params.slice) , [0, 2.5]); title ('FISTA-OS-HigherOrderDiffusion reconstruction');
+
+%%
+fprintf('%s\n', 'Reconstruction using FISTA-PB...');
+% !GPU version!
+clear params
+params.proj_geom = proj_geom; % pass geometry to the function
+params.vol_geom = vol_geom;
+params.sino = Sino3D(:,:,1);
+params.iterFISTA  = 25;
+params.Regul_LambdaPatchBased_GPU = 3; % PB regularization parameter 
+params.Regul_PB_h = 0.04; % threhsold parameter
+params.Regul_PB_SearchW = 3;
+params.Regul_PB_SimilW  = 1;
+params.Ring_LambdaR_L1 = 0.002; % Soft-Thresh L1 ring variable parameter
+params.Ring_Alpha = 21; % to boost ring removal procedure
+params.weights = Weights3D(:,:,1);
+params.show = 1;
+params.maxvalplot = 2.5; params.slice = 1;
+ 
+tic; [X_fista_GH_PB, outputPB] = FISTA_REC(params); toc;
+figure; imshow(X_fista_GH_PB(:,:,params.slice) , [0, 2.5]); title ('FISTA-OS-PB reconstruction');
+%%
+fprintf('%s\n', 'Reconstruction using FISTA-OS-GH-TGV...');
+% still testing...
+clear params
+params.proj_geom = proj_geom; % pass geometry to the function
+params.vol_geom = vol_geom;
+params.sino = Sino3D;
+params.iterFISTA  = 12;
+params.Regul_LambdaTGV = 0.5; % TGV regularization parameter 
+params.Regul_Iterations = 5;
+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.subsets = 16; % the number of ordered subsets 
+params.show = 1;
+params.maxvalplot = 2.5; params.slice = 1;
+
+tic; [X_fista_GH_TGV, outputTGV] = FISTA_REC(params); toc;
+figure; imshow(X_fista_GH_TGV(:,:,params.slice) , [0, 2.5]); title ('FISTA-OS-GH-TGV 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 = Sino3D(:,:,10);
+% params.iterFISTA  = 50;
+% params.L_const = 0.01; % Lipshitz constant
+% params.Regul_LambdaTV = 0.008; % TV regularization parameter for FISTA-TV
+% params.fidelity = 'student'; % choosing Student t penalty
+% params.weights = Weights3D(:,:,10);
+% params.show = 0;
+% params.initialize = 1;
+% 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');
+%%
diff --git a/Wrappers/Matlab/demos/exportDemoRD2Data.m b/Wrappers/Matlab/demos/exportDemoRD2Data.m
new file mode 100644
index 0000000..028353b
--- /dev/null
+++ b/Wrappers/Matlab/demos/exportDemoRD2Data.m
@@ -0,0 +1,35 @@
+clear all
+close all
+%%
+% % adding paths
+addpath('../data/');
+addpath('../main_func/'); addpath('../main_func/regularizers_CPU/');
+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
+
+hdf5write('DendrData.h5', '/Weights3D', Weights3D)
+hdf5write('DendrData.h5', '/Sino3D', Sino3D, 'WriteMode', 'append')
+hdf5write('DendrData.h5', '/angles_rad', angles_rad,  'WriteMode', 'append')
+hdf5write('DendrData.h5', '/size_det', size_det,  'WriteMode', 'append')
+hdf5write('DendrData.h5', '/angSize', angSize,  'WriteMode', 'append')
+hdf5write('DendrData.h5', '/slices_tot', slices_tot,  'WriteMode', 'append')
+hdf5write('DendrData.h5', '/recon_size', recon_size,  'WriteMode', 'append')
\ No newline at end of file