From 29636540aca6354b0f319765b0af9bf768593565 Mon Sep 17 00:00:00 2001
From: algol <dkazanc@hotmail.com>
Date: Wed, 5 Jul 2017 22:15:00 +0100
Subject: ring removal bug fixed
---
main_func/FISTA_REC.m | 157 +++++++++++++++++++++++++-------------------------
1 file changed, 77 insertions(+), 80 deletions(-)
(limited to 'main_func')
diff --git a/main_func/FISTA_REC.m b/main_func/FISTA_REC.m
index fcd19d8..db2b581 100644
--- a/main_func/FISTA_REC.m
+++ b/main_func/FISTA_REC.m
@@ -152,11 +152,6 @@ if (isfield(params,'Ring_Alpha'))
else
alpha_ring = 1;
end
-if (isfield(params,'fidelity'))
- fidelity = params.fidelity;
-else
- fidelity = 'LS';
-end
if (isfield(params,'show'))
show = params.show;
else
@@ -173,7 +168,7 @@ else
slice = 1;
end
if (isfield(params,'initialize'))
- % a 'warm start' with SIRT method
+ % a 'warm start' with SIRT method
% Create a data object for the reconstruction
rec_id = astra_mex_data3d('create', '-vol', vol_geom);
@@ -203,91 +198,93 @@ end
Resid_error = zeros(iterFISTA,1); % error vector
objective = zeros(iterFISTA,1); % obhective vector
- t = 1;
- X_t = X;
+t = 1;
+X_t = X;
+
+% add_ring = zeros(size(sino),'single'); % size of sinogram array
+r = zeros(Detectors,SlicesZ, 'single'); % 2D array (for 3D data) of sparse "ring" vectors
+r_x = r; % another ring variable
+residual = zeros(size(sino),'single');
+
+% Outer iterations loop
+for i = 1:iterFISTA
- % add_ring = zeros(size(sino),'single'); % size of sinogram array
- r = zeros(Detectors,SlicesZ, 'single'); % 2D array (for 3D data) of sparse "ring" vectors
- r_x = r; % another ring variable
+ X_old = X;
+ t_old = t;
+ r_old = r;
+ [sino_id, sino_updt] = astra_create_sino3d_cuda(X_t, proj_geom, vol_geom);
- % Outer iterations loop
- for i = 1:iterFISTA
-
- X_old = X;
- t_old = t;
- r_old = r;
-
- [sino_id, sino_updt] = astra_create_sino3d_cuda(X_t, proj_geom, vol_geom);
-
- if (lambdaR_L1 > 1)
+ if (lambdaR_L1 > 0)
% add ring removal part (Group-Huber fidelity)
- for kkk = 1:anglesNumb
+ for kkk = 1:anglesNumb
% add_ring(:,kkk,:) = squeeze(sino(:,kkk,:)) - alpha_ring.*r_x;
- residual(:,kkk,:) = weights(:,kkk,:).*(sino_updt(:,kkk,:) - (squeeze(sino(:,kkk,:)) - alpha_ring.*r_x));
- end
-
- vec = sum(residual,2);
- if (SlicesZ > 1)
- vec = squeeze(vec(:,1,:));
- end
- r = r_x - (1./L_const).*vec;
- else
- % no ring removal
- residual = weights.*(sino_updt - sino);
- end
- % residual = weights.*(sino_updt - add_ring);
-
- [id, x_temp] = astra_create_backprojection3d_cuda(residual, proj_geom, vol_geom);
-
- X = X_t - (1/L_const).*x_temp;
- astra_mex_data3d('delete', sino_id);
- astra_mex_data3d('delete', id);
-
- if (lambdaFGP_TV > 0)
- % FGP-TV regularization
- [X, f_val] = FGP_TV(single(X), lambdaFGP_TV, IterationsRegul, tol, 'iso');
- objective(i) = 0.5.*norm(residual(:))^2 + f_val;
- elseif (lambdaSB_TV > 0)
- % Split Bregman regularization
- X = SplitBregman_TV(single(X), lambdaSB_TV, IterationsRegul, tol); % (more memory efficent)
- objective(i) = 0.5.*norm(residual(:))^2;
- elseif (lambdaL1 > 0)
- % L1 soft-threhsolding regularization
- X = max(abs(X)-lambdaL1, 0).*sign(X);
- objective(i) = 0.5.*norm(residual(:))^2;
- elseif (lambdaHO > 0)
- % Higher Order (LLT) regularization
- X2 = LLT_model(single(X), lambdaHO, tauHO, IterationsRegul, 3.0e-05, 0);
- X = 0.5.*(X + X2); % averaged combination of two solutions
- objective(i) = 0.5.*norm(residual(:))^2;
+ residual(:,kkk,:) = squeeze(weights(:,kkk,:)).*(squeeze(sino_updt(:,kkk,:)) - (squeeze(sino(:,kkk,:)) - alpha_ring.*r_x));
end
- if (lambdaR_L1 > 1)
- r = max(abs(r)-lambdaR_L1, 0).*sign(r); % soft-thresholding operator
+ vec = sum(residual,2);
+ if (SlicesZ > 1)
+ vec = squeeze(vec(:,1,:));
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 > 1)
+ r = r_x - (1./L_const).*vec;
+ else
+ % no ring removal
+ residual = weights.*(sino_updt - sino);
+ end
+ % residual = weights.*(sino_updt - add_ring);
+
+ [id, x_temp] = astra_create_backprojection3d_cuda(residual, proj_geom, vol_geom);
+
+ X = X_t - (1/L_const).*x_temp;
+ astra_mex_data3d('delete', sino_id);
+ astra_mex_data3d('delete', id);
+
+ if (lambdaFGP_TV > 0)
+ % FGP-TV regularization
+ [X, f_val] = FGP_TV(single(X), lambdaFGP_TV, IterationsRegul, tol, 'iso');
+ objective(i) = 0.5.*norm(residual(:))^2 + f_val;
+ end
+ if (lambdaSB_TV > 0)
+ % Split Bregman regularization
+ X = SplitBregman_TV(single(X), lambdaSB_TV, IterationsRegul, tol); % (more memory efficent)
+ objective(i) = 0.5.*norm(residual(:))^2;
+ end
+ if (lambdaL1 > 0)
+ % L1 soft-threhsolding regularization
+ X = max(abs(X)-lambdaL1, 0).*sign(X);
+ objective(i) = 0.5.*norm(residual(:))^2;
+ end
+ if (lambdaHO > 0)
+ % Higher Order (LLT) regularization
+ X2 = LLT_model(single(X), lambdaHO, tauHO, iterHO, 3.0e-05, 0);
+ X = 0.5.*(X + X2); % averaged combination of two solutions
+ objective(i) = 0.5.*norm(residual(:))^2;
+ end
+
+ if (lambdaR_L1 > 0)
+ r = max(abs(r)-lambdaR_L1, 0).*sign(r); % soft-thresholding operator
+ 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
-
- if (show == 1)
- figure(10); imshow(X(:,:,slice), [0 maxvalplot]);
- if (lambdaR_L1 > 1)
+ end
+
+ if (show == 1)
+ figure(10); imshow(X(:,:,slice), [0 maxvalplot]);
+ if (lambdaR_L1 > 0)
figure(11); plot(r); title('Rings offset vector')
- end
- pause(0.01);
end
- if (strcmp(X_ideal, 'none' ) == 0)
- Resid_error(i) = RMSE(X(ROI), X_ideal(ROI));
- fprintf('%s %i %s %s %.4f %s %s %.4f \n', 'Iteration Number:', i, '|', 'Error RMSE:', Resid_error(i), '|', 'Objective:', objective(i));
- else
- fprintf('%s %i %s %s %.4f \n', 'Iteration Number:', i, '|', 'Objective:', objective(i));
- end
- end
+ pause(0.01);
+ end
+ if (strcmp(X_ideal, 'none' ) == 0)
+ Resid_error(i) = RMSE(X(ROI), X_ideal(ROI));
+ fprintf('%s %i %s %s %.4f %s %s %.4f \n', 'Iteration Number:', i, '|', 'Error RMSE:', Resid_error(i), '|', 'Objective:', objective(i));
+ else
+ fprintf('%s %i %s %s %.4f \n', 'Iteration Number:', i, '|', 'Objective:', objective(i));
+ end
end
output.Resid_error = Resid_error;
output.objective = objective;
--
cgit v1.2.3