further development using abstract forward/backward projector

4 files changed, 84 insertions, 59 deletions

diff --git a/src/Python/ccpi/reconstruction/AstraDevice.py b/src/Python/ccpi/reconstruction/AstraDevice.py
index d854183..ba174b3 100644
--- a/src/Python/ccpi/reconstruction/AstraDevice.py
+++ b/src/Python/ccpi/reconstruction/AstraDevice.py
@@ -1,5 +1,5 @@
 import astra
-from DeviceModel import DeviceModel
+from ccpi.reconstruction.DeviceModel import DeviceModel
 
 class AstraDevice(DeviceModel):
     '''Concrete class for Astra Device'''
@@ -9,6 +9,9 @@ class AstraDevice(DeviceModel):
                  data_aquisition_geometry,
                  reconstructed_volume_geometry):
         
+        super(AstraDevice, self).__init__(device_type,
+                                          data_aquisition_geometry,
+                                          reconstructed_volume_geometry)
 
         self.proj_geom = astra.creators.create_proj_geom(
             device_type,
@@ -17,7 +20,6 @@ class AstraDevice(DeviceModel):
             self.acquisition_data_geometry['cameraX'],
             self.acquisition_data_geometry['cameraY'],
             self.acquisition_data_geometry['angles'],
-            angles_rad
             )
 
         self.vol_geom = astra.creators.create_vol_geom(
@@ -48,19 +50,27 @@ Uses Astra-toolbox
                    self.vol_geom)
         astra.matlab.data3d('delete', idx)
         return volume
-
     def createReducedDevice(self):
-        return AstraDevice(self.proj_geom['type'],
-                    {'detectorSpacingX' : self.proj_geom['DetectorSpacingX'] ,
-                     'detectorSpacingY' : self.proj_geom['DetectorSpacingY'] ,
-                     'cameraX' : self.proj_geom['DetectorColCount'] ,
-                     'cameraY' : 1 ,
-                     'angles' : self.proj_geom['ProjectionAngles'] } ,
-                    {
-                        'X' : self.vol_geom['GridColCount'],
-                        'Y' : self.vol_geom['GridRowCount']
-                        'Z' : 1} ) 
+        proj_geom = astra.creators.create_proj_geom(
+            device_type,
+            self.acquisition_data_geometry['detectorSpacingX'],
+            self.acquisition_data_geometry['detectorSpacingX'],
+            self.acquisition_data_geometry['cameraX'],
+            1,
+            self.acquisition_data_geometry['angles'],
+            )
 
+        vol_geom = astra.creators.create_vol_geom(
+            self.reconstructed_volume_geometry['X'],
+            self.reconstructed_volume_geometry['Y'],
+            1
+            )
+        return AstraDevice(proj_geom['type'] ,
+                           proj_geom,
+                           vol_geom)
+        
+
+        
 if __name__=="main":
     a = AstraDevice()
 
diff --git a/src/Python/ccpi/reconstruction/DeviceModel.py b/src/Python/ccpi/reconstruction/DeviceModel.py
index 1717321..6214437 100644
--- a/src/Python/ccpi/reconstruction/DeviceModel.py
+++ b/src/Python/ccpi/reconstruction/DeviceModel.py
@@ -44,7 +44,7 @@ reconstructed_volume_geometry: tuple (dimX,dimY,dimZ)
             'Z': reconstructed_volume_geometry[2] }
 
     @abstractmethod
-    def doFowardProject(self, volume):
+    def doForwardProject(self, volume):
         '''Forward projects the volume according to the device geometry'''
         return NotImplemented
 
@@ -53,5 +53,10 @@ reconstructed_volume_geometry: tuple (dimX,dimY,dimZ)
     def doBackwardProject(self, projections):
         '''Backward projects the projections according to the device geometry'''
         return NotImplemented
+
+    @abstractmethod
+    def createReducedDevice(self):
+        '''Create a Device to do forward/backward projections on 2D slices'''
+        return NotImplemented
     
 
diff --git a/src/Python/ccpi/reconstruction/FISTAReconstructor.py b/src/Python/ccpi/reconstruction/FISTAReconstructor.py
index 1e464a1..d2eefc0 100644
--- a/src/Python/ccpi/reconstruction/FISTAReconstructor.py
+++ b/src/Python/ccpi/reconstruction/FISTAReconstructor.py
@@ -27,7 +27,7 @@ import numpy
 from enum import Enum
 
 import astra
-from ccpi.reconstruction.AstraDevice import AstraDevice
+#from ccpi.reconstruction.AstraDevice import AstraDevice
 
    
     
@@ -74,7 +74,10 @@ class FISTAReconstructor():
     # 3. "A novel tomographic reconstruction method based on the robust
     # Student's t function for suppressing data outliers" D. Kazantsev et.al.
     # D. Kazantsev, 2016-17
-    def __init__(self, projector_geometry, output_geometry, input_sinogram,
+    def __init__(self, projector_geometry,
+                 output_geometry,
+                 input_sinogram,
+                 device, 
                  **kwargs):
         # handle parmeters:
         # obligatory parameters
@@ -87,16 +90,10 @@ class FISTAReconstructor():
         self.pars['number_of_angles'] = nangles
         self.pars['SlicesZ'] = sliceZ
         self.pars['output_volume'] = None
-
+        self.pars['device'] = device
 
         
-        device = createAstraDevice(projector_geometry, output_geometry)
-        self.setParameter(device_model=device)
-        proj_geomT = projector_geometry.copy();
-        proj_geomT['DetectorRowCount'] = 1;
-        vol_geomT = output_geometry.copy();
-        vol_geomT['GridSliceCount'] = 1;
-        reduced_device = createAstraDevice(proj_geomT, vol_geomT)
+        reduced_device = device.createReducedDevice()
         self.setParameter(reduced_device_model=reduced_device)
 
         self.use_device = True
@@ -187,21 +184,7 @@ class FISTAReconstructor():
             self.pars['initialize'] = False
 
         
-    def createAstraDevice(self, projector_geometry, output_geometry):
-        '''TODO remove'''
-        
-        device = AstraDevice(DeviceModel.PARALLEL3D,
-                    {'detectorSpacingX' : projector_geometry['DetectorSpacingX'] ,
-                     'detectorSpacingY' : projector_geometry['DetectorSpacingY'] ,
-                     'cameraX' : projector_geometry['DetectorColCount'] ,
-                     'cameraY' : projector_geometry['DetectorRowCount'] ,
-                     'angles' : projector_geometry['ProjectionAngles'] } ,
-                    {
-                        'X' : output_geometry['GridColCount'],
-                        'Y' : output_geometry['GridRowCount']
-                        'Z' : output_geometry['GridSliceCount']} )
-        return device
-            
+                
     def setParameter(self, **kwargs):
         '''set named parameter for the reconstructor engine
         
diff --git a/src/Python/test/test_reconstructor.py b/src/Python/test/test_reconstructor.py
index 3342301..343b9bb 100644
--- a/src/Python/test/test_reconstructor.py
+++ b/src/Python/test/test_reconstructor.py
@@ -12,6 +12,9 @@ import numpy
 from ccpi.reconstruction.FISTAReconstructor import FISTAReconstructor
 import astra
 import matplotlib.pyplot as plt
+from ccpi.imaging.Regularizer import Regularizer
+from ccpi.reconstruction.AstraDevice import AstraDevice
+from ccpi.reconstruction.DeviceModel import DeviceModel
 
 def RMSE(signal1, signal2):
     '''RMSE Root Mean Squared Error'''
@@ -22,7 +25,23 @@ def RMSE(signal1, signal2):
         return err
     else:
         raise Exception('Input signals must have the same shape')
-  
+
+def createAstraDevice(projector_geometry, output_geometry):
+    '''TODO remove'''
+    
+    device = AstraDevice(DeviceModel.DeviceType.PARALLEL3D.value,
+                [projector_geometry['DetectorSpacingX'] ,
+                 projector_geometry['DetectorSpacingY'] ,
+                 projector_geometry['DetectorColCount'] ,
+                 projector_geometry['DetectorRowCount'] ,
+                 projector_geometry['ProjectionAngles']
+                 ],
+                [
+                    output_geometry['GridColCount'],
+                    output_geometry['GridRowCount'], 
+                    output_geometry['GridSliceCount'] ] )
+    return device
+
 filename = r'/home/ofn77899/Reconstruction/CCPi-FISTA_Reconstruction/demos/DendrData.h5'
 nx = h5py.File(filename, "r")
 #getEntry(nx, '/')
@@ -65,23 +84,23 @@ vol_geom = astra.creators.create_vol_geom( image_size_x,
 ## 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)
-
-reg = Regularizer(Regularizer.Algorithm.LLT_model)
-reg.setParameter(regularization_parameter=25,
-                          time_step=0.0003,
-                          tolerance_constant=0.0001,
-                          number_of_iterations=300)
-fistaRecon.setParameter(regularizer = reg)
+##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)
+##
+##reg = Regularizer(Regularizer.Algorithm.LLT_model)
+##reg.setParameter(regularization_parameter=25,
+##                          time_step=0.0003,
+##                          tolerance_constant=0.0001,
+##                          number_of_iterations=300)
+##fistaRecon.setParameter(regularizer=reg)
 
 ## Ordered subset
 if False:
@@ -294,16 +313,24 @@ if False:
 ##            fprintf('%s %i  %s %s %f \n', 'Iteration Number:', i, '|', 'Objective:', objective(i));
 ##        end
 else:
+    
+    # create a device for forward/backprojection
+    astradevice = createAstraDevice(proj_geom, vol_geom)
     fistaRecon = FISTAReconstructor(proj_geom,
                                 vol_geom,
                                 Sino3D ,
-                                weights=Weights3D)
+                                device = astradevice,
+                                weights=Weights3D
+                                )
 
     print ("Lipschitz Constant {0}".format(fistaRecon.pars['Lipschitz_constant']))
-    fistaRecon.setParameter(number_of_iterations = 12)
+    fistaRecon.setParameter(number_of_iterations = 3)
     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"))
     numpy.save("X_out.npy", X)