From 1f9e7441198a6f088dfed311a067e8acf05d4d3b Mon Sep 17 00:00:00 2001 From: Edoardo Pasca Date: Tue, 17 Apr 2018 18:06:41 +0100 Subject: removed LinearOperatorMatrix and added check of geometries in ccpi operator --- Wrappers/Python/ccpi/plugins/ops.py | 72 +++++++++++++++++++++---------------- 1 file changed, 42 insertions(+), 30 deletions(-) (limited to 'Wrappers') 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. -- cgit v1.2.3 From 1d7313fe6f66d4b87efab95c88cd510ff44e26e1 Mon Sep 17 00:00:00 2001 From: Edoardo Pasca Date: Tue, 17 Apr 2018 18:07:25 +0100 Subject: added setupCCPiGeometries --- Wrappers/Python/ccpi/plugins/processors.py | 45 ++++++++++++++++++++++++++++++ 1 file changed, 45 insertions(+) (limited to 'Wrappers') 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 -- cgit v1.2.3 From 0fc626543be8c0827ca1696ecae5f57ff05b0fbc Mon Sep 17 00:00:00 2001 From: Edoardo Pasca Date: Tue, 17 Apr 2018 18:08:22 +0100 Subject: started cleaning demo --- Wrappers/Python/wip/simple_demo_ccpi.py | 78 +++++++-------------------------- 1 file changed, 17 insertions(+), 61 deletions(-) (limited to 'Wrappers') 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') -- cgit v1.2.3