From 29636540aca6354b0f319765b0af9bf768593565 Mon Sep 17 00:00:00 2001 From: algol 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