diff options
| -rw-r--r-- | src/Python/test_reconstructor.py | 144 |
1 files changed, 6 insertions, 138 deletions
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 |