diff options
Diffstat (limited to 'Wrappers/Python')
| -rw-r--r-- | Wrappers/Python/ccpi/io/reader.py | 32 | ||||
| -rwxr-xr-x | Wrappers/Python/ccpi/processors.py | 2 | ||||
| -rwxr-xr-x | Wrappers/Python/wip/test_reader_reconstr.py | 181 |
3 files changed, 203 insertions, 12 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
'''
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
diff --git a/Wrappers/Python/wip/test_reader_reconstr.py b/Wrappers/Python/wip/test_reader_reconstr.py new file mode 100755 index 0000000..325437c --- /dev/null +++ b/Wrappers/Python/wip/test_reader_reconstr.py @@ -0,0 +1,181 @@ +# -*- 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
+import pickle
+
+
+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
+
+
+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())
+
+norm = Normalizer(flat_field=flat, dark_field=dark)
+
+norm.set_input(reader.get_acquisition_data())
+
+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)
+pickle.dump(x_fista0, open("fista0.pkl", "wb"))
+
+
+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)
+pickle.dump(x_fista1, open("fista1.pkl", "wb"))
+
+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)
+pickle.dump(x_fbpd1, open("fbpd1.pkl", "wb"))
+
+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)
+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.semilogy(criter_CGLS)
+plt.title('CGLS criterion')
+#plt.show()
+
+
+cols = 4
+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())
+
+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 |