Merge branch 'master' of https://github.com/vais-ral/CCPi-FISTA_Reconstruction into pythonize

4 files changed, 179 insertions, 265 deletions

diff --git a/demos/Demo1.m b/demos/Demo1.m
deleted file mode 100644
index 15e2e5b..0000000
--- a/demos/Demo1.m
+++ /dev/null
@@ -1,174 +0,0 @@
-% 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('../main_func/regularizers_CPU/');

-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(1,2,1, [0.05 0.05]); imshow(X_FISTA,[0 0.6]); title('FISTA-PWLS reconstruction'); colorbar;

-subplot(1,2,2, [0.05 0.05]); imshow((phantom - X_FISTA).^2,[0 0.1]);  title('residual'); colorbar;

-colormap(cmapnew); 

-figure(3); clf

-subplot(1,2,1, [0.05 0.05]); plot(error_FISTA);  title('RMSE plot'); colorbar;

-subplot(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(1,2,1, [0.05 0.05]); imshow(X_FISTA_TV,[0 0.6]); title('FISTA-PWLS-TV reconstruction'); colorbar;

-subplot(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(1,2,1, [0.05 0.05]); plot(error_FISTA_TV);  title('RMSE plot'); colorbar;

-subplot(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(1,2,1, [0.05 0.05]); imshow(X_FISTA_GH_TV,[0 0.6]); title('FISTA-GH-TV reconstruction'); colorbar;

-subplot(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(1,2,1, [0.05 0.05]);  plot(error_FISTA_GH_TV);  title('RMSE plot'); colorbar;

-subplot(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(1,2,1, [0.05 0.05]); imshow(X_FISTA_student_TV,[0 0.6]); title('FISTA-Student-TV reconstruction'); colorbar;

-subplot(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(1,2,1, [0.05 0.05]); plot(error_FISTA_student_TV);  title('RMSE plot'); colorbar;

-subplot(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/Demo_Phantom3D_Cone.m b/demos/Demo_Phantom3D_Cone.m
new file mode 100644
index 0000000..6419386
--- /dev/null
+++ b/demos/Demo_Phantom3D_Cone.m
@@ -0,0 +1,66 @@
+% 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
+% in order to run functions you have to go to the directory:
+cd /home/algol/Documents/MATLAB/TomoPhantom/functions/
+TomoPhantom = buildPhantom3D(modelNo,N); % 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
new file mode 100644
index 0000000..fd8096a
--- /dev/null
+++ b/demos/Demo_Phantom3D_Parallel.m
@@ -0,0 +1,49 @@
+% 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 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 and generate projection data

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

+cd /home/algol/Documents/MATLAB/TomoPhantom/functions/

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

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

+%%

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

+params.iterFISTA = 5; %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-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'); colorbar;

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

+%%
\ No newline at end of file
diff --git a/main_func/FISTA_REC.m b/main_func/FISTA_REC.m
index 6987dca..dde0e73 100644
--- a/main_func/FISTA_REC.m
+++ b/main_func/FISTA_REC.m
@@ -15,7 +15,7 @@ function [X,  output] = FISTA_REC(params)
 %----------------General Parameters------------------------
 %       - .proj_geom (geometry of the projector) [required]
 %       - .vol_geom (geometry of the reconstructed object) [required]
-%       - .sino (vectorized in 2D or 3D sinogram) [required]
+%       - .sino (2D or 3D sinogram) [required]
 %       - .iterFISTA (iterations for the main loop, default 40)
 %       - .L_const (Lipschitz constant, default Power method)                                                                                                    )
 %       - .X_ideal (ideal image, if given)
@@ -94,45 +94,36 @@ else
     weights = ones(size(sino));
 end
 if (isfield(params,'fidelity'))
-	studentt = 0;
+    studentt = 0;
     if (strcmp(params.fidelity,'studentt') == 1)
-    studentt = 1;
-    lambdaR_L1 = 0;
-    end    
+        studentt = 1;
+    end
 else
-    studentt = 0;    
+    studentt = 0;
 end
 if (isfield(params,'L_const'))
     L_const = params.L_const;
 else
     % using Power method (PM) to establish L constant
-    if (strcmp(proj_geom.type,'parallel') || strcmp(proj_geom.type,'parallel3d'))
-        % for parallel geometry we can do just one slice
-        fprintf('%s \n', 'Calculating Lipshitz constant for parallel beam geometry...');
+    fprintf('%s %s %s \n', 'Calculating Lipshitz constant for',proj_geom.type, 'beam geometry...');
+    if (strcmp(proj_geom.type,'parallel') || strcmp(proj_geom.type,'fanflat') || strcmp(proj_geom.type,'fanflat_vec'))       
+        % for 2D geometry we can do just one selected slice        
         niter = 15; % number of iteration for the PM
         x1 = rand(N,N,1);
         sqweight = sqrt(weights(:,:,1));
-        proj_geomT = proj_geom;
-        proj_geomT.DetectorRowCount = 1;
-        vol_geomT = vol_geom;
-        vol_geomT.GridSliceCount = 1;
-        [sino_id, y] = astra_create_sino3d_cuda(x1, proj_geomT, vol_geomT);
-        y = sqweight.*y;
-        astra_mex_data3d('delete', sino_id);
-        
+        [sino_id, y] = astra_create_sino_cuda(x1, proj_geom, vol_geom);
+        y = sqweight.*y';
+        astra_mex_data2d('delete', sino_id);        
         for i = 1:niter
-            [id,x1] = astra_create_backprojection3d_cuda(sqweight.*y, proj_geomT, vol_geomT);
+            [x1] = astra_create_backprojection_cuda((sqweight.*y)', proj_geom, vol_geom);
             s = norm(x1(:));
             x1 = x1./s;
-            [sino_id, y] = astra_create_sino3d_cuda(x1, proj_geomT, vol_geomT);
-            y = sqweight.*y;
-            astra_mex_data3d('delete', sino_id);
-            astra_mex_data3d('delete', id);
-        end
-        %clear proj_geomT vol_geomT
-    else
-        % divergen beam geometry
-        fprintf('%s \n', 'Calculating Lipshitz constant for divergen beam geometry... will take some time!');
+            [sino_id, y] = astra_create_sino_cuda(x1, proj_geom, vol_geom);
+            y = sqweight.*y';
+            astra_mex_data2d('delete', sino_id);
+        end        
+    elseif (strcmp(proj_geom.type,'cone') || strcmp(proj_geom.type,'parallel3d') || strcmp(proj_geom.type,'parallel3d_vec') || strcmp(proj_geom.type,'cone_vec'))
+        % 3D geometry        
         niter = 8; % number of iteration for PM
         x1 = rand(N,N,SlicesZ);
         sqweight = sqrt(weights);
@@ -150,6 +141,8 @@ else
             astra_mex_data3d('delete', id);
         end
         clear x1
+    else
+        error('%s \n', 'No suitable geometry has been found!');
     end
     L_const = s;
 end
@@ -272,28 +265,10 @@ else
     slice = 1;
 end
 if (isfield(params,'initialize'))
-    % a 'warm start' with SIRT method
-    % Create a data object for the reconstruction
-    rec_id = astra_mex_data3d('create', '-vol', vol_geom);
-    
-    sinogram_id = astra_mex_data3d('create', '-proj3d', proj_geom, sino);
-    
-    % Set up the parameters for a reconstruction algorithm using the GPU
-    cfg = astra_struct('SIRT3D_CUDA');
-    cfg.ReconstructionDataId = rec_id;
-    cfg.ProjectionDataId = sinogram_id;
-    
-    % Create the algorithm object from the configuration structure
-    alg_id = astra_mex_algorithm('create', cfg);
-    astra_mex_algorithm('iterate', alg_id, 35);
-    % Get the result
-    X = astra_mex_data3d('get', rec_id);
-    
-    % Clean up. Note that GPU memory is tied up in the algorithm object,
-    % and main RAM in the data objects.
-    astra_mex_algorithm('delete', alg_id);
-    astra_mex_data3d('delete', rec_id);
-    astra_mex_data3d('delete', sinogram_id);
+    X  = params.initialize;
+    if ((size(X,1) ~= N) || (size(X,2) ~= N) || (size(X,3) ~= SlicesZ))
+        error('%s \n', 'The initialized volume has different dimensions!');
+    end        
 else
     X = zeros(N,N,SlicesZ, 'single'); % storage for the solution
 end
@@ -345,16 +320,18 @@ if (subsets == 0)
         t_old = t;
         r_old = r;
         
-        % if the geometry is parallel use slice-by-slice projection-backprojection routine
-        if (strcmp(proj_geom.type,'parallel') || strcmp(proj_geom.type,'parallel3d'))
+        % if the geometry is 2D use slice-by-slice projection-backprojection routine
+        if (strcmp(proj_geom.type,'parallel') || strcmp(proj_geom.type,'fanflat') || strcmp(proj_geom.type,'fanflat_vec'))    
             sino_updt = zeros(size(sino),'single');
-            for kkk = 1:SlicesZ
-                [sino_id, sino_updt(:,:,kkk)] = astra_create_sino3d_cuda(X_t(:,:,kkk), proj_geomT, vol_geomT);
-                astra_mex_data3d('delete', sino_id);
+            for kkk = 1:SlicesZ                
+                [sino_id, sinoT] = astra_create_sino_cuda(X_t(:,:,kkk), proj_geom, vol_geom);
+                sino_updt(:,:,kkk) = sinoT';
+                astra_mex_data2d('delete', sino_id);
             end
         else
-            % for divergent 3D geometry (watch the GPU memory overflow in earlier ASTRA versions < 1.8)
+            % for 3D geometry (watch the GPU memory overflow in earlier ASTRA versions < 1.8)
             [sino_id, sino_updt] = astra_create_sino3d_cuda(X_t, proj_geom, vol_geom);
+            astra_mex_data3d('delete', sino_id);
         end
         
         if (lambdaR_L1 > 0)
@@ -368,40 +345,37 @@ if (subsets == 0)
             end
             r = r_x - (1./L_const).*vec;
             objective(i) = (0.5*sum(residual(:).^2)); % for the objective function output
-        else            
-			if (studentt == 1)
-				% artifacts removal with Students t penalty
-				residual = weights.*(sino_updt - sino);
-				for kkk = 1:SlicesZ
-				res_vec = reshape(residual(:,:,kkk), Detectors*anglesNumb, 1); % 1D vectorized sinogram
+        elseif (studentt > 0)
+            % artifacts removal with Students t penalty
+            residual = weights.*(sino_updt - sino);
+            for kkk = 1:SlicesZ
+                res_vec = reshape(residual(:,:,kkk), Detectors*anglesNumb, 1); % 1D vectorized sinogram
                 %s = 100;
                 %gr = (2)*res_vec./(s*2 + conj(res_vec).*res_vec);
                 [ff, gr] = studentst(res_vec, 1);
                 residual(:,:,kkk) = reshape(gr, Detectors, anglesNumb);
-				end            
-				objective(i) = ff; % for the objective function output
-            else
+            end
+            objective(i) = ff; % for the objective function output
+        else
             % no ring removal (LS model)
             residual = weights.*(sino_updt - sino);
             objective(i) = (0.5*sum(residual(:).^2)); % for the objective function output
-            end	
-        end       
+        end
+        
         
-        % if the geometry is parallel use slice-by-slice projection-backprojection routine
-        if (strcmp(proj_geom.type,'parallel') || strcmp(proj_geom.type,'parallel3d'))
+        % if the geometry is 2D use slice-by-slice projection-backprojection routine
+        if (strcmp(proj_geom.type,'parallel') || strcmp(proj_geom.type,'fanflat') || strcmp(proj_geom.type,'fanflat_vec'))    
             x_temp = zeros(size(X),'single');
             for kkk = 1:SlicesZ
-                [id, x_temp(:,:,kkk)] = astra_create_backprojection3d_cuda(squeeze(residual(:,:,kkk)), proj_geomT, vol_geomT);
-                astra_mex_data3d('delete', id);
+                [x_temp(:,:,kkk)] = astra_create_backprojection_cuda(squeeze(residual(:,:,kkk))', proj_geom, vol_geom);
             end
         else
             [id, x_temp] = astra_create_backprojection3d_cuda(residual, proj_geom, vol_geom);
+            astra_mex_data3d('delete', id);
         end
-        X = X_t - (1/L_const).*x_temp;
-        astra_mex_data3d('delete', sino_id);
-        astra_mex_data3d('delete', id);
+        X = X_t - (1/L_const).*x_temp;          
         
-        % regularization
+        % ----------------Regularization part------------------------
         if (lambdaFGP_TV > 0)
             % FGP-TV regularization
             if ((strcmp('2D', Dimension) == 1))
@@ -572,28 +546,27 @@ else
                 end
                 r = r_x - (1./L_const).*vec;
             else
-                [sino_id, sino_updt] = astra_create_sino3d_cuda(X_t, proj_geomSUB, vol_geom);                
+                [sino_id, sino_updt] = astra_create_sino3d_cuda(X_t, proj_geomSUB, vol_geom);
                 
                 if (studentt == 1)
-				% artifacts removal with Students t penalty
-				residual = squeeze(weights(:,CurrSubIndeces,:)).*(sino_updt - squeeze(sino(:,CurrSubIndeces,:)));
-				
-				for kkk = 1:SlicesZ
-				res_vec = reshape(residual(:,:,kkk), Detectors*numProjSub, 1); % 1D vectorized sinogram
-                %s = 100;
-                %gr = (2)*res_vec./(s*2 + conj(res_vec).*res_vec);
-                [ff, gr] = studentst(res_vec, 1);
-                residual(:,:,kkk) = reshape(gr, Detectors, numProjSub);
-				end            
-				objective(i) = ff; % for the objective function output
-				else
-				% no ring removal (LS model)
-				residual = squeeze(weights(:,CurrSubIndeces,:)).*(sino_updt - squeeze(sino(:,CurrSubIndeces,:)));				
-				end	               
+                    % artifacts removal with Students t penalty
+                    residual = squeeze(weights(:,CurrSubIndeces,:)).*(sino_updt - squeeze(sino(:,CurrSubIndeces,:)));
+                    
+                    for kkk = 1:SlicesZ
+                        res_vec = reshape(residual(:,:,kkk), Detectors*numProjSub, 1); % 1D vectorized sinogram
+                        %s = 100;
+                        %gr = (2)*res_vec./(s*2 + conj(res_vec).*res_vec);
+                        [ff, gr] = studentst(res_vec, 1);
+                        residual(:,:,kkk) = reshape(gr, Detectors, numProjSub);
+                    end
+                    objective(i) = ff; % for the objective function output
+                else
+                    % no ring removal (LS model)
+                    residual = squeeze(weights(:,CurrSubIndeces,:)).*(sino_updt - squeeze(sino(:,CurrSubIndeces,:)));
+                end
             end
             
             [id, x_temp] = astra_create_backprojection3d_cuda(residual, proj_geomSUB, vol_geom);
-            
             X = X_t - (1/L_const).*x_temp;
             astra_mex_data3d('delete', sino_id);
             astra_mex_data3d('delete', id);