diff options
| -rw-r--r-- | src/Python/ccpi/fista/FISTAReconstructor.py | 164 | ||||
| -rw-r--r-- | src/Python/test/readhd5.py | 2 | ||||
| -rw-r--r-- | src/Python/test_reconstructor.py | 144 |
3 files changed, 41 insertions, 269 deletions
diff --git a/src/Python/ccpi/fista/FISTAReconstructor.py b/src/Python/ccpi/fista/FISTAReconstructor.py index 1e76815..cbd27da 100644 --- a/src/Python/ccpi/fista/FISTAReconstructor.py +++ b/src/Python/ccpi/fista/FISTAReconstructor.py @@ -73,7 +73,8 @@ class FISTAReconstructor(): # 3. "A novel tomographic reconstruction method based on the robust # Student's t function for suppressing data outliers" D. Kazantsev et.al. # D. Kazantsev, 2016-17 - def __init__(self, projector_geometry, output_geometry, input_sinogram, **kwargs): + def __init__(self, projector_geometry, output_geometry, input_sinogram, + **kwargs): # handle parmeters: # obligatory parameters self.pars = dict() @@ -98,6 +99,7 @@ class FISTAReconstructor(): 'regularizer' , 'ring_lambda_R_L1', 'ring_alpha') + self.acceptedInputKeywords = kw # handle keyworded parameters if kwargs is not None: @@ -114,11 +116,14 @@ class FISTAReconstructor(): if 'weights' in kwargs.keys(): self.pars['weights'] = kwargs['weights'] else: - self.pars['weights'] = numpy.ones(numpy.shape(self.pars['input_sinogram'])) + self.pars['weights'] = \ + numpy.ones(numpy.shape( + self.pars['input_sinogram'])) if 'Lipschitz_constant' in kwargs.keys(): self.pars['Lipschitz_constant'] = kwargs['Lipschitz_constant'] else: - self.pars['Lipschitz_constant'] = self.calculateLipschitzConstantWithPowerMethod() + self.pars['Lipschitz_constant'] = \ + self.calculateLipschitzConstantWithPowerMethod() if not 'ideal_image' in kwargs.keys(): self.pars['ideal_image'] = None @@ -127,7 +132,8 @@ class FISTAReconstructor(): if self.pars['ideal_image'] == None: pass else: - self.pars['region_of_interest'] = numpy.nonzero(self.pars['ideal_image']>0.0) + self.pars['region_of_interest'] = numpy.nonzero( + self.pars['ideal_image']>0.0) if not 'regularizer' in kwargs.keys() : self.pars['regularizer'] = None @@ -140,7 +146,29 @@ class FISTAReconstructor(): + def setParameter(self, **kwargs): + '''set named parameter for the regularization engine + raises Exception if the named parameter is not recognized + Typical usage is: + + reg = Regularizer(Regularizer.Algorithm.SplitBregman_TV) + reg.setParameter(input=u0) + reg.setParameter(regularization_parameter=10.) + + it can be also used as + reg = Regularizer(Regularizer.Algorithm.SplitBregman_TV) + reg.setParameter(input=u0 , regularization_parameter=10.) + ''' + + for key , value in kwargs.items(): + if key in self.acceptedInputKeywords.keys(): + self.pars[key] = value + else: + raise Exception('Wrong parameter {0} for '.format(key) + + 'Reconstruction algorithm') + # setParameter + def calculateLipschitzConstantWithPowerMethod(self): ''' using Power method (PM) to establish L constant''' @@ -152,7 +180,8 @@ class FISTAReconstructor(): - if (proj_geom['type'] == 'parallel') or (proj_geom['type'] == 'parallel3d'): + if (proj_geom['type'] == 'parallel') or \ + (proj_geom['type'] == 'parallel3d'): #% for parallel geometry we can do just one slice #print('Calculating Lipshitz constant for parallel beam geometry...') niter = 5;# % number of iteration for the PM @@ -262,128 +291,3 @@ class FISTAReconstructor(): - - -def getEntry(location, nx): - for item in nx[location].keys(): - print (item) - - -print ("Loading Data") - -##fname = "D:\\Documents\\Dataset\\IMAT\\20170419_crabtomo\\crabtomo\\Sample\\IMAT00005153_crabstomo_Sample_000.tif" -####ind = [i * 1049 for i in range(360)] -#### use only 360 images -##images = 200 -##ind = [int(i * 1049 / images) for i in range(images)] -##stack_image = dxchange.reader.read_tiff_stack(fname, ind, digit=None, slc=None) - -#fname = "D:\\Documents\\Dataset\\CGLS\\24737_fd.nxs" -#fname = "C:\\Users\\ofn77899\\Documents\\CCPi\\CGLS\\24737_fd_2.nxs" -##fname = "/home/ofn77899/Reconstruction/CCPi-FISTA_Reconstruction/data/dendr.h5" -##nx = h5py.File(fname, "r") -## -### the data are stored in a particular location in the hdf5 -##for item in nx['entry1/tomo_entry/data'].keys(): -## print (item) -## -##data = nx.get('entry1/tomo_entry/data/rotation_angle') -##angles = numpy.zeros(data.shape) -##data.read_direct(angles) -##print (angles) -### angles should be in degrees -## -##data = nx.get('entry1/tomo_entry/data/data') -##stack = numpy.zeros(data.shape) -##data.read_direct(stack) -##print (data.shape) -## -##print ("Data Loaded") -## -## -### Normalize -##data = nx.get('entry1/tomo_entry/instrument/detector/image_key') -##itype = numpy.zeros(data.shape) -##data.read_direct(itype) -### 2 is dark field -##darks = [stack[i] for i in range(len(itype)) if itype[i] == 2 ] -##dark = darks[0] -##for i in range(1, len(darks)): -## dark += darks[i] -##dark = dark / len(darks) -###dark[0][0] = dark[0][1] -## -### 1 is flat field -##flats = [stack[i] for i in range(len(itype)) if itype[i] == 1 ] -##flat = flats[0] -##for i in range(1, len(flats)): -## flat += flats[i] -##flat = flat / len(flats) -###flat[0][0] = dark[0][1] -## -## -### 0 is projection data -##proj = [stack[i] for i in range(len(itype)) if itype[i] == 0 ] -##angle_proj = [angles[i] for i in range(len(itype)) if itype[i] == 0 ] -##angle_proj = numpy.asarray (angle_proj) -##angle_proj = angle_proj.astype(numpy.float32) -## -### normalized data are -### norm = (projection - dark)/(flat-dark) -## -##def normalize(projection, dark, flat, def_val=0.1): -## a = (projection - dark) -## b = (flat-dark) -## with numpy.errstate(divide='ignore', invalid='ignore'): -## c = numpy.true_divide( a, b ) -## c[ ~ numpy.isfinite( c )] = def_val # set to not zero if 0/0 -## return c -## -## -##norm = [normalize(projection, dark, flat) for projection in proj] -##norm = numpy.asarray (norm) -##norm = norm.astype(numpy.float32) - - -##niterations = 15 -##threads = 3 -## -##img_cgls = alg.cgls(norm, angle_proj, numpy.double(86.2), 1 , niterations, threads, False) -##img_mlem = alg.mlem(norm, angle_proj, numpy.double(86.2), 1 , niterations, threads, False) -##img_sirt = alg.sirt(norm, angle_proj, numpy.double(86.2), 1 , niterations, threads, False) -## -##iteration_values = numpy.zeros((niterations,)) -##img_cgls_conv = alg.cgls_conv(norm, angle_proj, numpy.double(86.2), 1 , niterations, threads, -## iteration_values, False) -##print ("iteration values %s" % str(iteration_values)) -## -##iteration_values = numpy.zeros((niterations,)) -##img_cgls_tikhonov = alg.cgls_tikhonov(norm, angle_proj, numpy.double(86.2), 1 , niterations, threads, -## numpy.double(1e-5), iteration_values , False) -##print ("iteration values %s" % str(iteration_values)) -##iteration_values = numpy.zeros((niterations,)) -##img_cgls_TVreg = alg.cgls_TVreg(norm, angle_proj, numpy.double(86.2), 1 , niterations, threads, -## numpy.double(1e-5), iteration_values , False) -##print ("iteration values %s" % str(iteration_values)) -## -## -####numpy.save("cgls_recon.npy", img_data) -##import matplotlib.pyplot as plt -##fig, ax = plt.subplots(1,6,sharey=True) -##ax[0].imshow(img_cgls[80]) -##ax[0].axis('off') # clear x- and y-axes -##ax[1].imshow(img_sirt[80]) -##ax[1].axis('off') # clear x- and y-axes -##ax[2].imshow(img_mlem[80]) -##ax[2].axis('off') # clear x- and y-axesplt.show() -##ax[3].imshow(img_cgls_conv[80]) -##ax[3].axis('off') # clear x- and y-axesplt.show() -##ax[4].imshow(img_cgls_tikhonov[80]) -##ax[4].axis('off') # clear x- and y-axesplt.show() -##ax[5].imshow(img_cgls_TVreg[80]) -##ax[5].axis('off') # clear x- and y-axesplt.show() -## -## -##plt.show() -## - diff --git a/src/Python/test/readhd5.py b/src/Python/test/readhd5.py index 10f25ed..eff6c43 100644 --- a/src/Python/test/readhd5.py +++ b/src/Python/test/readhd5.py @@ -12,7 +12,7 @@ def getEntry(nx, location): for item in nx[location].keys(): print (item) -filename = r'C:\Users\ofn77899\Documents\GitHub\CCPi-FISTA_reconstruction\Demos\DendrData.h5' +filename = r'/home/ofn77899/Reconstruction/CCPi-FISTA_Reconstruction/demos/DendrData.h5' nx = h5py.File(filename, "r") #getEntry(nx, '/') # I have exported the entries as children of / diff --git a/src/Python/test_reconstructor.py b/src/Python/test_reconstructor.py index a4a622b..a338d34 100644 --- a/src/Python/test_reconstructor.py +++ b/src/Python/test_reconstructor.py @@ -58,143 +58,11 @@ vol_geom = astra.creators.create_vol_geom( image_size_x, ## First pass the arguments to the FISTAReconstructor and test the ## Lipschitz constant -fistaRecon = FISTAReconstructor(proj_geom, vol_geom, Sino3D , weights=Weights3D) -print ("Lipschitz Constant {0}".format(fistaRecon.pars['Lipschitz_constant'])) - #N = params.vol_geom.GridColCount - -pars = dict() -pars['projector_geometry'] = proj_geom.copy() -pars['output_geometry'] = vol_geom.copy() -pars['input_sinogram'] = Sino3D.copy() -sliceZ , nangles , detectors = numpy.shape(Sino3D) -pars['detectors'] = detectors -pars['number_of_angles'] = nangles -pars['SlicesZ'] = sliceZ - - -#pars['weights'] = numpy.ones(numpy.shape(pars['input_sinogram'])) -pars['weights'] = Weights3D.copy() +fistaRecon = FISTAReconstructor(proj_geom, + vol_geom, + Sino3D , + weights=Weights3D) -N = pars['output_geometry']['GridColCount'] -proj_geom = pars['projector_geometry'] -vol_geom = pars['output_geometry'] -weights = pars['weights'] -SlicesZ = pars['SlicesZ'] - -if (proj_geom['type'] == 'parallel') or (proj_geom['type'] == 'parallel3d'): - #% for parallel geometry we can do just one slice - print('Calculating Lipshitz constant for parallel beam geometry...') - niter = 5;# % number of iteration for the PM - #N = params.vol_geom.GridColCount; - #x1 = rand(N,N,1); - x1 = numpy.random.rand(1,N,N) - #sqweight = sqrt(weights(:,:,1)); - sqweight = numpy.sqrt(weights[0]) - proj_geomT = proj_geom.copy(); - proj_geomT['DetectorRowCount'] = 1; - vol_geomT = vol_geom.copy(); - vol_geomT['GridSliceCount'] = 1; - - #[sino_id, y] = astra_create_sino3d_cuda(x1, proj_geomT, vol_geomT); - - import matplotlib.pyplot as plt - fig = [] - props = dict(boxstyle='round', facecolor='wheat', alpha=0.5) - - #a.set_title('Lipschitz') - for i in range(niter): -# [id,x1] = astra_create_backprojection3d_cuda(sqweight.*y, proj_geomT, vol_geomT); -# s = norm(x1(:)); -# x1 = x1/s; -# [sino_id, y] = astra_create_sino3d_cuda(x1, proj_geomT, vol_geomT); -# y = sqweight.*y; -# astra_mex_data3d('delete', sino_id); -# astra_mex_data3d('delete', id); - #print ("iteration {0}".format(i)) - fig.append(plt.figure()) +print ("Lipschitz Constant {0}".format(fistaRecon.pars['Lipschitz_constant'])) - a=fig[-1].add_subplot(1,2,1) - a.text(0.05, 0.95, "iteration {0}, x1".format(i), transform=a.transAxes, - fontsize=14,verticalalignment='top', bbox=props) - - imgplot = plt.imshow(x1[0].copy()) - - - sino_id, y = astra.creators.create_sino3d_gpu(x1, - proj_geomT, - vol_geomT) - a=fig[-1].add_subplot(1,2,2) - a.text(0.05, 0.95, "iteration {0}, y".format(i), - transform=a.transAxes, fontsize=14,verticalalignment='top', - bbox=props) - - imgplot = plt.imshow(y[0].copy()) - - y = (sqweight * y) # element wise multiplication - - #b=fig.add_subplot(2,1,2) - #imgplot = plt.imshow(x1[0]) - #plt.show() - - #astra_mex_data3d('delete', sino_id); - astra.matlab.data3d('delete', sino_id) - del x1 - - idx,x1 = astra.creators.create_backprojection3d_gpu((sqweight*y), - proj_geomT, - vol_geomT) - del y - - - s = numpy.linalg.norm(x1) - ### this line? - x1 = (x1/s) - -# ### this line? -# sino_id, y = astra.creators.create_sino3d_gpu(x1, -# proj_geomT, -# vol_geomT); -# y = sqweight * y; - astra.matlab.data3d('delete', sino_id); - astra.matlab.data3d('delete', idx) - print ("iteration {0} s= {1}".format(i,s)) - - #end - del proj_geomT - del vol_geomT - #plt.show() -else: - #% divergen beam geometry - print('Calculating Lipshitz constant for divergen beam geometry...') - niter = 8; #% number of iteration for PM - x1 = numpy.random.rand(SlicesZ , N , N); - #sqweight = sqrt(weights); - sqweight = numpy.sqrt(weights[0]) - - sino_id, y = astra.creators.create_sino3d_gpu(x1, proj_geom, vol_geom); - y = sqweight*y; - #astra_mex_data3d('delete', sino_id); - astra.matlab.data3d('delete', sino_id); - - for i in range(niter): - #[id,x1] = astra_create_backprojection3d_cuda(sqweight.*y, proj_geom, vol_geom); - idx,x1 = astra.creators.create_backprojection3d_gpu(sqweight*y, - proj_geom, - vol_geom) - s = numpy.linalg.norm(x1) - ### this line? - x1 = x1/s; - ### this line? - #[sino_id, y] = astra_create_sino3d_gpu(x1, proj_geom, vol_geom); - sino_id, y = astra.creators.create_sino3d_gpu(x1, - proj_geom, - vol_geom); - - y = sqweight*y; - #astra_mex_data3d('delete', sino_id); - #astra_mex_data3d('delete', id); - astra.matlab.data3d('delete', sino_id); - astra.matlab.data3d('delete', idx); - #end - #clear x1 - del x1 +## the calculation of the lipschitz constant should not start by itself |