summaryrefslogtreecommitdiffstats
path: root/main_func
diff options
context:
space:
mode:
Diffstat (limited to 'main_func')
-rw-r--r--main_func/FISTA_REC.m155
1 files changed, 64 insertions, 91 deletions
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);
generated by cgit v1.2.3 (git 2.25.1) at 2024-11-08 09:10:26 +0000