Merge branch 'pythonize' of https://github.com/vais-ral/CCPi-FISTA_Reconstruction into pythonize

3 files changed, 594 insertions, 89 deletions

diff --git a/src/Python/ccpi/fista/FISTAReconstructor.py b/src/Python/ccpi/fista/FISTAReconstructor.py
index 87dd2c0..fda9cf0 100644
--- a/src/Python/ccpi/fista/FISTAReconstructor.py
+++ b/src/Python/ccpi/fista/FISTAReconstructor.py
@@ -85,6 +85,7 @@ class FISTAReconstructor():
         self.pars['detectors'] = detectors
         self.pars['number_of_angles'] = nangles
         self.pars['SlicesZ'] = sliceZ
+        self.pars['output_volume'] = None
 
         print (self.pars)
         # handle optional input parameters (at instantiation)
@@ -109,8 +110,10 @@ class FISTAReconstructor():
               'ring_lambda_R_L1',
               'ring_alpha',
               'subsets',
-              'use_studentt_fidelity',
-              'studentt')
+              'output_volume',
+              'os_subsets',
+              'os_indices',
+              'os_bins')
         self.acceptedInputKeywords = list(kw)
         
         # handle keyworded parameters
@@ -141,10 +144,12 @@ class FISTAReconstructor():
         
         if not 'region_of_interest'in kwargs.keys() :
             if self.pars['ideal_image'] == None:
-                pass
+                self.pars['region_of_interest'] = None
             else:
-                self.pars['region_of_interest'] = numpy.nonzero(
-                    self.pars['ideal_image']>0.0)
+                ## nonzero if the image is larger than m
+                fsm = numpy.frompyfunc(lambda x,m: 1 if x>m else 0, 2,1)
+                
+                self.pars['region_of_interest'] = fsm(self.pars['ideal_image'], 0)
                 
         # the regularizer must be a correctly instantiated object    
         if not 'regularizer' in kwargs.keys() :
@@ -165,14 +170,7 @@ class FISTAReconstructor():
         if not 'initialize' in kwargs.keys():
             self.pars['initialize'] = False
 
-        if not 'use_studentt_fidelity' in kwargs.keys():
-            self.setParameter(studentt=False)
-        else:
-            print ("studentt {0}".format(kwargs['use_studentt_fidelity']))
-            if kwargs['use_studentt_fidelity']:
-                raise Exception('Not implemented')
-            
-            self.setParameter(studentt=kwargs['use_studentt_fidelity'])
+        
             
             
     def setParameter(self, **kwargs):
@@ -183,8 +181,6 @@ class FISTAReconstructor():
         '''
         for key , value in kwargs.items():
             if key in self.acceptedInputKeywords:
-                if key == 'use_studentt_fidelity':
-                    raise Exception('use_studentt_fidelity Not implemented')
                 self.pars[key] = value
             else:
                 raise Exception('Wrong parameter {0} for '.format(key) +
@@ -389,11 +385,15 @@ class FISTAReconstructor():
                     
                 counter = counter + binsDiscr[jj] - 1    
                 
-            
-        return IndicesReorg
+        # store the OS in parameters
+        self.setParameter(os_subsets=subsets,
+                          os_bins=binsDiscr,
+                          os_indices=IndicesReorg)
             
 
     def prepareForIteration(self):
+        print ("FISTA Reconstructor: prepare for iteration")
+        
         self.residual_error = numpy.zeros((self.pars['number_of_iterations']))
         self.objective = numpy.zeros((self.pars['number_of_iterations']))
 
@@ -408,19 +408,17 @@ class FISTAReconstructor():
         if self.getParameter('Lipschitz_constant') is None:
             self.pars['Lipschitz_constant'] = \
                             self.calculateLipschitzConstantWithPowerMethod()
+        # errors vector (if the ground truth is given)
+        self.Resid_error = numpy.zeros((self.getParameter('number_of_iterations')));
+        # objective function values vector
+        self.objective = numpy.zeros((self.getParameter('number_of_iterations')));      
         
 
     # prepareForIteration
 
     def iterate(self, Xin=None):
-        # convenience variable storage
-        proj_geom , vol_geom, sino , \
-                  SlicesZ = self.getParameter([ 'projector_geometry' ,
-                                                'output_geometry',
-                                                'input_sinogram',
-                                                'SlicesZ'])
-                        
-        t = 1
+        print ("FISTA Reconstructor: iterate")
+        
         if Xin is None:    
             if self.getParameter('initialize'):
                 X = self.initialize()
@@ -430,15 +428,25 @@ class FISTAReconstructor():
         else:
             # copy by reference
             X = Xin
-
+        # store the output volume in the parameters
+        self.setParameter(output_volume=X)
         X_t = X.copy()
+        # convenience variable storage
+        proj_geom , vol_geom, sino , \
+          SlicesZ  = self.getParameter([ 'projector_geometry' ,
+                                                'output_geometry',
+                                                'input_sinogram',
+                                                'SlicesZ' ])
+                   
+        t = 1
         
         for i in range(self.getParameter('number_of_iterations')):
             X_old = X.copy()
             t_old = t
             r_old = self.r.copy()
             if self.getParameter('projector_geometry')['type'] == 'parallel' or \
-               self.getParameter('projector_geometry')['type'] == 'parallel3d':
+               self.getParameter('projector_geometry')['type'] == 'fanflat' or \
+               self.getParameter('projector_geometry')['type'] == 'fanflat_vec':
                 # if the geometry is parallel use slice-by-slice
                 # projection-backprojection routine
                 #sino_updt = zeros(size(sino),'single');
@@ -446,10 +454,9 @@ class FISTAReconstructor():
                 proj_geomT['DetectorRowCount'] = 1
                 vol_geomT = vol_geom.copy()
                 vol_geomT['GridSliceCount'] = 1;
-                sino_updt = numpy.zeros(numpy.shape(sino), dtype=numpy.float)
+                self.sino_updt = numpy.zeros(numpy.shape(sino), dtype=numpy.float)
                 for kkk in range(SlicesZ):
-                    print (kkk)
-                    sino_id, sino_updt[kkk] = \
+                    sino_id, self.sino_updt[kkk] = \
                              astra.creators.create_sino3d_gpu(
                                  X_t[kkk:kkk+1], proj_geomT, vol_geomT)
                     astra.matlab.data3d('delete', sino_id)
@@ -457,11 +464,122 @@ class FISTAReconstructor():
                 # for divergent 3D geometry (watch the GPU memory overflow in
                 # ASTRA versions < 1.8)
                 #[sino_id, sino_updt] = astra_create_sino3d_cuda(X_t, proj_geom, vol_geom);
-                sino_id, sino_updt = astra.matlab.create_sino3d_gpu(
+                sino_id, self.sino_updt = astra.creators.create_sino3d_gpu(
                     X_t, proj_geom, vol_geom)
 
 
             ## RING REMOVAL
-
+            self.ringRemoval(i)
+            ## Projection/Backprojection Routine
+            self.projectionBackprojection(X, X_t)
+            astra.matlab.data3d('delete', sino_id)
             ## REGULARIZATION
+            X = self.regularize(X)
+            ## Update Loop
+            X , X_t, t = self.updateLoop(i, X, X_old, r_old, t, t_old)
+            self.setParameter(output_volume=X)
+        return X
+    ## iterate
+    
+    def ringRemoval(self, i):
+        print ("FISTA Reconstructor: ring removal")
+        residual = self.residual
+        lambdaR_L1 , alpha_ring , weights , L_const , sino= \
+                   self.getParameter(['ring_lambda_R_L1',
+                                      'ring_alpha' , 'weights',
+                                      'Lipschitz_constant',
+                                      'input_sinogram'])
+        r_x = self.r_x
+        sino_updt = self.sino_updt
+        
+        SlicesZ, anglesNumb, Detectors = \
+                    numpy.shape(self.getParameter('input_sinogram'))
+        if lambdaR_L1 > 0 :
+             for kkk in range(anglesNumb):
+                 
+                 residual[:,kkk,:] = (weights[:,kkk,:]).squeeze() * \
+                                       ((sino_updt[:,kkk,:]).squeeze() - \
+                                        (sino[:,kkk,:]).squeeze() -\
+                                        (alpha_ring * r_x)
+                                        )
+             vec = residual.sum(axis = 1)
+             #if SlicesZ > 1:
+             #    vec = vec[:,1,:].squeeze()
+             self.r = (r_x - (1./L_const) * vec).copy()
+             self.objective[i] = (0.5 * (residual ** 2).sum())
+
+    def projectionBackprojection(self, X, X_t):
+        print ("FISTA Reconstructor: projection-backprojection routine")
+        
+        # a few useful variables
+        SlicesZ, anglesNumb, Detectors = \
+                    numpy.shape(self.getParameter('input_sinogram'))
+        residual = self.residual
+        proj_geom , vol_geom , L_const = \
+                  self.getParameter(['projector_geometry' ,
+                                                  'output_geometry',
+                                                  'Lipschitz_constant'])
+        
+        
+        if self.getParameter('projector_geometry')['type'] == 'parallel' or \
+           self.getParameter('projector_geometry')['type'] == 'fanflat' or \
+           self.getParameter('projector_geometry')['type'] == 'fanflat_vec':
+            # if the geometry is parallel use slice-by-slice
+            # projection-backprojection routine
+            #sino_updt = zeros(size(sino),'single');
+            proj_geomT = proj_geom.copy()
+            proj_geomT['DetectorRowCount'] = 1
+            vol_geomT = vol_geom.copy()
+            vol_geomT['GridSliceCount'] = 1;
+            x_temp = numpy.zeros(numpy.shape(X),dtype=numpy.float32)
+            
+            for kkk in range(SlicesZ):
+                
+                x_id, x_temp[kkk] = \
+                         astra.creators.create_backprojection3d_gpu(
+                             residual[kkk:kkk+1],
+                             proj_geomT, vol_geomT)
+                astra.matlab.data3d('delete', x_id)
+        else:
+            x_id, x_temp = \
+                  astra.creators.create_backprojection3d_gpu(
+                      residual, proj_geom, vol_geom)            
+
+        X = X_t - (1/L_const) * x_temp
+        #astra.matlab.data3d('delete', sino_id)
+        astra.matlab.data3d('delete', x_id)
+
+    def regularize(self, X):
+        print ("FISTA Reconstructor: regularize")
+        
+        regularizer = self.getParameter('regularizer')
+        if regularizer is not None:
+            return regularizer(input=X)
+        else:
+            return X
+
+    def updateLoop(self, i, X, X_old, r_old, t, t_old):
+        print ("FISTA Reconstructor: update loop")
+        lambdaR_L1 = self.getParameter('ring_lambda_R_L1')
+        if lambdaR_L1 > 0:
+            self.r = numpy.max(
+                numpy.abs(self.r) - lambdaR_L1 , 0) * \
+                numpy.sign(self.r)
+        t = (1 + numpy.sqrt(1 + 4 * t**2))/2
+        X_t = X + (((t_old -1)/t) * (X - X_old))
+
+        if lambdaR_L1 > 0:
+            self.r_x = self.r + \
+                             (((t_old-1)/t) * (self.r - r_old))
+
+        if self.getParameter('region_of_interest') is None:
+            string = 'Iteration Number {0} | Objective {1} \n'
+            print (string.format( i, self.objective[i]))
+        else:
+            ROI , X_ideal = fistaRecon.getParameter('region_of_interest',
+                                                    'ideal_image')
             
+            Resid_error[i] = RMSE(X*ROI, X_ideal*ROI)
+            string = 'Iteration Number {0} | RMS Error {1} | Objective {2} \n'
+            print (string.format(i,Resid_error[i], self.objective[i]))
+        return (X , X_t, t)
diff --git a/src/Python/test_reconstructor-os.py b/src/Python/test_reconstructor-os.py
new file mode 100644
index 0000000..6f3721f
--- /dev/null
+++ b/src/Python/test_reconstructor-os.py
@@ -0,0 +1,379 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Wed Aug 23 16:34:49 2017
+
+@author: ofn77899
+Based on DemoRD2.m
+"""
+
+import h5py
+import numpy
+
+from ccpi.fista.FISTAReconstructor import FISTAReconstructor
+import astra
+import matplotlib.pyplot as plt
+
+def RMSE(signal1, signal2):
+    '''RMSE Root Mean Squared Error'''
+    if numpy.shape(signal1) == numpy.shape(signal2):
+        err = (signal1 - signal2)
+        err = numpy.sum( err * err )/numpy.size(signal1);  # MSE
+        err = sqrt(err);                                   # RMSE
+        return err
+    else:
+        raise Exception('Input signals must have the same shape')
+  
+filename = r'/home/ofn77899/Reconstruction/CCPi-FISTA_Reconstruction/demos/DendrData.h5'
+nx = h5py.File(filename, "r")
+#getEntry(nx, '/')
+# I have exported the entries as children of /
+entries = [entry for entry in nx['/'].keys()]
+print (entries)
+
+Sino3D = numpy.asarray(nx.get('/Sino3D'), dtype="float32")
+Weights3D = numpy.asarray(nx.get('/Weights3D'), dtype="float32")
+angSize = numpy.asarray(nx.get('/angSize'), dtype=int)[0]
+angles_rad = numpy.asarray(nx.get('/angles_rad'), dtype="float32")
+recon_size = numpy.asarray(nx.get('/recon_size'), dtype=int)[0]
+size_det = numpy.asarray(nx.get('/size_det'), dtype=int)[0]
+slices_tot = numpy.asarray(nx.get('/slices_tot'), dtype=int)[0]
+
+Z_slices = 20
+det_row_count = Z_slices
+# next definition is just for consistency of naming
+det_col_count = size_det
+
+detectorSpacingX = 1.0
+detectorSpacingY = detectorSpacingX
+
+
+proj_geom = astra.creators.create_proj_geom('parallel3d',
+                                            detectorSpacingX,
+                                            detectorSpacingY,
+                                            det_row_count,
+                                            det_col_count,
+                                            angles_rad)
+
+#vol_geom = astra_create_vol_geom(recon_size,recon_size,Z_slices);
+image_size_x = recon_size
+image_size_y = recon_size
+image_size_z = Z_slices
+vol_geom = astra.creators.create_vol_geom( image_size_x,
+                                           image_size_y,
+                                           image_size_z)
+
+## First pass the arguments to the FISTAReconstructor and test the
+## Lipschitz constant
+
+fistaRecon = FISTAReconstructor(proj_geom,
+                                vol_geom,
+                                Sino3D ,
+                                weights=Weights3D)
+
+print ("Lipschitz Constant {0}".format(fistaRecon.pars['Lipschitz_constant']))
+fistaRecon.setParameter(number_of_iterations = 12)
+fistaRecon.setParameter(Lipschitz_constant = 767893952.0)
+fistaRecon.setParameter(ring_alpha = 21)
+fistaRecon.setParameter(ring_lambda_R_L1 = 0.002)
+
+## Ordered subset
+if True:
+    subsets = 16
+    angles = fistaRecon.getParameter('projector_geometry')['ProjectionAngles']
+    #binEdges = numpy.linspace(angles.min(),
+    #                          angles.max(),
+    #                          subsets + 1)
+    binsDiscr, binEdges = numpy.histogram(angles, bins=subsets)
+    # get rearranged subset indices
+    IndicesReorg = numpy.zeros((numpy.shape(angles)))
+    counterM = 0
+    for ii in range(binsDiscr.max()):
+        counter = 0
+        for jj in range(subsets):
+            curr_index = ii + jj  + counter
+            #print ("{0} {1} {2}".format(binsDiscr[jj] , ii, counterM))
+            if binsDiscr[jj] > ii:
+                if (counterM < numpy.size(IndicesReorg)):
+                    IndicesReorg[counterM] = curr_index
+                counterM = counterM + 1
+                
+            counter = counter + binsDiscr[jj] - 1
+
+
+if True:
+    print ("Lipschitz Constant {0}".format(fistaRecon.pars['Lipschitz_constant']))
+    print ("prepare for iteration")
+    fistaRecon.prepareForIteration()
+    
+    
+
+    print("initializing  ...")
+    if False:
+        # if X doesn't exist
+        #N = params.vol_geom.GridColCount
+        N = vol_geom['GridColCount']
+        print ("N " + str(N))
+        X = numpy.zeros((N,N,SlicesZ), dtype=numpy.float)
+    else:
+        #X = fistaRecon.initialize()
+        X = numpy.load("X.npy")
+    
+    print (numpy.shape(X))
+    X_t = X.copy()
+    print ("initialized")
+    proj_geom , vol_geom, sino , \
+                      SlicesZ = fistaRecon.getParameter(['projector_geometry' ,
+                                                            'output_geometry',
+                                                            'input_sinogram',
+                                                         'SlicesZ'])
+
+    #fistaRecon.setParameter(number_of_iterations = 3)
+    iterFISTA = fistaRecon.getParameter('number_of_iterations')
+    # errors vector (if the ground truth is given)
+    Resid_error = numpy.zeros((iterFISTA));
+    # objective function values vector
+    objective = numpy.zeros((iterFISTA)); 
+
+      
+    t = 1
+
+    ## additional for 
+    proj_geomSUB = proj_geom.copy()
+    fistaRecon.residual2 = numpy.zeros(numpy.shape(self.pars['input_sinogram']))
+    print ("starting iterations")
+##    % Outer FISTA iterations loop
+    for i in range(fistaRecon.getParameter('number_of_iterations')):
+##        % With OS approach it becomes trickier to correlate independent subsets, hence additional work is required
+##        % one solution is to work with a full sinogram at times
+##        if ((i >= 3) && (lambdaR_L1 > 0))
+##            [sino_id2, sino_updt2] = astra_create_sino3d_cuda(X, proj_geom, vol_geom);
+##            astra_mex_data3d('delete', sino_id2);
+##        end
+        # With OS approach it becomes trickier to correlate independent subsets,
+        # hence additional work is required one solution is to work with a full
+        # sinogram at times
+
+        ## https://github.com/vais-ral/CCPi-FISTA_Reconstruction/issues/4
+        if (lambdaR_L1 > 0) :
+            sino_id2, sino_updt2 = astra.creators.create_sino3d_gpu(
+                X, proj_geom, vol_geom)
+            astra.matlab.data3d('delete', sino_id2)
+
+        # subset loop
+        counterInd = 1
+        for ss in range(fistaRecon.getParameter('subsets')):
+            print ("Subset {0}".format(ss))
+            X_old = X.copy()
+            t_old = t
+            r_old = fistaRecon.r.copy()
+
+            # the number of projections per subset
+            numProjSub = fistaRecon.getParameter('os_bins')[ss]
+            CurrSubIndices = fistaRecon.getParameter('os_indices')\
+                             [counterInd:counterInd+numProjSub-1]
+            proj_geomSUB['ProjectionAngles'] = angles[CurrSubIndeces]
+            
+##        if fistaRecon.getParameter('projector_geometry')['type'] == 'parallel' or \
+##           fistaRecon.getParameter('projector_geometry')['type'] == 'fanflat' or \
+##           fistaRecon.getParameter('projector_geometry')['type'] == 'fanflat_vec' :
+##            # if the geometry is parallel use slice-by-slice
+##            # projection-backprojection routine
+##            #sino_updt = zeros(size(sino),'single');
+##            proj_geomT = proj_geom.copy()
+##            proj_geomT['DetectorRowCount'] = 1
+##            vol_geomT = vol_geom.copy()
+##            vol_geomT['GridSliceCount'] = 1;
+##            sino_updt = numpy.zeros(numpy.shape(sino), dtype=numpy.float)
+##            for kkk in range(SlicesZ):
+##                sino_id, sino_updt[kkk] = \
+##                         astra.creators.create_sino3d_gpu(
+##                             X_t[kkk:kkk+1], proj_geom, vol_geom)
+##                astra.matlab.data3d('delete', sino_id)
+##        else:
+##            # for divergent 3D geometry (watch the GPU memory overflow in
+##            # ASTRA versions < 1.8)
+##            #[sino_id, sino_updt] = astra_create_sino3d_cuda(X_t, proj_geom, vol_geom);
+##            sino_id, sino_updt = astra.creators.create_sino3d_gpu(
+##                X_t, proj_geom, vol_geom)
+
+        ## RING REMOVAL
+        residual = fistaRecon.residual
+        residual2 = fistaRecon.residual2
+        
+        lambdaR_L1 , alpha_ring , weights , L_const= \
+                   fistaRecon.getParameter(['ring_lambda_R_L1',
+                                           'ring_alpha' , 'weights',
+                                           'Lipschitz_constant'])
+        r_x = fistaRecon.r_x
+        SlicesZ, anglesNumb, Detectors = \
+                    numpy.shape(fistaRecon.getParameter('input_sinogram'))
+        if lambdaR_L1 > 0 :
+             print ("ring removal")
+##             % the ring removal part (Group-Huber fidelity)
+##                % first 2 iterations do additional work reconstructing whole dataset to ensure
+##                % the stablility
+##                if (i < 3)
+##                    [sino_id2, sino_updt2] = astra_create_sino3d_cuda(X_t, proj_geom, vol_geom);
+##                    astra_mex_data3d('delete', sino_id2);
+##                else
+##                    [sino_id, sino_updt] = astra_create_sino3d_cuda(X_t, proj_geomSUB, vol_geom);
+##                end
+
+## https://github.com/vais-ral/CCPi-FISTA_Reconstruction/issues/4
+             if i < 3:
+                 pass
+             else:
+                 sino_id, sino_updt = astra.creators.create_sino3d_gpu(
+                      X_t, proj_geomSUB, vol_geom)
+##                 sino_id, sino_updt = astra.creators.create_sino3d_gpu(
+##                    X, proj_geom, vol_geom)
+##                 astra.matlab.data3d('delete', sino_id)
+                
+             for kkk in range(anglesNumb):
+                 
+                 residual2[:,kkk,:] = (weights[:,kkk,:]).squeeze() * \
+                                       ((sino_updt2[:,kkk,:]).squeeze() - \
+                                        (sino[:,kkk,:]).squeeze() -\
+                                        (alpha_ring * r_x)
+                                        )
+             shape = list(numpy.shape(fistaRecon.getParameter('input_sinogram')))
+             shape[1] = numProjSub
+             fistaRecon.residual = numpy.zeros(shape)
+             if fistaRecon.residual.__hash__() != residual.__hash__():
+                 residual = fistaRecon.residual
+##             for kkk = 1:numProjSub
+##                    indC = CurrSubIndeces(kkk);
+##                    if (i < 3)
+##                        residual(:,kkk,:) =  squeeze(residual2(:,indC,:));
+##                    else
+##                        residual(:,kkk,:) =  squeeze(weights(:,indC,:)).*(squeeze(sino_updt(:,kkk,:)) - (squeeze(sino(:,indC,:)) - alpha_ring.*r_x));
+##                    end
+##                end    
+             for kk in range(numProjSub):
+                 indC = fistaRecon.getParameter('os_indices')[kkk]
+                 if i < 3:
+                     residual[:,kkk,:] = residual2[:,indC,:].squeeze()
+                 else:
+                     residual(:,kkk,:) = \
+                        weights[:,indC,:].squeeze() * sino_updt[:,kkk,:].squeeze() - \
+                        sino[:,indC,:].squeeze() - alpha_ring * fistaRecon.r_x
+                        #squeeze(weights(:,indC,:)).* \
+                        #            (squeeze(sino_updt(:,kkk,:)) - \
+                        #(squeeze(sino(:,indC,:)) - alpha_ring.*r_x));
+                     
+
+                
+             vec = residual.sum(axis = 1)
+             #if SlicesZ > 1:
+             #    vec = vec[:,1,:].squeeze()
+             fistaRecon.r = (r_x - (1./L_const) * vec).copy()
+             objective[i] = (0.5 * (residual ** 2).sum())
+##            % the ring removal part (Group-Huber fidelity)
+##            for kkk = 1:anglesNumb
+##                residual(:,kkk,:) =  squeeze(weights(:,kkk,:)).*
+##                 (squeeze(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;
+##            objective(i) = (0.5*sum(residual(:).^2)); % for the objective function output
+
+        
+
+        # Projection/Backprojection Routine
+        if fistaRecon.getParameter('projector_geometry')['type'] == 'parallel' or \
+           fistaRecon.getParameter('projector_geometry')['type'] == 'fanflat' or\
+           fistaRecon.getParameter('projector_geometry')['type'] == 'fanflat_vec':
+            x_temp = numpy.zeros(numpy.shape(X),dtype=numpy.float32)
+            print ("Projection/Backprojection Routine")
+            for kkk in range(SlicesZ):
+                
+                x_id, x_temp[kkk] = \
+                         astra.creators.create_backprojection3d_gpu(
+                             residual[kkk:kkk+1],
+                             proj_geomT, vol_geomT)
+                astra.matlab.data3d('delete', x_id)
+        else:
+            x_id, x_temp = \
+                  astra.creators.create_backprojection3d_gpu(
+                      residual, proj_geom, vol_geom)            
+
+        X = X_t - (1/L_const) * x_temp
+        astra.matlab.data3d('delete', sino_id)
+        astra.matlab.data3d('delete', x_id)
+        
+
+        ## REGULARIZATION
+        ## SKIPPING FOR NOW
+        ## Should be simpli
+        # regularizer = fistaRecon.getParameter('regularizer')
+        # for slices:
+        # out = regularizer(input=X)
+        print ("skipping regularizer")
+
+
+        ## FINAL
+        print ("final")
+        lambdaR_L1 = fistaRecon.getParameter('ring_lambda_R_L1')
+        if lambdaR_L1 > 0:
+            fistaRecon.r = numpy.max(
+                numpy.abs(fistaRecon.r) - lambdaR_L1 , 0) * \
+                numpy.sign(fistaRecon.r)
+        t = (1 + numpy.sqrt(1 + 4 * t**2))/2
+        X_t = X + (((t_old -1)/t) * (X - X_old))
+
+        if lambdaR_L1 > 0:
+            fistaRecon.r_x = fistaRecon.r + \
+                             (((t_old-1)/t) * (fistaRecon.r - r_old))
+
+        if fistaRecon.getParameter('region_of_interest') is None:
+            string = 'Iteration Number {0} | Objective {1} \n'
+            print (string.format( i, objective[i]))
+        else:
+            ROI , X_ideal = fistaRecon.getParameter('region_of_interest',
+                                                    'ideal_image')
+            
+            Resid_error[i] = RMSE(X*ROI, X_ideal*ROI)
+            string = 'Iteration Number {0} | RMS Error {1} | Objective {2} \n'
+            print (string.format(i,Resid_error[i], objective[i]))
+            
+##        if (lambdaR_L1 > 0)
+##            r =  max(abs(r)-lambdaR_L1, 0).*sign(r); % soft-thresholding operator for ring vector
+##        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 > 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 %f \n', 'Iteration Number:', i, '|', 'Error RMSE:', Resid_error(i), '|', 'Objective:', objective(i));
+##        else
+##            fprintf('%s %i  %s %s %f \n', 'Iteration Number:', i, '|', 'Objective:', objective(i));
+##        end
+else:
+    fistaRecon = FISTAReconstructor(proj_geom,
+                                vol_geom,
+                                Sino3D ,
+                                weights=Weights3D)
+
+    print ("Lipschitz Constant {0}".format(fistaRecon.pars['Lipschitz_constant']))
+    fistaRecon.setParameter(number_of_iterations = 12)
+    fistaRecon.setParameter(Lipschitz_constant = 767893952.0)
+    fistaRecon.setParameter(ring_alpha = 21)
+    fistaRecon.setParameter(ring_lambda_R_L1 = 0.002)
+    fistaRecon.prepareForIteration()
+    X = fistaRecon.iterate(numpy.load("X.npy"))
diff --git a/src/Python/test_reconstructor.py b/src/Python/test_reconstructor.py
index f8f6b3c..07668ba 100644
--- a/src/Python/test_reconstructor.py
+++ b/src/Python/test_reconstructor.py
@@ -11,10 +11,17 @@ import numpy
 
 from ccpi.fista.FISTAReconstructor import FISTAReconstructor
 import astra
+import matplotlib.pyplot as plt
 
-##def getEntry(nx, location):
-##    for item in nx[location].keys():
-##        print (item)
+def RMSE(signal1, signal2):
+    '''RMSE Root Mean Squared Error'''
+    if numpy.shape(signal1) == numpy.shape(signal2):
+        err = (signal1 - signal2)
+        err = numpy.sum( err * err )/numpy.size(signal1);  # MSE
+        err = sqrt(err);                                   # RMSE
+        return err
+    else:
+        raise Exception('Input signals must have the same shape')
   
 filename = r'/home/ofn77899/Reconstruction/CCPi-FISTA_Reconstruction/demos/DendrData.h5'
 nx = h5py.File(filename, "r")
@@ -68,7 +75,6 @@ fistaRecon.setParameter(number_of_iterations = 12)
 fistaRecon.setParameter(Lipschitz_constant = 767893952.0)
 fistaRecon.setParameter(ring_alpha = 21)
 fistaRecon.setParameter(ring_lambda_R_L1 = 0.002)
-#fistaRecon.setParameter(use_studentt_fidelity= True)
 
 ## Ordered subset
 if False:
@@ -94,19 +100,34 @@ if False:
             counter = counter + binsDiscr[jj] - 1
 
 
-if True:
-    fistaRecon.prepareForIteration()
+if False:
     print ("Lipschitz Constant {0}".format(fistaRecon.pars['Lipschitz_constant']))
-
+    print ("prepare for iteration")
+    fistaRecon.prepareForIteration()
+    
     
 
+    print("initializing  ...")
+    if False:
+        # if X doesn't exist
+        #N = params.vol_geom.GridColCount
+        N = vol_geom['GridColCount']
+        print ("N " + str(N))
+        X = numpy.zeros((N,N,SlicesZ), dtype=numpy.float)
+    else:
+        #X = fistaRecon.initialize()
+        X = numpy.load("X.npy")
+    
+    print (numpy.shape(X))
+    X_t = X.copy()
+    print ("initialized")
     proj_geom , vol_geom, sino , \
                       SlicesZ = fistaRecon.getParameter(['projector_geometry' ,
                                                             'output_geometry',
                                                             'input_sinogram',
                                                          'SlicesZ'])
 
-    fistaRecon.setParameter(number_of_iterations = 3)
+    #fistaRecon.setParameter(number_of_iterations = 3)
     iterFISTA = fistaRecon.getParameter('number_of_iterations')
     # errors vector (if the ground truth is given)
     Resid_error = numpy.zeros((iterFISTA));
@@ -114,30 +135,18 @@ if True:
     objective = numpy.zeros((iterFISTA)); 
 
       
-    print ("line")                        
     t = 1
-    print ("line")
 
-    if False:
-        # if X doesn't exist
-        #N = params.vol_geom.GridColCount
-        N = vol_geom['GridColCount']
-        print ("N " + str(N))
-        X = numpy.zeros((N,N,SlicesZ), dtype=numpy.float)
-    else:
-        #X = fistaRecon.initialize()
-        X = numpy.load("X.npy")
-    
-    print (numpy.shape(X))
-    X_t = X.copy()
-    print ("X_t copy")
+
+    print ("starting iterations")
 ##    % Outer FISTA iterations loop
     for i in range(fistaRecon.getParameter('number_of_iterations')):
         X_old = X.copy()
         t_old = t
         r_old = fistaRecon.r.copy()
         if fistaRecon.getParameter('projector_geometry')['type'] == 'parallel' or \
-           fistaRecon.getParameter('projector_geometry')['type'] == 'parallel3d':
+           fistaRecon.getParameter('projector_geometry')['type'] == 'fanflat' or \
+           fistaRecon.getParameter('projector_geometry')['type'] == 'fanflat_vec' :
             # if the geometry is parallel use slice-by-slice
             # projection-backprojection routine
             #sino_updt = zeros(size(sino),'single');
@@ -147,16 +156,15 @@ if True:
             vol_geomT['GridSliceCount'] = 1;
             sino_updt = numpy.zeros(numpy.shape(sino), dtype=numpy.float)
             for kkk in range(SlicesZ):
-                print (kkk)
                 sino_id, sino_updt[kkk] = \
                          astra.creators.create_sino3d_gpu(
-                             X_t[kkk:kkk+1], proj_geomT, vol_geomT)
+                             X_t[kkk:kkk+1], proj_geom, vol_geom)
                 astra.matlab.data3d('delete', sino_id)
         else:
             # for divergent 3D geometry (watch the GPU memory overflow in
             # ASTRA versions < 1.8)
             #[sino_id, sino_updt] = astra_create_sino3d_cuda(X_t, proj_geom, vol_geom);
-            sino_id, sino_updt = astra.matlab.create_sino3d_gpu(
+            sino_id, sino_updt = astra.creators.create_sino3d_gpu(
                 X_t, proj_geom, vol_geom)
 
         ## RING REMOVAL
@@ -169,8 +177,9 @@ if True:
         SlicesZ, anglesNumb, Detectors = \
                     numpy.shape(fistaRecon.getParameter('input_sinogram'))
         if lambdaR_L1 > 0 :
+             print ("ring removal")
              for kkk in range(anglesNumb):
-                 print ("angles {0}".format(kkk))
+                 
                  residual[:,kkk,:] = (weights[:,kkk,:]).squeeze() * \
                                        ((sino_updt[:,kkk,:]).squeeze() - \
                                         (sino[:,kkk,:]).squeeze() -\
@@ -194,39 +203,16 @@ if True:
 ##            r = r_x - (1./L_const).*vec;
 ##            objective(i) = (0.5*sum(residual(:).^2)); % for the objective function output
 
-        else:
-            if fistaRecon.getParameter('use_studentt_fidelity'):
-                residual = weights * (sino_updt - sino)
-                for kkk in range(SlicesZ):
-                    # reshape(residual(:,:,kkk), Detectors*anglesNumb, 1)
-                    # 1D
-                    res_vec = numpy.reshape(residual[kkk], (Detectors * anglesNumb,1))
-                
-##            else            
-##			if (studentt == 1)
-##				% artifacts removal with Students t penalty
-##				residual = weights.*(sino_updt - sino);
-##				for kkk = 1:SlicesZ
-##				   res_vec = reshape(residual(:,:,kkk), Detectors*anglesNumb, 1); % 1D vectorized sinogram
-##                                 %s = 100;
-##                                 %gr = (2)*res_vec./(s*2 + conj(res_vec).*res_vec);
-##                                 [ff, gr] = studentst(res_vec, 1);
-##                                 residual(:,:,kkk) = reshape(gr, Detectors, anglesNumb);
-##				end            
-##				objective(i) = ff; % for the objective function output
-##            else
-##            % no ring removal (LS model)
-##            residual = weights.*(sino_updt - sino);
-##            objective(i) = (0.5*sum(residual(:).^2)); % for the objective function output
-##            end	
-##        end       
+        
 
         # Projection/Backprojection Routine
         if fistaRecon.getParameter('projector_geometry')['type'] == 'parallel' or \
-           fistaRecon.getParameter('projector_geometry')['type'] == 'parallel3d':
+           fistaRecon.getParameter('projector_geometry')['type'] == 'fanflat' or\
+           fistaRecon.getParameter('projector_geometry')['type'] == 'fanflat_vec':
             x_temp = numpy.zeros(numpy.shape(X),dtype=numpy.float32)
+            print ("Projection/Backprojection Routine")
             for kkk in range(SlicesZ):
-                print ("Projection/Backprojection Routine {0}".format( kkk ))
+                
                 x_id, x_temp[kkk] = \
                          astra.creators.create_backprojection3d_gpu(
                              residual[kkk:kkk+1],
@@ -248,9 +234,11 @@ if True:
         # regularizer = fistaRecon.getParameter('regularizer')
         # for slices:
         # out = regularizer(input=X)
+        print ("skipping regularizer")
 
 
         ## FINAL
+        print ("final")
         lambdaR_L1 = fistaRecon.getParameter('ring_lambda_R_L1')
         if lambdaR_L1 > 0:
             fistaRecon.r = numpy.max(
@@ -263,9 +251,16 @@ if True:
             fistaRecon.r_x = fistaRecon.r + \
                              (((t_old-1)/t) * (fistaRecon.r - r_old))
 
-        if fistaRecon.getParameter('ideal_image') is None:
+        if fistaRecon.getParameter('region_of_interest') is None:
             string = 'Iteration Number {0} | Objective {1} \n'
             print (string.format( i, objective[i]))
+        else:
+            ROI , X_ideal = fistaRecon.getParameter('region_of_interest',
+                                                    'ideal_image')
+            
+            Resid_error[i] = RMSE(X*ROI, X_ideal*ROI)
+            string = 'Iteration Number {0} | RMS Error {1} | Objective {2} \n'
+            print (string.format(i,Resid_error[i], objective[i]))
             
 ##        if (lambdaR_L1 > 0)
 ##            r =  max(abs(r)-lambdaR_L1, 0).*sign(r); % soft-thresholding operator for ring vector
@@ -291,3 +286,16 @@ if True:
 ##        else
 ##            fprintf('%s %i  %s %s %f \n', 'Iteration Number:', i, '|', 'Objective:', objective(i));
 ##        end
+else:
+    fistaRecon = FISTAReconstructor(proj_geom,
+                                vol_geom,
+                                Sino3D ,
+                                weights=Weights3D)
+
+    print ("Lipschitz Constant {0}".format(fistaRecon.pars['Lipschitz_constant']))
+    fistaRecon.setParameter(number_of_iterations = 12)
+    fistaRecon.setParameter(Lipschitz_constant = 767893952.0)
+    fistaRecon.setParameter(ring_alpha = 21)
+    fistaRecon.setParameter(ring_lambda_R_L1 = 0.002)
+    fistaRecon.prepareForIteration()
+    X = fistaRecon.iterate(numpy.load("X.npy"))