Block Framework (#23)

* fix for new block framework

* refactoring

closes #22

* updated setup.py

5 files changed, 81 insertions, 60 deletions

diff --git a/Wrappers/Python/ccpi/plugins/ops.py b/Wrappers/Python/ccpi/plugins/operators/CCPiProjectorSimple.py
index f2e0dd3..eedcac3 100755
--- a/Wrappers/Python/ccpi/plugins/ops.py
+++ b/Wrappers/Python/ccpi/plugins/operators/CCPiProjectorSimple.py
@@ -18,7 +18,7 @@
 #   limitations under the License.

 

 import numpy

-from ccpi.optimisation.ops import Operator, PowerMethodNonsquare

+from ccpi.optimisation.operators import Operator, LinearOperator

 from ccpi.framework import ImageData, DataContainer , \

                            ImageGeometry, AcquisitionGeometry

 from ccpi.plugins.processors import CCPiBackwardProjector, \

@@ -39,8 +39,7 @@ class CCPiProjectorSimple(Operator):
             raise TypeError('Can only handle parallel beam')

         

         # set-up the geometries if compatible

-        geoms = setupCCPiGeometries(geomv.voxel_num_x,geomv.voxel_num_y, 

-                                    geomv.voxel_num_z, geomp.angles, 0)

+        geoms = setupCCPiGeometries(geomv, geomp, 0)

         

 

         vg = ImageGeometry(voxel_num_x=geoms['output_volume_x'],

@@ -74,7 +73,8 @@ class CCPiProjectorSimple(Operator):
         

         self.bp = CCPiBackwardProjector(image_geometry=vg,

                                     acquisition_geometry=pg,

-                                    output_axes_order=['horizontal_x','horizontal_y','vertical'])

+                                    output_axes_order=[ ImageGeometry.HORIZONTAL_X,

+                 ImageGeometry.HORIZONTAL_Y, ImageGeometry.VERTICAL])

                 

         # Initialise empty for singular value.

         self.s1 = None

@@ -85,17 +85,21 @@ class CCPiProjectorSimple(Operator):
         return True

 

     def direct(self, image_data, out=None):

-        self.fp.set_input(image_data)

+        self.fp.set_input(image_data.subset(dimensions=['horizontal_x','horizontal_y','vertical']))

         if out is None:

             out = self.fp.get_output()

+            #return out.subset(dimensions=[AcquisitionGeometry.ANGLE , AcquisitionGeometry.VERTICAL ,

+            #     AcquisitionGeometry.HORIZONTAL])

             return out

         else:

             out.fill(self.fp.get_output())

 

     def adjoint(self, acquisition_data, out=None):

-        self.bp.set_input(acquisition_data)

+        self.bp.set_input(acquisition_data.subset(dimensions=['angle','vertical','horizontal']))

         if out is None:

             out = self.bp.get_output()

+            #return out.subset(dimensions=[ImageGeometry.VERTICAL, ImageGeometry.HORIZONTAL_Y,

+            #     ImageGeometry.HORIZONTAL_X])

             return out

         else:

             out.fill(self.bp.get_output())

@@ -103,8 +107,13 @@ class CCPiProjectorSimple(Operator):
     #def delete(self):

     #    astra.data2d.delete(self.proj_id)

     

-    def get_max_sing_val(self):

-        a = PowerMethodNonsquare(self,10)

+    def norm(self):

+        x0 = self.domain_geometry().allocate(ImageGeometry.RANDOM,

+        #dimension_labels=['horizontal_x', 'horizontal_y','vertical']

+        )

+        print (x0.dimension_labels)

+

+        a = LinearOperator.PowerMethod(self,10, x0)

         self.s1 = a[0] 

         return self.s1

     

@@ -133,8 +142,8 @@ class CCPiProjectorSimple(Operator):
             self.vg.center_x,

             self.vg.center_y,

             self.vg.center_z,

-            self.vg.channels,

-            ['horizontal_x','horizontal_y','vertical'] )

+            self.vg.channels

+            )

 

     def range_geometry(self):

         return AcquisitionGeometry(self.ag.geom_type,

@@ -147,5 +156,5 @@ class CCPiProjectorSimple(Operator):
                                self.ag.dist_source_center,

                                self.ag.dist_center_detector,

                                self.ag.channels,

-                               ['angle','vertical','horizontal'])

+                               )

 

diff --git a/Wrappers/Python/ccpi/plugins/operators/__init__.py b/Wrappers/Python/ccpi/plugins/operators/__init__.py
new file mode 100644
index 0000000..f5ffd0c
--- /dev/null
+++ b/Wrappers/Python/ccpi/plugins/operators/__init__.py
@@ -0,0 +1 @@
+from .CCPiProjectorSimple import CCPiProjectorSimple
\ No newline at end of file
diff --git a/Wrappers/Python/ccpi/plugins/processors/__init__.py b/Wrappers/Python/ccpi/plugins/processors/__init__.py
new file mode 100644
index 0000000..e688671
--- /dev/null
+++ b/Wrappers/Python/ccpi/plugins/processors/__init__.py
@@ -0,0 +1,4 @@
+from .processors import setupCCPiGeometries

+from .processors import CCPiForwardProjector

+from .processors import CCPiBackwardProjector

+from .processors import AcquisitionDataPadder

diff --git a/Wrappers/Python/ccpi/plugins/processors.py b/Wrappers/Python/ccpi/plugins/processors/processors.py
index d31c03b..84de79f 100755
--- a/Wrappers/Python/ccpi/plugins/processors.py
+++ b/Wrappers/Python/ccpi/plugins/processors/processors.py
@@ -22,46 +22,50 @@ from ccpi.framework import DataProcessor, AcquisitionData,\
 from ccpi.reconstruction.parallelbeam import alg as pbalg

 import numpy

 

-def setupCCPiGeometries(voxel_num_x, voxel_num_y, voxel_num_z, angles, counter):

+def setupCCPiGeometries(ig, ag, counter):

+        Phantom_ccpi = ig.allocate(dimension_labels=[ImageGeometry.HORIZONTAL_X, 

+                                                 ImageGeometry.HORIZONTAL_Y,

+                                                 ImageGeometry.VERTICAL])    

     

-    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

+        angles = ag.angles

+        geoms = pbalg.pb_setup_geometry_from_image(Phantom_ccpi.as_array(),

+                                                    angles,

+                                                    voxel_per_pixel )

         

-    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'], 1.0,

-                              geoms['n_v'], 1.0 #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 )

-    

-    counter+=1

-    

-    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)

+        pg = AcquisitionGeometry('parallel',

+                                  '3D',

+                                  angles,

+                                  geoms['n_h'], 1.0,

+                                  geoms['n_v'], 1.0 #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'])

+        ad = pg.allocate(dimension_labels=[AcquisitionGeometry.ANGLE, 

+                                           AcquisitionGeometry.VERTICAL,

+                                           AcquisitionGeometry.HORIZONTAL])

+        geoms_i = pbalg.pb_setup_geometry_from_acquisition(ad.as_array(),

+                                                    angles,

+                                                    center_of_rotation,

+                                                    voxel_per_pixel )

+        

+        counter+=1

+        

+        if counter < 4:

+            print (geoms, geoms_i)

+            if (not ( geoms_i == geoms )):

+                print ("not equal and {} {} {}".format(counter, geoms['output_volume_z'], geoms_i['output_volume_z']))

+                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 ("happy now {} {} {}".format(counter, geoms['output_volume_z'], geoms_i['output_volume_z']))

+                

+                return geoms

         else:

-            return geoms

-    else:

-        return geoms_i

+            return geoms_i

 

 

 class CCPiForwardProjector(DataProcessor):

@@ -102,7 +106,7 @@ class CCPiForwardProjector(DataProcessor):
             raise ValueError("Expected input dimensions is 2 or 3, got {0}"\

                              .format(dataset.number_of_dimensions))

 

-    def process(self):

+    def process(self, out=None):

         

         volume = self.get_input()

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

@@ -115,12 +119,15 @@ class CCPiForwardProjector(DataProcessor):
             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 = self.acquisition_geometry.allocate(

+                dimension_labels=self.output_axes_order)

             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)

+                print ("transpose")

+                pixels = numpy.transpose(pixels, out_axes)

+            out.fill(pixels)

+            

             return out

         else:

             raise ValueError('Cannot process cone beam')

@@ -146,7 +153,8 @@ class CCPiBackwardProjector(DataProcessor):
                  output_axes_order=None):

         if output_axes_order is None:

             # default ccpi projector image storing order

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

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

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

         kwargs = {

                   'image_geometry'       : image_geometry, 

                   'acquisition_geometry' : acquisition_geometry,

@@ -165,7 +173,7 @@ class CCPiBackwardProjector(DataProcessor):
             raise ValueError("Expected input dimensions is 2 or 3, got {0}"\

                              .format(dataset.number_of_dimensions))

 

-    def process(self):

+    def process(self, out=None):

         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]:

@@ -185,9 +193,7 @@ class CCPiBackwardProjector(DataProcessor):
                          center_of_rotation, 

                          pixel_per_voxel

                          )

-            out = ImageData(geometry=self.image_geometry, 

-                            dimension_labels=self.default_image_axes_order)

-            

+            out = self.image_geometry.allocate()

             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)

@@ -231,7 +237,7 @@ class AcquisitionDataPadder(DataProcessor):
             raise ValueError("Expected input dimensions is 2 or 3, got {0}"\

                              .format(dataset.number_of_dimensions))

 

-    def process(self):

+    def process(self, out=None):

         projections = self.get_input()

         w = projections.get_dimension_size('horizontal')

         delta = w - 2 * self.center_of_rotation

diff --git a/Wrappers/Python/setup.py b/Wrappers/Python/setup.py
index f4c42e8..ad12d1c 100644
--- a/Wrappers/Python/setup.py
+++ b/Wrappers/Python/setup.py
@@ -32,7 +32,8 @@ if  cil_version == '':
 setup(
     name="ccpi-plugins",
     version=cil_version,
-    packages=['ccpi' , 'ccpi.plugins'],
+    packages=['ccpi' , 'ccpi.plugins', 'ccpi.plugins.operators',
+        'ccpi.plugins.processors'],
 
     # Project uses reStructuredText, so ensure that the docutils get
     # installed or upgraded on the target machine