summaryrefslogtreecommitdiffstats
diff options
context:
space:
mode:
Diffstat
-rw-r--r--demos/Demo_Phantom3D_Parallel.m16
-rw-r--r--demos/Demo_RealData3D_Parallel.m26
-rw-r--r--main_func/FISTA_REC.m83
3 files changed, 55 insertions, 70 deletions
diff --git a/demos/Demo_Phantom3D_Parallel.m b/demos/Demo_Phantom3D_Parallel.m
index 402bdd2..4219bd1 100644
--- a/demos/Demo_Phantom3D_Parallel.m
+++ b/demos/Demo_Phantom3D_Parallel.m
@@ -2,7 +2,7 @@
% PARALLEL BEAM geometry
% requirements: ASTRA-toolbox and TomoPhantom toolbox
-close all;clc;clear all;
+close all;clc;clear;
% 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/');
@@ -71,12 +71,12 @@ clear params
params.proj_geom = proj_geom; % pass geometry to the function
params.vol_geom = vol_geom;
params.sino = single(sino3D_log); % sinogram
-params.iterFISTA = 12; %max number of outer iterations
+params.iterFISTA = 15; %max number of outer iterations
params.X_ideal = TomoPhantom; % ideal phantom
params.weights = dataRaw./max(dataRaw(:)); % statistical weight for PWLS
params.subsets = 12; % the number of subsets
params.show = 1; % visualize reconstruction on each iteration
-params.slice = 1; params.maxvalplot = 1.3;
+params.slice = 128; params.maxvalplot = 1.3;
tic; [X_FISTA, output] = FISTA_REC(params); toc;
error_FISTA = output.Resid_error; obj_FISTA = output.objective;
@@ -91,7 +91,7 @@ subplot(1,2,1); plot(error_FISTA); title('RMSE plot');
subplot(1,2,2); plot(obj_FISTA); title('Objective plot');
%%
%%
-fprintf('%s\n', 'Reconstruction using OS-FISTA-GH without FGP-TV regularization...');
+fprintf('%s\n', 'Reconstruction using OS-FISTA-GH with FGP-TV regularization...');
clear params
% define parameters
params.proj_geom = proj_geom; % pass geometry to the function
@@ -99,13 +99,13 @@ params.vol_geom = vol_geom;
params.sino = single(sino3D_log); % sinogram
params.iterFISTA = 15; %max number of outer iterations
params.X_ideal = TomoPhantom; % ideal phantom
-params.weights = dataRaw./max(dataRaw(:)); % statistical weight for PWLS
-params.subsets = 8; % the number of subsets
-params.Regul_Lambda_FGPTV = 0.003; % TV regularization parameter for FGP-TV
+params.weights = dataRaw./max(dataRaw(:)); % statistical weights for PWLS
+params.subsets = 12; % the number of subsets
+params.Regul_Lambda_FGPTV = 100; % TV regularization parameter for FGP-TV
params.Ring_LambdaR_L1 = 0.02; % Soft-Thresh L1 ring variable parameter
params.Ring_Alpha = 21; % to boost ring removal procedure
params.show = 1; % visualize reconstruction on each iteration
-params.slice = 1; params.maxvalplot = 1.3;
+params.slice = 128; params.maxvalplot = 1.3;
tic; [X_FISTA_GH_TV, output] = FISTA_REC(params); toc;
error_FISTA_GH_TV = output.Resid_error; obj_FISTA_GH_TV = output.objective;
diff --git a/demos/Demo_RealData3D_Parallel.m b/demos/Demo_RealData3D_Parallel.m
index e4c9eb0..f82e0b0 100644
--- a/demos/Demo_RealData3D_Parallel.m
+++ b/demos/Demo_RealData3D_Parallel.m
@@ -1,7 +1,7 @@
% Demonstration of tomographic 3D reconstruction from X-ray synchrotron
% dataset (dendrites) using various data fidelities
% ! It is advisable not to run the whole script, it will take lots of time to reconstruct the whole 3D data using many algorithms !
-clear all
+clear
close all
%%
% % adding paths
@@ -44,9 +44,9 @@ clear params
params.proj_geom = proj_geom; % pass geometry to the function
params.vol_geom = vol_geom;
params.sino = Sino3D;
-params.iterFISTA = 12;
+params.iterFISTA = 18;
params.weights = Weights3D;
-params.subsets = 16; % the number of ordered subsets
+params.subsets = 8; % the number of ordered subsets
params.show = 1;
params.maxvalplot = 2.5; params.slice = 1;
@@ -58,12 +58,12 @@ clear params
params.proj_geom = proj_geom; % pass geometry to the function
params.vol_geom = vol_geom;
params.sino = Sino3D;
-params.iterFISTA = 12;
-params.Regul_Lambda_FGPTV = 0.005; % TV regularization parameter for FGP-TV
+params.iterFISTA = 18;
+params.Regul_Lambda_FGPTV = 5.0000e+6; % TV regularization parameter for FGP-TV
params.weights = Weights3D;
-params.subsets = 16; % the number of ordered subsets
+params.subsets = 8; % the number of ordered subsets
params.show = 1;
-params.maxvalplot = 2.5; params.slice = 2;
+params.maxvalplot = 2.5; params.slice = 10;
tic; [X_fista_TV, outputTV] = FISTA_REC(params); toc;
figure; imshow(X_fista_TV(:,:,params.slice) , [0, 2.5]); title ('FISTA-OS-PWLS-TV reconstruction');
@@ -72,13 +72,13 @@ fprintf('%s\n', 'Reconstruction using FISTA-OS-GH-TV...');
clear params
params.proj_geom = proj_geom; % pass geometry to the function
params.vol_geom = vol_geom;
-params.sino = Sino3D;
-params.iterFISTA = 12;
-% params.Regul_Lambda_FGPTV = 0.005; % TV regularization parameter for FGP-TV
+params.sino = Sino3D(:,:,10);
+params.iterFISTA = 18;
+params.Regul_Lambda_FGPTV = 5.0000e+6; % TV regularization parameter for FGP-TV
params.Ring_LambdaR_L1 = 0.002; % Soft-Thresh L1 ring variable parameter
params.Ring_Alpha = 21; % to boost ring removal procedure
-params.weights = Weights3D;
-params.subsets = 16; % the number of ordered subsets
+params.weights = Weights3D(:,:,10);
+params.subsets = 8; % the number of ordered subsets
params.show = 1;
params.maxvalplot = 2.5; params.slice = 1;
@@ -91,7 +91,7 @@ params.proj_geom = proj_geom; % pass geometry to the function
params.vol_geom = vol_geom;
params.sino = Sino3D;
params.iterFISTA = 12;
-params.Regul_Lambda_FGPTV = 0.005; % TV regularization parameter for FGP-TV
+params.Regul_Lambda_FGPTV = 5.0000e+6; % TV regularization parameter for FGP-TV
params.Regul_LambdaLLT = 100; % regularization parameter for LLT problem
params.Regul_tauLLT = 0.0005; % time-step parameter for the explicit scheme
params.Ring_LambdaR_L1 = 0.002; % Soft-Thresh L1 ring variable parameter
diff --git a/main_func/FISTA_REC.m b/main_func/FISTA_REC.m
index 3d22b97..d717a03 100644
--- a/main_func/FISTA_REC.m
+++ b/main_func/FISTA_REC.m
@@ -169,12 +169,12 @@ end
if (isfield(params,'Regul_tol'))
tol = params.Regul_tol;
else
- tol = 1.0e-04;
+ tol = 1.0e-05;
end
if (isfield(params,'Regul_Iterations'))
IterationsRegul = params.Regul_Iterations;
else
- IterationsRegul = 25;
+ IterationsRegul = 45;
end
if (isfield(params,'Regul_LambdaLLT'))
lambdaHO = params.Regul_LambdaLLT;
@@ -381,11 +381,11 @@ if (subsets == 0)
if ((strcmp('2D', Dimension) == 1))
% 2D regularization
for kkk = 1:SlicesZ
- [X(:,:,kkk), f_val] = FGP_TV(single(X(:,:,kkk)), lambdaFGP_TV, IterationsRegul, tol, 'iso');
+ [X(:,:,kkk), f_val] = FGP_TV(single(X(:,:,kkk)), lambdaFGP_TV/L_const, IterationsRegul, tol, 'iso');
end
else
% 3D regularization
- [X, f_val] = FGP_TV(single(X), lambdaFGP_TV, IterationsRegul, tol, 'iso');
+ [X, f_val] = FGP_TV(single(X), lambdaFGP_TV/L_const, IterationsRegul, tol, 'iso');
end
objective(i) = (objective(i) + f_val)./(Detectors*anglesNumb*SlicesZ);
end
@@ -394,11 +394,11 @@ if (subsets == 0)
if ((strcmp('2D', Dimension) == 1))
% 2D regularization
for kkk = 1:SlicesZ
- X(:,:,kkk) = SplitBregman_TV(single(X(:,:,kkk)), lambdaSB_TV, IterationsRegul, tol); % (more memory efficent)
+ X(:,:,kkk) = SplitBregman_TV(single(X(:,:,kkk)), lambdaSB_TV/L_const, IterationsRegul, tol); % (more memory efficent)
end
else
% 3D regularization
- X = SplitBregman_TV(single(X), lambdaSB_TV, IterationsRegul, tol); % (more memory efficent)
+ X = SplitBregman_TV(single(X), lambdaSB_TV/L_const, IterationsRegul, tol); % (more memory efficent)
end
end
if (lambdaHO > 0)
@@ -407,11 +407,11 @@ if (subsets == 0)
if ((strcmp('2D', Dimension) == 1))
% 2D regularization
for kkk = 1:SlicesZ
- X2(:,:,kkk) = LLT_model(single(X(:,:,kkk)), lambdaHO, tauHO, iterHO, 3.0e-05, 0);
+ X2(:,:,kkk) = LLT_model(single(X(:,:,kkk)), lambdaHO/L_const, tauHO, iterHO, 3.0e-05, 0);
end
else
% 3D regularization
- X2 = LLT_model(single(X), lambdaHO, tauHO, iterHO, 3.0e-05, 0);
+ X2 = LLT_model(single(X), lambdaHO/L_const, tauHO, iterHO, 3.0e-05, 0);
end
X = 0.5.*(X + X2); % averaged combination of two solutions
@@ -421,10 +421,10 @@ if (subsets == 0)
if ((strcmp('2D', Dimension) == 1))
% 2D regularization
for kkk = 1:SlicesZ
- X(:,:,kkk) = PatchBased_Regul(single(X(:,:,kkk)), SearchW, SimilW, h_PB, lambdaPB);
+ X(:,:,kkk) = PatchBased_Regul(single(X(:,:,kkk)), SearchW, SimilW, h_PB, lambdaPB/L_const);
end
else
- X = PatchBased_Regul(single(X), SearchW, SimilW, h_PB, lambdaPB);
+ X = PatchBased_Regul(single(X), SearchW, SimilW, h_PB, lambdaPB/L_const);
end
end
if (lambdaPB_GPU > 0)
@@ -432,10 +432,10 @@ if (subsets == 0)
if ((strcmp('2D', Dimension) == 1))
% 2D regularization
for kkk = 1:SlicesZ
- X(:,:,kkk) = NLM_GPU(single(X(:,:,kkk)), SearchW, SimilW, h_PB, lambdaPB_GPU);
+ X(:,:,kkk) = NLM_GPU(single(X(:,:,kkk)), SearchW, SimilW, h_PB, lambdaPB_GPU/L_const);
end
else
- X = NLM_GPU(single(X), SearchW, SimilW, h_PB, lambdaPB_GPU);
+ X = NLM_GPU(single(X), SearchW, SimilW, h_PB, lambdaPB_GPU/L_const);
end
end
if (LambdaDiff_HO > 0)
@@ -443,10 +443,10 @@ if (subsets == 0)
if ((strcmp('2D', Dimension) == 1))
% 2D regularization
for kkk = 1:SlicesZ
- X(:,:,kkk) = Diff4thHajiaboli_GPU(single(X(:,:,kkk)), LambdaDiff_HO_EdgePar, LambdaDiff_HO, IterationsRegul);
+ X(:,:,kkk) = Diff4thHajiaboli_GPU(single(X(:,:,kkk)), LambdaDiff_HO_EdgePar, LambdaDiff_HO/L_const, IterationsRegul);
end
else
- X = Diff4thHajiaboli_GPU(X, LambdaDiff_HO_EdgePar, LambdaDiff_HO, IterationsRegul);
+ X = Diff4thHajiaboli_GPU(X, LambdaDiff_HO_EdgePar, LambdaDiff_HO/L_const, IterationsRegul);
end
end
if (LambdaTGV > 0)
@@ -454,7 +454,7 @@ if (subsets == 0)
lamTGV1 = 1.1; % smoothing trade-off parameters, see Pock's paper
lamTGV2 = 0.8; % second-order term
for kkk = 1:SlicesZ
- X(:,:,kkk) = TGV_PD(single(X(:,:,kkk)), LambdaTGV, lamTGV1, lamTGV2, IterationsRegul);
+ X(:,:,kkk) = TGV_PD(single(X(:,:,kkk)), LambdaTGV/L_const, lamTGV1, lamTGV2, IterationsRegul);
end
end
@@ -494,6 +494,7 @@ else
residual2 = zeros(size(sino),'single');
sino_updt_FULL = zeros(size(sino),'single');
+
% Outer FISTA iterations loop
for i = 1:iterFISTA
@@ -514,21 +515,9 @@ else
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
+
+ % subsets loop
+ counterInd = 1;
for ss = 1:subsets
X_old = X;
t_old = t;
@@ -553,8 +542,6 @@ else
if (lambdaR_L1 > 0)
% Group-Huber fidelity (ring removal)
-
-
residualSub = zeros(Detectors, numProjSub, SlicesZ,'single'); % residual for a chosen subset
for kkk = 1:numProjSub
indC = CurrSubIndeces(kkk);
@@ -564,7 +551,7 @@ else
elseif (studentt > 0)
% student t data fidelity
-
+
% artifacts removal with Students t penalty
residualSub = squeeze(weights(:,CurrSubIndeces,:)).*(sino_updt_Sub - squeeze(sino(:,CurrSubIndeces,:)));
@@ -577,12 +564,10 @@ else
end
objective(i) = ff; % for the objective function output
else
- % PWLS model
-
+ % PWLS model
residualSub = squeeze(weights(:,CurrSubIndeces,:)).*(sino_updt_Sub - squeeze(sino(:,CurrSubIndeces,:)));
objective(i) = 0.5*norm(residualSub(:)); % for the objective function output
end
-
% perform backprojection of a subset
if (strcmp(proj_geom.type,'parallel') || strcmp(proj_geom.type,'fanflat') || strcmp(proj_geom.type,'fanflat_vec'))
@@ -604,11 +589,11 @@ else
if ((strcmp('2D', Dimension) == 1))
% 2D regularization
for kkk = 1:SlicesZ
- [X(:,:,kkk), f_val] = FGP_TV(single(X(:,:,kkk)), lambdaFGP_TV/subsets, IterationsRegul, tol, 'iso');
+ [X(:,:,kkk), f_val] = FGP_TV(single(X(:,:,kkk)), lambdaFGP_TV/(subsets*L_const), IterationsRegul, tol, 'iso');
end
else
% 3D regularization
- [X, f_val] = FGP_TV(single(X), lambdaFGP_TV/subsets, IterationsRegul, tol, 'iso');
+ [X, f_val] = FGP_TV(single(X), lambdaFGP_TV/(subsets*L_const), IterationsRegul, tol, 'iso');
end
objective(i) = objective(i) + f_val;
end
@@ -617,11 +602,11 @@ else
if ((strcmp('2D', Dimension) == 1))
% 2D regularization
for kkk = 1:SlicesZ
- X(:,:,kkk) = SplitBregman_TV(single(X(:,:,kkk)), lambdaSB_TV/subsets, IterationsRegul, tol); % (more memory efficent)
+ X(:,:,kkk) = SplitBregman_TV(single(X(:,:,kkk)), lambdaSB_TV/(subsets*L_const), IterationsRegul, tol); % (more memory efficent)
end
else
% 3D regularization
- X = SplitBregman_TV(single(X), lambdaSB_TV/subsets, IterationsRegul, tol); % (more memory efficent)
+ X = SplitBregman_TV(single(X), lambdaSB_TV/(subsets*L_const), IterationsRegul, tol); % (more memory efficent)
end
end
if (lambdaHO > 0)
@@ -630,11 +615,11 @@ else
if ((strcmp('2D', Dimension) == 1))
% 2D regularization
for kkk = 1:SlicesZ
- X2(:,:,kkk) = LLT_model(single(X(:,:,kkk)), lambdaHO/subsets, tauHO/subsets, iterHO, 2.0e-05, 0);
+ X2(:,:,kkk) = LLT_model(single(X(:,:,kkk)), lambdaHO/(subsets*L_const), tauHO/subsets, iterHO, 2.0e-05, 0);
end
else
% 3D regularization
- X2 = LLT_model(single(X), lambdaHO/subsets, tauHO/subsets, iterHO, 2.0e-05, 0);
+ X2 = LLT_model(single(X), lambdaHO/(subsets*L_const), tauHO/subsets, iterHO, 2.0e-05, 0);
end
X = 0.5.*(X + X2); % the averaged combination of two solutions
end
@@ -643,10 +628,10 @@ else
if ((strcmp('2D', Dimension) == 1))
% 2D regularization
for kkk = 1:SlicesZ
- X(:,:,kkk) = PatchBased_Regul(single(X(:,:,kkk)), SearchW, SimilW, h_PB, lambdaPB/subsets);
+ X(:,:,kkk) = PatchBased_Regul(single(X(:,:,kkk)), SearchW, SimilW, h_PB, lambdaPB/(subsets*L_const));
end
else
- X = PatchBased_Regul(single(X), SearchW, SimilW, h_PB, lambdaPB/subsets);
+ X = PatchBased_Regul(single(X), SearchW, SimilW, h_PB, lambdaPB/(subsets*L_const));
end
end
if (lambdaPB_GPU > 0)
@@ -654,10 +639,10 @@ else
if ((strcmp('2D', Dimension) == 1))
% 2D regularization
for kkk = 1:SlicesZ
- X(:,:,kkk) = NLM_GPU(single(X(:,:,kkk)), SearchW, SimilW, h_PB, lambdaPB_GPU);
+ X(:,:,kkk) = NLM_GPU(single(X(:,:,kkk)), SearchW, SimilW, h_PB, lambdaPB_GPU/(subsets*L_const));
end
else
- X = NLM_GPU(single(X), SearchW, SimilW, h_PB, lambdaPB_GPU);
+ X = NLM_GPU(single(X), SearchW, SimilW, h_PB, lambdaPB_GPU/(subsets*L_const));
end
end
if (LambdaDiff_HO > 0)
@@ -665,10 +650,10 @@ else
if ((strcmp('2D', Dimension) == 1))
% 2D regularization
for kkk = 1:SlicesZ
- X(:,:,kkk) = Diff4thHajiaboli_GPU(single(X(:,:,kkk)), LambdaDiff_HO_EdgePar, LambdaDiff_HO, round(IterationsRegul/subsets));
+ X(:,:,kkk) = Diff4thHajiaboli_GPU(single(X(:,:,kkk)), LambdaDiff_HO_EdgePar, LambdaDiff_HO/(subsets*L_const), round(IterationsRegul/subsets));
end
else
- X = Diff4thHajiaboli_GPU(X, LambdaDiff_HO_EdgePar, LambdaDiff_HO, round(IterationsRegul/subsets));
+ X = Diff4thHajiaboli_GPU(X, LambdaDiff_HO_EdgePar, LambdaDiff_HO/(subsets*L_const), round(IterationsRegul/subsets));
end
end
if (LambdaTGV > 0)
@@ -676,7 +661,7 @@ else
lamTGV1 = 1.1; % smoothing trade-off parameters, see Pock's paper
lamTGV2 = 0.5; % second-order term
for kkk = 1:SlicesZ
- X(:,:,kkk) = TGV_PD(single(X(:,:,kkk)), LambdaTGV/subsets, lamTGV1, lamTGV2, IterationsRegul);
+ X(:,:,kkk) = TGV_PD(single(X(:,:,kkk)), LambdaTGV/(subsets*L_const), lamTGV1, lamTGV2, IterationsRegul);
end
end
generated by cgit v1.2.3 (git 2.25.1) at 2024-11-07 21:40:55 +0000