Lipshitz constant attached to the regularization parameter

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