Merge branch 'multichannel_adapt' of https://github.com/vais-ral/CCPi-ReconstructionFramework into dataset_for_multichannel

6 files changed, 582 insertions, 122 deletions

diff --git a/Wrappers/Python/ccpi/astra/astra_processors.py b/Wrappers/Python/ccpi/astra/astra_processors.py
new file mode 100644
index 0000000..650e11b
--- /dev/null
+++ b/Wrappers/Python/ccpi/astra/astra_processors.py
@@ -0,0 +1,273 @@
+from ccpi.framework import DataSetProcessor, DataSet, VolumeData, SinogramData
+from ccpi.astra.astra_utils import convert_geometry_to_astra
+import astra
+import numpy
+
+
+class AstraForwardProjector(DataSetProcessor):
+    '''AstraForwardProjector
+    
+    Forward project VolumeDataSet to SinogramDataSet using ASTRA proj_id.
+    
+    Input: VolumeDataSet
+    Parameter: proj_id
+    Output: SinogramDataSet
+    '''
+    
+    def __init__(self,
+                 volume_geometry=None,
+                 sinogram_geometry=None,
+                 proj_id=None,
+                 device='cpu'):
+        kwargs = {
+                  'volume_geometry'  : volume_geometry,
+                  'sinogram_geometry'  : sinogram_geometry,
+                  'proj_id'  : proj_id,
+                  'device'  : device
+                  }
+        
+        #DataSetProcessor.__init__(self, **kwargs)
+        super(AstraForwardProjector, self).__init__(**kwargs)
+        
+        self.setVolumeGeometry(volume_geometry)
+        self.setSinogramGeometry(sinogram_geometry)
+        
+        # Set up ASTRA Volume and projection geometry, not to be stored in self
+        vol_geom, proj_geom = convert_geometry_to_astra(self.volume_geometry,
+                                                        self.sinogram_geometry)
+        
+        # ASTRA projector, to be stored
+        if device == 'cpu':
+            # Note that 'line' is only for parallel (2D) and only one option
+            self.setProjector(astra.create_projector('line', proj_geom, vol_geom) )
+        elif device == 'gpu':
+            self.setProjector(astra.create_projector('cuda', proj_geom, vol_geom) )
+        else:
+            NotImplemented
+    
+    def checkInput(self, dataset):
+        if dataset.number_of_dimensions == 3 or dataset.number_of_dimensions == 2:
+               return True
+        else:
+            raise ValueError("Expected input dimensions is 2 or 3, got {0}"\
+                             .format(dataset.number_of_dimensions))
+    
+    def setProjector(self, proj_id):
+        self.proj_id = proj_id
+    
+    def setVolumeGeometry(self, volume_geometry):
+        self.volume_geometry = volume_geometry
+    
+    def setSinogramGeometry(self, sinogram_geometry):
+        self.sinogram_geometry = sinogram_geometry
+    
+    def process(self):
+        IM = self.getInput()
+        sinogram_id, DATA = astra.create_sino(IM.as_array(), self.proj_id)
+        astra.data2d.delete(sinogram_id)
+        return SinogramData(DATA,geometry=self.sinogram_geometry)
+
+class AstraBackProjector(DataSetProcessor):
+    '''AstraBackProjector
+    
+    Back project SinogramDataSet to VolumeDataSet using ASTRA proj_id.
+    
+    Input: SinogramDataSet
+    Parameter: proj_id
+    Output: VolumeDataSet 
+    '''
+    
+    def __init__(self,
+                 volume_geometry=None,
+                 sinogram_geometry=None,
+                 proj_id=None,
+                 device='cpu'):
+        kwargs = {
+                  'volume_geometry'  : volume_geometry,
+                  'sinogram_geometry'  : sinogram_geometry,
+                  'proj_id'  : proj_id,
+                  'device'  : device
+                  }
+        
+        #DataSetProcessor.__init__(self, **kwargs)
+        super(AstraBackProjector, self).__init__(**kwargs)
+        
+        self.setVolumeGeometry(volume_geometry)
+        self.setSinogramGeometry(sinogram_geometry)
+        
+                # Set up ASTRA Volume and projection geometry, not to be stored in self
+        vol_geom, proj_geom = convert_geometry_to_astra(self.volume_geometry,
+                                                        self.sinogram_geometry)
+        
+        # ASTRA projector, to be stored
+        if device == 'cpu':
+            # Note that 'line' is only for parallel (2D) and only one option
+            self.setProjector(astra.create_projector('line', proj_geom, vol_geom) )
+        elif device == 'gpu':
+            self.setProjector(astra.create_projector('cuda', proj_geom, vol_geom) )
+        else:
+            NotImplemented
+    
+    def checkInput(self, dataset):
+        if dataset.number_of_dimensions == 3 or dataset.number_of_dimensions == 2:
+               return True
+        else:
+            raise ValueError("Expected input dimensions is 2 or 3, got {0}"\
+                             .format(dataset.number_of_dimensions))
+    
+    def setProjector(self, proj_id):
+        self.proj_id = proj_id
+        
+    def setVolumeGeometry(self, volume_geometry):
+        self.volume_geometry = volume_geometry
+        
+    def setSinogramGeometry(self, sinogram_geometry):
+        self.sinogram_geometry = sinogram_geometry
+    
+    def process(self):
+        DATA = self.getInput()
+        rec_id, IM = astra.create_backprojection(DATA.as_array(), self.proj_id)
+        astra.data2d.delete(rec_id)
+        return VolumeData(IM,geometry=self.volume_geometry)
+
+class AstraForwardProjectorMC(DataSetProcessor):
+    '''AstraForwardProjector
+    
+    Forward project VolumeDataSet to SinogramDataSet using ASTRA proj_id.
+    
+    Input: VolumeDataSet
+    Parameter: proj_id
+    Output: SinogramDataSet
+    '''
+    
+    def __init__(self,
+                 volume_geometry=None,
+                 sinogram_geometry=None,
+                 proj_id=None,
+                 device='cpu'):
+        kwargs = {
+                  'volume_geometry'  : volume_geometry,
+                  'sinogram_geometry'  : sinogram_geometry,
+                  'proj_id'  : proj_id,
+                  'device'  : device
+                  }
+        
+        #DataSetProcessor.__init__(self, **kwargs)
+        super(AstraForwardProjectorMC, self).__init__(**kwargs)
+        
+        self.setVolumeGeometry(volume_geometry)
+        self.setSinogramGeometry(sinogram_geometry)
+        
+        # Set up ASTRA Volume and projection geometry, not to be stored in self
+        vol_geom, proj_geom = convert_geometry_to_astra(self.volume_geometry,
+                                                        self.sinogram_geometry)
+        
+        # ASTRA projector, to be stored
+        if device == 'cpu':
+            # Note that 'line' is only for parallel (2D) and only one option
+            self.setProjector(astra.create_projector('line', proj_geom, vol_geom) )
+        elif device == 'gpu':
+            self.setProjector(astra.create_projector('cuda', proj_geom, vol_geom) )
+        else:
+            NotImplemented
+    
+    def checkInput(self, dataset):
+        if dataset.number_of_dimensions == 3 or dataset.number_of_dimensions == 2:
+               return True
+        else:
+            return True
+            #raise ValueError("Expected input dimensions is 2 or 3, got {0}"\
+            #                 .format(dataset.number_of_dimensions))
+    
+    def setProjector(self, proj_id):
+        self.proj_id = proj_id
+    
+    def setVolumeGeometry(self, volume_geometry):
+        self.volume_geometry = volume_geometry
+    
+    def setSinogramGeometry(self, sinogram_geometry):
+        self.sinogram_geometry = sinogram_geometry
+    
+    def process(self):
+        IM = self.getInput()
+        DATA = numpy.zeros((self.sinogram_geometry.angles.size,
+                                self.sinogram_geometry.pixel_num_h,
+                                1,
+                                self.sinogram_geometry.channels),
+                                'float32')
+        for k in range(self.sinogram_geometry.channels):
+            sinogram_id, DATA[:,:,0,k] = astra.create_sino(IM.as_array()[:,:,0,k], 
+                                                           self.proj_id)
+            astra.data2d.delete(sinogram_id)
+        return SinogramData(DATA,geometry=self.sinogram_geometry)
+
+class AstraBackProjectorMC(DataSetProcessor):
+    '''AstraBackProjector
+    
+    Back project SinogramDataSet to VolumeDataSet using ASTRA proj_id.
+    
+    Input: SinogramDataSet
+    Parameter: proj_id
+    Output: VolumeDataSet 
+    '''
+    
+    def __init__(self,
+                 volume_geometry=None,
+                 sinogram_geometry=None,
+                 proj_id=None,
+                 device='cpu'):
+        kwargs = {
+                  'volume_geometry'  : volume_geometry,
+                  'sinogram_geometry'  : sinogram_geometry,
+                  'proj_id'  : proj_id,
+                  'device'  : device
+                  }
+        
+        #DataSetProcessor.__init__(self, **kwargs)
+        super(AstraBackProjectorMC, self).__init__(**kwargs)
+        
+        self.setVolumeGeometry(volume_geometry)
+        self.setSinogramGeometry(sinogram_geometry)
+        
+                # Set up ASTRA Volume and projection geometry, not to be stored in self
+        vol_geom, proj_geom = convert_geometry_to_astra(self.volume_geometry,
+                                                        self.sinogram_geometry)
+        
+        # ASTRA projector, to be stored
+        if device == 'cpu':
+            # Note that 'line' is only for parallel (2D) and only one option
+            self.setProjector(astra.create_projector('line', proj_geom, vol_geom) )
+        elif device == 'gpu':
+            self.setProjector(astra.create_projector('cuda', proj_geom, vol_geom) )
+        else:
+            NotImplemented
+    
+    def checkInput(self, dataset):
+        if dataset.number_of_dimensions == 3 or dataset.number_of_dimensions == 2:
+               return True
+        else:
+            return True
+            #raise ValueError("Expected input dimensions is 2 or 3, got {0}"\
+            #                 .format(dataset.number_of_dimensions))
+    
+    def setProjector(self, proj_id):
+        self.proj_id = proj_id
+        
+    def setVolumeGeometry(self, volume_geometry):
+        self.volume_geometry = volume_geometry
+        
+    def setSinogramGeometry(self, sinogram_geometry):
+        self.sinogram_geometry = sinogram_geometry
+    
+    def process(self):
+        DATA = self.getInput()
+        IM = numpy.zeros((self.volume_geometry.voxel_num_x,
+                                self.volume_geometry.voxel_num_y,
+                                1,
+                                self.volume_geometry.channels),
+                                'float32')
+        for k in range(self.volume_geometry.channels):
+            rec_id, IM[:,:,0,k] = astra.create_backprojection(DATA.as_array()[:,:,0,k], 
+                                  self.proj_id)
+            astra.data2d.delete(rec_id)
+        return VolumeData(IM,geometry=self.volume_geometry)
\ No newline at end of file
diff --git a/Wrappers/Python/ccpi/astra/astra_utils.py b/Wrappers/Python/ccpi/astra/astra_utils.py
new file mode 100644
index 0000000..8b5043a
--- /dev/null
+++ b/Wrappers/Python/ccpi/astra/astra_utils.py
@@ -0,0 +1,61 @@
+import astra
+
+def convert_geometry_to_astra(volume_geometry, sinogram_geometry):
+    # Set up ASTRA Volume and projection geometry, not stored
+    if sinogram_geometry.dimension == '2D':
+        vol_geom = astra.create_vol_geom(volume_geometry.voxel_num_x, 
+                                         volume_geometry.voxel_num_y, 
+                                         volume_geometry.getMinX(), 
+                                         volume_geometry.getMaxX(), 
+                                         volume_geometry.getMinY(), 
+                                         volume_geometry.getMaxY())
+        
+        if sinogram_geometry.geom_type == 'parallel':
+            proj_geom = astra.create_proj_geom('parallel',
+                                               sinogram_geometry.pixel_size_h,
+                                               sinogram_geometry.pixel_num_h,
+                                               sinogram_geometry.angles)
+        elif sinogram_geometry.geom_type == 'cone':
+            proj_geom = astra.create_proj_geom('fanflat',
+                                               sinogram_geometry.pixel_size_h,
+                                               sinogram_geometry.pixel_num_h,
+                                               sinogram_geometry.angles,
+                                               sinogram_geometry.dist_source_center,
+                                               sinogram_geometry.dist_center_detector)
+        else:
+            NotImplemented
+            
+    elif sinogram_geometry.dimension == '3D':
+        vol_geom = astra.create_vol_geom(volume_geometry.voxel_num_x, 
+                                         volume_geometry.voxel_num_y, 
+                                         volume_geometry.voxel_num_z, 
+                                         volume_geometry.getMinX(), 
+                                         volume_geometry.getMaxX(), 
+                                         volume_geometry.getMinY(), 
+                                         volume_geometry.getMaxY(), 
+                                         volume_geometry.getMinZ(), 
+                                         volume_geometry.getMaxZ())
+        
+        if sinogram_geometry.proj_geom == 'parallel':
+            proj_geom = astra.create_proj_geom('parallel3d',
+                                               sinogram_geometry.pixel_size_h,
+                                               sinogram_geometry.pixel_size_v,
+                                               sinogram_geometry.pixel_num_v,
+                                               sinogram_geometry.pixel_num_h,
+                                               sinogram_geometry.angles)
+        elif sinogram_geometry.geom_type == 'cone':
+            proj_geom = astra.create_proj_geom('cone',
+                                               sinogram_geometry.pixel_size_h,
+                                               sinogram_geometry.pixel_size_v,
+                                               sinogram_geometry.pixel_num_v,
+                                               sinogram_geometry.pixel_num_h,
+                                               sinogram_geometry.angles,
+                                               sinogram_geometry.dist_source_center,
+                                               sinogram_geometry.dist_center_detector)
+        else:
+            NotImplemented
+            
+    else:
+        NotImplemented
+    
+    return vol_geom, proj_geom
\ No newline at end of file
diff --git a/Wrappers/Python/ccpi/framework.py b/Wrappers/Python/ccpi/framework.py
index 754e0ee..82ad65c 100644
--- a/Wrappers/Python/ccpi/framework.py
+++ b/Wrappers/Python/ccpi/framework.py
@@ -437,20 +437,20 @@ class VolumeData(DataSet):
         
         if type(array) == DataSet:
             # if the array is a DataSet get the info from there
-            if not ( array.number_of_dimensions == 2 or \
-                     array.number_of_dimensions == 3 ):
-                raise ValueError('Number of dimensions are not 2 or 3: {0}'\
-                                 .format(array.number_of_dimensions))
+            #if not ( array.number_of_dimensions == 2 or \
+            #         array.number_of_dimensions == 3 ):
+            #    raise ValueError('Number of dimensions are not 2 or 3: {0}'\
+            #                     .format(array.number_of_dimensions))
             
             #DataSet.__init__(self, array.as_array(), deep_copy,
             #                 array.dimension_labels, **kwargs)
             super(VolumeData, self).__init__(array.as_array(), deep_copy,
                              array.dimension_labels, **kwargs)
         elif type(array) == numpy.ndarray:
-            if not ( array.ndim == 3 or array.ndim == 2 ):
-                raise ValueError(
-                        'Number of dimensions are not 3 or 2 : {0}'\
-                        .format(array.ndim))
+            #if not ( array.ndim == 3 or array.ndim == 2 ):
+            #    raise ValueError(
+            #            'Number of dimensions are not 3 or 2 : {0}'\
+            #            .format(array.ndim))
                 
             if dimension_labels is None:
                 if array.ndim == 3:
@@ -485,17 +485,17 @@ class SinogramData(DataSet):
         
         if type(array) == DataSet:
             # if the array is a DataSet get the info from there
-            if not ( array.number_of_dimensions == 2 or \
-                     array.number_of_dimensions == 3 ):
-                raise ValueError('Number of dimensions are not 2 or 3: {0}'\
-                                 .format(array.number_of_dimensions))
+            #if not ( array.number_of_dimensions == 2 or \
+            #         array.number_of_dimensions == 3 ):
+            #    raise ValueError('Number of dimensions are not 2 or 3: {0}'\
+            #                     .format(array.number_of_dimensions))
             
             DataSet.__init__(self, array.as_array(), deep_copy,
                              array.dimension_labels, **kwargs)
         elif type(array) == numpy.ndarray:
-            if not ( array.ndim == 3 or array.ndim == 2 ):
-                raise ValueError('Number of dimensions are != 3: {0}'\
-                                 .format(array.ndim))
+            #if not ( array.ndim == 3 or array.ndim == 2 ):
+            #    raise ValueError('Number of dimensions are != 3: {0}'\
+            #                     .format(array.ndim))
             if dimension_labels is None:
                 if array.ndim == 3:
                     dimension_labels = ['angle' , 
diff --git a/Wrappers/Python/ccpi/reconstruction/astra_ops.py b/Wrappers/Python/ccpi/reconstruction/astra_ops.py
index 452c86a..1346f22 100644
--- a/Wrappers/Python/ccpi/reconstruction/astra_ops.py
+++ b/Wrappers/Python/ccpi/reconstruction/astra_ops.py
@@ -16,118 +16,105 @@
 #    along with this program.  If not, see <http://www.gnu.org/licenses/>.
 
 from ccpi.reconstruction.ops import Operator
-import astra
 import numpy
 from ccpi.framework import SinogramData, VolumeData
 from ccpi.reconstruction.ops import PowerMethodNonsquare
+from ccpi.astra.astra_processors import *
 
 class AstraProjectorSimple(Operator):
     """ASTRA projector modified to use DataSet and geometry."""
     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
+        # Store volume and sinogram geometries.
         self.sinogram_geometry = geomp
         self.volume_geometry = geomv
         
-        # ASTRA Volume geometry
-        if geomp.dimension == '2D':
-            self.vol_geom = astra.create_vol_geom(geomv.voxel_num_x, 
-                                                  geomv.voxel_num_y, 
-                                                  geomv.getMinX(), 
-                                                  geomv.getMaxX(), 
-                                                  geomv.getMinY(), 
-                                                  geomv.getMaxY())
-        elif geomp.dimension == '3D':
-            self.vol_geom = astra.create_vol_geom(geomv.voxel_num_x, 
-                                                  geomv.voxel_num_y, 
-                                                  geomv.voxel_num_z, 
-                                                  geomv.getMinX(), 
-                                                  geomv.getMaxX(), 
-                                                  geomv.getMinY(), 
-                                                  geomv.getMaxY(), 
-                                                  geomv.getMinZ(), 
-                                                  geomv.getMaxZ())
-        else:
-            NotImplemented
-            
-        
-        # ASTRA Projections geometry
-        if geomp.dimension == '2D':
-            if geomp.geom_type == 'parallel':
-                self.proj_geom = astra.create_proj_geom('parallel',
-                                                        geomp.pixel_size_h,
-                                                        geomp.pixel_num_h,
-                                                        geomp.angles)
-            elif geomp.geom_type == 'cone':
-                self.proj_geom = astra.create_proj_geom('fanflat',
-                                                        geomp.pixel_size_h,
-                                                        geomp.pixel_num_h,
-                                                        geomp.angles,
-                                                        geomp.dist_source_center,
-                                                        geomp.dist_center_detector)
-            else:
-                NotImplemented
-        elif geomp.dimension == '3D':
-            if geomp.proj_geom == 'parallel':
-                self.proj_geom = astra.create_proj_geom('parallel3d',
-                                                        geomp.pixel_size_h,
-                                                        geomp.pixel_size_v,
-                                                        geomp.pixel_num_v,
-                                                        geomp.pixel_num_h,
-                                                        geomp.angles)
-            elif geomp.geom_type == 'cone':
-                self.proj_geom = astra.create_proj_geom('cone',
-                                                        geomp.pixel_size_h,
-                                                        geomp.pixel_size_v,
-                                                        geomp.pixel_num_v,
-                                                        geomp.pixel_num_h,
-                                                        geomp.angles,
-                                                        geomp.dist_source_center,
-                                                        geomp.dist_center_detector)
-            else:
-                NotImplemented
-        else:
-            NotImplemented
-        
-        # ASTRA projector
-        if device == 'cpu':
-            # Note that 'line' is only for parallel (2D) and only one option
-            self.proj_id = astra.create_projector('line', self.proj_geom, 
-                                                  self.vol_geom) # for CPU
-        elif device == 'gpu':
-            self.proj_id = astra.create_projector('cuda', self.proj_geom, 
-                                                  self.vol_geom) # for GPU
-        else:
-            NotImplemented
+        self.fp = AstraForwardProjector(volume_geometry=geomv,
+                                        sinogram_geometry=geomp,
+                                        proj_id=None,
+                                        device=device)
         
+        self.bp = AstraBackProjector(volume_geometry=geomv,
+                                        sinogram_geometry=geomp,
+                                        proj_id=None,
+                                        device=device)
+                
+        # Initialise empty for singular value.
         self.s1 = None
     
     def direct(self, IM):
-        
-        sinogram_id, DATA = astra.create_sino(IM.as_array(), self.proj_id)
-        astra.data2d.delete(sinogram_id)
-        return SinogramData(DATA,geometry=self.sinogram_geometry)
+        self.fp.setInput(IM)
+        out = self.fp.getOutput()
+        return out
     
     def adjoint(self, DATA):
-        rec_id, IM = astra.create_backprojection(DATA.as_array(), self.proj_id)
-        astra.data2d.delete(rec_id)
-        return VolumeData(IM,geometry=self.volume_geometry)
+        self.bp.setInput(DATA)
+        out = self.bp.getOutput()
+        return out
     
-    def delete(self):
-        astra.data2d.delete(self.proj_id)
+    #def delete(self):
+    #    astra.data2d.delete(self.proj_id)
     
     def get_max_sing_val(self):
         self.s1, sall, svec = PowerMethodNonsquare(self,10)
         return self.s1
     
     def size(self):
-        return ( (self.proj_geom['ProjectionAngles'].size, \
-                  self.proj_geom['DetectorCount']), \
-                 (self.vol_geom['GridColCount'], \
-                  self.vol_geom['GridRowCount']) )
+        # Only implemented for 2D
+        return ( (self.sinogram_geometry.angles.size, \
+                  self.sinogram_geometry.pixel_num_h), \
+                 (self.volume_geometry.voxel_num_x, \
+                  self.volume_geometry.voxel_num_y) )
+
+class AstraProjectorMC(Operator):
+    """ASTRA Multichannel projector"""
+    def __init__(self, geomv, geomp, device):
+        super(AstraProjectorMC, self).__init__()
+        
+        # Store volume and sinogram geometries.
+        self.sinogram_geometry = geomp
+        self.volume_geometry = geomv
+        
+        self.fp = AstraForwardProjectorMC(volume_geometry=geomv,
+                                        sinogram_geometry=geomp,
+                                        proj_id=None,
+                                        device=device)
+        
+        self.bp = AstraBackProjectorMC(volume_geometry=geomv,
+                                        sinogram_geometry=geomp,
+                                        proj_id=None,
+                                        device=device)
+                
+        # Initialise empty for singular value.
+        self.s1 = 50
+    
+    def direct(self, IM):
+        self.fp.setInput(IM)
+        out = self.fp.getOutput()
+        return out
+    
+    def adjoint(self, DATA):
+        self.bp.setInput(DATA)
+        out = self.bp.getOutput()
+        return out
+    
+    #def delete(self):
+    #    astra.data2d.delete(self.proj_id)
+    
+    def get_max_sing_val(self):
+        if self.s1 is None:
+            self.s1, sall, svec = PowerMethodNonsquare(self,10)
+            return self.s1
+        else:
+            return self.s1
+    
+    def size(self):
+        # Only implemented for 2D
+        return ( (self.sinogram_geometry.angles.size, \
+                  self.sinogram_geometry.pixel_num_h), \
+                 (self.volume_geometry.voxel_num_x, \
+                  self.volume_geometry.voxel_num_y) )
     
 
 class AstraProjector(Operator):
diff --git a/Wrappers/Python/ccpi/reconstruction/geoms.py b/Wrappers/Python/ccpi/reconstruction/geoms.py
index edce3b3..9f81fa0 100644
--- a/Wrappers/Python/ccpi/reconstruction/geoms.py
+++ b/Wrappers/Python/ccpi/reconstruction/geoms.py
@@ -1,16 +1,17 @@
 
 class VolumeGeometry:
     
-    def __init__(self, \
-                 voxel_num_x=None, \
-                 voxel_num_y=None, \
-                 voxel_num_z=None, \
-                 voxel_size_x=None, \
-                 voxel_size_y=None, \
-                 voxel_size_z=None, \
-                 center_x=0, \
-                 center_y=0, \
-                 center_z=0):
+    def __init__(self, 
+                 voxel_num_x=None, 
+                 voxel_num_y=None, 
+                 voxel_num_z=None, 
+                 voxel_size_x=None, 
+                 voxel_size_y=None, 
+                 voxel_size_z=None, 
+                 center_x=0, 
+                 center_y=0, 
+                 center_z=0, 
+                 channels=1):
         
         self.voxel_num_x = voxel_num_x
         self.voxel_num_y = voxel_num_y
@@ -21,6 +22,7 @@ class VolumeGeometry:
         self.center_x = center_x
         self.center_y = center_y
         self.center_z = center_z  
+        self.channels = channels
         
     def getMinX(self):
         return self.center_x - 0.5*self.voxel_num_x*self.voxel_size_x
@@ -43,16 +45,17 @@ class VolumeGeometry:
     
 class SinogramGeometry:
     
-    def __init__(self, \
-                 geom_type, \
-                 dimension, \
-                 angles, \
-                 pixel_num_h=None, \
-                 pixel_size_h=1, \
-                 pixel_num_v=None, \
-                 pixel_size_v=1, \
-                 dist_source_center=None, \
-                 dist_center_detector=None, \
+    def __init__(self, 
+                 geom_type, 
+                 dimension, 
+                 angles, 
+                 pixel_num_h=None, 
+                 pixel_size_h=1, 
+                 pixel_num_v=None, 
+                 pixel_size_v=1, 
+                 dist_source_center=None, 
+                 dist_center_detector=None, 
+                 channels=1
                  ):
         """
         General inputs for standard type projection geometries
@@ -91,6 +94,8 @@ class SinogramGeometry:
         self.pixel_size_h = pixel_size_h
         self.pixel_num_v = pixel_num_v
         self.pixel_size_v = pixel_size_v
+        
+        self.channels = channels
 
         
                 
diff --git a/Wrappers/Python/test/simple_mc_demo.py b/Wrappers/Python/test/simple_mc_demo.py
new file mode 100644
index 0000000..a6959f9
--- /dev/null
+++ b/Wrappers/Python/test/simple_mc_demo.py
@@ -0,0 +1,134 @@
+
+
+import sys
+
+sys.path.append("..")
+
+from ccpi.framework import *
+from ccpi.reconstruction.algs import *
+from ccpi.reconstruction.funcs import *
+from ccpi.reconstruction.ops import *
+from ccpi.reconstruction.astra_ops import *
+from ccpi.reconstruction.geoms import *
+
+import numpy as np
+import matplotlib.pyplot as plt
+
+test_case = 2   # 1=parallel2D, 2=cone2D
+
+# Set up phantom
+N = 128
+
+numchannels = 3
+
+x = np.zeros((N,N,1,numchannels))
+
+x[round(N/4):round(3*N/4),round(N/4):round(3*N/4),:,0] = 1.0
+x[round(N/8):round(7*N/8),round(3*N/8):round(5*N/8),:,0] = 2.0
+
+x[round(N/4):round(3*N/4),round(N/4):round(3*N/4),:,1] = 0.7
+x[round(N/8):round(7*N/8),round(3*N/8):round(5*N/8),:,1] = 1.2
+
+x[round(N/4):round(3*N/4),round(N/4):round(3*N/4),:,2] = 1.5
+x[round(N/8):round(7*N/8),round(3*N/8):round(5*N/8),:,2] = 2.2
+
+f, axarr = plt.subplots(1,numchannels)
+for k in numpy.arange(3):
+    axarr[k].imshow(x[:,:,0,k],vmin=0,vmax=2.5)
+plt.show()
+
+vg = VolumeGeometry(N,N,None, 1,1,None,channels=numchannels)
+
+
+Phantom = VolumeData(x,geometry=vg)
+
+# Set up measurement geometry
+angles_num = 20; # angles number
+
+if test_case==1:
+    angles = np.linspace(0,np.pi,angles_num,endpoint=False)
+elif test_case==2:
+    angles = np.linspace(0,2*np.pi,angles_num,endpoint=False)
+else:
+    NotImplemented
+
+det_w = 1.0
+det_num = N
+SourceOrig = 200
+OrigDetec = 0
+
+# Parallelbeam geometry test
+if test_case==1:
+    pg = SinogramGeometry('parallel',
+                          '2D',
+                          angles,
+                          det_num,
+                          det_w,
+                          channels=numchannels)
+elif test_case==2:
+    pg = SinogramGeometry('cone',
+                          '2D',
+                          angles,
+                          det_num,
+                          det_w,
+                          dist_source_center=SourceOrig, 
+                          dist_center_detector=OrigDetec,
+                          channels=numchannels)
+
+# ASTRA operator using volume and sinogram geometries
+Aop = AstraProjectorMC(vg, pg, 'gpu')
+
+
+# Try forward and backprojection
+b = Aop.direct(Phantom)
+
+fb, axarrb = plt.subplots(1,numchannels)
+for k in numpy.arange(3):
+    axarrb[k].imshow(b.as_array()[:,:,0,k],vmin=0,vmax=250)
+plt.show()
+
+out2 = Aop.adjoint(b)
+
+fo, axarro = plt.subplots(1,numchannels)
+for k in range(3):
+    axarro[k].imshow(out2.as_array()[:,:,0,k],vmin=0,vmax=3500)
+plt.show()
+
+# Create least squares object instance with projector and data.
+f = Norm2sq(Aop,b,c=0.5)
+
+# Initial guess
+x_init = VolumeData(np.zeros(x.shape),geometry=vg)
+
+# FISTA options
+opt = {'tol': 1e-4, 'iter': 200}
+
+# Run FISTA for least squares without regularization
+x_fista0, it0, timing0, criter0 = FISTA(x_init, f, None, opt)
+
+
+ff0, axarrf0 = plt.subplots(1,numchannels)
+for k in numpy.arange(3):
+    axarrf0[k].imshow(x_fista0.as_array()[:,:,0,k],vmin=0,vmax=2.5)
+plt.show()
+
+plt.semilogy(criter0)
+plt.title('Criterion vs iterations, least squares')
+plt.show()
+
+# Now least squares plus 1-norm regularization
+lam = 0.1
+g0 = Norm1(lam)
+
+
+# Run FISTA for least squares plus 1-norm function.
+x_fista1, it1, timing1, criter1 = FISTA(x_init, f, g0, opt)
+
+ff1, axarrf1 = plt.subplots(1,numchannels)
+for k in numpy.arange(3):
+    axarrf1[k].imshow(x_fista1.as_array()[:,:,0,k],vmin=0,vmax=2.5)
+plt.show()
+
+plt.semilogy(criter1)
+plt.title('Criterion vs iterations, least squares plus 1-norm regu')
+plt.show()
\ No newline at end of file