OS loop now fixed

3 files changed, 44 insertions, 27 deletions

diff --git a/demos/Demo_Phantom3D_Cone.m b/demos/Demo_Phantom3D_Cone.m
index 6419386..3a9178b 100644
--- a/demos/Demo_Phantom3D_Cone.m
+++ b/demos/Demo_Phantom3D_Cone.m
@@ -17,8 +17,8 @@ angles = 0:1.5:360; % angles vector in degrees
 angles_rad = angles*(pi/180); % conversion to radians
 det_size = round(sqrt(2)*N); % detector size
 % in order to run functions you have to go to the directory:
-cd /home/algol/Documents/MATLAB/TomoPhantom/functions/
-TomoPhantom = buildPhantom3D(modelNo,N); % generate 3D phantom
+pathTP = '/home/algol/Documents/MATLAB/TomoPhantom/functions/models/Phantom3DLibrary.dat'; % path to TomoPhantom parameters file
+TomoPhantom = buildPhantom3D(modelNo,N,pathTP); % generate 3D phantom
 %%
 % using ASTRA-toolbox to set the projection geometry (cone beam)
 % eg: astra.create_proj_geom('cone', 1.0 (resol), 1.0 (resol), detectorRowCount, detectorColCount, angles, originToSource, originToDetector)
diff --git a/demos/Demo_Phantom3D_Parallel.m b/demos/Demo_Phantom3D_Parallel.m
index ac9827c..6a54450 100644
--- a/demos/Demo_Phantom3D_Parallel.m
+++ b/demos/Demo_Phantom3D_Parallel.m
@@ -10,19 +10,35 @@ addpath('../supp/');
 

 %%

 % build 3D phantom using TomoPhantom and generate projection data

-modelNo = 3; % see Phantom3DLibrary.dat file in TomoPhantom

+modelNo = 2; % see Phantom3DLibrary.dat file in TomoPhantom

 N = 256; % x-y-z size (cubic image)

 angles = 1:0.5:180; % angles vector in degrees

 angles_rad = angles*(pi/180); % conversion to radians

 det_size = round(sqrt(2)*N); % detector size

 % in order to run functions you have to go to the directory:

-cd /home/algol/Documents/MATLAB/TomoPhantom/functions/

-TomoPhantom = buildPhantom3D(modelNo,N); % generate 3D phantom

-sino_tomophan3D = buildSino3D(modelNo, N, det_size, single(angles)); % generate data

+pathTP = '/home/algol/Documents/MATLAB/TomoPhantom/functions/models/Phantom3DLibrary.dat'; % path to TomoPhantom parameters file

+TomoPhantom = buildPhantom3D(modelNo,N,pathTP); % generate 3D phantom

+sino_tomophan3D = buildSino3D(modelNo, N, det_size, single(angles),pathTP); % generate ideal data

+% Adding noise and distortions if required

+sino_artifacts = sino_add_artifacts(sino_tomophan3D,'rings');

+%

 %%

 % using ASTRA-toolbox to set the projection geometry (parallel beam)

 proj_geom = astra_create_proj_geom('parallel', 1, det_size, angles_rad);

 vol_geom = astra_create_vol_geom(N,N);

+%%

+fprintf('%s\n', 'Reconstructing with FBP using ASTRA-toolbox ...');

+for i = 1:k

+vol_id = astra_mex_data2d('create', '-vol', vol_geom, 0);

+proj_id = astra_mex_data2d('create', '-proj3d', proj_geom, sino_artifacts(:,:,k)); 

+cfg = astra_struct('FBP_CUDA');

+cfg.ProjectionDataId = proj_id;

+cfg.ReconstructionDataId = vol_id;

+cfg.option.MinConstraint = 0;

+alg_id = astra_mex_algorithm('create', cfg);

+astra_mex_algorithm('iterate', alg_id, 15);

+reconASTRA_3D = astra_mex_data2d('get', vol_id);

+end

 %% 

 fprintf('%s\n', 'Reconstruction using FISTA-LS without regularization...');

 clear params

diff --git a/main_func/FISTA_REC.m b/main_func/FISTA_REC.m
index d177748..bdaeb18 100644
--- a/main_func/FISTA_REC.m
+++ b/main_func/FISTA_REC.m
@@ -487,7 +487,7 @@ else
     % Ordered Subsets (OS) FISTA reconstruction routine (normally one order of magnitude faster than the classical version)
     t = 1;
     X_t = X;
-    proj_geomSUB = proj_geom;    
+    proj_geomSUB = proj_geom;
     
     r = zeros(Detectors,SlicesZ, 'single'); % 2D array (for 3D data) of sparse "ring" vectors
     r_x = r; % another ring variable
@@ -495,7 +495,7 @@ else
     sino_updt_FULL = zeros(size(sino),'single');
     
     % Outer FISTA iterations loop
-    for i = 1:iterFISTA        
+    for i = 1:iterFISTA
         
         if ((i > 1) && (lambdaR_L1 > 0))
             % in order to make Group-Huber fidelity work with ordered subsets
@@ -517,31 +517,31 @@ else
         
         % 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
         for ss = 1:subsets
             X_old = X;
-            t_old = t;         
+            t_old = t;
             
             numProjSub = binsDiscr(ss); % the number of projections per subset
+            sino_updt_Sub = zeros(Detectors, numProjSub, SlicesZ,'single');
             CurrSubIndeces = IndicesReorg(counterInd:(counterInd + numProjSub - 1)); % extract indeces attached to the subset
             proj_geomSUB.ProjectionAngles = angles(CurrSubIndeces);
-            sino_updt_Sub = zeros(Detectors, numProjSub, SlicesZ,'single');
+            
+            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
             
             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);
@@ -551,8 +551,6 @@ else
                 
             elseif (studentt > 0)
                 % student t data fidelity
-                
-                
                 % artifacts removal with Students t penalty
                 residualSub = squeeze(weights(:,CurrSubIndeces,:)).*(sino_updt_Sub - squeeze(sino(:,CurrSubIndeces,:)));
                 
@@ -566,8 +564,11 @@ else
                 objective(i) = ff; % for the objective function output
             else
                 % 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'))
                 % if geometry is 2D use slice-by-slice projection-backprojection routine
                 x_temp = zeros(size(X),'single');
@@ -579,7 +580,7 @@ else
                 astra_mex_data3d('delete', id);
             end
             
-            X = X_t - (1/L_const).*x_temp;            
+            X = X_t - (1/L_const).*x_temp;
             
             % ----------------Regularization part------------------------%
             if (lambdaFGP_TV > 0)