Geoms into datasets #15, #16 (#34)

* Added vol and sino geoms into vol and sino data sets and in astra operator

* Moved geometry attribute up to DataSet

* Added copying of geometry in DataSet arithmetic operations

3 files changed, 55 insertions, 38 deletions

diff --git a/Wrappers/Python/ccpi/framework.py b/Wrappers/Python/ccpi/framework.py
index 2b0ba76..5ac17ee 100644
--- a/Wrappers/Python/ccpi/framework.py
+++ b/Wrappers/Python/ccpi/framework.py
@@ -105,7 +105,8 @@ class DataSet(object):
         self.shape = numpy.shape(array)
         self.number_of_dimensions = len (self.shape)
         self.dimension_labels = {}
-                
+        self.geometry = None # Only relevant for SinogramData and VolumeData
+        
         if dimension_labels is not None and \
            len (dimension_labels) == self.number_of_dimensions:
             for i in range(self.number_of_dimensions):
@@ -123,6 +124,13 @@ class DataSet(object):
             raise TypeError('Array must be NumpyArray, passed {0}'\
                             .format(type(array)))
         
+        # finally copy the geometry
+        if 'geometry' in kwargs.keys():
+            self.geometry = kwargs['geometry']
+        else:
+            # assume it is parallel beam
+            pass
+        
 
     def as_array(self, dimensions=None):
         '''Returns the DataSet as Numpy Array
@@ -212,7 +220,8 @@ class DataSet(object):
                 out = self.as_array() + other.as_array()
                 return type(self)(out, 
                                deep_copy=True, 
-                               dimension_labels=self.dimension_labels)
+                               dimension_labels=self.dimension_labels,
+                               geometry=self.geometry)
             else:
                 raise ValueError('Wrong shape: {0} and {1}'.format(self.shape, 
                                  other.shape))
@@ -220,7 +229,8 @@ class DataSet(object):
             return type(self)(
                     self.as_array() + other, 
                                deep_copy=True, 
-                               dimension_labels=self.dimension_labels)
+                               dimension_labels=self.dimension_labels,
+                               geometry=self.geometry)
         else:
             raise TypeError('Cannot {0} DataSet with {1}'.format("add" ,
                             type(other)))
@@ -232,14 +242,16 @@ class DataSet(object):
                 out = self.as_array() - other.as_array()
                 return type(self)(out, 
                                deep_copy=True, 
-                               dimension_labels=self.dimension_labels)
+                               dimension_labels=self.dimension_labels,
+                               geometry=self.geometry)
             else:
                 raise ValueError('Wrong shape: {0} and {1}'.format(self.shape, 
                                  other.shape))
         elif isinstance(other, (int, float, complex)):
             return type(self)(self.as_array() - other, 
                                deep_copy=True, 
-                               dimension_labels=self.dimension_labels)
+                               dimension_labels=self.dimension_labels,
+                               geometry=self.geometry)
         else:
             raise TypeError('Cannot {0} DataSet with {1}'.format("subtract" ,
                             type(other)))
@@ -254,14 +266,16 @@ class DataSet(object):
                 out = self.as_array() / other.as_array()
                 return type(self)(out, 
                                deep_copy=True, 
-                               dimension_labels=self.dimension_labels)
+                               dimension_labels=self.dimension_labels,
+                               geometry=self.geometry)
             else:
                 raise ValueError('Wrong shape: {0} and {1}'.format(self.shape, 
                                  other.shape))
         elif isinstance(other, (int, float, complex)):
             return type(self)(self.as_array() / other, 
                                deep_copy=True, 
-                               dimension_labels=self.dimension_labels)
+                               dimension_labels=self.dimension_labels,
+                               geometry=self.geometry)
         else:
             raise TypeError('Cannot {0} DataSet with {1}'.format("divide" ,
                             type(other)))
@@ -273,14 +287,16 @@ class DataSet(object):
                 out = self.as_array() ** other.as_array()
                 return type(self)(out, 
                                deep_copy=True, 
-                               dimension_labels=self.dimension_labels)
+                               dimension_labels=self.dimension_labels,
+                               geometry=self.geometry)
             else:
                 raise ValueError('Wrong shape: {0} and {1}'.format(self.shape, 
                                  other.shape))
         elif isinstance(other, (int, float, complex)):
             return type(self)(self.as_array() ** other, 
                                deep_copy=True, 
-                               dimension_labels=self.dimension_labels)
+                               dimension_labels=self.dimension_labels,
+                               geometry=self.geometry)
         else:
             raise TypeError('Cannot {0} DataSet with {1}'.format("power" ,
                             type(other)))
@@ -292,14 +308,16 @@ class DataSet(object):
                 out = self.as_array() * other.as_array()
                 return type(self)(out, 
                                deep_copy=True, 
-                               dimension_labels=self.dimension_labels)
+                               dimension_labels=self.dimension_labels,
+                               geometry=self.geometry)
             else:
                 raise ValueError('Wrong shape: {0} and {1}'.format(self.shape, 
                                  other.shape))
         elif isinstance(other, (int, float, complex)):
             return type(self)(self.as_array() * other, 
                                deep_copy=True, 
-                               dimension_labels=self.dimension_labels)
+                               dimension_labels=self.dimension_labels,
+                               geometry=self.geometry)
         else:
             raise TypeError('Cannot {0} DataSet with {1}'.format("multiply" ,
                             type(other)))
@@ -315,19 +333,22 @@ class DataSet(object):
         out = numpy.abs(self.as_array() )
         return type(self)(out,
                        deep_copy=True, 
-                       dimension_labels=self.dimension_labels)
+                       dimension_labels=self.dimension_labels,
+                       geometry=self.geometry)
     
     def maximum(self,otherscalar):
         out = numpy.maximum(self.as_array(),otherscalar)
         return type(self)(out,
                        deep_copy=True, 
-                       dimension_labels=self.dimension_labels)
+                       dimension_labels=self.dimension_labels,
+                       geometry=self.geometry)
     
     def sign(self):
         out = numpy.sign(self.as_array() )
         return type(self)(out,
                        deep_copy=True, 
-                       dimension_labels=self.dimension_labels)
+                       dimension_labels=self.dimension_labels,
+                       geometry=self.geometry)
     
     # reverse operand
     def __radd__(self, other):
@@ -353,11 +374,13 @@ class DataSet(object):
         if isinstance(other, (int, float)) :
             fother = numpy.ones(numpy.shape(self.array)) * other
             return type(self)(fother ** self.array , 
-                           dimension_labels=self.dimension_labels)
+                           dimension_labels=self.dimension_labels,
+                           geometry=self.geometry)
         elif issubclass(other, DataSet):
             if self.checkDimensions(other):
                 return type(self)(other.as_array() ** self.array , 
-                           dimension_labels=self.dimension_labels)
+                           dimension_labels=self.dimension_labels,
+                           geometry=self.geometry)
             else:
                 raise ValueError('Dimensions do not match')
     # __rpow__
@@ -473,13 +496,7 @@ class SinogramData(DataSet):
                     dimension_labels = ['angle' , 
                                         'horizontal']
             DataSet.__init__(self, array, deep_copy, dimension_labels, **kwargs)
-        
-        # finally copy the instrument geometry
-        if 'instrument_geometry' in kwargs.keys():
-            self.instrument_geometry = kwargs['instrument_geometry']
-        else:
-            # assume it is parallel beam
-            pass
+            
             
 class DataSetProcessor(object):
     '''Defines a generic DataSet processor
diff --git a/Wrappers/Python/ccpi/reconstruction/astra_ops.py b/Wrappers/Python/ccpi/reconstruction/astra_ops.py
index bc6f5e1..6d527ba 100644
--- a/Wrappers/Python/ccpi/reconstruction/astra_ops.py
+++ b/Wrappers/Python/ccpi/reconstruction/astra_ops.py
@@ -26,6 +26,12 @@ class AstraProjectorSimple(Operator):
     def __init__(self, geomv, geomp, device):
         super(AstraProjectorSimple, self).__init__()
         
+        # Store our volume and sinogram geometries. Redundant with also 
+        # storing in ASTRA format below but needed to assign to 
+        # SinogramData in "direct" method and VolumeData in "adjoint" method
+        self.sinogram_geometry = geomp
+        self.volume_geometry = geomv
+        
         # ASTRA Volume geometry
         self.vol_geom = astra.create_vol_geom(geomv.voxel_num_x, \
                                               geomv.voxel_num_y, \
@@ -63,12 +69,12 @@ class AstraProjectorSimple(Operator):
         
         sinogram_id, DATA = astra.create_sino(IM.as_array(), self.proj_id)
         astra.data2d.delete(sinogram_id)
-        return SinogramData(DATA)
+        return SinogramData(DATA,geometry=self.sinogram_geometry)
     
     def adjoint(self, DATA):
         rec_id, IM = astra.create_backprojection(DATA.as_array(), self.proj_id)
         astra.data2d.delete(rec_id)
-        return VolumeData(IM)
+        return VolumeData(IM,geometry=self.volume_geometry)
     
     def delete(self):
         astra.data2d.delete(self.proj_id)
diff --git a/Wrappers/Python/test/simple_demo.py b/Wrappers/Python/test/simple_demo.py
index 7a6b049..7a28ae9 100644
--- a/Wrappers/Python/test/simple_demo.py
+++ b/Wrappers/Python/test/simple_demo.py
@@ -25,32 +25,26 @@ plt.show()
 
 vg = VolumeGeometry(N,N,None, 1,1,None)
 
-Phantom = VolumeData(x)
-#Phantom = VolumeData(x)
+Phantom = VolumeData(x,geometry=vg)
 
 # Set up measurement geometry
 angles_num = 20; # angles number
 angles = np.linspace(0,np.pi,angles_num,endpoint=False)
-
 det_w = 1.0
 det_num = N
-
 SourceOrig = 500
 OrigDetec = 0
 
 # Parallelbeam geometry test
 pg = SinogramGeometry('parallel','2D',angles,det_num,det_w)
 
-# Set up ASTRA projector
-#Aop = AstraProjector(vg, angles, N,'gpu')
-#Aop = AstraProjector(det_w, det_num, angles, projtype='parallel',device='gpu')
-
-Aop_old = AstraProjector(det_w, det_num, SourceOrig, 
-                     OrigDetec, angles, 
-                     N,'fanbeam','gpu') 
-
+# ASTRA operator using volume and sinogram geometries
 Aop = AstraProjectorSimple(vg, pg, 'gpu')
 
+# Unused old astra projector without geometry
+# Aop_old = AstraProjector(det_w, det_num, SourceOrig, 
+#                      OrigDetec, angles, 
+#                      N,'fanbeam','gpu') 
 
 # Try forward and backprojection
 b = Aop.direct(Phantom)
@@ -66,7 +60,7 @@ plt.show()
 f = Norm2sq(Aop,b,c=0.5)
 
 # Initial guess
-x_init = VolumeData(np.zeros(x.shape))
+x_init = VolumeData(np.zeros(x.shape),geometry=vg)
 
 # Run FISTA for least squares without regularization
 x_fista0, it0, timing0, criter0 = FISTA(x_init, f, None)