bugfixes for AstraDevice use

4 files changed, 71 insertions, 35 deletions

diff --git a/src/Python/ccpi/reconstruction/AstraDevice.py b/src/Python/ccpi/reconstruction/AstraDevice.py
index 999533b..ac9206e 100644
--- a/src/Python/ccpi/reconstruction/AstraDevice.py
+++ b/src/Python/ccpi/reconstruction/AstraDevice.py
@@ -1,5 +1,6 @@
 import astra
 from ccpi.reconstruction.DeviceModel import DeviceModel
+import numpy
 
 class AstraDevice(DeviceModel):
     '''Concrete class for Astra Device'''
@@ -17,7 +18,7 @@ class AstraDevice(DeviceModel):
         self.proj_geom = astra.creators.create_proj_geom(
             device_type,
             self.acquisition_data_geometry['detectorSpacingX'],
-            self.acquisition_data_geometry['detectorSpacingX'],
+            self.acquisition_data_geometry['detectorSpacingY'],
             self.acquisition_data_geometry['cameraX'],
             self.acquisition_data_geometry['cameraY'],
             self.acquisition_data_geometry['angles'],
@@ -34,10 +35,15 @@ class AstraDevice(DeviceModel):
 
 Uses Astra-toolbox
 '''
-        sino_id, y = astra.creators.create_sino3d_gpu(
-            volume, self.proj_geom, self.vol_geom)
-        astra.matlab.data3d('delete', sino_id)
-        return y
+              
+        try:
+            sino_id, y = astra.creators.create_sino3d_gpu(
+                volume, self.proj_geom, self.vol_geom)
+            astra.matlab.data3d('delete', sino_id)
+            return y
+        except Exception(e):
+            print(e)
+            print("Value Error: ", self.proj_geom, self.vol_geom)
 
     def doBackwardProject(self, projections):
         '''Backward projects the projections according to the device geometry
@@ -49,27 +55,25 @@ Uses Astra-toolbox
                    projections,
                    self.proj_geom,
                    self.vol_geom)
+        
         astra.matlab.data3d('delete', idx)
         return volume
     
     def createReducedDevice(self):
-        proj_geom = astra.creators.create_proj_geom(
-            device_type,
-            self.acquisition_data_geometry['detectorSpacingX'],
-            self.acquisition_data_geometry['detectorSpacingX'],
+        proj_geom =  [ 
             self.acquisition_data_geometry['cameraX'],
             1,
-            self.acquisition_data_geometry['angles'],
-            )
+            self.acquisition_data_geometry['detectorSpacingX'],
+            self.acquisition_data_geometry['detectorSpacingY'],
+            self.acquisition_data_geometry['angles']
+            ]
 
-        vol_geom = astra.creators.create_vol_geom(
+        vol_geom = [
             self.reconstructed_volume_geometry['X'],
             self.reconstructed_volume_geometry['Y'],
             1
-            )
-        return AstraDevice(proj_geom['type'] ,
-                           proj_geom,
-                           vol_geom)
+            ]
+        return AstraDevice(self.type, proj_geom, vol_geom)
         
 
         
diff --git a/src/Python/ccpi/reconstruction/DeviceModel.py b/src/Python/ccpi/reconstruction/DeviceModel.py
index 6214437..eeb9a34 100644
--- a/src/Python/ccpi/reconstruction/DeviceModel.py
+++ b/src/Python/ccpi/reconstruction/DeviceModel.py
@@ -28,7 +28,8 @@ HELICAL'''
 
 Mandatory parameters are:
 device_type from DeviceType Enum
-data_acquisition_geometry: tuple (camera_X, camera_Y)
+data_acquisition_geometry: tuple (camera_X, camera_Y, detectorSpacingX,
+                                  detectorSpacingY, angles)
 reconstructed_volume_geometry: tuple (dimX,dimY,dimZ)
 '''
         self.device_geometry = device_type
diff --git a/src/Python/ccpi/reconstruction/FISTAReconstructor.py b/src/Python/ccpi/reconstruction/FISTAReconstructor.py
index d2eefc0..0607003 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
 
    
     
@@ -90,11 +90,7 @@ class FISTAReconstructor():
         self.pars['number_of_angles'] = nangles
         self.pars['SlicesZ'] = sliceZ
         self.pars['output_volume'] = None
-        self.pars['device'] = device
-
-        
-        reduced_device = device.createReducedDevice()
-        self.setParameter(reduced_device_model=reduced_device)
+        self.pars['device_model'] = device
 
         self.use_device = True
         
@@ -183,6 +179,9 @@ class FISTAReconstructor():
         if not 'initialize' in kwargs.keys():
             self.pars['initialize'] = False
 
+        reduced_device = device.createReducedDevice()
+        self.setParameter(reduced_device_model=reduced_device)
+
         
                 
     def setParameter(self, **kwargs):
@@ -455,22 +454,24 @@ class FISTAReconstructor():
                                                 'ring_lambda_R_L1'])
                    
         t = 1
-        
+
         device = self.getParameter('device_model')
         reduced_device = self.getParameter('reduced_device_model')
         
         for i in range(self.getParameter('number_of_iterations')):
+            print("iteration", i)
             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'] == 'fanflat' or \
-               self.getParameter('projector_geometry')['type'] == 'fanflat_vec':
+            pg = self.getParameter('projector_geometry')['type']
+            if pg == 'parallel' or \
+               pg == 'fanflat' or \
+               pg == 'fanflat_vec':
                 # if the geometry is parallel use slice-by-slice
                 # projection-backprojection routine
                 #sino_updt = zeros(size(sino),'single');
                 
-                if self.use_device:
+                if self.use_device :
                     self.sino_updt = numpy.zeros(numpy.shape(sino), dtype=numpy.float)
                     
                     for kkk in range(SlicesZ):
@@ -491,6 +492,7 @@ 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);
+                
                 if self.use_device:
                     self.sino_updt = device.doForwardProject(X_t)
                 else:
@@ -563,7 +565,7 @@ class FISTAReconstructor():
             #sino_updt = zeros(size(sino),'single');
             x_temp = numpy.zeros(numpy.shape(X),dtype=numpy.float32)
                 
-            if use_device:
+            if self.use_device:
                 proj_geomT = proj_geom.copy()
                 proj_geomT['DetectorRowCount'] = 1
                 vol_geomT = vol_geom.copy()
@@ -581,7 +583,7 @@ class FISTAReconstructor():
                     x_temp[kkk] = \
                         reduced_device.doBackwardProject(residual[kkk:kkk+1])
         else:
-            if use_device:
+            if self.use_device:
                 x_id, x_temp = \
                   astra.creators.create_backprojection3d_gpu(
                       residual, proj_geom, vol_geom)
diff --git a/src/Python/test/test_reconstructor.py b/src/Python/test/test_reconstructor.py
index 343b9bb..f4627d7 100644
--- a/src/Python/test/test_reconstructor.py
+++ b/src/Python/test/test_reconstructor.py
@@ -30,10 +30,10 @@ def createAstraDevice(projector_geometry, output_geometry):
     '''TODO remove'''
     
     device = AstraDevice(DeviceModel.DeviceType.PARALLEL3D.value,
-                [projector_geometry['DetectorSpacingX'] ,
-                 projector_geometry['DetectorSpacingY'] ,
+                [projector_geometry['DetectorRowCount'] ,
                  projector_geometry['DetectorColCount'] ,
-                 projector_geometry['DetectorRowCount'] ,
+                 projector_geometry['DetectorSpacingX'] ,
+                 projector_geometry['DetectorSpacingY'] ,
                  projector_geometry['ProjectionAngles']
                  ],
                 [
@@ -332,5 +332,34 @@ else:
     
     
     fistaRecon.prepareForIteration()
-    X = fistaRecon.iterate(numpy.load("X.npy"))
-    numpy.save("X_out.npy", X)
+    X = numpy.load("X.npy")
+##    rd = astradevice.createReducedDevice()
+##    print ("rd proj_geom" , rd.proj_geom)
+##    
+##        
+##    rd.doForwardProject(X[0:1])
+##    proj_geomT = proj_geom.copy()
+##    for ekey in rd.proj_geom.keys():
+##        if ekey == 'ProjectionAngles':
+##            valrd = numpy.shape(rd.proj_geom[ekey])
+##            valg  = numpy.shape(proj_geomT[ekey])
+##        else:
+##            valrd = rd.proj_geom[ekey]
+##            valg =  proj_geomT[ekey]
+##      
+##        print ("key {0}: RD {1} geomT {2}".format(ekey, valrd, valg))
+##        
+##        
+##    proj_geomT['DetectorRowCount'] = 1
+##    vol_geomT = vol_geom.copy()
+##    vol_geomT['GridSliceCount'] = 1;
+##    rd.proj_geom = proj_geomT.copy()
+##    rd.vol_geom = vol_geomT.copy()
+##
+##
+##
+##    sino_id, y = astra.creators.create_sino3d_gpu(
+##                X[0:1], rd.proj_geom, rd.vol_geom)
+    
+    X = fistaRecon.iterate(X)
+    #numpy.save("X_out.npy", X)