Merge pull request #47 from vais-ral/dataset_for_multichannel

Dataset for multichannel

10 files changed, 501 insertions, 111 deletions

diff --git a/Wrappers/Python/ccpi/astra/astra_processors.py b/Wrappers/Python/ccpi/astra/astra_processors.py
index c3127eb..f51aec9 100644
--- a/Wrappers/Python/ccpi/astra/astra_processors.py
+++ b/Wrappers/Python/ccpi/astra/astra_processors.py
@@ -1,6 +1,7 @@
 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):
@@ -50,7 +51,8 @@ class AstraForwardProjector(DataSetProcessor):
             NotImplemented
     
     def checkInput(self, dataset):
-        if dataset.number_of_dimensions == 3 or dataset.number_of_dimensions == 2:
+        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}"\
@@ -67,9 +69,14 @@ class AstraForwardProjector(DataSetProcessor):
     
     def process(self):
         IM = self.getInput()
-        sinogram_id, DATA = astra.create_sino(IM.as_array(), self.proj_id)
+        DATA = SinogramData(geometry=self.sinogram_geometry)
+        #sinogram_id, DATA = astra.create_sino( IM.as_array(), 
+        #                           self.proj_id)
+        sinogram_id, DATA.array = astra.create_sino(IM.as_array(), 
+                                                           self.proj_id)
         astra.data2d.delete(sinogram_id)
-        return SinogramData(DATA,geometry=self.sinogram_geometry)
+        #return SinogramData(array=DATA, geometry=self.sinogram_geometry)
+        return DATA
 
 class AstraBackProjector(DataSetProcessor):
     '''AstraBackProjector
@@ -135,6 +142,64 @@ class AstraBackProjector(DataSetProcessor):
     
     def process(self):
         DATA = self.getInput()
-        rec_id, IM = astra.create_backprojection(DATA.as_array(), self.proj_id)
+        IM = VolumeData(geometry=self.volume_geometry)
+        rec_id, IM.array = astra.create_backprojection(DATA.as_array(),
+                            self.proj_id)
         astra.data2d.delete(rec_id)
-        return VolumeData(IM,geometry=self.volume_geometry)
\ No newline at end of file
+        return IM
+
+class AstraForwardProjectorMC(AstraForwardProjector):
+    '''AstraForwardProjector Multi channel
+    
+    Forward project VolumeDataSet to SinogramDataSet using ASTRA proj_id.
+    
+    Input: VolumeDataSet
+    Parameter: proj_id
+    Output: SinogramDataSet
+    '''
+    def checkInput(self, dataset):
+        if dataset.number_of_dimensions == 2 or \
+           dataset.number_of_dimensions == 3 or \
+           dataset.number_of_dimensions == 4:
+               return True
+        else:
+            raise ValueError("Expected input dimensions is 2 or 3, got {0}"\
+                             .format(dataset.number_of_dimensions))
+    def process(self):
+        IM = self.getInput()
+        #create the output Sinogram
+        DATA = SinogramData(geometry=self.sinogram_geometry)
+        
+        for k in range(DATA.geometry.channels):
+            sinogram_id, DATA.as_array()[k] = astra.create_sino(IM.as_array()[k], 
+                                                           self.proj_id)
+            astra.data2d.delete(sinogram_id)
+        return DATA
+
+class AstraBackProjectorMC(AstraBackProjector):
+    '''AstraBackProjector Multi channel
+    
+    Back project SinogramDataSet to VolumeDataSet using ASTRA proj_id.
+    
+    Input: SinogramDataSet
+    Parameter: proj_id
+    Output: VolumeDataSet 
+    '''
+    def checkInput(self, dataset):
+        if dataset.number_of_dimensions == 2 or \
+           dataset.number_of_dimensions == 3 or \
+           dataset.number_of_dimensions == 4:
+               return True
+        else:
+            raise ValueError("Expected input dimensions is 2 or 3, got {0}"\
+                             .format(dataset.number_of_dimensions))
+    def process(self):
+        DATA = self.getInput()
+        
+        IM = VolumeData(geometry=self.volume_geometry)
+        
+        for k in range(IM.geometry.channels):
+            rec_id, IM.as_array()[k] = astra.create_backprojection(DATA.as_array()[k], 
+                                  self.proj_id)
+            astra.data2d.delete(rec_id)
+        return IM
\ No newline at end of file
diff --git a/Wrappers/Python/ccpi/framework.py b/Wrappers/Python/ccpi/framework.py
index 5ac17ee..5a507d9 100644
--- a/Wrappers/Python/ccpi/framework.py
+++ b/Wrappers/Python/ccpi/framework.py
@@ -23,6 +23,7 @@ import numpy
 import sys
 from datetime import timedelta, datetime
 import warnings
+from ccpi.reconstruction import geoms
 
 if sys.version_info[0] >= 3 and sys.version_info[1] >= 4:
     ABC = abc.ABC
@@ -119,7 +120,7 @@ class DataSet(object):
             if deep_copy:
                 self.array = array[:]
             else:
-                self.array = array
+                self.array = array    
         else:
             raise TypeError('Array must be NumpyArray, passed {0}'\
                             .format(type(array)))
@@ -141,8 +142,9 @@ class DataSet(object):
         if dimensions is not None:
             return self.subset(dimensions).as_array()
         return self.array
-        
-    def subset(self, dimensions=None):
+    
+    
+    def subset(self, dimensions=None, **kw):
         '''Creates a DataSet containing a subset of self according to the 
         labels in dimensions'''
         if dimensions is None:
@@ -164,7 +166,7 @@ class DataSet(object):
                     else:
                         axis_order.append(find_key(self.dimension_labels, dl))
                 if not proceed:
-                    raise KeyError('Unknown key specified {0}'.format(dl))
+                    raise KeyError('Subset error: Unknown key specified {0}'.format(dl))
                     
                 # slice away the unwanted data from the array
                 unwanted_dimensions = self.dimension_labels.copy()
@@ -176,13 +178,21 @@ class DataSet(object):
                 #print ("left_dimensions {0}".format(left_dimensions))
                 #new_shape = [self.shape[ax] for ax in axis_order]
                 #print ("new_shape {0}".format(new_shape))
-                command = "self.array"
+                command = "self.array["
                 for i in range(self.number_of_dimensions):
                     if self.dimension_labels[i] in unwanted_dimensions.values():
-                        command = command + "[0]"
+                        value = 0
+                        for k,v in kw.items():
+                            if k == self.dimension_labels[i]:
+                                value = v
+                                
+                        command = command + str(value)
                     else:
-                        command = command + "[:]"
-                #print ("command {0}".format(command))
+                        command = command + ":"
+                    if i < self.number_of_dimensions -1:
+                        command = command + ','
+                command = command + ']'
+                
                 cleaned = eval(command)
                 # cleaned has collapsed dimensions in the same order of
                 # self.array, but we want it in the order stated in the 
@@ -413,7 +423,7 @@ class DataSet(object):
         repres += "Number of dimensions: {0}\n".format(self.number_of_dimensions)
         repres += "Shape: {0}\n".format(self.shape)
         repres += "Axis labels: {0}\n".format(self.dimension_labels)
-        repres += "Representation: {0}\n".format(self.array)
+        repres += "Representation: \n{0}\n".format(self.array)
         return repres
                 
                     
@@ -421,40 +431,81 @@ class DataSet(object):
 class VolumeData(DataSet):
     '''DataSet for holding 2D or 3D dataset'''
     def __init__(self, 
-                 array, 
+                 array = None, 
                  deep_copy=True, 
                  dimension_labels=None, 
                  **kwargs):
         
-        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))
-            
-            #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))
+        self.geometry = None
+        if array is None:
+            if 'geometry' in kwargs.keys():
+                geometry  = kwargs['geometry']
+                self.geometry = geometry
+                channels  = geometry.channels
+                horiz_x   = geometry.voxel_num_x
+                horiz_y   = geometry.voxel_num_y
+                vert      = 1 if geometry.voxel_num_z is None\
+                              else geometry.voxel_num_z # this should be 1 for 2D
                 
-            if dimension_labels is None:
-                if array.ndim == 3:
-                    dimension_labels = ['horizontal_x' , 
-                                        'horizontal_y' , 
-                                        'vertical']
+                if channels > 1:
+                    if vert > 1:
+                        shape = (channels, vert, horiz_y, horiz_x)
+                        dim_labels = ['channel' ,'vertical' , 'horizontal_y' , 
+                                      'horizontal_x']
+                    else:
+                        shape = (channels , horiz_y, horiz_x)
+                        dim_labels = ['channel' , 'horizontal_y' , 
+                                      'horizontal_x']
                 else:
-                    dimension_labels = ['horizontal' , 
-                                        'vertical']   
-            
-            #DataSet.__init__(self, array, deep_copy, dimension_labels, **kwargs)
-            super(VolumeData, self).__init__(array, deep_copy, 
-                 dimension_labels, **kwargs)
+                    if vert > 1:
+                        shape = (vert, horiz_y, horiz_x)
+                        dim_labels = ['vertical' , 'horizontal_y' , 
+                                      'horizontal_x']
+                    else:
+                        shape = (horiz_y, horiz_x)
+                        dim_labels = ['horizontal_y' , 
+                                      'horizontal_x']
+                        
+                array = numpy.zeros( shape , dtype=numpy.float32) 
+                super(VolumeData, self).__init__(array, deep_copy,
+                                 dim_labels, **kwargs)
+                
+            else:
+                raise ValueError('Please pass either a DataSet, ' +\
+                                 'a numpy array or a geometry')
+        else:
+            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 or \
+                         array.number_of_dimensions == 4):
+                    raise ValueError('Number of dimensions are not 2 or 3 or 4: {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 == 2 or array.ndim == 3 or array.ndim == 4 ):
+                    raise ValueError(
+                            'Number of dimensions are not 2 or 3 or 4 : {0}'\
+                            .format(array.ndim))
+                    
+                if dimension_labels is None:
+                    if array.ndim == 4:
+                        dimension_labels = ['channel' ,'vertical' , 'horizontal_y' , 
+                                      'horizontal_x']
+                    elif array.ndim == 3:
+                        dimension_labels = ['vertical' , 'horizontal_y' , 
+                                      'horizontal_x']
+                    else:
+                        dimension_labels = ['horizontal_y' , 
+                                      'horizontal_x']   
+                
+                #DataSet.__init__(self, array, deep_copy, dimension_labels, **kwargs)
+                super(VolumeData, self).__init__(array, deep_copy, 
+                     dimension_labels, **kwargs)
        
         # load metadata from kwargs if present
         for key, value in kwargs.items():
@@ -469,34 +520,79 @@ class VolumeData(DataSet):
 class SinogramData(DataSet):
     '''DataSet for holding 2D or 3D sinogram'''
     def __init__(self, 
-                 array, 
+                 array = None, 
                  deep_copy=True, 
                  dimension_labels=None, 
                  **kwargs):
-        
-        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))
-            
-            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 dimension_labels is None:
-                if array.ndim == 3:
-                    dimension_labels = ['angle' , 
-                                        'horizontal' , 
-                                        'vertical']
+        self.geometry = None
+        if array is None:
+            if 'geometry' in kwargs.keys():
+                geometry  = kwargs['geometry']
+                self.geometry = geometry
+                channels  = geometry.channels
+                horiz   = geometry.pixel_num_h
+                vert   = geometry.pixel_num_v
+                angles = geometry.angles
+                num_of_angles = numpy.shape(angles)[0]
+                
+                
+                if channels > 1:
+                    if vert > 1:
+                        shape = (channels, num_of_angles , vert, horiz)
+                        dim_labels = ['channel' , ' angle' ,
+                                      'vertical' , 'horizontal']
+                    else:
+                        shape = (channels , num_of_angles, horiz)
+                        dim_labels = ['channel' , 'angle' , 
+                                      'horizontal']
                 else:
-                    dimension_labels = ['angle' , 
-                                        'horizontal']
-            DataSet.__init__(self, array, deep_copy, dimension_labels, **kwargs)
+                    if vert > 1:
+                        shape = (num_of_angles, vert, horiz)
+                        dim_labels = ['angles' , 'vertical' , 
+                                      'horizontal']
+                    else:
+                        shape = (num_of_angles, horiz)
+                        dim_labels = ['angles' , 
+                                      'horizontal']
+                
+                array = numpy.zeros( shape , dtype=numpy.float32) 
+                super(SinogramData, self).__init__(array, deep_copy,
+                                 dim_labels, **kwargs)
+        else:
             
+            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 or \
+                         array.number_of_dimensions == 4):
+                    raise ValueError('Number of dimensions are not 2 or 3 or 4: {0}'\
+                                     .format(array.number_of_dimensions))
+                
+                #DataSet.__init__(self, array.as_array(), deep_copy,
+                #                 array.dimension_labels, **kwargs)
+                super(SinogramData, self).__init__(array.as_array(), deep_copy,
+                                 array.dimension_labels, **kwargs)
+            elif type(array) == numpy.ndarray:
+                if not ( array.ndim == 2 or array.ndim == 3 or array.ndim == 4 ):
+                    raise ValueError(
+                            'Number of dimensions are not 2 or 3 or 4 : {0}'\
+                            .format(array.ndim))
+                    
+                if dimension_labels is None:
+                    if array.ndim == 4:
+                        dimension_labels = ['channel' ,'angle' , 'vertical' , 
+                                      'horizontal']
+                    elif array.ndim == 3:
+                        dimension_labels = ['angle' , 'vertical' , 
+                                      'horizontal']
+                    else:
+                        dimension_labels = ['angle' , 
+                                      'horizontal']   
+                
+                #DataSet.__init__(self, array, deep_copy, dimension_labels, **kwargs)
+                super(SinogramData, self).__init__(array, deep_copy, 
+                     dimension_labels, **kwargs)
+                
             
 class DataSetProcessor(object):
     '''Defines a generic DataSet processor
@@ -772,5 +868,31 @@ if __name__ == '__main__':
     s = numpy.reshape(numpy.asarray(s), (3,4,4))
     sino = SinogramData( s )
     
-    
+    shape = (4,3,2)
+    a = [i for i in range(2*3*4)]
+    a = numpy.asarray(a)
+    a = numpy.reshape(a, shape)
+    print (numpy.shape(a))
+    ds = DataSet(a, True, ['X', 'Y','Z'])
+    # this means that I expect the X to be of length 2 ,
+    # y of length 3 and z of length 4
+    subset = ['Y' ,'Z']
+    b0 = ds.subset( subset )
+    print ("shape b 3,2? {0}".format(numpy.shape(b0.as_array())))
+    # expectation on b is that it is 
+    # 3x2 cut at z = 0
+    
+    subset = ['X' ,'Y']
+    b1 = ds.subset( subset , Z=1)
+    print ("shape b 2,3? {0}".format(numpy.shape(b1.as_array())))
+    
+    
+    # create VolumeData from geometry
+    vgeometry = geoms.VolumeGeometry(voxel_num_x=2, voxel_num_y=3, channels=2)
+    vol = VolumeData(geometry=vgeometry)
+    
+    sgeometry = geoms.SinogramGeometry(dimension=2, angles=numpy.linspace(0, 180, num=20), 
+                                       geom_type='parallel', pixel_num_v=3,
+                                       pixel_num_h=5 , channels=2)
+    sino = SinogramData(geometry=sgeometry)
     
\ No newline at end of file
diff --git a/Wrappers/Python/ccpi/processors.py b/Wrappers/Python/ccpi/processors.py
index 09e7229..c509a4e 100755
--- a/Wrappers/Python/ccpi/processors.py
+++ b/Wrappers/Python/ccpi/processors.py
@@ -32,19 +32,23 @@ class Normalizer(DataSetProcessor):
     Input: SinogramDataSet

     Parameter: 2D projection with flat field (or stack)

                2D projection with dark field (or stack)

-    Output: SinogramDataSet

+    Output: SinogramDataSetn

     '''

     

-    def __init__(self):

+    def __init__(self, flat_field = None, dark_field = None, tolerance = 1e-5):

         kwargs = {

-                  'flat_field'  :None, 

-                  'dark_field'  :None,

+                  'flat_field'  : None, 

+                  'dark_field'  : None,

                   # very small number. Used when there is a division by zero

-                  'tolerance'   : 1e-5

+                  'tolerance'   : tolerance

                   }

         

         #DataSetProcessor.__init__(self, **kwargs)

         super(Normalizer, self).__init__(**kwargs)

+        if not flat_field is None:

+            self.setFlatField(flat_field)

+        if not dark_field is None:

+            self.setDarkField(dark_field)

     

     def checkInput(self, dataset):

         if dataset.number_of_dimensions == 3:

diff --git a/Wrappers/Python/ccpi/reconstruction/astra_ops.py b/Wrappers/Python/ccpi/reconstruction/astra_ops.py
index 6af7af7..1346f22 100644
--- a/Wrappers/Python/ccpi/reconstruction/astra_ops.py
+++ b/Wrappers/Python/ccpi/reconstruction/astra_ops.py
@@ -19,7 +19,7 @@ from ccpi.reconstruction.ops import Operator
 import numpy
 from ccpi.framework import SinogramData, VolumeData
 from ccpi.reconstruction.ops import PowerMethodNonsquare
-from ccpi.astra.astra_processors import AstraForwardProjector, AstraBackProjector
+from ccpi.astra.astra_processors import *
 
 class AstraProjectorSimple(Operator):
     """ASTRA projector modified to use DataSet and geometry."""
@@ -66,6 +66,55 @@ class AstraProjectorSimple(Operator):
                   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..a887c0c 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=0, 
+                 voxel_num_y=0, 
+                 voxel_num_z=0, 
+                 voxel_size_x=1, 
+                 voxel_size_y=1, 
+                 voxel_size_z=1, 
+                 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
@@ -35,24 +37,31 @@ class VolumeGeometry:
         return self.center_y + 0.5*self.voxel_num_y*self.voxel_size_y
         
     def getMinZ(self):
-        return self.center_z - 0.5*self.voxel_num_z*self.voxel_size_z
+        if not voxel_num_z == 0:
+            return self.center_z - 0.5*self.voxel_num_z*self.voxel_size_z
+        else:
+            return 0
         
     def getMaxZ(self):
-        return self.center_z + 0.5*self.voxel_num_z*self.voxel_size_z
+        if not voxel_num_z == 0:
+            return self.center_z + 0.5*self.voxel_num_z*self.voxel_size_z
+        else:
+            return 0
         
     
 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=0, 
+                 pixel_size_h=1, 
+                 pixel_num_v=0, 
+                 pixel_size_v=1, 
+                 dist_source_center=None, 
+                 dist_center_detector=None, 
+                 channels=1
                  ):
         """
         General inputs for standard type projection geometries
@@ -91,6 +100,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/setup.py b/Wrappers/Python/setup.py
index bbd210f..3f83081 100644
--- a/Wrappers/Python/setup.py
+++ b/Wrappers/Python/setup.py
@@ -35,7 +35,7 @@ if  cil_version == '':
 setup(
     name="ccpi-common",
     version=cil_version,
-    packages=['ccpi' , 'ccpi.reconstruction'],
+    packages=['ccpi' , 'ccpi.reconstruction' , 'ccpi.astra'],
 
     # Project uses reStructuredText, so ensure that the docutils get
     # installed or upgraded on the target machine
diff --git a/Wrappers/Python/test/DemoRecIP.py b/Wrappers/Python/wip/DemoRecIP.py
index 483a202..483a202 100644
--- a/Wrappers/Python/test/DemoRecIP.py
+++ b/Wrappers/Python/wip/DemoRecIP.py
diff --git a/Wrappers/Python/test/regularizers.py b/Wrappers/Python/wip/regularizers.py
index 04ac3aa..04ac3aa 100644
--- a/Wrappers/Python/test/regularizers.py
+++ b/Wrappers/Python/wip/regularizers.py
diff --git a/Wrappers/Python/test/simple_demo.py b/Wrappers/Python/wip/simple_demo.py
index 0bbb687..7c5f43e 100644
--- a/Wrappers/Python/test/simple_demo.py
+++ b/Wrappers/Python/wip/simple_demo.py
@@ -1,34 +1,30 @@
+#import sys
+#sys.path.append("..")
 
-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 *
+from ccpi.framework import VolumeData
+from ccpi.reconstruction.algs import FISTA
+from ccpi.reconstruction.funcs import Norm2sq, Norm1
+from ccpi.reconstruction.astra_ops import AstraProjectorSimple
+from ccpi.reconstruction.geoms import VolumeGeometry, SinogramGeometry
 
 import numpy as np
 import matplotlib.pyplot as plt
 
-test_case = 2   # 1=parallel2D, 2=cone2D
+test_case = 1   # 1=parallel2D, 2=cone2D
 
 # Set up phantom
 N = 128
 
-x = np.zeros((N,N))
+vg = VolumeGeometry(voxel_num_x=N,voxel_num_y=N)
+Phantom = VolumeData(geometry=vg)
+
+x = Phantom.as_array()
 x[round(N/4):round(3*N/4),round(N/4):round(3*N/4)] = 1.0
 x[round(N/8):round(7*N/8),round(3*N/8):round(5*N/8)] = 2.0
 
 plt.imshow(x)
 plt.show()
 
-vg = VolumeGeometry(N,N,None, 1,1,None)
-
-Phantom = VolumeData(x,geometry=vg)
-
 # Set up measurement geometry
 angles_num = 20; # angles number
 
diff --git a/Wrappers/Python/wip/simple_mc_demo.py b/Wrappers/Python/wip/simple_mc_demo.py
new file mode 100755
index 0000000..0d976d7
--- /dev/null
+++ b/Wrappers/Python/wip/simple_mc_demo.py
@@ -0,0 +1,143 @@
+#import sys
+#sys.path.append("..")
+
+from ccpi.framework import VolumeData, SinogramData
+from ccpi.reconstruction.algs import FISTA
+from ccpi.reconstruction.funcs import Norm2sq, Norm1
+from ccpi.reconstruction.astra_ops import AstraProjectorMC
+from ccpi.reconstruction.geoms import VolumeGeometry, SinogramGeometry
+
+import numpy
+import matplotlib.pyplot as plt
+
+test_case = 1   # 1=parallel2D, 2=cone2D
+
+# Set up phantom
+N = 128
+M = 100
+numchannels = 3
+
+vg = VolumeGeometry(voxel_num_x=N,voxel_num_y=N,channels=numchannels)
+Phantom = VolumeData(geometry=vg)
+
+x = Phantom.as_array()
+x[0 , round(N/4):round(3*N/4) , round(N/4):round(3*N/4)  ] = 1.0
+x[0 , round(N/8):round(7*N/8) , round(3*N/8):round(5*N/8)] = 2.0
+
+x[1 , round(N/4):round(3*N/4) , round(N/4):round(3*N/4)  ] = 0.7
+x[1 , round(N/8):round(7*N/8) , round(3*N/8):round(5*N/8)] = 1.2
+
+x[2 , round(N/4):round(3*N/4) , round(N/4):round(3*N/4)  ] = 1.5
+x[2 , round(N/8):round(7*N/8) , round(3*N/8):round(5*N/8)] = 2.2
+
+#x = numpy.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[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 = numpy.linspace(0,numpy.pi,angles_num,endpoint=False)
+elif test_case==2:
+    angles = numpy.linspace(0,2*numpy.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, 'cpu')
+
+
+# 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()[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()[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(numpy.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()[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()[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