Ccpi forward backprojector (#77)

* added simple test for SinogramData

* initial version of processors

This contains a version of a simple normalization algorithm.

* added default pixel size = 1 (mm)

* Added FindCenter of rotation for parallel beam geometry

closes #5

* WIP for ccpi forward/backprojectors

* WIP

* Added Forward and Backprojectors from ccpi

added processor to pad acquisition data to center on center of rotation.
requires update of ccpi-reconstruction.

* Added simple_demo_ccpi.py

Currently there is a problem with the order of the axis after a projection/backprojection with the ccpi operator.
In the calculation of the PowerMethod there is a wrong order of the axis.
Possibly the creation of the data should be handled by the operator.

* fixed calculation of power method

* further work

* fixed the x_init for the FISTA example

* fixed ccpi FISTA demo

* Fix missing subset call for display

* Added create_image_data method to AstraOperatorSimple

* Stack direction 0, not 2, in TV

* Added finite different forward processor

* wip

* add io to setup.py

* wip

* comment out FiniteDifferentiator

* uniformed algorithm interface

Added some docstring.
Added parsing of opt.

* removed ccpi.io

* little changes

* fix setupCCPiGeometries

* added simple astra demo

* use minimal size for example

7 files changed, 1183 insertions, 91 deletions

diff --git a/Wrappers/Python/ccpi/astra/astra_ops.py b/Wrappers/Python/ccpi/astra/astra_ops.py
index 5f25eec..d8c2f4b 100644
--- a/Wrappers/Python/ccpi/astra/astra_ops.py
+++ b/Wrappers/Python/ccpi/astra/astra_ops.py
@@ -18,7 +18,7 @@
 from ccpi.reconstruction.ops import Operator
 import numpy
 import astra
-from ccpi.framework import AcquisitionData, ImageData
+from ccpi.framework import AcquisitionData, ImageData, DataContainer
 from ccpi.reconstruction.ops import PowerMethodNonsquare
 from ccpi.astra.astra_processors import AstraForwardProjector, AstraBackProjector, \
      AstraForwardProjectorMC, AstraBackProjectorMC
@@ -68,6 +68,10 @@ class AstraProjectorSimple(Operator):
                   self.sinogram_geometry.pixel_num_h), \
                  (self.volume_geometry.voxel_num_x, \
                   self.volume_geometry.voxel_num_y) )
+    
+    def create_image_data(self):
+        inputsize = self.size()[1]
+        return DataContainer(numpy.random.randn(inputsize[0],inputsize[1]))
 
 class AstraProjectorMC(Operator):
     """ASTRA Multichannel projector"""
diff --git a/Wrappers/Python/ccpi/framework.py b/Wrappers/Python/ccpi/framework.py
index 877becf..15dbe30 100644
--- a/Wrappers/Python/ccpi/framework.py
+++ b/Wrappers/Python/ccpi/framework.py
@@ -83,6 +83,20 @@ class ImageGeometry:
         else:
             return 0
         
+    def clone(self):
+        '''returns a copy of ImageGeometry'''
+        return ImageGeometry(
+                            self.voxel_num_x, 
+                            self.voxel_num_y, 
+                            self.voxel_num_z, 
+                            self.voxel_size_x, 
+                            self.voxel_size_y, 
+                            self.voxel_size_z, 
+                            self.center_x, 
+                            self.center_y, 
+                            self.center_z, 
+                            self.channels)
+        
     
 class AcquisitionGeometry:
     
@@ -137,6 +151,31 @@ class AcquisitionGeometry:
         self.pixel_size_v = pixel_size_v
         
         self.channels = channels
+        
+    def clone(self):
+        '''returns a copy of the AcquisitionGeometry'''
+        return AcquisitionGeometry(self.geom_type,
+                                   self.dimension, 
+                                   self.angles, 
+                                   self.pixel_num_h, 
+                                   self.pixel_size_h, 
+                                   self.pixel_num_v, 
+                                   self.pixel_size_v, 
+                                   self.dist_source_center, 
+                                   self.dist_center_detector, 
+                                   self.channels)
+        
+    def __str__ (self):
+        repres = ""
+        repres += "Number of dimensions: {0}\n".format(self.dimension)
+        repres += "angles: {0}\n".format(len(self.angles))
+        repres += "voxel_num : h{0},v{1}\n".format(self.pixel_num_h, self.pixel_num_v)
+        repres += "voxel size: h{0},v{1}\n".format(self.pixel_size_h, self.pixel_size_v)
+        repres += "geometry type: {0}\n".format(self.geom_type)
+        repres += "distance source-detector: {0}\n".format(self.dist_source_center)
+        repres += "distance center-detector: {0}\n".format(self.dist_source_center)
+        repres += "number of channels: {0}\n".format(self.channels)
+        return repres
 
 
             
@@ -178,6 +217,25 @@ class DataContainer(object):
             # assume it is parallel beam
             pass
         
+    def get_dimension_size(self, dimension_label):
+        if dimension_label in self.dimension_labels.values():
+            acq_size = -1
+            for k,v in self.dimension_labels.items():
+                if v == dimension_label:
+                    acq_size = self.shape[k]
+            return acq_size
+        else:
+            raise ValueError('Unknown dimension {0}. Should be one of'.format(dimension_label,
+                             self.dimension_labels))
+    def get_dimension_axis(self, dimension_label):
+        if dimension_label in self.dimension_labels.values():
+            for k,v in self.dimension_labels.items():
+                if v == dimension_label:
+                    return k
+        else:
+            raise ValueError('Unknown dimension {0}. Should be one of'.format(dimension_label,
+                             self.dimension_labels.values()))
+                        
 
     def as_array(self, dimensions=None):
         '''Returns the DataContainer as Numpy Array
@@ -194,7 +252,15 @@ class DataContainer(object):
         '''Creates a DataContainer containing a subset of self according to the 
         labels in dimensions'''
         if dimensions is None:
-            return self.array.copy()
+            if kw == {}:
+                return self.array.copy()
+            else:
+                reduced_dims = [v for k,v in self.dimension_labels.items()]
+                for dim_l, dim_v in kw.items():
+                    for k,v in self.dimension_labels.items():
+                        if v == dim_l:
+                            reduced_dims.pop(k)
+                return self.subset(dimensions=reduced_dims, **kw)
         else:
             # check that all the requested dimensions are in the array
             # this is done by checking the dimension_labels
@@ -258,14 +324,31 @@ class DataContainer(object):
                 
                 return type(self)(cleaned , True, dimensions)
     
-    def fill(self, array):
+    def fill(self, array, **dimension):
         '''fills the internal numpy array with the one provided'''
-        if numpy.shape(array) != numpy.shape(self.array):
-            raise ValueError('Cannot fill with the provided array.' + \
-                             'Expecting {0} got {1}'.format(
-                                     numpy.shape(self.array),
-                                     numpy.shape(array)))
-        self.array = array[:]
+        if dimension == {}:
+            if numpy.shape(array) != numpy.shape(self.array):
+                raise ValueError('Cannot fill with the provided array.' + \
+                                 'Expecting {0} got {1}'.format(
+                                         numpy.shape(self.array),
+                                         numpy.shape(array)))
+            self.array = array[:]
+        else:
+            
+            command = 'self.array['
+            i = 0
+            for k,v in self.dimension_labels.items():
+                for dim_label, dim_value in dimension.items():    
+                    if dim_label == v:
+                        command = command + str(dim_value)
+                    else:
+                        command = command + ":"
+                if i < self.number_of_dimensions -1:
+                    command = command + ','
+                i += 1
+            command = command + "] = array[:]" 
+            exec(command)
+            
         
     def check_dimensions(self, other):
         return self.shape == other.shape
@@ -301,7 +384,7 @@ class DataContainer(object):
                                dimension_labels=self.dimension_labels,
                                geometry=self.geometry)
             else:
-                raise ValueError('Wrong shape: {0} and {1}'.format(self.shape, 
+                raise ValueError('__sub__ Wrong shape: {0} and {1}'.format(self.shape, 
                                  other.shape))
         elif isinstance(other, (int, float, complex)):
             return type(self)(self.as_array() - other, 
@@ -325,7 +408,7 @@ class DataContainer(object):
                                dimension_labels=self.dimension_labels,
                                geometry=self.geometry)
             else:
-                raise ValueError('Wrong shape: {0} and {1}'.format(self.shape, 
+                raise ValueError('__div__ Wrong shape: {0} and {1}'.format(self.shape, 
                                  other.shape))
         elif isinstance(other, (int, float, complex)):
             return type(self)(self.as_array() / other, 
@@ -346,7 +429,7 @@ class DataContainer(object):
                                dimension_labels=self.dimension_labels,
                                geometry=self.geometry)
             else:
-                raise ValueError('Wrong shape: {0} and {1}'.format(self.shape, 
+                raise ValueError('__pow__ Wrong shape: {0} and {1}'.format(self.shape, 
                                  other.shape))
         elif isinstance(other, (int, float, complex)):
             return type(self)(self.as_array() ** other, 
@@ -354,7 +437,7 @@ class DataContainer(object):
                                dimension_labels=self.dimension_labels,
                                geometry=self.geometry)
         else:
-            raise TypeError('Cannot {0} DataContainer with {1}'.format("power" ,
+            raise TypeError('pow: Cannot {0} DataContainer with {1}'.format("power" ,
                             type(other)))
     # __pow__
     
@@ -367,7 +450,7 @@ class DataContainer(object):
                                dimension_labels=self.dimension_labels,
                                geometry=self.geometry)
             else:
-                raise ValueError('Wrong shape: {0} and {1}'.format(self.shape, 
+                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, 
@@ -464,12 +547,13 @@ class DataContainer(object):
         return self / other
     # __idiv__
     
-    def __str__ (self):
+    def __str__ (self, representation=False):
         repres = ""
         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: \n{0}\n".format(self.array)
+        if representation:
+            repres += "Representation: \n{0}\n".format(self.array)
         return repres
     
     def clone(self):
@@ -479,6 +563,40 @@ class DataContainer(object):
                           dimension_labels=self.dimension_labels,
                           deep_copy=True,
                           geometry=self.geometry )
+    
+    def get_data_axes_order(self,new_order=None):
+        '''returns the axes label of self as a list
+        
+        if new_order is None returns the labels of the axes as a sorted-by-key list
+        if new_order is a list of length number_of_dimensions, returns a list
+        with the indices of the axes in new_order with respect to those in 
+        self.dimension_labels: i.e.
+          self.dimension_labels = {0:'horizontal',1:'vertical'}
+          new_order = ['vertical','horizontal']
+          returns [1,0]
+        '''
+        if new_order is None:
+            
+            axes_order = [i for i in range(len(self.shape))]
+            for k,v in self.dimension_labels.items():
+                axes_order[k] = v
+            return axes_order
+        else:
+            if len(new_order) == self.number_of_dimensions:
+                axes_order = [i for i in range(self.number_of_dimensions)]
+                
+                for i in range(len(self.shape)):
+                    found = False
+                    for k,v in self.dimension_labels.items():
+                        if new_order[i] == v:
+                            axes_order[i] = k
+                            found = True
+                    if not found:
+                        raise ValueError('Axis label {0} not found.'.format(new_order[i]))
+                return axes_order
+            else:
+                raise ValueError('Expecting {0} axes, got {2}'\
+                                 .format(len(self.shape),len(new_order)))
         
                 
                     
@@ -501,29 +619,45 @@ class ImageData(DataContainer):
                 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 channels > 1:
-                    if vert > 1:
-                        shape = (channels, vert, horiz_y, horiz_x)
-                        dim_labels = ['channel' ,'vertical' , 'horizontal_y' , 
-                                      'horizontal_x']
+                if dimension_labels is None:
+                    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:
-                        shape = (channels , horiz_y, horiz_x)
-                        dim_labels = ['channel' , 'horizontal_y' , 
-                                      'horizontal_x']
+                        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']
+                    dimension_labels = dim_labels
                 else:
-                    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']
-                        
+                    shape = []
+                    for dim in dimension_labels:
+                        if dim == 'channel':
+                            shape.append(channels)
+                        elif dim == 'horizontal_y':
+                            shape.append(horiz_y)
+                        elif dim == 'vertical':
+                            shape.append(vert)
+                        elif dim == 'horizontal_x':
+                            shape.append(horiz_x)
+                    if len(shape) != len(dimension_labels):
+                        raise ValueError('Missing {0} axes'.format(
+                                len(dimension_labels) - len(shape)))
+                    shape = tuple(shape)
+                    
                 array = numpy.zeros( shape , dtype=numpy.float32) 
                 super(ImageData, self).__init__(array, deep_copy,
-                                 dim_labels, **kwargs)
+                                 dimension_labels, **kwargs)
                 
             else:
                 raise ValueError('Please pass either a DataContainer, ' +\
@@ -571,6 +705,12 @@ class ImageData(DataContainer):
                     if key == 'spacing' :
                         self.spacing = value
                         
+        def subset(self, dimensions=None, **kw):
+            out = super(ImageData, self).subset(dimensions, **kw)
+            #out.geometry = self.recalculate_geometry(dimensions , **kw)
+            out.geometry = self.geometry
+            return out
+                        
 
 class AcquisitionData(DataContainer):
     '''DataContainer for holding 2D or 3D sinogram'''
@@ -582,37 +722,54 @@ class AcquisitionData(DataContainer):
         self.geometry = None
         if array is None:
             if 'geometry' in kwargs.keys():
-                geometry  = kwargs['geometry']
+                geometry      = kwargs['geometry']
                 self.geometry = geometry
-                channels  = geometry.channels
-                horiz   = geometry.pixel_num_h
-                vert   = geometry.pixel_num_v
-                angles = geometry.angles
+                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']
+                if dimension_labels is None:
+                    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:
-                        shape = (channels , num_of_angles, horiz)
-                        dim_labels = ['channel' , 'angle' , 
-                                      'horizontal']
+                        if vert > 1:
+                            shape = (num_of_angles, vert, horiz)
+                            dim_labels = ['angle' , 'vertical' , 
+                                          'horizontal']
+                        else:
+                            shape = (num_of_angles, horiz)
+                            dim_labels = ['angle' , 
+                                          'horizontal']
+                    
+                    dimension_labels = dim_labels
                 else:
-                    if vert > 1:
-                        shape = (num_of_angles, vert, horiz)
-                        dim_labels = ['angles' , 'vertical' , 
-                                      'horizontal']
-                    else:
-                        shape = (num_of_angles, horiz)
-                        dim_labels = ['angles' , 
-                                      'horizontal']
-                
+                    shape = []
+                    for dim in dimension_labels:
+                        if dim == 'channel':
+                            shape.append(channels)
+                        elif dim == 'angle':
+                            shape.append(num_of_angles)
+                        elif dim == 'vertical':
+                            shape.append(vert)
+                        elif dim == 'horizontal':
+                            shape.append(horiz)
+                    if len(shape) != len(dimension_labels):
+                        raise ValueError('Missing {0} axes'.format(
+                                len(dimension_labels) - len(shape)))
+                    shape = tuple(shape)
+                    
                 array = numpy.zeros( shape , dtype=numpy.float32) 
                 super(AcquisitionData, self).__init__(array, deep_copy,
-                                 dim_labels, **kwargs)
+                                 dimension_labels, **kwargs)
         else:
             
             if type(array) == DataContainer:
@@ -732,6 +889,9 @@ class DataSetProcessor(object):
         
     def process(self):
         raise NotImplementedError('process must be implemented')
+        
+    
+    
 
 class DataSetProcessor23D(DataSetProcessor):
     '''Regularizers DataSetProcessor
@@ -945,7 +1105,7 @@ if __name__ == '__main__':
 
     # create VolumeData from geometry
     vgeometry = ImageGeometry(voxel_num_x=2, voxel_num_y=3, channels=2)
-    vol = VolumeData(geometry=vgeometry)
+    vol = ImageData(geometry=vgeometry)
     
     sgeometry = AcquisitionGeometry(dimension=2, angles=numpy.linspace(0, 180, num=20), 
                                        geom_type='parallel', pixel_num_v=3,
diff --git a/Wrappers/Python/ccpi/processors.py b/Wrappers/Python/ccpi/processors.py
index 3010053..d98ef12 100755
--- a/Wrappers/Python/ccpi/processors.py
+++ b/Wrappers/Python/ccpi/processors.py
@@ -17,12 +17,16 @@
 #   See the License for the specific language governing permissions and

 #   limitations under the License

 

-from ccpi.framework import DataSetProcessor, DataSet, AcquisitionData,\

- AcquisitionGeometry

+from ccpi.framework import DataSetProcessor, DataContainer, AcquisitionData,\

+ AcquisitionGeometry, ImageGeometry, ImageData

+from ccpi.reconstruction.parallelbeam import alg as pbalg

 import numpy

 import h5py

 from scipy import ndimage

 

+import matplotlib.pyplot as plt

+

+

 class Normalizer(DataSetProcessor):

     '''Normalization based on flat and dark

     

@@ -130,7 +134,8 @@ class CenterOfRotationFinder(DataSetProcessor):
             if dataset.geometry.geom_type == 'parallel':

                 return True

             else:

-                raise ValueError('This algorithm is suitable only for parallel beam geometry')

+                raise ValueError('{0} is suitable only for parallel beam geometry'\

+                                 .format(self.__class__.__name__))

         else:

             raise ValueError("Expected input dimensions is 3, got {0}"\

                              .format(dataset.number_of_dimensions))

@@ -385,6 +390,275 @@ class CenterOfRotationFinder(DataSetProcessor):
         

         return cor

 

+class CCPiForwardProjector(DataSetProcessor):

+    '''Normalization based on flat and dark

+    

+    This processor read in a AcquisitionData and normalises it based on 

+    the instrument reading with and without incident photons or neutrons.

+    

+    Input: AcquisitionData

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

+               2D projection with dark field (or stack)

+    Output: AcquisitionDataSetn

+    '''

+    

+    def __init__(self,

+                 image_geometry       = None, 

+                 acquisition_geometry = None,

+                 output_axes_order    = None):

+        if output_axes_order is None:

+            # default ccpi projector image storing order

+            output_axes_order = ['angle','vertical','horizontal']

+        

+        kwargs = {

+                  'image_geometry'       : image_geometry, 

+                  'acquisition_geometry' : acquisition_geometry,

+                  'output_axes_order'    : output_axes_order,

+                  'default_image_axes_order' : ['horizontal_x','horizontal_y','vertical'],

+                  'default_acquisition_axes_order' : ['angle','vertical','horizontal'] 

+                  }

+        

+        super(CCPiForwardProjector, self).__init__(**kwargs)

+        

+    def check_input(self, dataset):

+        if dataset.number_of_dimensions == 3 or dataset.number_of_dimensions == 2:

+            # sort in the order that this projector needs it

+            return True

+        else:

+            raise ValueError("Expected input dimensions is 2 or 3, got {0}"\

+                             .format(dataset.number_of_dimensions))

+

+    def process(self):

+        

+        volume = self.get_input()

+        volume_axes = volume.get_data_axes_order(new_order=self.default_image_axes_order)

+        if not volume_axes == [0,1,2]:

+            volume.array = numpy.transpose(volume.array, volume_axes)

+        pixel_per_voxel = 1 # should be estimated from image_geometry and 

+                            # acquisition_geometry

+        if self.acquisition_geometry.geom_type == 'parallel':

+            #int msize = ndarray_volume.shape(0) > ndarray_volume.shape(1) ? ndarray_volume.shape(0) : ndarray_volume.shape(1);

+        	  #int detector_width = msize;

+            # detector_width is the max between the shape[0] and shape[1]

+            

+            

+            #double rotation_center = (double)detector_width/2.;

+        	  #int detector_height = ndarray_volume.shape(2);

+        	  

+            #int number_of_projections = ndarray_angles.shape(0);

+        	

+            ##numpy_3d pixels(reinterpret_cast<float*>(ndarray_volume.get_data()),

+		     #boost::extents[number_of_projections][detector_height][detector_width]);

+

+            pixels = pbalg.pb_forward_project(volume.as_array(), 

+                                                  self.acquisition_geometry.angles, 

+                                                  pixel_per_voxel)

+            out = AcquisitionData(geometry=self.acquisition_geometry, 

+                                  label_dimensions=self.default_acquisition_axes_order)

+            out.fill(pixels)

+            out_axes = out.get_data_axes_order(new_order=self.output_axes_order)

+            if not out_axes == [0,1,2]:

+                out.array = numpy.transpose(out.array, out_axes)

+            return out

+        else:

+            raise ValueError('Cannot process cone beam')

+

+class CCPiBackwardProjector(DataSetProcessor):

+    '''Backward projector

+    

+    This processor reads in a AcquisitionData and performs a backward projection, 

+    i.e. project to reconstruction space.

+    Notice that it assumes that the center of rotation is in the middle

+    of the horizontal axis: in case when that's not the case it can be chained 

+    with the AcquisitionDataPadder.

+    

+    Input: AcquisitionData

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

+               2D projection with dark field (or stack)

+    Output: AcquisitionDataSetn

+    '''

+    

+    def __init__(self, 

+                 image_geometry = None, 

+                 acquisition_geometry = None,

+                 output_axes_order=None):

+        if output_axes_order is None:

+            # default ccpi projector image storing order

+            output_axes_order = ['horizontal_x','horizontal_y','vertical']

+        kwargs = {

+                  'image_geometry'       : image_geometry, 

+                  'acquisition_geometry' : acquisition_geometry,

+                  'output_axes_order'    : output_axes_order,

+                  'default_image_axes_order' : ['horizontal_x','horizontal_y','vertical'],

+                  'default_acquisition_axes_order' : ['angle','vertical','horizontal'] 

+                  }

+        

+        super(CCPiBackwardProjector, self).__init__(**kwargs)

+        

+    def check_input(self, dataset):

+        if dataset.number_of_dimensions == 3 or dataset.number_of_dimensions == 2:

+            #number_of_projections][detector_height][detector_width

+            

+            return True

+        else:

+            raise ValueError("Expected input dimensions is 2 or 3, got {0}"\

+                             .format(dataset.number_of_dimensions))

+

+    def process(self):

+        projections = self.get_input()

+        projections_axes = projections.get_data_axes_order(new_order=self.default_acquisition_axes_order)

+        if not projections_axes == [0,1,2]:

+            projections.array = numpy.transpose(projections.array, projections_axes)

+        

+        pixel_per_voxel = 1 # should be estimated from image_geometry and acquisition_geometry

+        image_geometry = ImageGeometry(voxel_num_x = self.acquisition_geometry.pixel_num_h,

+                                       voxel_num_y = self.acquisition_geometry.pixel_num_h,

+                                       voxel_num_z = self.acquisition_geometry.pixel_num_v)

+        # input centered/padded acquisitiondata

+        center_of_rotation = projections.get_dimension_size('horizontal') / 2

+        #print (center_of_rotation)

+        if self.acquisition_geometry.geom_type == 'parallel':

+            back = pbalg.pb_backward_project(

+                         projections.as_array(), 

+                         self.acquisition_geometry.angles, 

+                         center_of_rotation, 

+                         pixel_per_voxel

+                         )

+            out = ImageData(geometry=self.image_geometry, 

+                            dimension_labels=self.default_image_axes_order)

+            

+            out_axes = out.get_data_axes_order(new_order=self.output_axes_order)

+            if not out_axes == [0,1,2]:

+                back = numpy.transpose(back, out_axes)

+            out.fill(back)

+            

+            return out

+            

+        else:

+            raise ValueError('Cannot process cone beam')

+            

+class AcquisitionDataPadder(DataSetProcessor):

+    '''Normalization based on flat and dark

+    

+    This processor read in a AcquisitionData and normalises it based on 

+    the instrument reading with and without incident photons or neutrons.

+    

+    Input: AcquisitionData

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

+               2D projection with dark field (or stack)

+    Output: AcquisitionDataSetn

+    '''

+    

+    def __init__(self, 

+                 center_of_rotation   = None,

+                 acquisition_geometry = None,

+                 pad_value            = 1e-5):

+        kwargs = {

+                  'acquisition_geometry' : acquisition_geometry, 

+                  'center_of_rotation'   : center_of_rotation,

+                  'pad_value'            : pad_value

+                  }

+        

+        super(AcquisitionDataPadder, self).__init__(**kwargs)

+        

+    def check_input(self, dataset):

+        if self.acquisition_geometry is None:

+            self.acquisition_geometry = dataset.geometry

+        if dataset.number_of_dimensions == 3:

+               return True

+        else:

+            raise ValueError("Expected input dimensions is 2 or 3, got {0}"\

+                             .format(dataset.number_of_dimensions))

+

+    def process(self):

+        projections = self.get_input()

+        w = projections.get_dimension_size('horizontal')

+        delta = w - 2 * cor

+               

+        padded_width = int (

+                numpy.ceil(abs(delta)) + w

+                )

+        delta_pix = padded_width - w

+        

+        voxel_per_pixel = 1

+        geom = pbalg.pb_setup_geometry_from_acquisition(projections.as_array(),

+                                            self.acquisition_geometry.angles,

+                                            self.center_of_rotation,

+                                            voxel_per_pixel )

+        

+        padded_geometry = self.acquisition_geometry.clone()

+        

+        padded_geometry.pixel_num_h = geom['n_h']

+        padded_geometry.pixel_num_v = geom['n_v']

+        

+        delta_pix_h = padded_geometry.pixel_num_h - self.acquisition_geometry.pixel_num_h

+        delta_pix_v = padded_geometry.pixel_num_v - self.acquisition_geometry.pixel_num_v

+        

+        if delta_pix_h == 0:

+            delta_pix_h = delta_pix

+            padded_geometry.pixel_num_h = padded_width

+        #initialize a new AcquisitionData with values close to 0

+        out = AcquisitionData(geometry=padded_geometry)

+        out = out + self.pad_value

+        

+        

+        #pad in the horizontal-vertical plane -> slice on angles

+        if delta > 0:

+            #pad left of middle

+            command = "out.array["

+            for i in range(out.number_of_dimensions):

+                if out.dimension_labels[i] == 'horizontal':

+                    value = '{0}:{1}'.format(delta_pix_h, delta_pix_h+w)

+                    command = command + str(value)

+                else:

+                    if out.dimension_labels[i] == 'vertical' :

+                        value = '{0}:'.format(delta_pix_v)

+                        command = command + str(value)

+                    else:

+                        command = command + ":"

+                if i < out.number_of_dimensions -1:

+                    command = command + ','

+            command = command + '] = projections.array'

+            #print (command)    

+        else:

+            #pad right of middle

+            command = "out.array["

+            for i in range(out.number_of_dimensions):

+                if out.dimension_labels[i] == 'horizontal':

+                    value = '{0}:{1}'.format(0, w)

+                    command = command + str(value)

+                else:

+                    if out.dimension_labels[i] == 'vertical' :

+                        value = '{0}:'.format(delta_pix_v)

+                        command = command + str(value)

+                    else:

+                        command = command + ":"

+                if i < out.number_of_dimensions -1:

+                    command = command + ','

+            command = command + '] = projections.array'

+            #print (command)    

+            #cleaned = eval(command)

+        exec(command)

+        return out

+        

+#class FiniteDifferentiator(DataSetProcessor):

+#    def __init__(self):

+#        kwargs = {

+#                  

+#                  }

+#        

+#        super(FiniteDifferentiator, self).__init__(**kwargs)

+#        

+#    def check_input(self, dataset):

+#        return True

+#    

+#    def process(self):

+#        axis = 0

+#        d1 = numpy.diff(x,n=1,axis=axis)

+#        d1 = numpy.resize(d1, numpy.shape(x))

+

+        

+

 def loadNexus(filename):

     '''Load a dataset stored in a NeXuS file (HDF5)'''

     ###########################################################################

@@ -437,35 +711,94 @@ if __name__ == '__main__':
                                  pixel_num_h=numpy.shape(proj)[2],

                                  pixel_num_v=numpy.shape(proj)[1],

                                  )    

-    sino = AcquisitionData( proj , geometry=parallelbeam)

+    

+    dim_labels = ['angles' , 'vertical' , 'horizontal']

+    sino = AcquisitionData( proj , geometry=parallelbeam, 

+                            dimension_labels=dim_labels)

     

     normalizer = Normalizer()

     normalizer.set_input(sino)

     normalizer.set_flat_field(flat)

     normalizer.set_dark_field(dark)

     norm = normalizer.get_output()

-    print ("Processor min {0} max {1}".format(norm.as_array().min(), norm.as_array().max()))

+    #print ("Processor min {0} max {1}".format(norm.as_array().min(), norm.as_array().max()))

     

-    norm1 = numpy.asarray(

-            [Normalizer.normalize_projection( p, flat, dark, 1e-5 ) 

-            for p in proj]

-            )

+    #norm1 = numpy.asarray(

+    #        [Normalizer.normalize_projection( p, flat, dark, 1e-5 ) 

+    #       for p in proj]

+    #        )

     

-    print ("Numpy min {0} max {1}".format(norm1.min(), norm1.max()))

+    #print ("Numpy min {0} max {1}".format(norm1.min(), norm1.max()))

     

     cor_finder = CenterOfRotationFinder()

     cor_finder.set_input(sino)

     cor = cor_finder.get_output()

     print ("center of rotation {0} == 86.25?".format(cor))

     

+        

+    padder = AcquisitionDataPadder(center_of_rotation=cor, 

+                                   pad_value=0.75,

+                                   acquisition_geometry=normalizer.get_output().geometry)

+    #padder.set_input(normalizer.get_output())

+    padder.set_input_processor(normalizer)

+    

+    #print ("padder ", padder.get_output())

+    

+    volume_geometry = ImageGeometry()

+    

+    back = CCPiBackwardProjector(acquisition_geometry=parallelbeam, 

+                                 image_geometry=volume_geometry)

+    back.set_input_processor(padder)

+    #back.set_input(padder.get_output())

+    #print (back.image_geometry)

+    forw = CCPiForwardProjector(acquisition_geometry=parallelbeam, 

+                                 image_geometry=volume_geometry)

+    forw.set_input_processor(back)

+    

+

+    out  = padder.get_output()

+    fig = plt.figure()

+    # projections row

+    a = fig.add_subplot(1,4,1)    

+    a.imshow(norm.array[80])

+    a.set_title('orig')

+    a = fig.add_subplot(1,4,2)    

+    a.imshow(out.array[80])

+    a.set_title('padded')

+    

+    a = fig.add_subplot(1,4,3)    

+    a.imshow(back.get_output().as_array()[15])

+    a.set_title('back')

+    

+    a = fig.add_subplot(1,4,4)    

+    a.imshow(forw.get_output().as_array()[80])

+    a.set_title('forw')

+    

+    

+    plt.show()

+

+#    back = CCPiBackwardProjector(acquisition_geometry=parallelbeam, 

+#                                 image_geometry=volume_geometry)

+    

+#    int msize = ndarray_volume.shape(0) > ndarray_volume.shape(1) ? ndarray_volume.shape(0) : ndarray_volume.shape(1);#

+#	int detector_width = msize;

+

+#	double rotation_center = (double)detector_width/2.;

+#	int detector_height = ndarray_volume.shape(2);

+#	int number_of_projections = ndarray_angles.shape(0);

+#	std::cout << "pb_forward_project rotation_center " << rotation_center << std::endl;

+#		boost::extents[number_of_projections][detector_height][detector_width]);

+

+    

     conebeam = AcquisitionGeometry('cone', '3D' , 

                                  angles=angles, 

                                  pixel_num_h=numpy.shape(proj)[2],

                                  pixel_num_v=numpy.shape(proj)[1],

                                  )    

-    sino = AcquisitionData( proj , geometry=conebeam)

+    

     try:

-        cor_finder.set_input(sino)

+        sino2 = AcquisitionData( proj , geometry=conebeam)

+        cor_finder.set_input(sino2)

         cor = cor_finder.get_output()

     except ValueError as err:

         print (err)
\ No newline at end of file
diff --git a/Wrappers/Python/ccpi/reconstruction/algs.py b/Wrappers/Python/ccpi/reconstruction/algs.py
index ec52fee..46865a8 100644
--- a/Wrappers/Python/ccpi/reconstruction/algs.py
+++ b/Wrappers/Python/ccpi/reconstruction/algs.py
@@ -23,14 +23,36 @@ import time
 from ccpi.reconstruction.funcs import BaseFunction
 
 def FISTA(x_init, f=None, g=None, opt=None):
-
+    '''Fast Iterative Shrinkage-Thresholding Algorithm
+    
+    Beck, A. and Teboulle, M., 2009. A fast iterative shrinkage-thresholding 
+    algorithm for linear inverse problems. 
+    SIAM journal on imaging sciences,2(1), pp.183-202.
+    
+    Parameters:
+      x_init: initial guess
+      f: data fidelity
+      g: regularizer
+      h:
+      opt: additional algorithm 
+    '''
     # default inputs
     if f   is None: f = BaseFunction()
     if g   is None: g = BaseFunction()
-    if opt is None: opt = {'tol': 1e-4, 'iter': 1000}
-
+    
     # algorithmic parameters
-    tol      = opt['tol']
+    if opt is None: 
+        opt = {'tol': 1e-4, 'iter': 1000}
+    else:
+        try:
+            max_iter = opt['iter']
+        except KeyError as ke:
+            opt[ke] = 1000
+        try:
+            opt['tol'] = 1000
+        except KeyError as ke:
+            opt[ke] = 1e-4
+    tol = opt['tol']
     max_iter = opt['iter']
     
     # initialization
@@ -76,15 +98,33 @@ def FISTA(x_init, f=None, g=None, opt=None):
     return x, it, timing, criter
 
 def FBPD(x_init, f=None, g=None, h=None, opt=None):
-
+    '''FBPD Algorithm
+    
+    Parameters:
+      x_init: initial guess
+      f: constraint
+      g: data fidelity
+      h: regularizer
+      opt: additional algorithm 
+    '''
     # default inputs
     if f   is None: f = BaseFunction()
     if g   is None: g = BaseFunction()
     if h   is None: h = BaseFunction()
-    if opt is None: opt = {'tol': 1e-4, 'iter': 1000}
-
+    
     # algorithmic parameters
-    tol      = opt['tol']
+    if opt is None: 
+        opt = {'tol': 1e-4, 'iter': 1000}
+    else:
+        try:
+            max_iter = opt['iter']
+        except KeyError as ke:
+            opt[ke] = 1000
+        try:
+            opt['tol'] = 1000
+        except KeyError as ke:
+            opt[ke] = 1e-4
+    tol = opt['tol']
     max_iter = opt['iter']
     
     # step-sizes
@@ -127,13 +167,34 @@ def FBPD(x_init, f=None, g=None, h=None, opt=None):
     
     return x, it, timing, criter
 
-def CGLS(A,b,max_iter,x_init):
-    '''Conjugate Gradient Least Squares algorithm'''
+def CGLS(x_init, operator , data , opt=None):
+    '''Conjugate Gradient Least Squares algorithm
+    
+    Parameters:
+      x_init: initial guess
+      operator: operator for forward/backward projections
+      data: data to operate on
+      opt: additional algorithm 
+    '''
+    
+    if opt is None: 
+        opt = {'tol': 1e-4, 'iter': 1000}
+    else:
+        try:
+            max_iter = opt['iter']
+        except KeyError as ke:
+            opt[ke] = 1000
+        try:
+            opt['tol'] = 1000
+        except KeyError as ke:
+            opt[ke] = 1e-4
+    tol = opt['tol']
+    max_iter = opt['iter']
     
-    r = b.clone()
+    r = data.clone()
     x = x_init.clone()
     
-    d = A.adjoint(r)
+    d = operator.adjoint(r)
     
     normr2 = (d**2).sum()
     
@@ -145,11 +206,11 @@ def CGLS(A,b,max_iter,x_init):
     
         t = time.time()
         
-        Ad = A.direct(d)
+        Ad = operator.direct(d)
         alpha = normr2/( (Ad**2).sum() )
         x  = x + alpha*d
         r  = r - alpha*Ad
-        s  = A.adjoint(r)
+        s  = operator.adjoint(r)
         
         normr2_new = (s**2).sum()
         beta = normr2_new/normr2
diff --git a/Wrappers/Python/ccpi/reconstruction/ops.py b/Wrappers/Python/ccpi/reconstruction/ops.py
index d6f31eb..fa30b0c 100644
--- a/Wrappers/Python/ccpi/reconstruction/ops.py
+++ b/Wrappers/Python/ccpi/reconstruction/ops.py
@@ -19,7 +19,8 @@
 
 import numpy
 from scipy.sparse.linalg import svds
-from ccpi.framework import DataContainer
+from ccpi.framework import DataContainer, ImageGeometry , ImageData
+from ccpi.processors import CCPiBackwardProjector, CCPiForwardProjector
 
 # Maybe operators need to know what types they take as inputs/outputs
 # to not just use generic DataContainer
@@ -110,7 +111,13 @@ class FiniteDiff2D(Operator):
     
     def adjoint(self,x):
         '''Backward differences, Newumann BC.'''
-        Nrows, Ncols, Nchannels = x.as_array().shape
+        #Nrows, Ncols, Nchannels = x.as_array().shape
+        print (x)
+        Nrows = x.get_dimension_size('horizontal_x')
+        Ncols = x.get_dimension_size('horizontal_x')
+        Nchannels = 1
+        if len(x.shape) == 4:
+            Nchannels = x.get_dimension_size('channel')
         zer = numpy.zeros((Nrows,1))
         xxx = x.as_array()[:,:-1,0]
         h = numpy.concatenate((zer,xxx), 1) - numpy.concatenate((xxx,zer), 1)
@@ -129,13 +136,27 @@ class FiniteDiff2D(Operator):
 
 def PowerMethodNonsquare(op,numiters):
     # Initialise random
-    inputsize = op.size()[1]
-    x0 = DataContainer(numpy.random.randn(inputsize[0],inputsize[1]))
+    # Jakob's
+    #inputsize = op.size()[1]
+    #x0 = ImageContainer(numpy.random.randn(*inputsize)
+    # Edo's
+    #vg = ImageGeometry(voxel_num_x=inputsize[0],
+    #                   voxel_num_y=inputsize[1], 
+    #                   voxel_num_z=inputsize[2])
+    #
+    #x0 = ImageData(geometry = vg, dimension_labels=['vertical','horizontal_y','horizontal_x'])
+    #print (x0)
+    #x0.fill(numpy.random.randn(*x0.shape))
+    
+    x0 = op.create_image_data()
+    
     s = numpy.zeros(numiters)
     # Loop
     for it in numpy.arange(numiters):
         x1 = op.adjoint(op.direct(x0))
         x1norm = numpy.sqrt((x1**2).sum())
+        #print ("x0 **********" ,x0)
+        #print ("x1 **********" ,x1)
         s[it] = (x1*x0).sum() / (x0*x0).sum()
         x0 = (1.0/x1norm)*x1
     return numpy.sqrt(s[-1]), numpy.sqrt(s), x0
@@ -150,4 +171,58 @@ def PowerMethodNonsquare(op,numiters):
 #        x1norm = np.sqrt(np.sum(np.dot(x1,x1)))
 #        s[it] = np.dot(x1,x0) / np.dot(x1,x0)
 #        x0 = (1.0/x1norm)*x1
-#    return s, x0
\ No newline at end of file
+#    return s, x0
+    
+class CCPiProjectorSimple(Operator):
+    """ASTRA projector modified to use DataSet and geometry."""
+    def __init__(self, geomv, geomp):
+        super(CCPiProjectorSimple, self).__init__()
+        
+        # Store volume and sinogram geometries.
+        self.acquisition_geometry = geomp
+        self.volume_geometry = geomv
+        
+        self.fp = CCPiForwardProjector(image_geometry=geomv,
+                                       acquisition_geometry=geomp,
+                                       output_axes_order=['angle','vertical','horizontal'])
+        
+        self.bp = CCPiBackwardProjector(image_geometry=geomv,
+                                    acquisition_geometry=geomp,
+                                    output_axes_order=['horizontal_x','horizontal_y','vertical'])
+                
+        # Initialise empty for singular value.
+        self.s1 = None
+    
+    def direct(self, image_data):
+        self.fp.set_input(image_data)
+        out = self.fp.get_output()
+        return out
+    
+    def adjoint(self, acquisition_data):
+        self.bp.set_input(acquisition_data)
+        out = self.bp.get_output()
+        return out
+    
+    #def delete(self):
+    #    astra.data2d.delete(self.proj_id)
+    
+    def get_max_sing_val(self):
+        a = PowerMethodNonsquare(self,10)
+        self.s1 = a[0] 
+        return self.s1
+    
+    def size(self):
+        # Only implemented for 3D
+        return ( (self.acquisition_geometry.angles.size, \
+                  self.acquisition_geometry.pixel_num_v,
+                  self.acquisition_geometry.pixel_num_h), \
+                 (self.volume_geometry.voxel_num_x, \
+                  self.volume_geometry.voxel_num_y,
+                  self.volume_geometry.voxel_num_z) )
+    def create_image_data(self):
+        x0 = ImageData(geometry = self.volume_geometry, 
+                       dimension_labels=self.bp.output_axes_order)#\
+                       #.subset(['horizontal_x','horizontal_y','vertical'])
+        print (x0)
+        x0.fill(numpy.random.randn(*x0.shape))
+        return x0
diff --git a/Wrappers/Python/wip/simple_demo_astra.py b/Wrappers/Python/wip/simple_demo_astra.py
new file mode 100755
index 0000000..3365504
--- /dev/null
+++ b/Wrappers/Python/wip/simple_demo_astra.py
@@ -0,0 +1,189 @@
+#import sys
+#sys.path.append("..")
+
+from ccpi.framework import ImageData , ImageGeometry, AcquisitionGeometry
+from ccpi.reconstruction.algs import FISTA, FBPD, CGLS
+from ccpi.reconstruction.funcs import Norm2sq, Norm1 , TV2D
+from ccpi.astra.astra_ops import AstraProjectorSimple
+
+import numpy as np
+import matplotlib.pyplot as plt
+
+test_case = 1   # 1=parallel2D, 2=cone2D
+
+# Set up phantom
+N = 128
+
+
+vg = ImageGeometry(voxel_num_x=N,voxel_num_y=N)
+Phantom = ImageData(geometry=vg)
+
+x = Phantom.as_array()
+x[round(N/4):round(3*N/4),round(N/4):round(3*N/4)] = 0.5
+x[round(N/8):round(7*N/8),round(3*N/8):round(5*N/8)] = 1
+
+plt.imshow(x)
+plt.show()
+
+# 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 = AcquisitionGeometry('parallel',
+                          '2D',
+                          angles,
+                          det_num,det_w)
+elif test_case==2:
+    pg = AcquisitionGeometry('cone',
+                          '2D',
+                          angles,
+                          det_num,
+                          det_w,
+                          dist_source_center=SourceOrig, 
+                          dist_center_detector=OrigDetec)
+
+# ASTRA operator using volume and sinogram geometries
+Aop = AstraProjectorSimple(vg, pg, 'cpu')
+
+# 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)
+out2 = Aop.adjoint(b)
+
+plt.imshow(b.array)
+plt.show()
+
+plt.imshow(out2.array)
+plt.show()
+
+# Create least squares object instance with projector and data.
+f = Norm2sq(Aop,b,c=0.5)
+
+# Initial guess
+x_init = ImageData(np.zeros(x.shape),geometry=vg)
+
+# Run FISTA for least squares without regularization
+x_fista0, it0, timing0, criter0 = FISTA(x_init, f, None)
+
+plt.imshow(x_fista0.array)
+plt.title('FISTA0')
+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)
+
+plt.imshow(x_fista1.array)
+plt.title('FISTA1')
+plt.show()
+
+plt.semilogy(criter1)
+plt.show()
+
+# Run FBPD=Forward Backward Primal Dual method on least squares plus 1-norm
+opt = {'tol': 1e-4, 'iter': 100}
+x_fbpd1, it_fbpd1, timing_fbpd1, criter_fbpd1 = FBPD(x_init,None,f,g0,opt=opt)
+
+plt.imshow(x_fbpd1.array)
+plt.title('FBPD1')
+plt.show()
+
+plt.semilogy(criter_fbpd1)
+plt.show()
+
+# Now FBPD for least squares plus TV
+#lamtv = 1
+#gtv = TV2D(lamtv)
+
+#x_fbpdtv, it_fbpdtv, timing_fbpdtv, criter_fbpdtv = FBPD(x_init,None,f,gtv,opt=opt)
+
+#plt.imshow(x_fbpdtv.array)
+#plt.show()
+
+#plt.semilogy(criter_fbpdtv)
+#plt.show()
+
+
+# Run CGLS, which should agree with the FISTA0
+x_CGLS, it_CGLS, timing_CGLS, criter_CGLS = CGLS(x_init, Aop, b, opt )
+
+plt.imshow(x_CGLS.array)
+plt.title('CGLS')
+#plt.title('CGLS recon, compare FISTA0')
+plt.show()
+
+plt.semilogy(criter_CGLS)
+plt.title('CGLS criterion')
+plt.show()
+
+
+#%%
+
+clims = (0,1)
+cols = 3
+rows = 2
+current = 1
+fig = plt.figure()
+# projections row
+a=fig.add_subplot(rows,cols,current)
+a.set_title('phantom {0}'.format(np.shape(Phantom.as_array())))
+
+imgplot = plt.imshow(Phantom.as_array(),vmin=clims[0],vmax=clims[1])
+
+current = current + 1
+a=fig.add_subplot(rows,cols,current)
+a.set_title('FISTA0')
+imgplot = plt.imshow(x_fista0.as_array(),vmin=clims[0],vmax=clims[1])
+
+current = current + 1
+a=fig.add_subplot(rows,cols,current)
+a.set_title('FISTA1')
+imgplot = plt.imshow(x_fista1.as_array(),vmin=clims[0],vmax=clims[1])
+
+current = current + 1
+a=fig.add_subplot(rows,cols,current)
+a.set_title('FBPD1')
+imgplot = plt.imshow(x_fbpd1.as_array(),vmin=clims[0],vmax=clims[1])
+
+current = current + 1
+a=fig.add_subplot(rows,cols,current)
+a.set_title('CGLS')
+imgplot = plt.imshow(x_CGLS.as_array(),vmin=clims[0],vmax=clims[1])
+
+#current = current + 1
+#a=fig.add_subplot(rows,cols,current)
+#a.set_title('FBPD TV')
+#imgplot = plt.imshow(x_fbpdtv.as_array(),vmin=clims[0],vmax=clims[1])
+
+fig = plt.figure()
+# projections row
+b=fig.add_subplot(1,1,1)
+b.set_title('criteria')
+imgplot = plt.loglog(criter0 , label='FISTA0')
+imgplot = plt.loglog(criter1 , label='FISTA1')
+imgplot = plt.loglog(criter_fbpd1, label='FBPD1')
+imgplot = plt.loglog(criter_CGLS, label='CGLS')
+#imgplot = plt.loglog(criter_fbpdtv, label='FBPD TV')
+b.legend(loc='right')
+plt.show()
+#%%
\ No newline at end of file
diff --git a/Wrappers/Python/wip/simple_demo_ccpi.py b/Wrappers/Python/wip/simple_demo_ccpi.py
new file mode 100755
index 0000000..4eeaaae
--- /dev/null
+++ b/Wrappers/Python/wip/simple_demo_ccpi.py
@@ -0,0 +1,270 @@
+#import sys
+#sys.path.append("..")
+
+from ccpi.framework import ImageData , AcquisitionData, ImageGeometry, AcquisitionGeometry
+from ccpi.reconstruction.algs import FISTA, FBPD, CGLS
+from ccpi.reconstruction.funcs import Norm2sq, Norm1 , TV2D
+from ccpi.astra.astra_ops import AstraProjectorSimple
+from ccpi.reconstruction.ops import CCPiProjectorSimple
+from ccpi.reconstruction.parallelbeam import alg as pbalg
+from ccpi.processors import CCPiForwardProjector, CCPiBackwardProjector 
+
+import numpy as np
+import matplotlib.pyplot as plt
+
+test_case = 1   # 1=parallel2D, 2=cone2D, 3=parallel3D
+
+# Set up phantom
+N = 128
+vert = 1
+# Set up measurement geometry
+angles_num = 20; # angles number
+det_w = 1.0
+det_num = N
+SourceOrig = 200
+OrigDetec = 0
+
+if test_case==1:
+    angles = np.linspace(0,np.pi,angles_num,endpoint=False,dtype=np.float32)*180/np.pi
+    #nangles = angles_num
+    #angles = np.linspace(0,360, nangles, dtype=np.float32)
+
+elif test_case==2:
+    angles = np.linspace(0,2*np.pi,angles_num,endpoint=False)
+elif test_case == 3:
+    angles = np.linspace(0,np.pi,angles_num,endpoint=False)
+else:
+    NotImplemented
+
+
+def setupCCPiGeometries(voxel_num_x, voxel_num_y, voxel_num_z, angles, counter):
+    vg = ImageGeometry(voxel_num_x=voxel_num_x,voxel_num_y=voxel_num_y, voxel_num_z=voxel_num_z)
+    Phantom_ccpi = ImageData(geometry=vg,
+                        dimension_labels=['horizontal_x','horizontal_y','vertical'])
+    #.subset(['horizontal_x','horizontal_y','vertical'])
+    # ask the ccpi code what dimensions it would like
+        
+    voxel_per_pixel = 1
+    geoms = pbalg.pb_setup_geometry_from_image(Phantom_ccpi.as_array(),
+                                                angles,
+                                                voxel_per_pixel )
+    
+    pg = AcquisitionGeometry('parallel',
+                              '3D',
+                              angles,
+                              geoms['n_h'],det_w,
+                              geoms['n_v'], det_w #2D in 3D is a slice 1 pixel thick
+                              )
+    
+    center_of_rotation = Phantom_ccpi.get_dimension_size('horizontal_x') / 2
+    ad = AcquisitionData(geometry=pg,dimension_labels=['angle','vertical','horizontal'])
+    geoms_i = pbalg.pb_setup_geometry_from_acquisition(ad.as_array(),
+                                                angles,
+                                                center_of_rotation,
+                                                voxel_per_pixel )
+
+    #print (counter)
+    counter+=1
+    #print (geoms , geoms_i)
+    if counter < 4:
+        if (not ( geoms_i == geoms )):
+            print ("not equal and {0}".format(counter))
+            X = max(geoms['output_volume_x'], geoms_i['output_volume_x'])
+            Y = max(geoms['output_volume_y'], geoms_i['output_volume_y'])
+            Z = max(geoms['output_volume_z'], geoms_i['output_volume_z'])
+            return setupCCPiGeometries(X,Y,Z,angles, counter)
+        else:
+            print ("return geoms {0}".format(geoms))
+            return geoms
+    else:
+        print ("return geoms_i {0}".format(geoms_i))
+        return geoms_i
+
+geoms = setupCCPiGeometries(N,N,vert,angles,0)
+#%%
+#geoms = {'n_v': 12, 'output_volume_y': 128, 'n_angles': 20, 
+# 'output_volume_x': 128, 'output_volume_z': 12, 'n_h': 128}
+vg = ImageGeometry(voxel_num_x=geoms['output_volume_x'],
+                   voxel_num_y=geoms['output_volume_y'], 
+                   voxel_num_z=geoms['output_volume_z'])
+Phantom = ImageData(geometry=vg,dimension_labels=['horizontal_x','horizontal_y','vertical']) + 0.1
+    
+
+#x = Phantom.as_array()
+i = 0
+while i < geoms['n_v']:
+    #x = Phantom.subset(vertical=i, dimensions=['horizontal_x','horizontal_y']).array
+    x = Phantom.subset(vertical=i).array
+    x[round(N/4):round(3*N/4),round(N/4):round(3*N/4)] = 0.5
+    x[round(N/8):round(7*N/8),round(3*N/8):round(5*N/8)] = 0.98
+    Phantom.fill(x, vertical=i)
+    i += 1
+
+plt.imshow(Phantom.subset(vertical=0).as_array())
+plt.show()
+
+
+
+# Parallelbeam geometry test
+if test_case==1:
+    #Phantom_ccpi = Phantom.subset(['horizontal_x','horizontal_y','vertical'])
+    #Phantom_ccpi.geometry = vg.clone()
+    center_of_rotation = Phantom.get_dimension_size('horizontal_x') / 2
+        
+    pg = AcquisitionGeometry('parallel',
+                          '3D',
+                          angles,
+                          geoms['n_h'],det_w,
+                          geoms['n_v'], det_w #2D in 3D is a slice 1 pixel thick
+                          )
+elif test_case==2:
+    raise NotImplemented('cone beam projector not yet available')
+    pg = AcquisitionGeometry('cone',
+                          '2D',
+                          angles,
+                          det_num,
+                          det_w,
+                          vert, det_w, #2D in 3D is a slice 1 pixel thick 
+                          dist_source_center=SourceOrig, 
+                          dist_center_detector=OrigDetec)
+
+# ASTRA operator using volume and sinogram geometries
+#Aop = AstraProjectorSimple(vg, pg, 'cpu')
+Cop = CCPiProjectorSimple(vg, pg)
+
+# Try forward and backprojection
+b = Cop.direct(Phantom)
+out2 = Cop.adjoint(b)
+
+#%%
+for i in range(b.get_dimension_size('vertical')):
+    plt.imshow(b.subset(vertical=i).array)
+    #plt.imshow(Phantom.subset( vertical=i).array)
+    #plt.imshow(b.array[:,i,:])
+    plt.show()
+#%%
+
+plt.imshow(out2.subset( vertical=0).array)
+plt.show()
+
+# Create least squares object instance with projector and data.
+f = Norm2sq(Cop,b,c=0.5)
+
+# Initial guess
+x_init = ImageData(geometry=vg, dimension_labels=['horizontal_x','horizontal_y','vertical'])
+#invL = 0.5
+#g = f.grad(x_init)
+#print (g)
+#u = x_init - invL*f.grad(x_init)
+        
+#%%
+# Run FISTA for least squares without regularization
+opt = {'tol': 1e-4, 'iter': 100}
+x_fista0, it0, timing0, criter0 = FISTA(x_init, f, None, opt=opt)
+
+plt.imshow(x_fista0.subset(vertical=0).array)
+plt.title('FISTA0')
+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=opt)
+
+plt.imshow(x_fista0.subset(vertical=0).array)
+plt.title('FISTA1')
+plt.show()
+
+plt.semilogy(criter1)
+plt.show()
+
+# Run FBPD=Forward Backward Primal Dual method on least squares plus 1-norm
+x_fbpd1, it_fbpd1, timing_fbpd1, criter_fbpd1 = FBPD(x_init,None,f,g0,opt=opt)
+
+plt.imshow(x_fbpd1.subset(vertical=0).array)
+plt.title('FBPD1')
+plt.show()
+
+plt.semilogy(criter_fbpd1)
+plt.show()
+
+# Now FBPD for least squares plus TV
+#lamtv = 1
+#gtv = TV2D(lamtv)
+
+#x_fbpdtv, it_fbpdtv, timing_fbpdtv, criter_fbpdtv = FBPD(x_init,None,f,gtv,opt=opt)
+
+#plt.imshow(x_fbpdtv.subset(vertical=0).array)
+#plt.show()
+
+#plt.semilogy(criter_fbpdtv)
+#plt.show()  
+
+
+# Run CGLS, which should agree with the FISTA0
+x_CGLS, it_CGLS, timing_CGLS, criter_CGLS = CGLS(x_init, Cop, b, opt=opt)
+
+plt.imshow(x_CGLS.subset(vertical=0).array)
+plt.title('CGLS')
+plt.title('CGLS recon, compare FISTA0')
+plt.show()
+
+plt.semilogy(criter_CGLS)
+plt.title('CGLS criterion')
+plt.show()
+
+
+#%%
+
+clims = (0,1)
+cols = 3
+rows = 2
+current = 1
+fig = plt.figure()
+# projections row
+a=fig.add_subplot(rows,cols,current)
+a.set_title('phantom {0}'.format(np.shape(Phantom.as_array())))
+
+imgplot = plt.imshow(Phantom.subset(vertical=0).as_array(),vmin=clims[0],vmax=clims[1])
+
+current = current + 1
+a=fig.add_subplot(rows,cols,current)
+a.set_title('FISTA0')
+imgplot = plt.imshow(x_fista0.subset(vertical=0).as_array(),vmin=clims[0],vmax=clims[1])
+
+current = current + 1
+a=fig.add_subplot(rows,cols,current)
+a.set_title('FISTA1')
+imgplot = plt.imshow(x_fista1.subset(vertical=0).as_array(),vmin=clims[0],vmax=clims[1])
+
+current = current + 1
+a=fig.add_subplot(rows,cols,current)
+a.set_title('FBPD1')
+imgplot = plt.imshow(x_fbpd1.subset(vertical=0).as_array(),vmin=clims[0],vmax=clims[1])
+
+current = current + 1
+a=fig.add_subplot(rows,cols,current)
+a.set_title('CGLS')
+imgplot = plt.imshow(x_CGLS.subset(vertical=0).as_array(),vmin=clims[0],vmax=clims[1])
+
+plt.show()
+#%%
+#current = current + 1
+#a=fig.add_subplot(rows,cols,current)
+#a.set_title('FBPD TV')
+#imgplot = plt.imshow(x_fbpdtv.subset(vertical=0).as_array(),vmin=clims[0],vmax=clims[1])
+
+fig = plt.figure()
+# projections row
+b=fig.add_subplot(1,1,1)
+b.set_title('criteria')
+imgplot = plt.loglog(criter0 , label='FISTA0')
+imgplot = plt.loglog(criter1 , label='FISTA1')
+imgplot = plt.loglog(criter_fbpd1, label='FBPD1')
+imgplot = plt.loglog(criter_CGLS, label='CGLS')
+#imgplot = plt.loglog(criter_fbpdtv, label='FBPD TV')
+b.legend(loc='right')
+plt.show()
+#%%
\ No newline at end of file