all Matlab related stuff have been moved to wrappers

5 files changed, 0 insertions, 409 deletions

diff --git a/demos/Demo_Phantom3D_Cone.m b/demos/Demo_Phantom3D_Cone.m
deleted file mode 100644
index a8f2c92..0000000
--- a/demos/Demo_Phantom3D_Cone.m
+++ /dev/null
@@ -1,67 +0,0 @@
-% 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/demos/Demo_Phantom3D_Parallel.m b/demos/Demo_Phantom3D_Parallel.m
deleted file mode 100644
index 4219bd1..0000000
--- a/demos/Demo_Phantom3D_Parallel.m
+++ /dev/null
@@ -1,121 +0,0 @@
-% 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/demos/Demo_RealData3D_Parallel.m b/demos/Demo_RealData3D_Parallel.m
deleted file mode 100644
index f82e0b0..0000000
--- a/demos/Demo_RealData3D_Parallel.m
+++ /dev/null
@@ -1,186 +0,0 @@
-% 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/demos/DendrData.h5 b/demos/DendrData.h5
deleted file mode 100644
index f048268..0000000
--- a/demos/DendrData.h5
+++ /dev/null
diff --git a/demos/exportDemoRD2Data.m b/demos/exportDemoRD2Data.m
deleted file mode 100644
index 028353b..0000000
--- a/demos/exportDemoRD2Data.m
+++ /dev/null
@@ -1,35 +0,0 @@
-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