Lipshitz constant attached to the regularization parameter

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'))