diff options
Diffstat (limited to 'src/Python/test')
| -rw-r--r-- | src/Python/test/astra_test.py | 85 | ||||
| -rw-r--r-- | src/Python/test/create_phantom_projections.py | 49 | ||||
| -rw-r--r-- | src/Python/test/readhd5.py | 42 | ||||
| -rw-r--r-- | src/Python/test/simple_astra_test.py | 25 | ||||
| -rw-r--r-- | src/Python/test/test_reconstructor-os.py | 403 | ||||
| -rw-r--r-- | src/Python/test/test_reconstructor-os_phantom.py | 480 | ||||
| -rw-r--r-- | src/Python/test/test_reconstructor.py | 359 | ||||
| -rw-r--r-- | src/Python/test/test_regularizers.py | 412 | ||||
| -rw-r--r-- | src/Python/test/test_regularizers_3d.py | 425 |
9 files changed, 0 insertions, 2280 deletions
diff --git a/src/Python/test/astra_test.py b/src/Python/test/astra_test.py deleted file mode 100644 index 42c375a..0000000 --- a/src/Python/test/astra_test.py +++ /dev/null @@ -1,85 +0,0 @@ -import astra -import numpy -import filefun - - -# read in the same data as the DemoRD2 -angles = filefun.dlmread("DemoRD2/angles.csv") -darks_ar = filefun.dlmread("DemoRD2/darks_ar.csv", separator=",") -flats_ar = filefun.dlmread("DemoRD2/flats_ar.csv", separator=",") - -if True: - Sino3D = numpy.load("DemoRD2/Sino3D.npy") -else: - sino = filefun.dlmread("DemoRD2/sino_01.csv", separator=",") - a = map (lambda x:x, numpy.shape(sino)) - a.append(20) - - Sino3D = numpy.zeros(tuple(a), dtype="float") - - for i in range(1,numpy.shape(Sino3D)[2]+1): - print("Read file DemoRD2/sino_%02d.csv" % i) - sino = filefun.dlmread("DemoRD2/sino_%02d.csv" % i, separator=",") - Sino3D.T[i-1] = sino.T - -Weights3D = numpy.asarray(Sino3D, dtype="float") - -##angles_rad = angles*(pi/180); % conversion to radians -##size_det = size(data_raw3D,1); % detectors dim -##angSize = size(data_raw3D, 2); % angles dim -##slices_tot = size(data_raw3D, 3); % no of slices -##recon_size = 950; % reconstruction size - - -angles_rad = angles * numpy.pi /180. -size_det, angSize, slices_tot = numpy.shape(Sino3D) -size_det, angSize, slices_tot = [int(i) for i in numpy.shape(Sino3D)] -recon_size = 950 -Z_slices = 3; -det_row_count = Z_slices; - -#proj_geom = astra_create_proj_geom('parallel3d', 1, 1, -# det_row_count, size_det, angles_rad); - -detectorSpacingX = 1.0 -detectorSpacingY = detectorSpacingX -proj_geom = astra.create_proj_geom('parallel3d', - detectorSpacingX, - detectorSpacingY, - det_row_count, - size_det, - angles_rad) - -#vol_geom = astra_create_vol_geom(recon_size,recon_size,Z_slices); -vol_geom = astra.create_vol_geom(recon_size,recon_size,Z_slices); - -sino = numpy.zeros((size_det, angSize, slices_tot), dtype="float") - -#weights = ones(size(sino)); -weights = numpy.ones(numpy.shape(sino)) - -##################################################################### -## PowerMethod for Lipschitz constant - -N = vol_geom['GridColCount'] -x1 = numpy.random.rand(1,N,N) -#sqweight = sqrt(weights(:,:,1)); -sqweight = numpy.sqrt(weights.T[0]).T -##proj_geomT = proj_geom; -proj_geomT = proj_geom.copy() -##proj_geomT.DetectorRowCount = 1; -proj_geomT['DetectorRowCount'] = 1 -##vol_geomT = vol_geom; -vol_geomT = vol_geom.copy() -##vol_geomT.GridSliceCount = 1; -vol_geomT['GridSliceCount'] = 1 - -##[sino_id, y] = astra_create_sino3d_cuda(x1, proj_geomT, vol_geomT); - -#sino_id, y = astra.create_sino3d_gpu(x1, proj_geomT, vol_geomT); -sino_id, y = astra.create_sino(x1, proj_geomT, vol_geomT); - -##y = sqweight.*y; -##astra_mex_data3d('delete', sino_id); - - diff --git a/src/Python/test/create_phantom_projections.py b/src/Python/test/create_phantom_projections.py deleted file mode 100644 index 20a9278..0000000 --- a/src/Python/test/create_phantom_projections.py +++ /dev/null @@ -1,49 +0,0 @@ -from ccpi.reconstruction.AstraDevice import AstraDevice -from ccpi.reconstruction.DeviceModel import DeviceModel -import h5py -import numpy -import matplotlib.pyplot as plt - -nx = h5py.File('phant3D_256.h5', "r") -phantom = numpy.asarray(nx.get('/dataset1')) -pX,pY,pZ = numpy.shape(phantom) - -filename = r'/home/ofn77899/Reconstruction/CCPi-FISTA_Reconstruction/demos/DendrData.h5' -nxa = h5py.File(filename, "r") -#getEntry(nx, '/') -# I have exported the entries as children of / -entries = [entry for entry in nxa['/'].keys()] -print (entries) - -angles_rad = numpy.asarray(nxa.get('/angles_rad'), dtype="float32") - - -device = AstraDevice( - DeviceModel.DeviceType.PARALLEL3D.value, - [ pX , pY , 1., 1., angles_rad], - [ pX, pY, pZ ] ) - - -proj = device.doForwardProject(phantom) -stack = [proj[:,i,:] for i in range(len(angles_rad))] -stack = numpy.asarray(stack) - - -fig = plt.figure() -a=fig.add_subplot(1,2,1) -a.set_title('proj') -imgplot = plt.imshow(proj[:,100,:]) -a=fig.add_subplot(1,2,2) -a.set_title('stack') -imgplot = plt.imshow(stack[100]) -plt.show() - -pf = h5py.File("phantom3D256_projections.h5" , "w") -pf.create_dataset("/projections", data=stack) -pf.create_dataset("/sinogram", data=proj) -pf.create_dataset("/angles", data=angles_rad) -pf.create_dataset("/reconstruction_volume" , data=numpy.asarray([pX, pY, pZ])) -pf.create_dataset("/camera/size" , data=numpy.asarray([pX , pY ])) -pf.create_dataset("/camera/spacing" , data=numpy.asarray([1.,1.])) -pf.flush() -pf.close() diff --git a/src/Python/test/readhd5.py b/src/Python/test/readhd5.py deleted file mode 100644 index eff6c43..0000000 --- a/src/Python/test/readhd5.py +++ /dev/null @@ -1,42 +0,0 @@ -# -*- coding: utf-8 -*- -""" -Created on Wed Aug 23 16:34:49 2017 - -@author: ofn77899 -""" - -import h5py -import numpy - -def getEntry(nx, location): - for item in nx[location].keys(): - print (item) - -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 / -entries = [entry for entry in nx['/'].keys()] -print (entries) - -Sino3D = numpy.asarray(nx.get('/Sino3D')) -Weights3D = numpy.asarray(nx.get('/Weights3D')) -angSize = numpy.asarray(nx.get('/angSize'), dtype=int)[0] -angles_rad = numpy.asarray(nx.get('/angles_rad')) -recon_size = numpy.asarray(nx.get('/recon_size'), dtype=int)[0] -size_det = numpy.asarray(nx.get('/size_det'), dtype=int)[0] - -slices_tot = numpy.asarray(nx.get('/slices_tot'), dtype=int)[0] - -#from ccpi.viewer.CILViewer2D import CILViewer2D -#v = CILViewer2D() -#v.setInputAsNumpy(Weights3D) -#v.startRenderLoop() - -import matplotlib.pyplot as plt -fig = plt.figure() - -a=fig.add_subplot(1,1,1) -a.set_title('noise') -imgplot = plt.imshow(Weights3D[0].T) -plt.show() diff --git a/src/Python/test/simple_astra_test.py b/src/Python/test/simple_astra_test.py deleted file mode 100644 index 905eeea..0000000 --- a/src/Python/test/simple_astra_test.py +++ /dev/null @@ -1,25 +0,0 @@ -import astra -import numpy - -detectorSpacingX = 1.0 -detectorSpacingY = 1.0 -det_row_count = 128 -det_col_count = 128 - -angles_rad = numpy.asarray([i for i in range(360)], dtype=float) / 180. * numpy.pi - -proj_geom = astra.creators.create_proj_geom('parallel3d', - detectorSpacingX, - detectorSpacingY, - det_row_count, - det_col_count, - angles_rad) - -image_size_x = 64 -image_size_y = 64 -image_size_z = 32 - -vol_geom = astra.creators.create_vol_geom(image_size_x,image_size_y,image_size_z) - -x1 = numpy.random.rand(image_size_z,image_size_y,image_size_x) -sino_id, y = astra.creators.create_sino3d_gpu(x1, proj_geom, vol_geom) diff --git a/src/Python/test/test_reconstructor-os.py b/src/Python/test/test_reconstructor-os.py deleted file mode 100644 index 21b7ecd..0000000 --- a/src/Python/test/test_reconstructor-os.py +++ /dev/null @@ -1,403 +0,0 @@ -# -*- coding: utf-8 -*- -""" -Created on Wed Aug 23 16:34:49 2017 - -@author: ofn77899 -Based on DemoRD2.m -""" - -import h5py -import numpy - -from ccpi.reconstruction.FISTAReconstructor import FISTAReconstructor -import astra -import matplotlib.pyplot as plt -from ccpi.imaging.Regularizer import Regularizer -from ccpi.reconstruction.AstraDevice import AstraDevice -from ccpi.reconstruction.DeviceModel import DeviceModel - -def RMSE(signal1, signal2): - '''RMSE Root Mean Squared Error''' - if numpy.shape(signal1) == numpy.shape(signal2): - err = (signal1 - signal2) - err = numpy.sum( err * err )/numpy.size(signal1); # MSE - err = sqrt(err); # RMSE - return err - else: - raise Exception('Input signals must have the same shape') - -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 / -entries = [entry for entry in nx['/'].keys()] -print (entries) - -Sino3D = numpy.asarray(nx.get('/Sino3D'), dtype="float32") -Weights3D = numpy.asarray(nx.get('/Weights3D'), dtype="float32") -angSize = numpy.asarray(nx.get('/angSize'), dtype=int)[0] -angles_rad = numpy.asarray(nx.get('/angles_rad'), dtype="float32") -recon_size = numpy.asarray(nx.get('/recon_size'), dtype=int)[0] -size_det = numpy.asarray(nx.get('/size_det'), dtype=int)[0] -slices_tot = numpy.asarray(nx.get('/slices_tot'), dtype=int)[0] - -Z_slices = 20 -det_row_count = Z_slices -# next definition is just for consistency of naming -det_col_count = size_det - -detectorSpacingX = 1.0 -detectorSpacingY = detectorSpacingX - - -proj_geom = astra.creators.create_proj_geom('parallel3d', - detectorSpacingX, - detectorSpacingY, - det_row_count, - det_col_count, - angles_rad) - -#vol_geom = astra_create_vol_geom(recon_size,recon_size,Z_slices); -image_size_x = recon_size -image_size_y = recon_size -image_size_z = Z_slices -vol_geom = astra.creators.create_vol_geom( image_size_x, - image_size_y, - image_size_z) - -## First pass the arguments to the FISTAReconstructor and test the -## Lipschitz constant -astradevice = AstraDevice(DeviceModel.DeviceType.PARALLEL3D.value, - [proj_geom['DetectorRowCount'] , - proj_geom['DetectorColCount'] , - proj_geom['DetectorSpacingX'] , - proj_geom['DetectorSpacingY'] , - proj_geom['ProjectionAngles'] - ], - [ - vol_geom['GridColCount'], - vol_geom['GridRowCount'], - vol_geom['GridSliceCount'] ] ) -fistaRecon = FISTAReconstructor(proj_geom, - vol_geom, - Sino3D , - weights=Weights3D, - device=astradevice) - -print ("Lipschitz Constant {0}".format(fistaRecon.pars['Lipschitz_constant'])) -fistaRecon.setParameter(number_of_iterations = 2) -fistaRecon.setParameter(Lipschitz_constant = 767893952.0) -fistaRecon.setParameter(ring_alpha = 21) -fistaRecon.setParameter(ring_lambda_R_L1 = 0.002) -fistaRecon.setParameter(ring_lambda_R_L1 = 0) -subsets = 8 -fistaRecon.setParameter(subsets=subsets) - - -#reg = Regularizer(Regularizer.Algorithm.FGP_TV) -#reg.setParameter(regularization_parameter=0.005, -# number_of_iterations=50) -reg = Regularizer(Regularizer.Algorithm.FGP_TV) -reg.setParameter(regularization_parameter=5e6, - tolerance_constant=0.0001, - number_of_iterations=50) - -fistaRecon.setParameter(regularizer=reg) -lc = fistaRecon.getParameter('Lipschitz_constant') -reg.setParameter(regularization_parameter=5e6/lc) - -## Ordered subset -if True: - subsets = 8 - fistaRecon.setParameter(subsets=subsets) - fistaRecon.createOrderedSubsets() -else: - angles = fistaRecon.getParameter('projector_geometry')['ProjectionAngles'] - #binEdges = numpy.linspace(angles.min(), - # angles.max(), - # subsets + 1) - binsDiscr, binEdges = numpy.histogram(angles, bins=subsets) - # get rearranged subset indices - IndicesReorg = numpy.zeros((numpy.shape(angles))) - counterM = 0 - for ii in range(binsDiscr.max()): - counter = 0 - for jj in range(subsets): - curr_index = ii + jj + counter - #print ("{0} {1} {2}".format(binsDiscr[jj] , ii, counterM)) - if binsDiscr[jj] > ii: - if (counterM < numpy.size(IndicesReorg)): - IndicesReorg[counterM] = curr_index - counterM = counterM + 1 - - counter = counter + binsDiscr[jj] - 1 - - -if False: - print ("Lipschitz Constant {0}".format(fistaRecon.pars['Lipschitz_constant'])) - print ("prepare for iteration") - fistaRecon.prepareForIteration() - - - - print("initializing ...") - if False: - # if X doesn't exist - #N = params.vol_geom.GridColCount - N = vol_geom['GridColCount'] - print ("N " + str(N)) - X = numpy.zeros((N,N,SlicesZ), dtype=numpy.float) - else: - #X = fistaRecon.initialize() - X = numpy.load("X.npy") - - print (numpy.shape(X)) - X_t = X.copy() - print ("initialized") - proj_geom , vol_geom, sino , \ - SlicesZ, weights , alpha_ring = fistaRecon.getParameter( - ['projector_geometry' , 'output_geometry', - 'input_sinogram', 'SlicesZ' , 'weights', 'ring_alpha']) - lambdaR_L1 , alpha_ring , weights , L_const= \ - fistaRecon.getParameter(['ring_lambda_R_L1', - 'ring_alpha' , 'weights', - 'Lipschitz_constant']) - - #fistaRecon.setParameter(number_of_iterations = 3) - iterFISTA = fistaRecon.getParameter('number_of_iterations') - # errors vector (if the ground truth is given) - Resid_error = numpy.zeros((iterFISTA)); - # objective function values vector - objective = numpy.zeros((iterFISTA)); - - - t = 1 - - - ## additional for - proj_geomSUB = proj_geom.copy() - fistaRecon.residual2 = numpy.zeros(numpy.shape(fistaRecon.pars['input_sinogram'])) - residual2 = fistaRecon.residual2 - sino_updt_FULL = fistaRecon.residual.copy() - r_x = fistaRecon.r.copy() - - print ("starting iterations") -## % Outer FISTA iterations loop - for i in range(fistaRecon.getParameter('number_of_iterations')): -## % With OS approach it becomes trickier to correlate independent subsets, hence additional work is required -## % one solution is to work with a full sinogram at times -## if ((i >= 3) && (lambdaR_L1 > 0)) -## [sino_id2, sino_updt2] = astra_create_sino3d_cuda(X, proj_geom, vol_geom); -## astra_mex_data3d('delete', sino_id2); -## end - # With OS approach it becomes trickier to correlate independent subsets, - # hence additional work is required one solution is to work with a full - # sinogram at times - - r_old = fistaRecon.r.copy() - t_old = t - SlicesZ, anglesNumb, Detectors = \ - numpy.shape(fistaRecon.getParameter('input_sinogram')) ## https://github.com/vais-ral/CCPi-FISTA_Reconstruction/issues/4 - if (i > 1 and lambdaR_L1 > 0) : - for kkk in range(anglesNumb): - - residual2[:,kkk,:] = (weights[:,kkk,:]).squeeze() * \ - ((sino_updt_FULL[:,kkk,:]).squeeze() - \ - (sino[:,kkk,:]).squeeze() -\ - (alpha_ring * r_x) - ) - - vec = fistaRecon.residual.sum(axis = 1) - #if SlicesZ > 1: - # vec = vec[:,1,:] # 1 or 0? - r_x = fistaRecon.r_x - # update ring variable - fistaRecon.r = (r_x - (1./L_const) * vec).copy() - - # subset loop - counterInd = 1 - geometry_type = fistaRecon.getParameter('projector_geometry')['type'] - angles = fistaRecon.getParameter('projector_geometry')['ProjectionAngles'] - -## if geometry_type == 'parallel' or \ -## geometry_type == 'fanflat' or \ -## geometry_type == 'fanflat_vec' : -## -## for kkk in range(SlicesZ): -## sino_id, sinoT[kkk] = \ -## astra.creators.create_sino3d_gpu( -## X_t[kkk:kkk+1], proj_geomSUB, vol_geom) -## sino_updt_Sub[kkk] = sinoT.T.copy() -## -## else: -## sino_id, sino_updt_Sub = \ -## astra.creators.create_sino3d_gpu(X_t, proj_geomSUB, vol_geom) -## -## astra.matlab.data3d('delete', sino_id) - - for ss in range(fistaRecon.getParameter('subsets')): - print ("Subset {0}".format(ss)) - X_old = X.copy() - t_old = t - - # the number of projections per subset - numProjSub = fistaRecon.getParameter('os_bins')[ss] - CurrSubIndices = fistaRecon.getParameter('os_indices')\ - [counterInd:counterInd+numProjSub] - #print ("Len CurrSubIndices {0}".format(numProjSub)) - mask = numpy.zeros(numpy.shape(angles), dtype=bool) - cc = 0 - for j in range(len(CurrSubIndices)): - mask[int(CurrSubIndices[j])] = True - proj_geomSUB['ProjectionAngles'] = angles[mask] - - shape = list(numpy.shape(fistaRecon.getParameter('input_sinogram'))) - shape[1] = numProjSub - sino_updt_Sub = numpy.zeros(shape) - - if geometry_type == 'parallel' or \ - geometry_type == 'fanflat' or \ - geometry_type == 'fanflat_vec' : - - for kkk in range(SlicesZ): - sino_id, sinoT = astra.creators.create_sino3d_gpu ( - X_t[kkk:kkk+1] , proj_geomSUB, vol_geom) - sino_updt_Sub[kkk] = sinoT.T.copy() - - else: - # for 3D geometry (watch the GPU memory overflow in ASTRA < 1.8) - sino_id, sino_updt_Sub = \ - astra.creators.create_sino3d_gpu (X_t, proj_geomSUB, vol_geom) - - astra.matlab.data3d('delete', sino_id) - - - - - ## RING REMOVAL - residual = fistaRecon.residual - - - if lambdaR_L1 > 0 : - print ("ring removal") - residualSub = numpy.zeros(shape) - ## for a chosen subset - ## for kkk = 1:numProjSub - ## indC = CurrSubIndeces(kkk); - ## residualSub(:,kkk,:) = squeeze(weights(:,indC,:)).*(squeeze(sino_updt_Sub(:,kkk,:)) - (squeeze(sino(:,indC,:)) - alpha_ring.*r_x)); - ## sino_updt_FULL(:,indC,:) = squeeze(sino_updt_Sub(:,kkk,:)); % filling the full sinogram - ## end - for kkk in range(numProjSub): - #print ("ring removal indC ... {0}".format(kkk)) - indC = int(CurrSubIndices[kkk]) - residualSub[:,kkk,:] = weights[:,indC,:].squeeze() * \ - (sino_updt_Sub[:,kkk,:].squeeze() - \ - sino[:,indC,:].squeeze() - alpha_ring * r_x) - # filling the full sinogram - sino_updt_FULL[:,indC,:] = sino_updt_Sub[:,kkk,:].squeeze() - - else: - #PWLS model - # I guess we need to use mask here instead - residualSub = weights[:,CurrSubIndices,:] * \ - ( sino_updt_Sub - sino[:,CurrSubIndices,:].squeeze() ) - objective[i] = 0.5 * numpy.linalg.norm(residualSub) - - if geometry_type == 'parallel' or \ - geometry_type == 'fanflat' or \ - geometry_type == 'fanflat_vec' : - # if geometry is 2D use slice-by-slice projection-backprojection - # routine - x_temp = numpy.zeros(numpy.shape(X), dtype=numpy.float32) - for kkk in range(SlicesZ): - - x_id, x_temp[kkk] = \ - astra.creators.create_backprojection3d_gpu( - residualSub[kkk:kkk+1], - proj_geomSUB, vol_geom) - - else: - x_id, x_temp = \ - astra.creators.create_backprojection3d_gpu( - residualSub, proj_geomSUB, vol_geom) - - astra.matlab.data3d('delete', x_id) - X = X_t - (1/L_const) * x_temp - - - - ## REGULARIZATION - ## SKIPPING FOR NOW - ## Should be simpli - # regularizer = fistaRecon.getParameter('regularizer') - # for slices: - # out = regularizer(input=X) - print ("regularizer") - reg = fistaRecon.getParameter('regularizer') - - X = reg(input=X, - output_all=False) - - - ## FINAL - print ("final") - lambdaR_L1 = fistaRecon.getParameter('ring_lambda_R_L1') - if lambdaR_L1 > 0: - fistaRecon.r = numpy.max( - numpy.abs(fistaRecon.r) - lambdaR_L1 , 0) * \ - numpy.sign(fistaRecon.r) - # updating r - r_x = fistaRecon.r + ((t_old-1)/t) * (fistaRecon.r - r_old) - - - if fistaRecon.getParameter('region_of_interest') is None: - string = 'Iteration Number {0} | Objective {1} \n' - print (string.format( i, objective[i])) - else: - ROI , X_ideal = fistaRecon.getParameter('region_of_interest', - 'ideal_image') - - Resid_error[i] = RMSE(X*ROI, X_ideal*ROI) - string = 'Iteration Number {0} | RMS Error {1} | Objective {2} \n' - print (string.format(i,Resid_error[i], objective[i])) - - numpy.save("X_out_os.npy", X) - -else: - astradevice = AstraDevice(DeviceModel.DeviceType.PARALLEL3D.value, - [proj_geom['DetectorRowCount'] , - proj_geom['DetectorColCount'] , - proj_geom['DetectorSpacingX'] , - proj_geom['DetectorSpacingY'] , - proj_geom['ProjectionAngles'] - ], - [ - vol_geom['GridColCount'], - vol_geom['GridRowCount'], - vol_geom['GridSliceCount'] ] ) - regul = Regularizer(Regularizer.Algorithm.FGP_TV) - regul.setParameter(regularization_parameter=5e6, - number_of_iterations=50, - tolerance_constant=1e-4, - TV_penalty=Regularizer.TotalVariationPenalty.isotropic) - - fistaRecon = FISTAReconstructor(proj_geom, - vol_geom, - Sino3D , - weights=Weights3D, - device=astradevice, - regularizer = regul, - subsets=8) - - print ("Lipschitz Constant {0}".format(fistaRecon.pars['Lipschitz_constant'])) - fistaRecon.setParameter(number_of_iterations = 2) - fistaRecon.setParameter(Lipschitz_constant = 767893952.0) - fistaRecon.setParameter(ring_alpha = 21) - fistaRecon.setParameter(ring_lambda_R_L1 = 0.002) - #fistaRecon.setParameter(subsets=8) - - lc = fistaRecon.getParameter('Lipschitz_constant') - fistaRecon.getParameter('regularizer').setParameter(regularization_parameter=5e6/lc) - - fistaRecon.prepareForIteration() - X = fistaRecon.iterate(numpy.load("X.npy")) diff --git a/src/Python/test/test_reconstructor-os_phantom.py b/src/Python/test/test_reconstructor-os_phantom.py deleted file mode 100644 index 01f1354..0000000 --- a/src/Python/test/test_reconstructor-os_phantom.py +++ /dev/null @@ -1,480 +0,0 @@ -# -*- coding: utf-8 -*- -""" -Created on Wed Aug 23 16:34:49 2017 - -@author: ofn77899 -Based on DemoRD2.m -""" - -import h5py -import numpy - -from ccpi.reconstruction.FISTAReconstructor import FISTAReconstructor -import astra -import matplotlib.pyplot as plt -from ccpi.imaging.Regularizer import Regularizer -from ccpi.reconstruction.AstraDevice import AstraDevice -from ccpi.reconstruction.DeviceModel import DeviceModel - -#from ccpi.viewer.CILViewer2D import * - - -def RMSE(signal1, signal2): - '''RMSE Root Mean Squared Error''' - if numpy.shape(signal1) == numpy.shape(signal2): - err = (signal1 - signal2) - err = numpy.sum( err * err )/numpy.size(signal1); # MSE - err = sqrt(err); # RMSE - return err - else: - raise Exception('Input signals must have the same shape') - -filename = r'/home/ofn77899/Reconstruction/CCPi-FISTA_Reconstruction/src/Python/test/phantom3D256_projections.h5' -nx = h5py.File(filename, "r") -#getEntry(nx, '/') -# I have exported the entries as children of / -entries = [entry for entry in nx['/'].keys()] -print (entries) - -projections = numpy.asarray(nx.get('/projections'), dtype="float32") -#Weights3D = numpy.asarray(nx.get('/Weights3D'), dtype="float32") -#angSize = numpy.asarray(nx.get('/angSize'), dtype=int)[0] -angles_rad = numpy.asarray(nx.get('/angles'), dtype="float32") -angSize = numpy.size(angles_rad) -image_size_x, image_size_y, image_size_z = \ - numpy.asarray(nx.get('/reconstruction_volume'), dtype=int) -det_col_count, det_row_count = \ - numpy.asarray(nx.get('/camera/size')) -#slices_tot = numpy.asarray(nx.get('/slices_tot'), dtype=int)[0] -detectorSpacingX, detectorSpacingY = numpy.asarray(nx.get('/camera/spacing'), dtype=int) - -Z_slices = 20 -#det_row_count = image_size_y -# next definition is just for consistency of naming -#det_col_count = image_size_x - -detectorSpacingX = 1.0 -detectorSpacingY = detectorSpacingX - - -proj_geom = astra.creators.create_proj_geom('parallel3d', - detectorSpacingX, - detectorSpacingY, - det_row_count, - det_col_count, - angles_rad) - -#vol_geom = astra_create_vol_geom(recon_size,recon_size,Z_slices); -##image_size_x = recon_size -##image_size_y = recon_size -##image_size_z = Z_slices -vol_geom = astra.creators.create_vol_geom( image_size_x, - image_size_y, - image_size_z) - -## First pass the arguments to the FISTAReconstructor and test the -## Lipschitz constant -astradevice = AstraDevice(DeviceModel.DeviceType.PARALLEL3D.value, - [proj_geom['DetectorRowCount'] , - proj_geom['DetectorColCount'] , - proj_geom['DetectorSpacingX'] , - proj_geom['DetectorSpacingY'] , - proj_geom['ProjectionAngles'] - ], - [ - vol_geom['GridColCount'], - vol_geom['GridRowCount'], - vol_geom['GridSliceCount'] ] ) -## create the sinogram -Sino3D = numpy.transpose(projections, axes=[1,0,2]) - -fistaRecon = FISTAReconstructor(proj_geom, - vol_geom, - Sino3D , - #weights=Weights3D, - device=astradevice) - -print ("Lipschitz Constant {0}".format(fistaRecon.pars['Lipschitz_constant'])) -fistaRecon.setParameter(number_of_iterations = 4) -#fistaRecon.setParameter(Lipschitz_constant = 767893952.0) -fistaRecon.setParameter(ring_alpha = 21) -fistaRecon.setParameter(ring_lambda_R_L1 = 0.002) -#fistaRecon.setParameter(ring_lambda_R_L1 = 0) -subsets = 8 -fistaRecon.setParameter(subsets=subsets) - - -#reg = Regularizer(Regularizer.Algorithm.FGP_TV) -#reg.setParameter(regularization_parameter=0.005, -# number_of_iterations=50) -reg = Regularizer(Regularizer.Algorithm.FGP_TV) -reg.setParameter(regularization_parameter=5e6, - tolerance_constant=0.0001, - number_of_iterations=50) - -#fistaRecon.setParameter(regularizer=reg) -#lc = fistaRecon.getParameter('Lipschitz_constant') -#reg.setParameter(regularization_parameter=5e6/lc) - -## Ordered subset -if True: - #subsets = 8 - fistaRecon.setParameter(subsets=subsets) - fistaRecon.createOrderedSubsets() -else: - angles = fistaRecon.getParameter('projector_geometry')['ProjectionAngles'] - #binEdges = numpy.linspace(angles.min(), - # angles.max(), - # subsets + 1) - binsDiscr, binEdges = numpy.histogram(angles, bins=subsets) - # get rearranged subset indices - IndicesReorg = numpy.zeros((numpy.shape(angles))) - counterM = 0 - for ii in range(binsDiscr.max()): - counter = 0 - for jj in range(subsets): - curr_index = ii + jj + counter - #print ("{0} {1} {2}".format(binsDiscr[jj] , ii, counterM)) - if binsDiscr[jj] > ii: - if (counterM < numpy.size(IndicesReorg)): - IndicesReorg[counterM] = curr_index - counterM = counterM + 1 - - counter = counter + binsDiscr[jj] - 1 - - -if True: - print ("Lipschitz Constant {0}".format(fistaRecon.pars['Lipschitz_constant'])) - print ("prepare for iteration") - fistaRecon.prepareForIteration() - - - - print("initializing ...") - if True: - # if X doesn't exist - #N = params.vol_geom.GridColCount - N = vol_geom['GridColCount'] - print ("N " + str(N)) - X = numpy.asarray(numpy.ones((image_size_x,image_size_y,image_size_z)), - dtype=numpy.float) * 0.001 - X = numpy.asarray(numpy.zeros((image_size_x,image_size_y,image_size_z)), - dtype=numpy.float) - else: - #X = fistaRecon.initialize() - X = numpy.load("X.npy") - - print (numpy.shape(X)) - X_t = X.copy() - print ("initialized") - proj_geom , vol_geom, sino , \ - SlicesZ, weights , alpha_ring = fistaRecon.getParameter( - ['projector_geometry' , 'output_geometry', - 'input_sinogram', 'SlicesZ' , 'weights', 'ring_alpha']) - lambdaR_L1 , alpha_ring , weights , L_const= \ - fistaRecon.getParameter(['ring_lambda_R_L1', - 'ring_alpha' , 'weights', - 'Lipschitz_constant']) - - #fistaRecon.setParameter(number_of_iterations = 3) - iterFISTA = fistaRecon.getParameter('number_of_iterations') - # errors vector (if the ground truth is given) - Resid_error = numpy.zeros((iterFISTA)); - # objective function values vector - objective = numpy.zeros((iterFISTA)); - - - t = 1 - - - ## additional for - proj_geomSUB = proj_geom.copy() - fistaRecon.residual2 = numpy.zeros(numpy.shape(fistaRecon.pars['input_sinogram'])) - residual2 = fistaRecon.residual2 - sino_updt_FULL = fistaRecon.residual.copy() - r_x = fistaRecon.r.copy() - - results = [] - print ("starting iterations") -## % Outer FISTA iterations loop - for i in range(fistaRecon.getParameter('number_of_iterations')): -## % With OS approach it becomes trickier to correlate independent subsets, hence additional work is required -## % one solution is to work with a full sinogram at times -## if ((i >= 3) && (lambdaR_L1 > 0)) -## [sino_id2, sino_updt2] = astra_create_sino3d_cuda(X, proj_geom, vol_geom); -## astra_mex_data3d('delete', sino_id2); -## end - # With OS approach it becomes trickier to correlate independent subsets, - # hence additional work is required one solution is to work with a full - # sinogram at times - - - #t_old = t - SlicesZ, anglesNumb, Detectors = \ - numpy.shape(fistaRecon.getParameter('input_sinogram')) - ## https://github.com/vais-ral/CCPi-FISTA_Reconstruction/issues/4 - r_old = fistaRecon.r.copy() - - if (i > 1 and lambdaR_L1 > 0) : - for kkk in range(anglesNumb): - - residual2[:,kkk,:] = (weights[:,kkk,:]).squeeze() * \ - ((sino_updt_FULL[:,kkk,:]).squeeze() - \ - (sino[:,kkk,:]).squeeze() -\ - (alpha_ring * r_x) - ) - #r_old = fistaRecon.r.copy() - vec = fistaRecon.residual.sum(axis = 1) - #if SlicesZ > 1: - # vec = vec[:,1,:] # 1 or 0? - r_x = fistaRecon.r_x - # update ring variable - fistaRecon.r = (r_x - (1./L_const) * vec) - - # subset loop - counterInd = 1 - geometry_type = fistaRecon.getParameter('projector_geometry')['type'] - angles = fistaRecon.getParameter('projector_geometry')['ProjectionAngles'] - -## if geometry_type == 'parallel' or \ -## geometry_type == 'fanflat' or \ -## geometry_type == 'fanflat_vec' : -## -## for kkk in range(SlicesZ): -## sino_id, sinoT[kkk] = \ -## astra.creators.create_sino3d_gpu( -## X_t[kkk:kkk+1], proj_geomSUB, vol_geom) -## sino_updt_Sub[kkk] = sinoT.T.copy() -## -## else: -## sino_id, sino_updt_Sub = \ -## astra.creators.create_sino3d_gpu(X_t, proj_geomSUB, vol_geom) -## -## astra.matlab.data3d('delete', sino_id) - - for ss in range(fistaRecon.getParameter('subsets')): - print ("Subset {0}".format(ss)) - X_old = X.copy() - t_old = t - print ("X[0][0][0] {0} t {1}".format(X[0][0][0], t)) - - # the number of projections per subset - numProjSub = fistaRecon.getParameter('os_bins')[ss] - CurrSubIndices = fistaRecon.getParameter('os_indices')\ - [counterInd:counterInd+numProjSub] - shape = list(numpy.shape(fistaRecon.getParameter('input_sinogram'))) - shape[1] = numProjSub - sino_updt_Sub = numpy.zeros(shape) - - #print ("Len CurrSubIndices {0}".format(numProjSub)) - mask = numpy.zeros(numpy.shape(angles), dtype=bool) - cc = 0 - for j in range(len(CurrSubIndices)): - mask[int(CurrSubIndices[j])] = True - - ## this is a reduced device - rdev = fistaRecon.getParameter('device_model')\ - .createReducedDevice(proj_par={'angles' : angles[mask]}, - vol_par={}) - proj_geomSUB['ProjectionAngles'] = angles[mask] - - - - if geometry_type == 'parallel' or \ - geometry_type == 'fanflat' or \ - geometry_type == 'fanflat_vec' : - - for kkk in range(SlicesZ): - sino_id, sinoT = astra.creators.create_sino3d_gpu ( - X_t[kkk:kkk+1] , proj_geomSUB, vol_geom) - sino_updt_Sub[kkk] = sinoT.T.copy() - astra.matlab.data3d('delete', sino_id) - else: - # for 3D geometry (watch the GPU memory overflow in ASTRA < 1.8) - sino_id, sino_updt_Sub = \ - astra.creators.create_sino3d_gpu (X_t, - proj_geomSUB, - vol_geom) - - astra.matlab.data3d('delete', sino_id) - - - - - ## RING REMOVAL - residual = fistaRecon.residual - - - if lambdaR_L1 > 0 : - print ("ring removal") - residualSub = numpy.zeros(shape) - ## for a chosen subset - ## for kkk = 1:numProjSub - ## indC = CurrSubIndeces(kkk); - ## residualSub(:,kkk,:) = squeeze(weights(:,indC,:)).*(squeeze(sino_updt_Sub(:,kkk,:)) - (squeeze(sino(:,indC,:)) - alpha_ring.*r_x)); - ## sino_updt_FULL(:,indC,:) = squeeze(sino_updt_Sub(:,kkk,:)); % filling the full sinogram - ## end - for kkk in range(numProjSub): - #print ("ring removal indC ... {0}".format(kkk)) - indC = int(CurrSubIndices[kkk]) - residualSub[:,kkk,:] = weights[:,indC,:].squeeze() * \ - (sino_updt_Sub[:,kkk,:].squeeze() - \ - sino[:,indC,:].squeeze() - alpha_ring * r_x) - # filling the full sinogram - sino_updt_FULL[:,indC,:] = sino_updt_Sub[:,kkk,:].squeeze() - - else: - #PWLS model - # I guess we need to use mask here instead - residualSub = weights[:,CurrSubIndices,:] * \ - ( sino_updt_Sub - \ - sino[:,CurrSubIndices,:].squeeze() ) - # it seems that in the original code the following like is not - # calculated in the case of ring removal - objective[i] = 0.5 * numpy.linalg.norm(residualSub) - - #backprojection - if geometry_type == 'parallel' or \ - geometry_type == 'fanflat' or \ - geometry_type == 'fanflat_vec' : - # if geometry is 2D use slice-by-slice projection-backprojection - # routine - x_temp = numpy.zeros(numpy.shape(X), dtype=numpy.float32) - for kkk in range(SlicesZ): - - x_id, x_temp[kkk] = \ - astra.creators.create_backprojection3d_gpu( - residualSub[kkk:kkk+1], - proj_geomSUB, vol_geom) - astra.matlab.data3d('delete', x_id) - - else: - x_id, x_temp = \ - astra.creators.create_backprojection3d_gpu( - residualSub, proj_geomSUB, vol_geom) - - astra.matlab.data3d('delete', x_id) - - X = X_t - (1/L_const) * x_temp - - - - ## REGULARIZATION - ## SKIPPING FOR NOW - ## Should be simpli - # regularizer = fistaRecon.getParameter('regularizer') - # for slices: - # out = regularizer(input=X) - print ("regularizer") - reg = fistaRecon.getParameter('regularizer') - - if reg is not None: - X = reg(input=X, - output_all=False) - - t = (1 + numpy.sqrt(1 + 4 * t **2))/2 - X_t = X + (((t_old -1)/t) * (X-X_old)) - counterInd = counterInd + numProjSub - 1 - if i == 1: - r_old = fistaRecon.r.copy() - - ## FINAL - print ("final") - lambdaR_L1 = fistaRecon.getParameter('ring_lambda_R_L1') - if lambdaR_L1 > 0: - fistaRecon.r = numpy.max( - numpy.abs(fistaRecon.r) - lambdaR_L1 , 0) * \ - numpy.sign(fistaRecon.r) - # updating r - r_x = fistaRecon.r + ((t_old-1)/t) * (fistaRecon.r - r_old) - - - if fistaRecon.getParameter('region_of_interest') is None: - string = 'Iteration Number {0} | Objective {1} \n' - print (string.format( i, objective[i])) - else: - ROI , X_ideal = fistaRecon.getParameter('region_of_interest', - 'ideal_image') - - Resid_error[i] = RMSE(X*ROI, X_ideal*ROI) - string = 'Iteration Number {0} | RMS Error {1} | Objective {2} \n' - print (string.format(i,Resid_error[i], objective[i])) - - results.append(X[10]) - numpy.save("X_out_os.npy", X) - -else: - - - - astradevice = AstraDevice(DeviceModel.DeviceType.PARALLEL3D.value, - [proj_geom['DetectorRowCount'] , - proj_geom['DetectorColCount'] , - proj_geom['DetectorSpacingX'] , - proj_geom['DetectorSpacingY'] , - proj_geom['ProjectionAngles'] - ], - [ - vol_geom['GridColCount'], - vol_geom['GridRowCount'], - vol_geom['GridSliceCount'] ] ) - regul = Regularizer(Regularizer.Algorithm.FGP_TV) - regul.setParameter(regularization_parameter=5e6, - number_of_iterations=50, - tolerance_constant=1e-4, - TV_penalty=Regularizer.TotalVariationPenalty.isotropic) - - fistaRecon = FISTAReconstructor(proj_geom, - vol_geom, - Sino3D , - weights=Weights3D, - device=astradevice, - #regularizer = regul, - subsets=8) - - print ("Lipschitz Constant {0}".format(fistaRecon.pars['Lipschitz_constant'])) - fistaRecon.setParameter(number_of_iterations = 1) - fistaRecon.setParameter(Lipschitz_constant = 767893952.0) - fistaRecon.setParameter(ring_alpha = 21) - fistaRecon.setParameter(ring_lambda_R_L1 = 0.002) - #fistaRecon.setParameter(subsets=8) - - #lc = fistaRecon.getParameter('Lipschitz_constant') - #fistaRecon.getParameter('regularizer').setParameter(regularization_parameter=5e6/lc) - - fistaRecon.prepareForIteration() - X = fistaRecon.iterate(numpy.load("X.npy")) - - -# plot -fig = plt.figure() -cols = 3 - -## add the difference -rd = [] -for i in range(1,len(results)): - rd.append(results[i-1]) - rd.append(results[i]) - rd.append(results[i] - results[i-1]) - -rows = (lambda x: int(numpy.floor(x/cols) + 1) if x%cols != 0 else int(x/cols)) \ - (len (rd)) -for i in range(len (results)): - a=fig.add_subplot(rows,cols,i+1) - imgplot = plt.imshow(results[i], vmin=0, vmax=1) - a.text(0.05, 0.95, "iteration {0}".format(i), - verticalalignment='top') -## i = i + 1 -## a=fig.add_subplot(rows,cols,i+1) -## imgplot = plt.imshow(results[i], vmin=0, vmax=10) -## a.text(0.05, 0.95, "iteration {0}".format(i), -## verticalalignment='top') - -## a=fig.add_subplot(rows,cols,i+2) -## imgplot = plt.imshow(results[i]-results[i-1], vmin=0, vmax=10) -## a.text(0.05, 0.95, "difference {0}-{1}".format(i, i-1), -## verticalalignment='top') - - - -plt.show() diff --git a/src/Python/test/test_reconstructor.py b/src/Python/test/test_reconstructor.py deleted file mode 100644 index 40065e7..0000000 --- a/src/Python/test/test_reconstructor.py +++ /dev/null @@ -1,359 +0,0 @@ -# -*- coding: utf-8 -*- -""" -Created on Wed Aug 23 16:34:49 2017 - -@author: ofn77899 -Based on DemoRD2.m -""" - -import h5py -import numpy - -from ccpi.reconstruction.FISTAReconstructor import FISTAReconstructor -import astra -import matplotlib.pyplot as plt -from ccpi.imaging.Regularizer import Regularizer -from ccpi.reconstruction.AstraDevice import AstraDevice -from ccpi.reconstruction.DeviceModel import DeviceModel - -def RMSE(signal1, signal2): - '''RMSE Root Mean Squared Error''' - if numpy.shape(signal1) == numpy.shape(signal2): - err = (signal1 - signal2) - err = numpy.sum( err * err )/numpy.size(signal1); # MSE - err = sqrt(err); # RMSE - return err - else: - raise Exception('Input signals must have the same shape') - -def createAstraDevice(projector_geometry, output_geometry): - '''TODO remove''' - - device = AstraDevice(DeviceModel.DeviceType.PARALLEL3D.value, - [projector_geometry['DetectorRowCount'] , - projector_geometry['DetectorColCount'] , - projector_geometry['DetectorSpacingX'] , - projector_geometry['DetectorSpacingY'] , - projector_geometry['ProjectionAngles'] - ], - [ - output_geometry['GridColCount'], - output_geometry['GridRowCount'], - output_geometry['GridSliceCount'] ] ) - return device - -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 / -entries = [entry for entry in nx['/'].keys()] -print (entries) - -Sino3D = numpy.asarray(nx.get('/Sino3D'), dtype="float32") -Weights3D = numpy.asarray(nx.get('/Weights3D'), dtype="float32") -angSize = numpy.asarray(nx.get('/angSize'), dtype=int)[0] -angles_rad = numpy.asarray(nx.get('/angles_rad'), dtype="float32") -recon_size = numpy.asarray(nx.get('/recon_size'), dtype=int)[0] -size_det = numpy.asarray(nx.get('/size_det'), dtype=int)[0] -slices_tot = numpy.asarray(nx.get('/slices_tot'), dtype=int)[0] - -Z_slices = 20 -det_row_count = Z_slices -# next definition is just for consistency of naming -det_col_count = size_det - -detectorSpacingX = 1.0 -detectorSpacingY = detectorSpacingX - - -proj_geom = astra.creators.create_proj_geom('parallel3d', - detectorSpacingX, - detectorSpacingY, - det_row_count, - det_col_count, - angles_rad) - -#vol_geom = astra_create_vol_geom(recon_size,recon_size,Z_slices); -image_size_x = recon_size -image_size_y = recon_size -image_size_z = Z_slices -vol_geom = astra.creators.create_vol_geom( image_size_x, - image_size_y, - image_size_z) - -## 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'])) -##fistaRecon.setParameter(number_of_iterations = 12) -##fistaRecon.setParameter(Lipschitz_constant = 767893952.0) -##fistaRecon.setParameter(ring_alpha = 21) -##fistaRecon.setParameter(ring_lambda_R_L1 = 0.002) -## -##reg = Regularizer(Regularizer.Algorithm.LLT_model) -##reg.setParameter(regularization_parameter=25, -## time_step=0.0003, -## tolerance_constant=0.0001, -## number_of_iterations=300) -##fistaRecon.setParameter(regularizer=reg) - -## Ordered subset -if False: - subsets = 16 - angles = fistaRecon.getParameter('projector_geometry')['ProjectionAngles'] - #binEdges = numpy.linspace(angles.min(), - # angles.max(), - # subsets + 1) - binsDiscr, binEdges = numpy.histogram(angles, bins=subsets) - # get rearranged subset indices - IndicesReorg = numpy.zeros((numpy.shape(angles))) - counterM = 0 - for ii in range(binsDiscr.max()): - counter = 0 - for jj in range(subsets): - curr_index = ii + jj + counter - #print ("{0} {1} {2}".format(binsDiscr[jj] , ii, counterM)) - if binsDiscr[jj] > ii: - if (counterM < numpy.size(IndicesReorg)): - IndicesReorg[counterM] = curr_index - counterM = counterM + 1 - - counter = counter + binsDiscr[jj] - 1 - - -if False: - print ("Lipschitz Constant {0}".format(fistaRecon.pars['Lipschitz_constant'])) - print ("prepare for iteration") - fistaRecon.prepareForIteration() - - - - print("initializing ...") - if False: - # if X doesn't exist - #N = params.vol_geom.GridColCount - N = vol_geom['GridColCount'] - print ("N " + str(N)) - X = numpy.zeros((N,N,SlicesZ), dtype=numpy.float) - else: - #X = fistaRecon.initialize() - X = numpy.load("X.npy") - - print (numpy.shape(X)) - X_t = X.copy() - print ("initialized") - proj_geom , vol_geom, sino , \ - SlicesZ = fistaRecon.getParameter(['projector_geometry' , - 'output_geometry', - 'input_sinogram', - 'SlicesZ']) - - #fistaRecon.setParameter(number_of_iterations = 3) - iterFISTA = fistaRecon.getParameter('number_of_iterations') - # errors vector (if the ground truth is given) - Resid_error = numpy.zeros((iterFISTA)); - # objective function values vector - objective = numpy.zeros((iterFISTA)); - - - t = 1 - - - print ("starting iterations") -## % Outer FISTA iterations loop - for i in range(fistaRecon.getParameter('number_of_iterations')): - X_old = X.copy() - t_old = t - r_old = fistaRecon.r.copy() - if fistaRecon.getParameter('projector_geometry')['type'] == 'parallel' or \ - fistaRecon.getParameter('projector_geometry')['type'] == 'fanflat' or \ - fistaRecon.getParameter('projector_geometry')['type'] == 'fanflat_vec' : - # if the geometry is parallel use slice-by-slice - # projection-backprojection routine - #sino_updt = zeros(size(sino),'single'); - proj_geomT = proj_geom.copy() - proj_geomT['DetectorRowCount'] = 1 - vol_geomT = vol_geom.copy() - vol_geomT['GridSliceCount'] = 1; - sino_updt = numpy.zeros(numpy.shape(sino), dtype=numpy.float) - for kkk in range(SlicesZ): - sino_id, sino_updt[kkk] = \ - astra.creators.create_sino3d_gpu( - X_t[kkk:kkk+1], proj_geom, vol_geom) - astra.matlab.data3d('delete', sino_id) - else: - # for divergent 3D geometry (watch the GPU memory overflow in - # ASTRA versions < 1.8) - #[sino_id, sino_updt] = astra_create_sino3d_cuda(X_t, proj_geom, vol_geom); - sino_id, sino_updt = astra.creators.create_sino3d_gpu( - X_t, proj_geom, vol_geom) - - ## RING REMOVAL - residual = fistaRecon.residual - lambdaR_L1 , alpha_ring , weights , L_const= \ - fistaRecon.getParameter(['ring_lambda_R_L1', - 'ring_alpha' , 'weights', - 'Lipschitz_constant']) - r_x = fistaRecon.r_x - SlicesZ, anglesNumb, Detectors = \ - numpy.shape(fistaRecon.getParameter('input_sinogram')) - if lambdaR_L1 > 0 : - print ("ring removal") - for kkk in range(anglesNumb): - - residual[:,kkk,:] = (weights[:,kkk,:]).squeeze() * \ - ((sino_updt[:,kkk,:]).squeeze() - \ - (sino[:,kkk,:]).squeeze() -\ - (alpha_ring * r_x) - ) - vec = residual.sum(axis = 1) - #if SlicesZ > 1: - # vec = vec[:,1,:].squeeze() - fistaRecon.r = (r_x - (1./L_const) * vec).copy() - objective[i] = (0.5 * (residual ** 2).sum()) -## % the ring removal part (Group-Huber fidelity) -## for kkk = 1:anglesNumb -## residual(:,kkk,:) = squeeze(weights(:,kkk,:)).* -## (squeeze(sino_updt(:,kkk,:)) - -## (squeeze(sino(:,kkk,:)) - alpha_ring.*r_x)); -## end -## vec = sum(residual,2); -## if (SlicesZ > 1) -## vec = squeeze(vec(:,1,:)); -## end -## r = r_x - (1./L_const).*vec; -## objective(i) = (0.5*sum(residual(:).^2)); % for the objective function output - - - - # Projection/Backprojection Routine - if fistaRecon.getParameter('projector_geometry')['type'] == 'parallel' or \ - fistaRecon.getParameter('projector_geometry')['type'] == 'fanflat' or\ - fistaRecon.getParameter('projector_geometry')['type'] == 'fanflat_vec': - x_temp = numpy.zeros(numpy.shape(X),dtype=numpy.float32) - print ("Projection/Backprojection Routine") - for kkk in range(SlicesZ): - - x_id, x_temp[kkk] = \ - astra.creators.create_backprojection3d_gpu( - residual[kkk:kkk+1], - proj_geomT, vol_geomT) - astra.matlab.data3d('delete', x_id) - else: - x_id, x_temp = \ - astra.creators.create_backprojection3d_gpu( - residual, proj_geom, vol_geom) - - X = X_t - (1/L_const) * x_temp - astra.matlab.data3d('delete', sino_id) - astra.matlab.data3d('delete', x_id) - - - ## REGULARIZATION - ## SKIPPING FOR NOW - ## Should be simpli - # regularizer = fistaRecon.getParameter('regularizer') - # for slices: - # out = regularizer(input=X) - print ("skipping regularizer") - - - ## FINAL - print ("final") - lambdaR_L1 = fistaRecon.getParameter('ring_lambda_R_L1') - if lambdaR_L1 > 0: - fistaRecon.r = numpy.max( - numpy.abs(fistaRecon.r) - lambdaR_L1 , 0) * \ - numpy.sign(fistaRecon.r) - t = (1 + numpy.sqrt(1 + 4 * t**2))/2 - X_t = X + (((t_old -1)/t) * (X - X_old)) - - if lambdaR_L1 > 0: - fistaRecon.r_x = fistaRecon.r + \ - (((t_old-1)/t) * (fistaRecon.r - r_old)) - - if fistaRecon.getParameter('region_of_interest') is None: - string = 'Iteration Number {0} | Objective {1} \n' - print (string.format( i, objective[i])) - else: - ROI , X_ideal = fistaRecon.getParameter('region_of_interest', - 'ideal_image') - - Resid_error[i] = RMSE(X*ROI, X_ideal*ROI) - string = 'Iteration Number {0} | RMS Error {1} | Objective {2} \n' - print (string.format(i,Resid_error[i], objective[i])) - -## if (lambdaR_L1 > 0) -## r = max(abs(r)-lambdaR_L1, 0).*sign(r); % soft-thresholding operator for ring vector -## end -## -## t = (1 + sqrt(1 + 4*t^2))/2; % updating t -## X_t = X + ((t_old-1)/t).*(X - X_old); % updating X -## -## if (lambdaR_L1 > 0) -## r_x = r + ((t_old-1)/t).*(r - r_old); % updating r -## end -## -## if (show == 1) -## figure(10); imshow(X(:,:,slice), [0 maxvalplot]); -## if (lambdaR_L1 > 0) -## figure(11); plot(r); title('Rings offset vector') -## end -## pause(0.01); -## end -## if (strcmp(X_ideal, 'none' ) == 0) -## Resid_error(i) = RMSE(X(ROI), X_ideal(ROI)); -## fprintf('%s %i %s %s %.4f %s %s %f \n', 'Iteration Number:', i, '|', 'Error RMSE:', Resid_error(i), '|', 'Objective:', objective(i)); -## else -## fprintf('%s %i %s %s %f \n', 'Iteration Number:', i, '|', 'Objective:', objective(i)); -## end -else: - - # create a device for forward/backprojection - #astradevice = createAstraDevice(proj_geom, vol_geom) - - astradevice = AstraDevice(DeviceModel.DeviceType.PARALLEL3D.value, - [proj_geom['DetectorRowCount'] , - proj_geom['DetectorColCount'] , - proj_geom['DetectorSpacingX'] , - proj_geom['DetectorSpacingY'] , - proj_geom['ProjectionAngles'] - ], - [ - vol_geom['GridColCount'], - vol_geom['GridRowCount'], - vol_geom['GridSliceCount'] ] ) - - regul = Regularizer(Regularizer.Algorithm.FGP_TV) - regul.setParameter(regularization_parameter=5e6, - number_of_iterations=50, - tolerance_constant=1e-4, - TV_penalty=Regularizer.TotalVariationPenalty.isotropic) - - fistaRecon = FISTAReconstructor(proj_geom, - vol_geom, - Sino3D , - device = astradevice, - weights=Weights3D, - regularizer = regul - ) - - print ("Lipschitz Constant {0}".format(fistaRecon.pars['Lipschitz_constant'])) - fistaRecon.setParameter(number_of_iterations = 18) - fistaRecon.setParameter(Lipschitz_constant = 767893952.0) - fistaRecon.setParameter(ring_alpha = 21) - fistaRecon.setParameter(ring_lambda_R_L1 = 0.002) - - - - fistaRecon.prepareForIteration() - X = numpy.load("X.npy") - - - X = fistaRecon.iterate(X) - #numpy.save("X_out.npy", X) diff --git a/src/Python/test/test_regularizers.py b/src/Python/test/test_regularizers.py deleted file mode 100644 index 27e4ed3..0000000 --- a/src/Python/test/test_regularizers.py +++ /dev/null @@ -1,412 +0,0 @@ -# -*- coding: utf-8 -*- -""" -Created on Fri Aug 4 11:10:05 2017 - -@author: ofn77899 -""" - -#from ccpi.viewer.CILViewer2D import Converter -#import vtk - -import matplotlib.pyplot as plt -import numpy as np -import os -from enum import Enum -import timeit -#from PIL import Image -#from Regularizer import Regularizer -from ccpi.imaging.Regularizer import Regularizer - -############################################################################### -#https://stackoverflow.com/questions/13875989/comparing-image-in-url-to-image-in-filesystem-in-python/13884956#13884956 -#NRMSE a normalization of the root of the mean squared error -#NRMSE is simply 1 - [RMSE / (maxval - minval)]. Where maxval is the maximum -# intensity from the two images being compared, and respectively the same for -# minval. RMSE is given by the square root of MSE: -# sqrt[(sum(A - B) ** 2) / |A|], -# where |A| means the number of elements in A. By doing this, the maximum value -# given by RMSE is maxval. - -def nrmse(im1, im2): - a, b = im1.shape - rmse = np.sqrt(np.sum((im2 - im1) ** 2) / float(a * b)) - max_val = max(np.max(im1), np.max(im2)) - min_val = min(np.min(im1), np.min(im2)) - return 1 - (rmse / (max_val - min_val)) -############################################################################### - -############################################################################### -# -# 2D Regularizers -# -############################################################################### -#Example: -# figure; -# Im = double(imread('lena_gray_256.tif'))/255; % loading image -# u0 = Im + .05*randn(size(Im)); u0(u0 < 0) = 0; -# u = SplitBregman_TV(single(u0), 10, 30, 1e-04); - - -#filename = r"C:\Users\ofn77899\Documents\GitHub\CCPi-FISTA_reconstruction\data\lena_gray_512.tif" -filename = r"/home/ofn77899/Reconstruction/CCPi-FISTA_Reconstruction/data/lena_gray_512.tif" -#filename = r'/home/algol/Documents/Python/STD_test_images/lena_gray_512.tif' - -#reader = vtk.vtkTIFFReader() -#reader.SetFileName(os.path.normpath(filename)) -#reader.Update() -Im = plt.imread(filename) -#Im = Image.open('/home/algol/Documents/Python/STD_test_images/lena_gray_512.tif')/255 -#img.show() -Im = np.asarray(Im, dtype='float32') - - - - -#imgplot = plt.imshow(Im) -perc = 0.05 -u0 = Im + (perc* np.random.normal(size=np.shape(Im))) -# map the u0 u0->u0>0 -f = np.frompyfunc(lambda x: 0 if x < 0 else x, 1,1) -u0 = f(u0).astype('float32') - -## plot -fig = plt.figure() -#a=fig.add_subplot(3,3,1) -#a.set_title('Original') -#imgplot = plt.imshow(Im) - -a=fig.add_subplot(2,3,1) -a.set_title('noise') -imgplot = plt.imshow(u0,cmap="gray") - -reg_output = [] -############################################################################## -# Call regularizer - -####################### SplitBregman_TV ##################################### -# u = SplitBregman_TV(single(u0), 10, 30, 1e-04); - -use_object = True -if use_object: - reg = Regularizer(Regularizer.Algorithm.SplitBregman_TV) - print (reg.pars) - reg.setParameter(input=u0) - reg.setParameter(regularization_parameter=10.) - # or - # reg.setParameter(input=u0, regularization_parameter=10., #number_of_iterations=30, - #tolerance_constant=1e-4, - #TV_Penalty=Regularizer.TotalVariationPenalty.l1) - plotme = reg() [0] - pars = reg.pars - textstr = reg.printParametersToString() - - #out = reg(input=u0, regularization_parameter=10., #number_of_iterations=30, - #tolerance_constant=1e-4, - # TV_Penalty=Regularizer.TotalVariationPenalty.l1) - -#out2 = Regularizer.SplitBregman_TV(input=u0, regularization_parameter=10., number_of_iterations=30, -# tolerance_constant=1e-4, -# TV_Penalty=Regularizer.TotalVariationPenalty.l1) - -else: - out2 = Regularizer.SplitBregman_TV(input=u0, regularization_parameter=10. ) - pars = out2[2] - reg_output.append(out2) - plotme = reg_output[-1][0] - textstr = out2[-1] - -a=fig.add_subplot(2,3,2) - - -# these are matplotlib.patch.Patch properties -props = dict(boxstyle='round', facecolor='wheat', alpha=0.5) -# place a text box in upper left in axes coords -a.text(0.05, 0.95, textstr, transform=a.transAxes, fontsize=14, - verticalalignment='top', bbox=props) -imgplot = plt.imshow(plotme,cmap="gray") - -###################### FGP_TV ######################################### -# u = FGP_TV(single(u0), 0.05, 100, 1e-04); -out2 = Regularizer.FGP_TV(input=u0, regularization_parameter=0.0005, - number_of_iterations=50) -pars = out2[-2] - -reg_output.append(out2) - -a=fig.add_subplot(2,3,3) - -textstr = out2[-1] - -# these are matplotlib.patch.Patch properties -props = dict(boxstyle='round', facecolor='wheat', alpha=0.5) -# place a text box in upper left in axes coords -a.text(0.05, 0.95, textstr, transform=a.transAxes, fontsize=14, - verticalalignment='top', bbox=props) -imgplot = plt.imshow(reg_output[-1][0]) -# place a text box in upper left in axes coords -a.text(0.05, 0.95, textstr, transform=a.transAxes, fontsize=14, - verticalalignment='top', bbox=props) -imgplot = plt.imshow(reg_output[-1][0],cmap="gray") - -###################### LLT_model ######################################### -# * u0 = Im + .03*randn(size(Im)); % adding noise -# [Den] = LLT_model(single(u0), 10, 0.1, 1); -#Den = LLT_model(single(u0), 25, 0.0003, 300, 0.0001, 0); -#input, regularization_parameter , time_step, number_of_iterations, -# tolerance_constant, restrictive_Z_smoothing=0 -out2 = Regularizer.LLT_model(input=u0, regularization_parameter=25, - time_step=0.0003, - tolerance_constant=0.0001, - number_of_iterations=300) -pars = out2[-2] - -reg_output.append(out2) - -a=fig.add_subplot(2,3,4) - -textstr = out2[-1] - -# these are matplotlib.patch.Patch properties -props = dict(boxstyle='round', facecolor='wheat', alpha=0.5) -# place a text box in upper left in axes coords -a.text(0.05, 0.95, textstr, transform=a.transAxes, fontsize=14, - verticalalignment='top', bbox=props) -imgplot = plt.imshow(reg_output[-1][0],cmap="gray") - - -# ###################### PatchBased_Regul ######################################### -# # Quick 2D denoising example in Matlab: -# # Im = double(imread('lena_gray_256.tif'))/255; % loading image -# # u0 = Im + .03*randn(size(Im)); u0(u0<0) = 0; % adding noise -# # ImDen = PB_Regul_CPU(single(u0), 3, 1, 0.08, 0.05); - -out2 = Regularizer.PatchBased_Regul(input=u0, regularization_parameter=0.05, - searching_window_ratio=3, - similarity_window_ratio=1, - PB_filtering_parameter=0.08) -pars = out2[-2] -reg_output.append(out2) - -a=fig.add_subplot(2,3,5) - - -textstr = out2[-1] - -# these are matplotlib.patch.Patch properties -props = dict(boxstyle='round', facecolor='wheat', alpha=0.5) -# place a text box in upper left in axes coords -a.text(0.05, 0.95, textstr, transform=a.transAxes, fontsize=14, - verticalalignment='top', bbox=props) -imgplot = plt.imshow(reg_output[-1][0],cmap="gray") - - -# ###################### TGV_PD ######################################### -# # Quick 2D denoising example in Matlab: -# # Im = double(imread('lena_gray_256.tif'))/255; % loading image -# # u0 = Im + .03*randn(size(Im)); u0(u0<0) = 0; % adding noise -# # u = PrimalDual_TGV(single(u0), 0.02, 1.3, 1, 550); - - -out2 = Regularizer.TGV_PD(input=u0, regularization_parameter=0.05, - first_order_term=1.3, - second_order_term=1, - number_of_iterations=550) -pars = out2[-2] -reg_output.append(out2) - -a=fig.add_subplot(2,3,6) - - -textstr = out2[-1] - - -# these are matplotlib.patch.Patch properties -props = dict(boxstyle='round', facecolor='wheat', alpha=0.5) -# place a text box in upper left in axes coords -a.text(0.05, 0.95, textstr, transform=a.transAxes, fontsize=14, - verticalalignment='top', bbox=props) -imgplot = plt.imshow(reg_output[-1][0],cmap="gray") - - -plt.show() - -################################################################################ -## -## 3D Regularizers -## -################################################################################ -##Example: -## figure; -## Im = double(imread('lena_gray_256.tif'))/255; % loading image -## u0 = Im + .05*randn(size(Im)); u0(u0 < 0) = 0; -## u = SplitBregman_TV(single(u0), 10, 30, 1e-04); -# -##filename = r"C:\Users\ofn77899\Documents\GitHub\CCPi-Reconstruction\python\test\reconstruction_example.mha" -#filename = r"C:\Users\ofn77899\Documents\GitHub\CCPi-Simpleflex\data\head.mha" -# -#reader = vtk.vtkMetaImageReader() -#reader.SetFileName(os.path.normpath(filename)) -#reader.Update() -##vtk returns 3D images, let's take just the one slice there is as 2D -#Im = Converter.vtk2numpy(reader.GetOutput()) -#Im = Im.astype('float32') -##imgplot = plt.imshow(Im) -#perc = 0.05 -#u0 = Im + (perc* np.random.normal(size=np.shape(Im))) -## map the u0 u0->u0>0 -#f = np.frompyfunc(lambda x: 0 if x < 0 else x, 1,1) -#u0 = f(u0).astype('float32') -#converter = Converter.numpy2vtkImporter(u0, reader.GetOutput().GetSpacing(), -# reader.GetOutput().GetOrigin()) -#converter.Update() -#writer = vtk.vtkMetaImageWriter() -#writer.SetInputData(converter.GetOutput()) -#writer.SetFileName(r"C:\Users\ofn77899\Documents\GitHub\CCPi-FISTA_reconstruction\data\noisy_head.mha") -##writer.Write() -# -# -### plot -#fig3D = plt.figure() -##a=fig.add_subplot(3,3,1) -##a.set_title('Original') -##imgplot = plt.imshow(Im) -#sliceNo = 32 -# -#a=fig3D.add_subplot(2,3,1) -#a.set_title('noise') -#imgplot = plt.imshow(u0.T[sliceNo]) -# -#reg_output3d = [] -# -############################################################################### -## Call regularizer -# -######################## SplitBregman_TV ##################################### -## u = SplitBregman_TV(single(u0), 10, 30, 1e-04); -# -##reg = Regularizer(Regularizer.Algorithm.SplitBregman_TV) -# -##out = reg(input=u0, regularization_parameter=10., #number_of_iterations=30, -## #tolerance_constant=1e-4, -## TV_Penalty=Regularizer.TotalVariationPenalty.l1) -# -#out2 = Regularizer.SplitBregman_TV(input=u0, regularization_parameter=10., number_of_iterations=30, -# tolerance_constant=1e-4, -# TV_Penalty=Regularizer.TotalVariationPenalty.l1) -# -# -#pars = out2[-2] -#reg_output3d.append(out2) -# -#a=fig3D.add_subplot(2,3,2) -# -# -#textstr = out2[-1] -# -# -## these are matplotlib.patch.Patch properties -#props = dict(boxstyle='round', facecolor='wheat', alpha=0.5) -## place a text box in upper left in axes coords -#a.text(0.05, 0.95, textstr, transform=a.transAxes, fontsize=14, -# verticalalignment='top', bbox=props) -#imgplot = plt.imshow(reg_output3d[-1][0].T[sliceNo]) -# -####################### FGP_TV ######################################### -## u = FGP_TV(single(u0), 0.05, 100, 1e-04); -#out2 = Regularizer.FGP_TV(input=u0, regularization_parameter=0.005, -# number_of_iterations=200) -#pars = out2[-2] -#reg_output3d.append(out2) -# -#a=fig3D.add_subplot(2,3,2) -# -# -#textstr = out2[-1] -# -# -## these are matplotlib.patch.Patch properties -#props = dict(boxstyle='round', facecolor='wheat', alpha=0.5) -## place a text box in upper left in axes coords -#a.text(0.05, 0.95, textstr, transform=a.transAxes, fontsize=14, -# verticalalignment='top', bbox=props) -#imgplot = plt.imshow(reg_output3d[-1][0].T[sliceNo]) -# -####################### LLT_model ######################################### -## * u0 = Im + .03*randn(size(Im)); % adding noise -## [Den] = LLT_model(single(u0), 10, 0.1, 1); -##Den = LLT_model(single(u0), 25, 0.0003, 300, 0.0001, 0); -##input, regularization_parameter , time_step, number_of_iterations, -## tolerance_constant, restrictive_Z_smoothing=0 -#out2 = Regularizer.LLT_model(input=u0, regularization_parameter=25, -# time_step=0.0003, -# tolerance_constant=0.0001, -# number_of_iterations=300) -#pars = out2[-2] -#reg_output3d.append(out2) -# -#a=fig3D.add_subplot(2,3,2) -# -# -#textstr = out2[-1] -# -# -## these are matplotlib.patch.Patch properties -#props = dict(boxstyle='round', facecolor='wheat', alpha=0.5) -## place a text box in upper left in axes coords -#a.text(0.05, 0.95, textstr, transform=a.transAxes, fontsize=14, -# verticalalignment='top', bbox=props) -#imgplot = plt.imshow(reg_output3d[-1][0].T[sliceNo]) -# -####################### PatchBased_Regul ######################################### -## Quick 2D denoising example in Matlab: -## Im = double(imread('lena_gray_256.tif'))/255; % loading image -## u0 = Im + .03*randn(size(Im)); u0(u0<0) = 0; % adding noise -## ImDen = PB_Regul_CPU(single(u0), 3, 1, 0.08, 0.05); -# -#out2 = Regularizer.PatchBased_Regul(input=u0, regularization_parameter=0.05, -# searching_window_ratio=3, -# similarity_window_ratio=1, -# PB_filtering_parameter=0.08) -#pars = out2[-2] -#reg_output3d.append(out2) -# -#a=fig3D.add_subplot(2,3,2) -# -# -#textstr = out2[-1] -# -# -## these are matplotlib.patch.Patch properties -#props = dict(boxstyle='round', facecolor='wheat', alpha=0.5) -## place a text box in upper left in axes coords -#a.text(0.05, 0.95, textstr, transform=a.transAxes, fontsize=14, -# verticalalignment='top', bbox=props) -#imgplot = plt.imshow(reg_output3d[-1][0].T[sliceNo]) -# - -###################### TGV_PD ######################################### -# Quick 2D denoising example in Matlab: -# Im = double(imread('lena_gray_256.tif'))/255; % loading image -# u0 = Im + .03*randn(size(Im)); u0(u0<0) = 0; % adding noise -# u = PrimalDual_TGV(single(u0), 0.02, 1.3, 1, 550); - - -#out2 = Regularizer.TGV_PD(input=u0, regularization_parameter=0.05, -# first_order_term=1.3, -# second_order_term=1, -# number_of_iterations=550) -#pars = out2[-2] -#reg_output3d.append(out2) -# -#a=fig3D.add_subplot(2,3,2) -# -# -#textstr = out2[-1] -# -# -## these are matplotlib.patch.Patch properties -#props = dict(boxstyle='round', facecolor='wheat', alpha=0.5) -## place a text box in upper left in axes coords -#a.text(0.05, 0.95, textstr, transform=a.transAxes, fontsize=14, -# verticalalignment='top', bbox=props) -#imgplot = plt.imshow(reg_output3d[-1][0].T[sliceNo]) diff --git a/src/Python/test/test_regularizers_3d.py b/src/Python/test/test_regularizers_3d.py deleted file mode 100644 index 2d11a7e..0000000 --- a/src/Python/test/test_regularizers_3d.py +++ /dev/null @@ -1,425 +0,0 @@ -# -*- coding: utf-8 -*- -""" -Created on Fri Aug 4 11:10:05 2017 - -@author: ofn77899 -""" - -#from ccpi.viewer.CILViewer2D import Converter -#import vtk - -import matplotlib.pyplot as plt -import numpy as np -import os -from enum import Enum -import timeit -#from PIL import Image -#from Regularizer import Regularizer -from ccpi.imaging.Regularizer import Regularizer - -############################################################################### -#https://stackoverflow.com/questions/13875989/comparing-image-in-url-to-image-in-filesystem-in-python/13884956#13884956 -#NRMSE a normalization of the root of the mean squared error -#NRMSE is simply 1 - [RMSE / (maxval - minval)]. Where maxval is the maximum -# intensity from the two images being compared, and respectively the same for -# minval. RMSE is given by the square root of MSE: -# sqrt[(sum(A - B) ** 2) / |A|], -# where |A| means the number of elements in A. By doing this, the maximum value -# given by RMSE is maxval. - -def nrmse(im1, im2): - a, b = im1.shape - rmse = np.sqrt(np.sum((im2 - im1) ** 2) / float(a * b)) - max_val = max(np.max(im1), np.max(im2)) - min_val = min(np.min(im1), np.min(im2)) - return 1 - (rmse / (max_val - min_val)) -############################################################################### - -############################################################################### -# -# 2D Regularizers -# -############################################################################### -#Example: -# figure; -# Im = double(imread('lena_gray_256.tif'))/255; % loading image -# u0 = Im + .05*randn(size(Im)); u0(u0 < 0) = 0; -# u = SplitBregman_TV(single(u0), 10, 30, 1e-04); - - -#filename = r"C:\Users\ofn77899\Documents\GitHub\CCPi-FISTA_reconstruction\data\lena_gray_512.tif" -filename = r"/home/ofn77899/Reconstruction/CCPi-FISTA_Reconstruction/data/lena_gray_512.tif" -#filename = r'/home/algol/Documents/Python/STD_test_images/lena_gray_512.tif' - -#reader = vtk.vtkTIFFReader() -#reader.SetFileName(os.path.normpath(filename)) -#reader.Update() -Im = plt.imread(filename) -#Im = Image.open('/home/algol/Documents/Python/STD_test_images/lena_gray_512.tif')/255 -#img.show() -Im = np.asarray(Im, dtype='float32') - -# create a 3D image by stacking N of this images - - -#imgplot = plt.imshow(Im) -perc = 0.05 -u_n = Im + (perc* np.random.normal(size=np.shape(Im))) -y,z = np.shape(u_n) -u_n = np.reshape(u_n , (1,y,z)) - -u0 = u_n.copy() -for i in range (19): - u_n = Im + (perc* np.random.normal(size=np.shape(Im))) - u_n = np.reshape(u_n , (1,y,z)) - - u0 = np.vstack ( (u0, u_n) ) - -# map the u0 u0->u0>0 -f = np.frompyfunc(lambda x: 0 if x < 0 else x, 1,1) -u0 = f(u0).astype('float32') - -print ("Passed image shape {0}".format(np.shape(u0))) - -## plot -fig = plt.figure() -#a=fig.add_subplot(3,3,1) -#a.set_title('Original') -#imgplot = plt.imshow(Im) -sliceno = 10 - -a=fig.add_subplot(2,3,1) -a.set_title('noise') -imgplot = plt.imshow(u0[sliceno],cmap="gray") - -reg_output = [] -############################################################################## -# Call regularizer - -####################### SplitBregman_TV ##################################### -# u = SplitBregman_TV(single(u0), 10, 30, 1e-04); - -use_object = True -if use_object: - reg = Regularizer(Regularizer.Algorithm.SplitBregman_TV) - print (reg.pars) - reg.setParameter(input=u0) - reg.setParameter(regularization_parameter=10.) - # or - # reg.setParameter(input=u0, regularization_parameter=10., #number_of_iterations=30, - #tolerance_constant=1e-4, - #TV_Penalty=Regularizer.TotalVariationPenalty.l1) - plotme = reg() [0] - pars = reg.pars - textstr = reg.printParametersToString() - - #out = reg(input=u0, regularization_parameter=10., #number_of_iterations=30, - #tolerance_constant=1e-4, - # TV_Penalty=Regularizer.TotalVariationPenalty.l1) - -#out2 = Regularizer.SplitBregman_TV(input=u0, regularization_parameter=10., number_of_iterations=30, -# tolerance_constant=1e-4, -# TV_Penalty=Regularizer.TotalVariationPenalty.l1) - -else: - out2 = Regularizer.SplitBregman_TV(input=u0, regularization_parameter=10. ) - pars = out2[2] - reg_output.append(out2) - plotme = reg_output[-1][0] - textstr = out2[-1] - -a=fig.add_subplot(2,3,2) - - -# these are matplotlib.patch.Patch properties -props = dict(boxstyle='round', facecolor='wheat', alpha=0.5) -# place a text box in upper left in axes coords -a.text(0.05, 0.95, textstr, transform=a.transAxes, fontsize=14, - verticalalignment='top', bbox=props) -imgplot = plt.imshow(plotme[sliceno],cmap="gray") - -###################### FGP_TV ######################################### -# u = FGP_TV(single(u0), 0.05, 100, 1e-04); -out2 = Regularizer.FGP_TV(input=u0, regularization_parameter=0.0005, - number_of_iterations=50) -pars = out2[-2] - -reg_output.append(out2) - -a=fig.add_subplot(2,3,3) - -textstr = out2[-1] - -# these are matplotlib.patch.Patch properties -props = dict(boxstyle='round', facecolor='wheat', alpha=0.5) -# place a text box in upper left in axes coords -a.text(0.05, 0.95, textstr, transform=a.transAxes, fontsize=14, - verticalalignment='top', bbox=props) -imgplot = plt.imshow(reg_output[-1][0][sliceno]) -# place a text box in upper left in axes coords -a.text(0.05, 0.95, textstr, transform=a.transAxes, fontsize=14, - verticalalignment='top', bbox=props) -imgplot = plt.imshow(reg_output[-1][0][sliceno],cmap="gray") - -###################### LLT_model ######################################### -# * u0 = Im + .03*randn(size(Im)); % adding noise -# [Den] = LLT_model(single(u0), 10, 0.1, 1); -#Den = LLT_model(single(u0), 25, 0.0003, 300, 0.0001, 0); -#input, regularization_parameter , time_step, number_of_iterations, -# tolerance_constant, restrictive_Z_smoothing=0 -out2 = Regularizer.LLT_model(input=u0, regularization_parameter=25, - time_step=0.0003, - tolerance_constant=0.0001, - number_of_iterations=300) -pars = out2[-2] - -reg_output.append(out2) - -a=fig.add_subplot(2,3,4) - -textstr = out2[-1] - -# these are matplotlib.patch.Patch properties -props = dict(boxstyle='round', facecolor='wheat', alpha=0.5) -# place a text box in upper left in axes coords -a.text(0.05, 0.95, textstr, transform=a.transAxes, fontsize=14, - verticalalignment='top', bbox=props) -imgplot = plt.imshow(reg_output[-1][0][sliceno],cmap="gray") - - -# ###################### PatchBased_Regul ######################################### -# # Quick 2D denoising example in Matlab: -# # Im = double(imread('lena_gray_256.tif'))/255; % loading image -# # u0 = Im + .03*randn(size(Im)); u0(u0<0) = 0; % adding noise -# # ImDen = PB_Regul_CPU(single(u0), 3, 1, 0.08, 0.05); - -out2 = Regularizer.PatchBased_Regul(input=u0, regularization_parameter=0.05, - searching_window_ratio=3, - similarity_window_ratio=1, - PB_filtering_parameter=0.08) -pars = out2[-2] -reg_output.append(out2) - -a=fig.add_subplot(2,3,5) - - -textstr = out2[-1] - -# these are matplotlib.patch.Patch properties -props = dict(boxstyle='round', facecolor='wheat', alpha=0.5) -# place a text box in upper left in axes coords -a.text(0.05, 0.95, textstr, transform=a.transAxes, fontsize=14, - verticalalignment='top', bbox=props) -imgplot = plt.imshow(reg_output[-1][0][sliceno],cmap="gray") - - -# ###################### TGV_PD ######################################### -# # Quick 2D denoising example in Matlab: -# # Im = double(imread('lena_gray_256.tif'))/255; % loading image -# # u0 = Im + .03*randn(size(Im)); u0(u0<0) = 0; % adding noise -# # u = PrimalDual_TGV(single(u0), 0.02, 1.3, 1, 550); - - -out2 = Regularizer.TGV_PD(input=u0, regularization_parameter=0.05, - first_order_term=1.3, - second_order_term=1, - number_of_iterations=550) -pars = out2[-2] -reg_output.append(out2) - -a=fig.add_subplot(2,3,6) - - -textstr = out2[-1] - - -# these are matplotlib.patch.Patch properties -props = dict(boxstyle='round', facecolor='wheat', alpha=0.5) -# place a text box in upper left in axes coords -a.text(0.05, 0.95, textstr, transform=a.transAxes, fontsize=14, - verticalalignment='top', bbox=props) -imgplot = plt.imshow(reg_output[-1][0][sliceno],cmap="gray") - - -plt.show() - -################################################################################ -## -## 3D Regularizers -## -################################################################################ -##Example: -## figure; -## Im = double(imread('lena_gray_256.tif'))/255; % loading image -## u0 = Im + .05*randn(size(Im)); u0(u0 < 0) = 0; -## u = SplitBregman_TV(single(u0), 10, 30, 1e-04); -# -##filename = r"C:\Users\ofn77899\Documents\GitHub\CCPi-Reconstruction\python\test\reconstruction_example.mha" -#filename = r"C:\Users\ofn77899\Documents\GitHub\CCPi-Simpleflex\data\head.mha" -# -#reader = vtk.vtkMetaImageReader() -#reader.SetFileName(os.path.normpath(filename)) -#reader.Update() -##vtk returns 3D images, let's take just the one slice there is as 2D -#Im = Converter.vtk2numpy(reader.GetOutput()) -#Im = Im.astype('float32') -##imgplot = plt.imshow(Im) -#perc = 0.05 -#u0 = Im + (perc* np.random.normal(size=np.shape(Im))) -## map the u0 u0->u0>0 -#f = np.frompyfunc(lambda x: 0 if x < 0 else x, 1,1) -#u0 = f(u0).astype('float32') -#converter = Converter.numpy2vtkImporter(u0, reader.GetOutput().GetSpacing(), -# reader.GetOutput().GetOrigin()) -#converter.Update() -#writer = vtk.vtkMetaImageWriter() -#writer.SetInputData(converter.GetOutput()) -#writer.SetFileName(r"C:\Users\ofn77899\Documents\GitHub\CCPi-FISTA_reconstruction\data\noisy_head.mha") -##writer.Write() -# -# -### plot -#fig3D = plt.figure() -##a=fig.add_subplot(3,3,1) -##a.set_title('Original') -##imgplot = plt.imshow(Im) -#sliceNo = 32 -# -#a=fig3D.add_subplot(2,3,1) -#a.set_title('noise') -#imgplot = plt.imshow(u0.T[sliceNo]) -# -#reg_output3d = [] -# -############################################################################### -## Call regularizer -# -######################## SplitBregman_TV ##################################### -## u = SplitBregman_TV(single(u0), 10, 30, 1e-04); -# -##reg = Regularizer(Regularizer.Algorithm.SplitBregman_TV) -# -##out = reg(input=u0, regularization_parameter=10., #number_of_iterations=30, -## #tolerance_constant=1e-4, -## TV_Penalty=Regularizer.TotalVariationPenalty.l1) -# -#out2 = Regularizer.SplitBregman_TV(input=u0, regularization_parameter=10., number_of_iterations=30, -# tolerance_constant=1e-4, -# TV_Penalty=Regularizer.TotalVariationPenalty.l1) -# -# -#pars = out2[-2] -#reg_output3d.append(out2) -# -#a=fig3D.add_subplot(2,3,2) -# -# -#textstr = out2[-1] -# -# -## these are matplotlib.patch.Patch properties -#props = dict(boxstyle='round', facecolor='wheat', alpha=0.5) -## place a text box in upper left in axes coords -#a.text(0.05, 0.95, textstr, transform=a.transAxes, fontsize=14, -# verticalalignment='top', bbox=props) -#imgplot = plt.imshow(reg_output3d[-1][0].T[sliceNo]) -# -####################### FGP_TV ######################################### -## u = FGP_TV(single(u0), 0.05, 100, 1e-04); -#out2 = Regularizer.FGP_TV(input=u0, regularization_parameter=0.005, -# number_of_iterations=200) -#pars = out2[-2] -#reg_output3d.append(out2) -# -#a=fig3D.add_subplot(2,3,2) -# -# -#textstr = out2[-1] -# -# -## these are matplotlib.patch.Patch properties -#props = dict(boxstyle='round', facecolor='wheat', alpha=0.5) -## place a text box in upper left in axes coords -#a.text(0.05, 0.95, textstr, transform=a.transAxes, fontsize=14, -# verticalalignment='top', bbox=props) -#imgplot = plt.imshow(reg_output3d[-1][0].T[sliceNo]) -# -####################### LLT_model ######################################### -## * u0 = Im + .03*randn(size(Im)); % adding noise -## [Den] = LLT_model(single(u0), 10, 0.1, 1); -##Den = LLT_model(single(u0), 25, 0.0003, 300, 0.0001, 0); -##input, regularization_parameter , time_step, number_of_iterations, -## tolerance_constant, restrictive_Z_smoothing=0 -#out2 = Regularizer.LLT_model(input=u0, regularization_parameter=25, -# time_step=0.0003, -# tolerance_constant=0.0001, -# number_of_iterations=300) -#pars = out2[-2] -#reg_output3d.append(out2) -# -#a=fig3D.add_subplot(2,3,2) -# -# -#textstr = out2[-1] -# -# -## these are matplotlib.patch.Patch properties -#props = dict(boxstyle='round', facecolor='wheat', alpha=0.5) -## place a text box in upper left in axes coords -#a.text(0.05, 0.95, textstr, transform=a.transAxes, fontsize=14, -# verticalalignment='top', bbox=props) -#imgplot = plt.imshow(reg_output3d[-1][0].T[sliceNo]) -# -####################### PatchBased_Regul ######################################### -## Quick 2D denoising example in Matlab: -## Im = double(imread('lena_gray_256.tif'))/255; % loading image -## u0 = Im + .03*randn(size(Im)); u0(u0<0) = 0; % adding noise -## ImDen = PB_Regul_CPU(single(u0), 3, 1, 0.08, 0.05); -# -#out2 = Regularizer.PatchBased_Regul(input=u0, regularization_parameter=0.05, -# searching_window_ratio=3, -# similarity_window_ratio=1, -# PB_filtering_parameter=0.08) -#pars = out2[-2] -#reg_output3d.append(out2) -# -#a=fig3D.add_subplot(2,3,2) -# -# -#textstr = out2[-1] -# -# -## these are matplotlib.patch.Patch properties -#props = dict(boxstyle='round', facecolor='wheat', alpha=0.5) -## place a text box in upper left in axes coords -#a.text(0.05, 0.95, textstr, transform=a.transAxes, fontsize=14, -# verticalalignment='top', bbox=props) -#imgplot = plt.imshow(reg_output3d[-1][0].T[sliceNo]) -# - -###################### TGV_PD ######################################### -# Quick 2D denoising example in Matlab: -# Im = double(imread('lena_gray_256.tif'))/255; % loading image -# u0 = Im + .03*randn(size(Im)); u0(u0<0) = 0; % adding noise -# u = PrimalDual_TGV(single(u0), 0.02, 1.3, 1, 550); - - -#out2 = Regularizer.TGV_PD(input=u0, regularization_parameter=0.05, -# first_order_term=1.3, -# second_order_term=1, -# number_of_iterations=550) -#pars = out2[-2] -#reg_output3d.append(out2) -# -#a=fig3D.add_subplot(2,3,2) -# -# -#textstr = out2[-1] -# -# -## these are matplotlib.patch.Patch properties -#props = dict(boxstyle='round', facecolor='wheat', alpha=0.5) -## place a text box in upper left in axes coords -#a.text(0.05, 0.95, textstr, transform=a.transAxes, fontsize=14, -# verticalalignment='top', bbox=props) -#imgplot = plt.imshow(reg_output3d[-1][0].T[sliceNo]) |