implemented non ordered subset FISTA in reconstructor

2 files changed, 165 insertions, 25 deletions

diff --git a/src/Python/ccpi/fista/FISTAReconstructor.py b/src/Python/ccpi/fista/FISTAReconstructor.py
index 33e67a3..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)
@@ -108,7 +109,11 @@ class FISTAReconstructor():
               'regularizer' , 
               'ring_lambda_R_L1',
               'ring_alpha',
-              'subsets')
+              'subsets',
+              'output_volume',
+              'os_subsets',
+              'os_indices',
+              'os_bins')
         self.acceptedInputKeywords = list(kw)
         
         # handle keyworded parameters
@@ -176,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) +
@@ -382,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']))
 
@@ -401,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()
@@ -423,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');
@@ -439,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)
@@ -450,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.py b/src/Python/test_reconstructor.py
index 2f188b4..07668ba 100644
--- a/src/Python/test_reconstructor.py
+++ b/src/Python/test_reconstructor.py
@@ -100,7 +100,7 @@ if False:
             counter = counter + binsDiscr[jj] - 1
 
 
-if True:
+if False:
     print ("Lipschitz Constant {0}".format(fistaRecon.pars['Lipschitz_constant']))
     print ("prepare for iteration")
     fistaRecon.prepareForIteration()
@@ -145,7 +145,8 @@ if True:
         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');
@@ -157,13 +158,13 @@ if True:
             for kkk in range(SlicesZ):
                 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
@@ -206,7 +207,8 @@ if True:
 
         # 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):
@@ -284,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"))