ring removal bug fixed

2 files changed, 80 insertions, 83 deletions

diff --git a/demos/DemoRD2.m b/demos/DemoRD2.m
index b991e70..0f829a8 100644
--- a/demos/DemoRD2.m
+++ b/demos/DemoRD2.m
@@ -63,7 +63,7 @@ params.L_const = 7.6789e+08;
 params.Regul_Lambda_FGPTV = 0.005; % TV regularization parameter for FGP-TV
 params.weights = Weights3D;
 params.show = 1;
-params.maxvalplot = 2.5; params.slice = 10;
+params.maxvalplot = 2.5; params.slice = 4;
 
 tic; [X_fista_TV] = FISTA_REC(params); toc;
 figure; imshow(X_fista_TV(:,:,1) , [0, 2.5]); title ('FISTA-PWLS-TV reconstruction');
@@ -81,7 +81,7 @@ 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.show = 1;
-params.maxvalplot = 2.5; params.slice = 10;
+params.maxvalplot = 2.5; params.slice = 4;
 
 tic; [X_fista_GH_TV] = FISTA_REC(params); toc;
 figure; imshow(X_fista_GH_TV(:,:,1) , [0, 2.5]); title ('FISTA-GH-TV reconstruction');
@@ -94,7 +94,7 @@ params.vol_geom = vol_geom;
 params.sino = Sino3D;
 params.iterFISTA  = 40;
 params.Regul_Lambda_FGPTV = 0.005; % TV regularization parameter for FGP-TV
-params.Regul_LambdaHO = 200;  % regularization parameter for LLT problem
+params.Regul_LambdaHO = 250;  % 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
 params.Ring_Alpha = 21; % to boost ring removal procedure
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;