Merge pull request #2 from vais-ral/demo_cleanup

Demo cleanup

3 files changed, 104 insertions, 91 deletions

diff --git a/Wrappers/Python/ccpi/plugins/ops.py b/Wrappers/Python/ccpi/plugins/ops.py
index 0028cf1..aeb51af 100755
--- a/Wrappers/Python/ccpi/plugins/ops.py
+++ b/Wrappers/Python/ccpi/plugins/ops.py
@@ -19,49 +19,61 @@
 

 import numpy

 from ccpi.optimisation.ops import Operator, PowerMethodNonsquare

-from ccpi.framework import ImageData, DataContainer

-from ccpi.plugins.processors import CCPiBackwardProjector, CCPiForwardProjector

-

-class LinearOperatorMatrix(Operator):

-    def __init__(self,A):

-        self.A = A

-        self.s1 = None   # Largest singular value, initially unknown

-        super(LinearOperatorMatrix, self).__init__()

-        

-    def direct(self,x):

-        return DataContainer(numpy.dot(self.A,x.as_array()))

-    

-    def adjoint(self,x):

-        return DataContainer(numpy.dot(self.A.transpose(),x.as_array()))

-    

-    def size(self):

-        return self.A.shape

-    

-    def get_max_sing_val(self):

-        # If unknown, compute and store. If known, simply return it.

-        if self.s1 is None:

-            self.s1 = svds(self.A,1,return_singular_vectors=False)[0]

-            return self.s1

-        else:

-            return self.s1

+from ccpi.framework import ImageData, DataContainer , \

+                           ImageGeometry, AcquisitionGeometry

+from ccpi.plugins.processors import CCPiBackwardProjector, \

+                                    CCPiForwardProjector , setupCCPiGeometries

         

 

 

 class CCPiProjectorSimple(Operator):

     """ASTRA projector modified to use DataSet and geometry."""

-    def __init__(self, geomv, geomp):

+    def __init__(self, geomv, geomp, default=False):

         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,

+        if geomp.geom_type == "cone":

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

+        

+

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

+                           voxel_num_y=geoms['output_volume_y'], 

+                           voxel_num_z=geoms['output_volume_z'])

+

+        pg = AcquisitionGeometry('parallel',

+                          '3D',

+                          geomp.angles,

+                          geoms['n_h'], geomp.pixel_size_h,

+                          geoms['n_v'], geomp.pixel_size_v #2D in 3D is a slice 1 pixel thick

+                          )

+        if not default:

+            # check if geometry is the same (on the voxels)

+            if not ( vg.voxel_num_x == geomv.voxel_num_x and \

+                     vg.voxel_num_y == geomv.voxel_num_y and \

+                     vg.voxel_num_z == geomv.voxel_num_z ):

+                msg = 'The required volume geometry will not work\nThe following would\n'

+                msg += vg.__str__()

+                raise ValueError(msg)

+            if not (pg.pixel_num_h == geomp.pixel_num_h and \

+                    pg.pixel_num_v == geomp.pixel_num_v and \

+                    len( pg.angles ) == len( geomp.angles ) ) :

+                msg = 'The required acquisition geometry will not work\nThe following would\n'

+                msg += pg.__str__()

+                raise ValueError(msg)

+        

+        self.fp = CCPiForwardProjector(image_geometry=vg,

+                                       acquisition_geometry=pg,

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

         

-        self.bp = CCPiBackwardProjector(image_geometry=geomv,

-                                    acquisition_geometry=geomp,

+        self.bp = CCPiBackwardProjector(image_geometry=vg,

+                                    acquisition_geometry=pg,

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

                 

         # Initialise empty for singular value.

diff --git a/Wrappers/Python/ccpi/plugins/processors.py b/Wrappers/Python/ccpi/plugins/processors.py
index e05c6ca..df580e0 100755
--- a/Wrappers/Python/ccpi/plugins/processors.py
+++ b/Wrappers/Python/ccpi/plugins/processors.py
@@ -27,6 +27,51 @@ from scipy import ndimage
 import matplotlib.pyplot as plt

 

 

+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'], 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 )

+

+    #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

+

+

 class CCPiForwardProjector(DataProcessor):

     '''Normalization based on flat and dark

     

diff --git a/Wrappers/Python/wip/simple_demo_ccpi.py b/Wrappers/Python/wip/simple_demo_ccpi.py
index 10bb75f..3fdc2d4 100755
--- a/Wrappers/Python/wip/simple_demo_ccpi.py
+++ b/Wrappers/Python/wip/simple_demo_ccpi.py
@@ -8,7 +8,6 @@ from ccpi.optimisation.funcs import Norm2sq, Norm1 , TV2D
 
 from ccpi.plugins.ops import CCPiProjectorSimple
 from ccpi.plugins.processors import CCPiForwardProjector, CCPiBackwardProjector 
-
 from ccpi.reconstruction.parallelbeam import alg as pbalg
 
 import numpy as np
@@ -18,7 +17,7 @@ test_case = 1   # 1=parallel2D, 2=cone2D, 3=parallel3D
 
 # Set up phantom
 N = 128
-vert = 1
+vert = 4
 # Set up measurement geometry
 angles_num = 20; # angles number
 det_w = 1.0
@@ -38,71 +37,28 @@ elif test_case == 3:
 else:
     NotImplemented
 
+vg = ImageGeometry(voxel_num_x=N,
+                   voxel_num_y=N, 
+                   voxel_num_z=vert)
 
-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
+while i < vert:
+    if vert > 1:
+        x = Phantom.subset(vertical=i).array
+    else:
+        x = Phantom.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)
+    if vert > 1 :
+        Phantom.fill(x, vertical=i)
     i += 1
 
-plt.imshow(Phantom.subset(vertical=0).as_array())
+if vert > 1:
+    plt.imshow(Phantom.subset(vertical=0).as_array())
+else:
+    plt.imshow(Phantom.as_array())
 plt.show()
 
 
@@ -116,8 +72,8 @@ if test_case==1:
     pg = AcquisitionGeometry('parallel',
                           '3D',
                           angles,
-                          geoms['n_h'],det_w,
-                          geoms['n_v'], det_w #2D in 3D is a slice 1 pixel thick
+                          N , det_w,
+                          vert , det_w #2D in 3D is a slice 1 pixel thick
                           )
 elif test_case==2:
     raise NotImplemented('cone beam projector not yet available')