GH-OS fixed, demos updated, add_artifacts function

7 files changed, 105 insertions, 126 deletions

diff --git a/demos/DemoRD1.m b/demos/DemoRD1.m
deleted file mode 100644
index 5bb5f6b..0000000
--- a/demos/DemoRD1.m
+++ /dev/null
@@ -1,100 +0,0 @@
-% 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('../main_func/regularizers_CPU/');
-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.initialize = 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/Demo_Phantom3D_Cone.m b/demos/Demo_Phantom3D_Cone.m
index 3a9178b..a8f2c92 100644
--- a/demos/Demo_Phantom3D_Cone.m
+++ b/demos/Demo_Phantom3D_Cone.m
@@ -8,7 +8,6 @@ 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
@@ -16,9 +15,11 @@ 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:
+
+%---------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)
diff --git a/demos/Demo_Phantom3D_Parallel.m b/demos/Demo_Phantom3D_Parallel.m
index 6a54450..9ef1f3c 100644
--- a/demos/Demo_Phantom3D_Parallel.m
+++ b/demos/Demo_Phantom3D_Parallel.m
@@ -9,18 +9,27 @@ addpath('../main_func/'); addpath('../main_func/regularizers_CPU/'); addpath('..
 addpath('../supp/');

 

 %%

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

+% 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

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

+

+%---------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_artifacts = sino_add_artifacts(sino_tomophan3D,'rings');

+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

 %

 %%

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

@@ -28,29 +37,34 @@ 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 ...');

-for i = 1:k

+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', '-proj3d', proj_geom, sino_artifacts(:,:,k)); 

+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, 15);

-reconASTRA_3D = astra_mex_data2d('get', vol_id);

+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 FISTA-LS without regularization...');

+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(sino_tomophan3D); % sinogram

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

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

+params.iterFISTA = 12; %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.subsets = 12; 

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

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

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

 

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

@@ -63,4 +77,33 @@ subplot(1,2,2); imshow(Resid3D(:,:,params.slice),[0 0.1]);  title('residual'); c
 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 without 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 weight for PWLS

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

+params.Regul_Lambda_FGPTV = 0.003; % 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 = 1; 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/DemoRD2.m b/demos/Demo_RealData3D_Parallel.m
index 717a55d..e4c9eb0 100644
--- a/demos/DemoRD2.m
+++ b/demos/Demo_RealData3D_Parallel.m
@@ -11,12 +11,12 @@ 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
+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(size_det, angSize, slices_tot, 'single'); % log-corrected sino
+Sino3D = zeros(det_size, angSize, slices_tot, 'single'); % log-corrected sino
 % normalizing the data
 for  jj = 1:slices_tot
     sino = data_raw3D(:,:,jj);
@@ -30,10 +30,8 @@ Weights3D = single(data_raw3D); % weights for PW model
 clear data_raw3D
 %%
 % set projection/reconstruction geometry here
-Z_slices = 5;
-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);
+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);
@@ -50,7 +48,7 @@ params.iterFISTA  = 12;
 params.weights = Weights3D;
 params.subsets = 16; % the number of ordered subsets 
 params.show = 1;
-params.maxvalplot = 2.5; params.slice = 2;
+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');
@@ -76,13 +74,13 @@ 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 = 0.005; % TV regularization parameter for FGP-TV
+% params.Regul_Lambda_FGPTV = 0.005; % 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;
 params.subsets = 16; % the number of ordered subsets 
 params.show = 1;
-params.maxvalplot = 2.5; params.slice = 2;
+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');
diff --git a/main_func/FISTA_REC.m b/main_func/FISTA_REC.m
index bdaeb18..fa98360 100644
--- a/main_func/FISTA_REC.m
+++ b/main_func/FISTA_REC.m
@@ -669,6 +669,10 @@ else
             counterInd = counterInd + numProjSub;
         end
         
+        if (i == 1)
+            r_old = r;
+        end
+        
         % working with a 'ring vector'
         if (lambdaR_L1 > 0)
             r =  max(abs(r)-lambdaR_L1, 0).*sign(r); % soft-thresholding operator for ring vector
diff --git a/main_func/regularizers_CPU/FGP_TV.c b/main_func/regularizers_CPU/FGP_TV.c
index cfe5b9e..e0dae57 100644
--- a/main_func/regularizers_CPU/FGP_TV.c
+++ b/main_func/regularizers_CPU/FGP_TV.c
@@ -65,7 +65,7 @@ void mexFunction(
     A  = (float *) mxGetData(prhs[0]); /*noisy image (2D/3D) */
     lambda =  (float) mxGetScalar(prhs[1]); /* regularization parameter */
     iter = 50; /* default iterations number */
-    epsil = 0.001; /* default tolerance constant */
+    epsil = 0.0001; /* default tolerance constant */
     methTV = 0;  /* default isotropic TV penalty */
     
     if ((nrhs == 3) || (nrhs == 4) || (nrhs == 5))  iter = (int) mxGetScalar(prhs[2]); /* iterations number */
@@ -88,7 +88,7 @@ void mexFunction(
     
     tk = 1.0f;
     tkp1=1.0f;
-    count = 1;
+    count = 0;
     re_old = 0.0f;
     
     if (number_of_dims == 2) {
@@ -127,7 +127,7 @@ void mexFunction(
             }
             re = sqrt(re)/sqrt(re1);
             if (re < epsil)  count++;
-            if (count > 3) {
+            if (count > 4) {
                 Obj_func_CALC2D(A, D, funcvalA, lambda, dimX, dimY); 
                 break; }
             
diff --git a/supp/sino_add_artifacts.m b/supp/sino_add_artifacts.m
new file mode 100644
index 0000000..f601914
--- /dev/null
+++ b/supp/sino_add_artifacts.m
@@ -0,0 +1,33 @@
+function sino_artifacts = sino_add_artifacts(sino,artifact_type)
+% function to add various distortions to the sinogram space, current
+% version includes: random rings and zingers (streaks)
+% Input: 
+% 1. sinogram
+% 2. artifact type: 'rings' or 'zingers' (streaks)
+
+
+[Detectors, anglesNumb, SlicesZ] = size(sino);
+fprintf('%s %i %s %i %s %i %s \n', 'Sinogram has a dimension of', Detectors, 'detectors;', anglesNumb, 'projections;', SlicesZ, 'vertical slices.');
+
+sino_artifacts = sino;
+
+if (strcmp(artifact_type,'rings'))
+    fprintf('%s \n', 'Adding rings...');    
+    NumRings = round(Detectors/20); % Number of rings relatively to the size of Detectors
+    IntenOff = linspace(0.05,0.5,NumRings); % the intensity of rings in the selected range
+    
+    for k = 1:SlicesZ
+        % generate random indices to propagate rings
+        RandInd = randperm(Detectors,Detectors);
+        for jj = 1:NumRings
+            ind_c = RandInd(jj);
+            sino_artifacts(ind_c,1:end,k) = sino_artifacts(ind_c,1:end,k) + IntenOff(jj).*sino_artifacts(ind_c,1:end,k); % generate a constant offset            
+        end
+        
+    end
+elseif (strcmp(artifact_type,'zingers'))
+    fprintf('%s \n', 'Adding zingers...');
+else
+    fprintf('%s \n', 'Nothing selected, the same sinogram returned...');
+end
+end
\ No newline at end of file