summaryrefslogtreecommitdiffstats
diff options
context:
space:
mode:
Diffstat
-rw-r--r--demos/Demo1.m174
-rw-r--r--demos/Demo_Phantom3D_Cone.m66
-rw-r--r--demos/Demo_Phantom3D_Parallel.m50
-rw-r--r--main_func/FISTA_REC.m256
4 files changed, 234 insertions, 312 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..ac9827c
--- /dev/null
+++ b/demos/Demo_Phantom3D_Parallel.m
@@ -0,0 +1,50 @@
+% 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.subsets = 12;
+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');
+subplot(1,2,2); plot(obj_FISTA); title('Objective plot');
+%% \ No newline at end of file
diff --git a/main_func/FISTA_REC.m b/main_func/FISTA_REC.m
index 6987dca..d177748 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);
+ [sino_id, y] = astra_create_sino_cuda(x1, proj_geom, vol_geom);
+ y = sqweight.*y';
+ astra_mex_data2d('delete', sino_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!');
+ 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,19 @@ 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 (strcmp(proj_geom.type,'parallel') || strcmp(proj_geom.type,'fanflat') || strcmp(proj_geom.type,'fanflat_vec'))
+ % if geometry is 2D use slice-by-slice projection-backprojection routine
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);
+ [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 +346,36 @@ 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
+ objective(i) = 0.5*norm(residual(:)); % for the objective function output
+ 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);
- % regularization
+ % ----------------Regularization part------------------------%
if (lambdaFGP_TV > 0)
% FGP-TV regularization
if ((strcmp('2D', Dimension) == 1))
@@ -510,95 +484,104 @@ if (subsets == 0)
end
end
else
- % Ordered Subsets (OS) FISTA reconstruction routine (normally one order of magnitude faster than classical)
+ % Ordered Subsets (OS) FISTA reconstruction routine (normally one order of magnitude faster than the classical version)
t = 1;
X_t = X;
- proj_geomSUB = proj_geom;
-
+ proj_geomSUB = proj_geom;
r = zeros(Detectors,SlicesZ, 'single'); % 2D array (for 3D data) of sparse "ring" vectors
r_x = r; % another ring variable
residual2 = zeros(size(sino),'single');
+ sino_updt_FULL = zeros(size(sino),'single');
% Outer FISTA iterations loop
- for i = 1:iterFISTA
+ for i = 1:iterFISTA
- % With OS approach it becomes trickier to correlate independent subsets, hence additional work is required
- % one solution is to work with a full sinogram at times
- if ((i >= 3) && (lambdaR_L1 > 0))
- [sino_id2, sino_updt2] = astra_create_sino3d_cuda(X, proj_geom, vol_geom);
- astra_mex_data3d('delete', sino_id2);
+ if ((i > 1) && (lambdaR_L1 > 0))
+ % in order to make Group-Huber fidelity work with ordered subsets
+ % we still need to work with full sinogram
+
+ % the offset variable must be calculated for the whole
+ % updated sinogram - sino_updt_FULL
+ for kkk = 1:anglesNumb
+ residual2(:,kkk,:) = squeeze(weights(:,kkk,:)).*(squeeze(sino_updt_FULL(:,kkk,:)) - (squeeze(sino(:,kkk,:)) - alpha_ring.*r_x));
+ end
+
+ r_old = r;
+ vec = sum(residual2,2);
+ if (SlicesZ > 1)
+ vec = squeeze(vec(:,1,:));
+ end
+ r = r_x - (1./L_const).*vec; % update ring variable
end
% subsets loop
counterInd = 1;
+ if (strcmp(proj_geom.type,'parallel') || strcmp(proj_geom.type,'fanflat') || strcmp(proj_geom.type,'fanflat_vec'))
+ % if geometry is 2D use slice-by-slice projection-backprojection routine
+ for kkk = 1:SlicesZ
+ [sino_id, sinoT] = astra_create_sino_cuda(X_t(:,:,kkk), proj_geomSUB, vol_geom);
+ sino_updt_Sub(:,:,kkk) = sinoT';
+ astra_mex_data2d('delete', sino_id);
+ end
+ else
+ % for 3D geometry (watch the GPU memory overflow in earlier ASTRA versions < 1.8)
+ [sino_id, sino_updt_Sub] = astra_create_sino3d_cuda(X_t, proj_geomSUB, vol_geom);
+ astra_mex_data3d('delete', sino_id);
+ end
for ss = 1:subsets
X_old = X;
- t_old = t;
- r_old = r;
+ t_old = t;
numProjSub = binsDiscr(ss); % the number of projections per subset
CurrSubIndeces = IndicesReorg(counterInd:(counterInd + numProjSub - 1)); % extract indeces attached to the subset
proj_geomSUB.ProjectionAngles = angles(CurrSubIndeces);
+ sino_updt_Sub = zeros(Detectors, numProjSub, SlicesZ,'single');
if (lambdaR_L1 > 0)
+ % Group-Huber fidelity (ring removal)
- % the ring removal part (Group-Huber fidelity)
- % first 2 iterations do additional work reconstructing whole dataset to ensure
- % the stablility
- if (i < 3)
- [sino_id2, sino_updt2] = astra_create_sino3d_cuda(X_t, proj_geom, vol_geom);
- astra_mex_data3d('delete', sino_id2);
- else
- [sino_id, sino_updt] = astra_create_sino3d_cuda(X_t, proj_geomSUB, vol_geom);
- end
- for kkk = 1:anglesNumb
- residual2(:,kkk,:) = squeeze(weights(:,kkk,:)).*(squeeze(sino_updt2(:,kkk,:)) - (squeeze(sino(:,kkk,:)) - alpha_ring.*r_x));
- end
-
- residual = zeros(Detectors, numProjSub, SlicesZ,'single');
+ residualSub = zeros(Detectors, numProjSub, SlicesZ,'single'); % residual for a chosen subset
for kkk = 1:numProjSub
indC = CurrSubIndeces(kkk);
- if (i < 3)
- residual(:,kkk,:) = squeeze(residual2(:,indC,:));
- else
- residual(:,kkk,:) = squeeze(weights(:,indC,:)).*(squeeze(sino_updt(:,kkk,:)) - (squeeze(sino(:,indC,:)) - alpha_ring.*r_x));
- end
+ residualSub(:,kkk,:) = squeeze(weights(:,indC,:)).*(squeeze(sino_updt_Sub(:,kkk,:)) - (squeeze(sino(:,indC,:)) - alpha_ring.*r_x));
+ sino_updt_FULL(:,indC,:) = squeeze(sino_updt_Sub(:,kkk,:)); % filling the full sinogram
end
- vec = sum(residual2,2);
- if (SlicesZ > 1)
- vec = squeeze(vec(:,1,:));
+
+ elseif (studentt > 0)
+ % student t data fidelity
+
+
+ % artifacts removal with Students t penalty
+ residualSub = squeeze(weights(:,CurrSubIndeces,:)).*(sino_updt_Sub - squeeze(sino(:,CurrSubIndeces,:)));
+
+ for kkk = 1:SlicesZ
+ res_vec = reshape(residualSub(:,:,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);
+ residualSub(:,:,kkk) = reshape(gr, Detectors, numProjSub);
end
- r = r_x - (1./L_const).*vec;
+ objective(i) = ff; % for the objective function output
else
- [sino_id, sino_updt] = astra_create_sino3d_cuda(X_t, proj_geomSUB, vol_geom);
+ % PWLS model
- 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
- end
- [id, x_temp] = astra_create_backprojection3d_cuda(residual, proj_geomSUB, vol_geom);
+ if (strcmp(proj_geom.type,'parallel') || strcmp(proj_geom.type,'fanflat') || strcmp(proj_geom.type,'fanflat_vec'))
+ % if geometry is 2D use slice-by-slice projection-backprojection routine
+ x_temp = zeros(size(X),'single');
+ for kkk = 1:SlicesZ
+ [x_temp(:,:,kkk)] = astra_create_backprojection_cuda(squeeze(residualSub(:,:,kkk))', proj_geomSUB, vol_geom);
+ end
+ else
+ [id, x_temp] = astra_create_backprojection3d_cuda(residualSub, proj_geomSUB, 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))
@@ -680,20 +663,17 @@ else
end
end
- if (lambdaR_L1 > 0)
- r = max(abs(r)-lambdaR_L1, 0).*sign(r); % soft-thresholding operator for ring vector
- end
-
t = (1 + sqrt(1 + 4*t^2))/2; % updating t
X_t = X + ((t_old-1)/t).*(X - X_old); % updating X
-
- if (lambdaR_L1 > 0)
- r_x = r + ((t_old-1)/t).*(r - r_old); % updating r
- end
-
counterInd = counterInd + numProjSub;
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
+ r_x = r + ((t_old-1)/t).*(r - r_old); % updating r
+ end
+
if (show == 1)
figure(10); imshow(X(:,:,slice), [0 maxvalplot]);
if (lambdaR_L1 > 0)
generated by cgit v1.2.3 (git 2.25.1) at 2024-11-26 18:31:26 +0000