From a4ee6ded8036fa1a10ef860478f49d2eb2cd11e0 Mon Sep 17 00:00:00 2001 From: algol Date: Thu, 26 Oct 2017 15:28:34 +0100 Subject: Lipshitz constant attached to the regularization parameter --- demos/Demo_Phantom3D_Parallel.m | 14 +++++----- demos/Demo_RealData3D_Parallel.m | 22 ++++++++-------- main_func/FISTA_REC.m | 56 ++++++++++++++++++++-------------------- 3 files changed, 46 insertions(+), 46 deletions(-) diff --git a/demos/Demo_Phantom3D_Parallel.m b/demos/Demo_Phantom3D_Parallel.m index 402bdd2..b0737f1 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/'); @@ -76,7 +76,7 @@ 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..d680aca 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'); @@ -73,12 +73,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.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.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 fa98360..658d2a9 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 @@ -588,11 +588,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 @@ -601,11 +601,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) @@ -614,11 +614,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 @@ -627,10 +627,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) @@ -638,10 +638,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) @@ -649,10 +649,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) @@ -660,7 +660,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 -- cgit v1.2.3 From 09f9bf9828c39bcdd870cfefbcb52e61451802eb Mon Sep 17 00:00:00 2001 From: algol Date: Thu, 26 Oct 2017 15:45:36 +0100 Subject: Lipshitz constant attached to the regularization parameter --- demos/Demo_Phantom3D_Parallel.m | 2 +- demos/Demo_RealData3D_Parallel.m | 4 ++-- main_func/FISTA_REC.m | 27 ++++++--------------------- 3 files changed, 9 insertions(+), 24 deletions(-) diff --git a/demos/Demo_Phantom3D_Parallel.m b/demos/Demo_Phantom3D_Parallel.m index b0737f1..4219bd1 100644 --- a/demos/Demo_Phantom3D_Parallel.m +++ b/demos/Demo_Phantom3D_Parallel.m @@ -71,7 +71,7 @@ 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 diff --git a/demos/Demo_RealData3D_Parallel.m b/demos/Demo_RealData3D_Parallel.m index d680aca..f82e0b0 100644 --- a/demos/Demo_RealData3D_Parallel.m +++ b/demos/Demo_RealData3D_Parallel.m @@ -72,12 +72,12 @@ 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.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.weights = Weights3D(:,:,10); params.subsets = 8; % the number of ordered subsets params.show = 1; params.maxvalplot = 2.5; params.slice = 1; diff --git a/main_func/FISTA_REC.m b/main_func/FISTA_REC.m index 202ebc2..d717a03 100644 --- a/main_func/FISTA_REC.m +++ b/main_func/FISTA_REC.m @@ -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')) -- cgit v1.2.3