uses AstraDevice

1 files changed, 92 insertions, 32 deletions

diff --git a/src/Python/ccpi/reconstruction/FISTAReconstructor.py b/src/Python/ccpi/reconstruction/FISTAReconstructor.py
index c903712..1e464a1 100644
--- a/src/Python/ccpi/reconstruction/FISTAReconstructor.py
+++ b/src/Python/ccpi/reconstruction/FISTAReconstructor.py
@@ -27,6 +27,7 @@ import numpy
 from enum import Enum
 
 import astra
+from ccpi.reconstruction.AstraDevice import AstraDevice
 
    
     
@@ -87,6 +88,19 @@ class FISTAReconstructor():
         self.pars['SlicesZ'] = sliceZ
         self.pars['output_volume'] = None
 
+
+        
+        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)
+        self.setParameter(reduced_device_model=reduced_device)
+
+        self.use_device = True
+        
         print (self.pars)
         # handle optional input parameters (at instantiation)
         
@@ -113,7 +127,9 @@ class FISTAReconstructor():
               'output_volume',
               'os_subsets',
               'os_indices',
-              'os_bins')
+              'os_bins',
+              'device_model',
+              'reduced_device_model')
         self.acceptedInputKeywords = list(kw)
         
         # handle keyworded parameters
@@ -171,7 +187,20 @@ 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
@@ -436,13 +465,17 @@ class FISTAReconstructor():
         X_t = X.copy()
         # convenience variable storage
         proj_geom , vol_geom, sino , \
-          SlicesZ  = self.getParameter([ 'projector_geometry' ,
+          SlicesZ , ring_lambda_R_L1 = self.getParameter([ 'projector_geometry' ,
                                                 'output_geometry',
                                                 'input_sinogram',
-                                                'SlicesZ' ])
+                                                'SlicesZ' ,
+                                                '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')):
             X_old = X.copy()
             t_old = t
@@ -453,29 +486,42 @@ class FISTAReconstructor():
                 # 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;
-                self.sino_updt = numpy.zeros(numpy.shape(sino), dtype=numpy.float)
-                for kkk in range(SlicesZ):
-                    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)
+                
+                if self.use_device:
+                    self.sino_updt = numpy.zeros(numpy.shape(sino), dtype=numpy.float)
+                    
+                    for kkk in range(SlicesZ):
+                        self.sino_updt[kkk] = \
+                            reduced_device.doForwardProject( X_t[kkk:kkk+1] )
+                else:
+                    proj_geomT = proj_geom.copy()
+                    proj_geomT['DetectorRowCount'] = 1
+                    vol_geomT = vol_geom.copy()
+                    vol_geomT['GridSliceCount'] = 1;
+                    self.sino_updt = numpy.zeros(numpy.shape(sino), dtype=numpy.float)
+                    for kkk in range(SlicesZ):
+                        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)
             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, self.sino_updt = astra.creators.create_sino3d_gpu(
-                    X_t, proj_geom, vol_geom)
+                if self.use_device:
+                    self.sino_updt = device.doForwardProject(X_t)
+                else:
+                    sino_id, self.sino_updt = astra.creators.create_sino3d_gpu(
+                        X_t, proj_geom, vol_geom)
+                    astra.matlab.data3d('delete', sino_id)
 
 
             ## RING REMOVAL
-            self.ringRemoval(i)
+            if ring_lambda_R_L1 != 0:
+                self.ringRemoval(i)
             ## Projection/Backprojection Routine
             self.projectionBackprojection(X, X_t)
-            astra.matlab.data3d('delete', sino_id)
+            
             ## REGULARIZATION
             X = self.regularize(X)
             ## Update Loop
@@ -523,6 +569,8 @@ class FISTAReconstructor():
                                                   'output_geometry',
                                                   'Lipschitz_constant'])
         
+        device, reduced_device = self.getParameter(['device_model',
+                                                    'reduced_device_model'])
         
         if self.getParameter('projector_geometry')['type'] == 'parallel' or \
            self.getParameter('projector_geometry')['type'] == 'fanflat' or \
@@ -530,27 +578,39 @@ class FISTAReconstructor():
             # 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)
+            if use_device:
+                proj_geomT = proj_geom.copy()
+                proj_geomT['DetectorRowCount'] = 1
+                vol_geomT = vol_geom.copy()
+                vol_geomT['GridSliceCount'] = 1;
+                
+                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:
+                for kkk in range(SliceZ):
+                    x_temp[kkk] = \
+                        reduced_device.doBackwardProject(residual[kkk:kkk+1])
         else:
-            x_id, x_temp = \
+            if use_device:
+                x_id, x_temp = \
                   astra.creators.create_backprojection3d_gpu(
-                      residual, proj_geom, vol_geom)            
+                      residual, proj_geom, vol_geom)
+                astra.matlab.data3d('delete', x_id)
+            else:
+                x_temp = \
+                    device.doBackwardProject(residual)
+                       
 
         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")