From e453b6b6c4dd4869fdcc6f2b76efe96b1178d864 Mon Sep 17 00:00:00 2001 From: Edoardo Pasca Date: Wed, 21 Mar 2018 17:26:54 +0000 Subject: fix and begin to change to lower_case --- Wrappers/Python/ccpi/io/reader.py | 32 +++++++++++++++++++++----------- 1 file changed, 21 insertions(+), 11 deletions(-) diff --git a/Wrappers/Python/ccpi/io/reader.py b/Wrappers/Python/ccpi/io/reader.py index ca5b380..c66b1e8 100644 --- a/Wrappers/Python/ccpi/io/reader.py +++ b/Wrappers/Python/ccpi/io/reader.py @@ -83,28 +83,28 @@ class NexusReader(object): print("Error reading nexus file") raise - def loadProjection(self, dimensions=None): + def load_projection(self, dimensions=None): ''' Loads the projection data from the nexus file. returns: numpy array with projection data ''' return self.load(dimensions, 0) - def loadFlat(self, dimensions=None): + def load_flat(self, dimensions=None): ''' Loads the flat field data from the nexus file. returns: numpy array with flat field data ''' return self.load(dimensions, 1) - def loadDark(self, dimensions=None): + def load_dark(self, dimensions=None): ''' Loads the Dark field data from the nexus file. returns: numpy array with dark field data ''' return self.load(dimensions, 2) - def getProjectionAngles(self): + def get_projection_angles(self): ''' This function returns the projection angles ''' @@ -122,7 +122,7 @@ class NexusReader(object): raise - def getSinogramDimensions(self): + def get_sinogram_dimensions(self): ''' Return the dimensions of the dataset ''' @@ -142,7 +142,7 @@ class NexusReader(object): print("Error reading nexus file") raise - def getProjectionDimensions(self): + def get_projection_dimensions(self): ''' Return the dimensions of the dataset ''' @@ -161,13 +161,23 @@ class NexusReader(object): print("Error reading nexus file") raise - def getAcquisitionData(self, dimensions=None): + def get_acquisition_data(self, dimensions=None): ''' This method load the acquisition data and given dimension and returns an AcquisitionData Object ''' - data = self.loadProjection(dimensions) - geometry = AcquisitionGeometry('parallel', '3D', self.getProjectionAngles()) - return AcquisitionData(data, geometry=geometry) + data = self.load_projection(dimensions) + dims = self.get_projection_dimensions() + geometry = AcquisitionGeometry('parallel', '3D', + self.get_projection_angles(), + pixel_num_h = dims[2], + pixel_size_h = 1 , + pixel_num_v = dims[1], + pixel_size_v = 1, + dist_source_center = None, + dist_center_detector = None, + channels = 1) + return AcquisitionData(data, geometry=geometry, + dimension_labels=['angle','vertical','horizontal']) class XTEKReader(object): @@ -263,7 +273,7 @@ class XTEKReader(object): raise RuntimeError("Can't find angles file") return angles - def loadProjection(self, dimensions=None): + def load_projection(self, dimensions=None): ''' This method reads the projection images from the directory and returns a numpy array ''' -- cgit v1.2.3 From 12e0e7104a2f9da2681e5a6a6ddd309fe259426f Mon Sep 17 00:00:00 2001 From: Edoardo Pasca Date: Wed, 21 Mar 2018 17:27:33 +0000 Subject: bugfix, cor undefined --- Wrappers/Python/ccpi/processors.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Wrappers/Python/ccpi/processors.py b/Wrappers/Python/ccpi/processors.py index d98ef12..8f60cf4 100755 --- a/Wrappers/Python/ccpi/processors.py +++ b/Wrappers/Python/ccpi/processors.py @@ -573,7 +573,7 @@ class AcquisitionDataPadder(DataSetProcessor): def process(self): projections = self.get_input() w = projections.get_dimension_size('horizontal') - delta = w - 2 * cor + delta = w - 2 * self.center_of_rotation padded_width = int ( numpy.ceil(abs(delta)) + w -- cgit v1.2.3 From f082d54dfdef5335524806b703d630985afb247c Mon Sep 17 00:00:00 2001 From: Edoardo Pasca Date: Wed, 21 Mar 2018 17:28:09 +0000 Subject: initial revision --- Wrappers/Python/wip/test_reader_reconstr.py | 198 ++++++++++++++++++++++++++++ 1 file changed, 198 insertions(+) create mode 100755 Wrappers/Python/wip/test_reader_reconstr.py diff --git a/Wrappers/Python/wip/test_reader_reconstr.py b/Wrappers/Python/wip/test_reader_reconstr.py new file mode 100755 index 0000000..33ab461 --- /dev/null +++ b/Wrappers/Python/wip/test_reader_reconstr.py @@ -0,0 +1,198 @@ +# -*- coding: utf-8 -*- +""" +Created on Wed Mar 21 14:26:21 2018 + +@author: ofn77899 +""" + +from ccpi.framework import ImageData , AcquisitionData, ImageGeometry, AcquisitionGeometry +from ccpi.reconstruction.algs import FISTA, FBPD, CGLS +from ccpi.reconstruction.funcs import Norm2sq, Norm1 +from ccpi.reconstruction.ops import CCPiProjectorSimple +from ccpi.reconstruction.parallelbeam import alg as pbalg +from ccpi.processors import CCPiForwardProjector, CCPiBackwardProjector , \ +Normalizer , CenterOfRotationFinder , AcquisitionDataPadder + +from ccpi.io.reader import NexusReader + +import numpy +import matplotlib.pyplot as plt + +import os + +def add_dimension(data, fill_with, axis, start=True): + delta = data.shape[data.get_dimension_axis(axis)] - fill_with.shape[fill_with.get_dimension_axis(axis)] + command = 'data.array[' + i = 0 + for k,v in data.dimension_labels.items(): + if axis == v: + if start: + command = command + str(delta) + ":" + else: + l = data.get_dimension_size(axis) - delta + command = command + "0:" + str(l) + else: + command = command + ":" + + if i < data.number_of_dimensions -1: + command = command + ',' + i += 1 + command = command + "] = fill_with.array[:]" + #print (command) + exec(command) + #return command + + +def avg_img(image): + shape = list(numpy.shape(image)) + l = shape.pop(0) + avg = numpy.zeros(shape) + for i in range(l): + avg += image[i] / l + return avg + +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, + geoms['n_v'], 1 #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 + +reader = NexusReader(os.path.join(".." ,".." ,".." , "data" , "24737_fd.nxs" )) + +dims = reader.get_projection_dimensions() +print (dims) + +flat = avg_img(reader.load_flat()) +dark = avg_img(reader.load_dark()) + +data = reader.getAcquisitionData() +norm = Normalizer(flat_field=flat, dark_field=dark) + +norm.set_input(reader.getAcquisitionData()) + +cor = CenterOfRotationFinder() +cor.set_input(norm.get_output()) +center_of_rotation = cor.get_output() +voxel_per_pixel = 1 + +padder = AcquisitionDataPadder(center_of_rotation=center_of_rotation) +padder.set_input(norm.get_output()) +padded_data = padder.get_output() + +pg = padded_data.geometry +geoms = pbalg.pb_setup_geometry_from_acquisition(padded_data.as_array(), + pg.angles, + center_of_rotation, + voxel_per_pixel ) +vg = ImageGeometry(voxel_num_x=geoms['output_volume_x'], + voxel_num_y=geoms['output_volume_y'], + voxel_num_z=geoms['output_volume_z']) +#data = numpy.reshape(reader.getAcquisitionData()) +print ("define projector") +Cop = CCPiProjectorSimple(vg, pg) +# Create least squares object instance with projector and data. +print ("Create least squares object instance with projector and data.") +f = Norm2sq(Cop,padded_data,c=0.5) +print ("Initial guess") +# 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) + +#%% +print ("run FISTA") +# Run FISTA for least squares without regularization +opt = {'tol': 1e-4, 'iter': 10} +x_fista0, it0, timing0, criter0 = FISTA(x_init, f, None, opt=opt) + +plt.imshow(x_fista0.subset(horizontal_x=80).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(horizontal_x=80).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(horizontal_x=80).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, padded_data, opt=opt) + +plt.imshow(x_CGLS.subset(horizontal_x=80).array) +plt.title('CGLS') +plt.title('CGLS recon, compare FISTA0') +plt.show() + +plt.semilogy(criter_CGLS) +plt.title('CGLS criterion') +plt.show() \ No newline at end of file -- cgit v1.2.3 From 7abea9ef222c1c63267945d210cd9a4ec39ee184 Mon Sep 17 00:00:00 2001 From: Edoardo Pasca Date: Wed, 21 Mar 2018 17:51:37 +0000 Subject: reasonable long run --- Wrappers/Python/wip/test_reader_reconstr.py | 139 ++++++++++++++-------------- 1 file changed, 67 insertions(+), 72 deletions(-) diff --git a/Wrappers/Python/wip/test_reader_reconstr.py b/Wrappers/Python/wip/test_reader_reconstr.py index 33ab461..79afd61 100755 --- a/Wrappers/Python/wip/test_reader_reconstr.py +++ b/Wrappers/Python/wip/test_reader_reconstr.py @@ -19,28 +19,7 @@ import numpy import matplotlib.pyplot as plt import os - -def add_dimension(data, fill_with, axis, start=True): - delta = data.shape[data.get_dimension_axis(axis)] - fill_with.shape[fill_with.get_dimension_axis(axis)] - command = 'data.array[' - i = 0 - for k,v in data.dimension_labels.items(): - if axis == v: - if start: - command = command + str(delta) + ":" - else: - l = data.get_dimension_size(axis) - delta - command = command + "0:" + str(l) - else: - command = command + ":" - - if i < data.number_of_dimensions -1: - command = command + ',' - i += 1 - command = command + "] = fill_with.array[:]" - #print (command) - exec(command) - #return command +import pickle def avg_img(image): @@ -51,48 +30,6 @@ def avg_img(image): avg += image[i] / l return avg -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, - geoms['n_v'], 1 #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 reader = NexusReader(os.path.join(".." ,".." ,".." , "data" , "24737_fd.nxs" )) @@ -141,12 +78,14 @@ x_init = ImageData(geometry=vg, dimension_labels=['horizontal_x','horizontal_y', #%% print ("run FISTA") # Run FISTA for least squares without regularization -opt = {'tol': 1e-4, 'iter': 10} +opt = {'tol': 1e-4, 'iter': 500} x_fista0, it0, timing0, criter0 = FISTA(x_init, f, None, opt=opt) +pickle.dump(x_fista0, open("fista0.pkl", "wb")) + plt.imshow(x_fista0.subset(horizontal_x=80).array) plt.title('FISTA0') -plt.show() +#plt.show() # Now least squares plus 1-norm regularization lam = 0.1 @@ -154,23 +93,25 @@ g0 = Norm1(lam) # Run FISTA for least squares plus 1-norm function. x_fista1, it1, timing1, criter1 = FISTA(x_init, f, g0,opt=opt) +pickle.dump(x_fista1, open("fista1.pkl", "wb")) plt.imshow(x_fista0.subset(horizontal_x=80).array) plt.title('FISTA1') -plt.show() +#plt.show() plt.semilogy(criter1) -plt.show() +#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) +pickle.dump(x_fbpd1, open("fbpd1.pkl", "wb")) plt.imshow(x_fbpd1.subset(horizontal_x=80).array) plt.title('FBPD1') -plt.show() +#plt.show() plt.semilogy(criter_fbpd1) -plt.show() +#plt.show() # Now FBPD for least squares plus TV #lamtv = 1 @@ -187,12 +128,66 @@ plt.show() # Run CGLS, which should agree with the FISTA0 x_CGLS, it_CGLS, timing_CGLS, criter_CGLS = CGLS(x_init, Cop, padded_data, opt=opt) - +pickle.dump(x_CGLS, open("cgls.pkl", "wb")) plt.imshow(x_CGLS.subset(horizontal_x=80).array) plt.title('CGLS') plt.title('CGLS recon, compare FISTA0') -plt.show() +#plt.show() plt.semilogy(criter_CGLS) plt.title('CGLS criterion') +#plt.show() + + +clims = (0,1) +cols = 5 +rows = 1 +current = 1 +fig = plt.figure() +# projections row + +current = current +a=fig.add_subplot(rows,cols,current) +a.set_title('FISTA0') +imgplot = plt.imshow(x_fista0.subset(horizontal_x=80).as_array()) + +current = current + 1 +a=fig.add_subplot(rows,cols,current) +a.set_title('FISTA1') +imgplot = plt.imshow(x_fista1.subset(horizontal_x=80).as_array()) + +current = current + 1 +a=fig.add_subplot(rows,cols,current) +a.set_title('FBPD1') +imgplot = plt.imshow(x_fbpd1.subset(horizontal_x=80).as_array()) + +current = current + 1 +a=fig.add_subplot(rows,cols,current) +a.set_title('CGLS') +imgplot = plt.imshow(x_CGLS.subset(horizontal_x=80).as_array()) + +current = current + 1 +a=fig.add_subplot(rows,cols,current) +a.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') +a.legend(loc='right') + +plt.show() + + +#%% +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 -- cgit v1.2.3 From ebbf875a6baeb0e0a77f2d90094538314c4fe91c Mon Sep 17 00:00:00 2001 From: Edoardo Pasca Date: Thu, 22 Mar 2018 14:33:57 +0000 Subject: bugfix: use get_acquisition_data --- Wrappers/Python/wip/test_reader_reconstr.py | 18 +++--------------- 1 file changed, 3 insertions(+), 15 deletions(-) diff --git a/Wrappers/Python/wip/test_reader_reconstr.py b/Wrappers/Python/wip/test_reader_reconstr.py index 79afd61..325437c 100755 --- a/Wrappers/Python/wip/test_reader_reconstr.py +++ b/Wrappers/Python/wip/test_reader_reconstr.py @@ -39,10 +39,9 @@ print (dims) flat = avg_img(reader.load_flat()) dark = avg_img(reader.load_dark()) -data = reader.getAcquisitionData() norm = Normalizer(flat_field=flat, dark_field=dark) -norm.set_input(reader.getAcquisitionData()) +norm.set_input(reader.get_acquisition_data()) cor = CenterOfRotationFinder() cor.set_input(norm.get_output()) @@ -78,7 +77,7 @@ x_init = ImageData(geometry=vg, dimension_labels=['horizontal_x','horizontal_y', #%% print ("run FISTA") # Run FISTA for least squares without regularization -opt = {'tol': 1e-4, 'iter': 500} +opt = {'tol': 1e-4, 'iter': 10} x_fista0, it0, timing0, criter0 = FISTA(x_init, f, None, opt=opt) pickle.dump(x_fista0, open("fista0.pkl", "wb")) @@ -139,8 +138,7 @@ plt.title('CGLS criterion') #plt.show() -clims = (0,1) -cols = 5 +cols = 4 rows = 1 current = 1 fig = plt.figure() @@ -166,16 +164,6 @@ a=fig.add_subplot(rows,cols,current) a.set_title('CGLS') imgplot = plt.imshow(x_CGLS.subset(horizontal_x=80).as_array()) -current = current + 1 -a=fig.add_subplot(rows,cols,current) -a.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') -a.legend(loc='right') - plt.show() -- cgit v1.2.3