summaryrefslogtreecommitdiffstats
path: root/main_func
diff options
context:
space:
mode:
Diffstat (limited to 'main_func')
-rw-r--r--main_func/FISTA_REC.m267
1 files changed, 132 insertions, 135 deletions
diff --git a/main_func/FISTA_REC.m b/main_func/FISTA_REC.m
index 6987dca..1e4228d 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,121 @@ 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;
-
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
- % 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;
numProjSub = binsDiscr(ss); % the number of projections per subset
+ sino_updt_Sub = zeros(Detectors, numProjSub, SlicesZ,'single');
CurrSubIndeces = IndicesReorg(counterInd:(counterInd + numProjSub - 1)); % extract indeces attached to the subset
proj_geomSUB.ProjectionAngles = angles(CurrSubIndeces);
- if (lambdaR_L1 > 0)
-
- % 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));
+ 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
+
+ if (lambdaR_L1 > 0)
+ % Group-Huber fidelity (ring removal)
+
- 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);
-
- 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
+ % PWLS model
+
+ residualSub = squeeze(weights(:,CurrSubIndeces,:)).*(sino_updt_Sub - squeeze(sino(:,CurrSubIndeces,:)));
+ objective(i) = 0.5*norm(residualSub(:)); % for the objective function output
end
+
- [id, x_temp] = astra_create_backprojection3d_cuda(residual, proj_geomSUB, vol_geom);
+ % perform backprojection of a subset
+ 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);
- % regularization
+ % ----------------Regularization part------------------------%
if (lambdaFGP_TV > 0)
% FGP-TV regularization
if ((strcmp('2D', Dimension) == 1))
@@ -680,20 +680,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 2025-01-07 17:23:42 +0000