diff options
Diffstat (limited to 'src')
| -rw-r--r-- | src/Python/ccpi/fista/FISTAReconstructor.py | 182 | ||||
| -rw-r--r-- | src/Python/test_reconstructor-os.py | 379 | ||||
| -rw-r--r-- | src/Python/test_reconstructor.py | 122 |
3 files changed, 594 insertions, 89 deletions
diff --git a/src/Python/ccpi/fista/FISTAReconstructor.py b/src/Python/ccpi/fista/FISTAReconstructor.py index 87dd2c0..fda9cf0 100644 --- a/src/Python/ccpi/fista/FISTAReconstructor.py +++ b/src/Python/ccpi/fista/FISTAReconstructor.py @@ -85,6 +85,7 @@ class FISTAReconstructor(): self.pars['detectors'] = detectors self.pars['number_of_angles'] = nangles self.pars['SlicesZ'] = sliceZ + self.pars['output_volume'] = None print (self.pars) # handle optional input parameters (at instantiation) @@ -109,8 +110,10 @@ class FISTAReconstructor(): 'ring_lambda_R_L1', 'ring_alpha', 'subsets', - 'use_studentt_fidelity', - 'studentt') + 'output_volume', + 'os_subsets', + 'os_indices', + 'os_bins') self.acceptedInputKeywords = list(kw) # handle keyworded parameters @@ -141,10 +144,12 @@ class FISTAReconstructor(): if not 'region_of_interest'in kwargs.keys() : if self.pars['ideal_image'] == None: - pass + self.pars['region_of_interest'] = None else: - self.pars['region_of_interest'] = numpy.nonzero( - self.pars['ideal_image']>0.0) + ## nonzero if the image is larger than m + fsm = numpy.frompyfunc(lambda x,m: 1 if x>m else 0, 2,1) + + self.pars['region_of_interest'] = fsm(self.pars['ideal_image'], 0) # the regularizer must be a correctly instantiated object if not 'regularizer' in kwargs.keys() : @@ -165,14 +170,7 @@ class FISTAReconstructor(): if not 'initialize' in kwargs.keys(): self.pars['initialize'] = False - if not 'use_studentt_fidelity' in kwargs.keys(): - self.setParameter(studentt=False) - else: - print ("studentt {0}".format(kwargs['use_studentt_fidelity'])) - if kwargs['use_studentt_fidelity']: - raise Exception('Not implemented') - - self.setParameter(studentt=kwargs['use_studentt_fidelity']) + def setParameter(self, **kwargs): @@ -183,8 +181,6 @@ class FISTAReconstructor(): ''' for key , value in kwargs.items(): if key in self.acceptedInputKeywords: - if key == 'use_studentt_fidelity': - raise Exception('use_studentt_fidelity Not implemented') self.pars[key] = value else: raise Exception('Wrong parameter {0} for '.format(key) + @@ -389,11 +385,15 @@ class FISTAReconstructor(): counter = counter + binsDiscr[jj] - 1 - - return IndicesReorg + # store the OS in parameters + self.setParameter(os_subsets=subsets, + os_bins=binsDiscr, + os_indices=IndicesReorg) def prepareForIteration(self): + print ("FISTA Reconstructor: prepare for iteration") + self.residual_error = numpy.zeros((self.pars['number_of_iterations'])) self.objective = numpy.zeros((self.pars['number_of_iterations'])) @@ -408,19 +408,17 @@ class FISTAReconstructor(): if self.getParameter('Lipschitz_constant') is None: self.pars['Lipschitz_constant'] = \ self.calculateLipschitzConstantWithPowerMethod() + # errors vector (if the ground truth is given) + self.Resid_error = numpy.zeros((self.getParameter('number_of_iterations'))); + # objective function values vector + self.objective = numpy.zeros((self.getParameter('number_of_iterations'))); # prepareForIteration def iterate(self, Xin=None): - # convenience variable storage - proj_geom , vol_geom, sino , \ - SlicesZ = self.getParameter([ 'projector_geometry' , - 'output_geometry', - 'input_sinogram', - 'SlicesZ']) - - t = 1 + print ("FISTA Reconstructor: iterate") + if Xin is None: if self.getParameter('initialize'): X = self.initialize() @@ -430,15 +428,25 @@ class FISTAReconstructor(): else: # copy by reference X = Xin - + # store the output volume in the parameters + self.setParameter(output_volume=X) X_t = X.copy() + # convenience variable storage + proj_geom , vol_geom, sino , \ + SlicesZ = self.getParameter([ 'projector_geometry' , + 'output_geometry', + 'input_sinogram', + 'SlicesZ' ]) + + t = 1 for i in range(self.getParameter('number_of_iterations')): X_old = X.copy() t_old = t r_old = self.r.copy() if self.getParameter('projector_geometry')['type'] == 'parallel' or \ - self.getParameter('projector_geometry')['type'] == 'parallel3d': + self.getParameter('projector_geometry')['type'] == 'fanflat' or \ + self.getParameter('projector_geometry')['type'] == 'fanflat_vec': # if the geometry is parallel use slice-by-slice # projection-backprojection routine #sino_updt = zeros(size(sino),'single'); @@ -446,10 +454,9 @@ class FISTAReconstructor(): proj_geomT['DetectorRowCount'] = 1 vol_geomT = vol_geom.copy() vol_geomT['GridSliceCount'] = 1; - sino_updt = numpy.zeros(numpy.shape(sino), dtype=numpy.float) + self.sino_updt = numpy.zeros(numpy.shape(sino), dtype=numpy.float) for kkk in range(SlicesZ): - print (kkk) - sino_id, sino_updt[kkk] = \ + sino_id, self.sino_updt[kkk] = \ astra.creators.create_sino3d_gpu( X_t[kkk:kkk+1], proj_geomT, vol_geomT) astra.matlab.data3d('delete', sino_id) @@ -457,11 +464,122 @@ class FISTAReconstructor(): # 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.matlab.create_sino3d_gpu( + sino_id, self.sino_updt = astra.creators.create_sino3d_gpu( X_t, proj_geom, vol_geom) ## RING REMOVAL - + self.ringRemoval(i) + ## Projection/Backprojection Routine + self.projectionBackprojection(X, X_t) + astra.matlab.data3d('delete', sino_id) ## REGULARIZATION + X = self.regularize(X) + ## Update Loop + X , X_t, t = self.updateLoop(i, X, X_old, r_old, t, t_old) + self.setParameter(output_volume=X) + return X + ## iterate + + def ringRemoval(self, i): + print ("FISTA Reconstructor: ring removal") + residual = self.residual + lambdaR_L1 , alpha_ring , weights , L_const , sino= \ + self.getParameter(['ring_lambda_R_L1', + 'ring_alpha' , 'weights', + 'Lipschitz_constant', + 'input_sinogram']) + r_x = self.r_x + sino_updt = self.sino_updt + + SlicesZ, anglesNumb, Detectors = \ + numpy.shape(self.getParameter('input_sinogram')) + if lambdaR_L1 > 0 : + 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() + self.r = (r_x - (1./L_const) * vec).copy() + self.objective[i] = (0.5 * (residual ** 2).sum()) + + def projectionBackprojection(self, X, X_t): + print ("FISTA Reconstructor: projection-backprojection routine") + + # a few useful variables + SlicesZ, anglesNumb, Detectors = \ + numpy.shape(self.getParameter('input_sinogram')) + residual = self.residual + proj_geom , vol_geom , L_const = \ + self.getParameter(['projector_geometry' , + 'output_geometry', + 'Lipschitz_constant']) + + + if self.getParameter('projector_geometry')['type'] == 'parallel' or \ + self.getParameter('projector_geometry')['type'] == 'fanflat' or \ + self.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; + 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( + 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) + + def regularize(self, X): + print ("FISTA Reconstructor: regularize") + + regularizer = self.getParameter('regularizer') + if regularizer is not None: + return regularizer(input=X) + else: + return X + + def updateLoop(self, i, X, X_old, r_old, t, t_old): + print ("FISTA Reconstructor: update loop") + lambdaR_L1 = self.getParameter('ring_lambda_R_L1') + if lambdaR_L1 > 0: + self.r = numpy.max( + numpy.abs(self.r) - lambdaR_L1 , 0) * \ + numpy.sign(self.r) + t = (1 + numpy.sqrt(1 + 4 * t**2))/2 + X_t = X + (((t_old -1)/t) * (X - X_old)) + + if lambdaR_L1 > 0: + self.r_x = self.r + \ + (((t_old-1)/t) * (self.r - r_old)) + + if self.getParameter('region_of_interest') is None: + string = 'Iteration Number {0} | Objective {1} \n' + print (string.format( i, self.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], self.objective[i])) + return (X , X_t, t) diff --git a/src/Python/test_reconstructor-os.py b/src/Python/test_reconstructor-os.py new file mode 100644 index 0000000..6f3721f --- /dev/null +++ b/src/Python/test_reconstructor-os.py @@ -0,0 +1,379 @@ +# -*- coding: utf-8 -*- +""" +Created on Wed Aug 23 16:34:49 2017 + +@author: ofn77899 +Based on DemoRD2.m +""" + +import h5py +import numpy + +from ccpi.fista.FISTAReconstructor import FISTAReconstructor +import astra +import matplotlib.pyplot as plt + +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 + +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) + +## Ordered subset +if True: + 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 True: + 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 + + ## additional for + proj_geomSUB = proj_geom.copy() + fistaRecon.residual2 = numpy.zeros(numpy.shape(self.pars['input_sinogram'])) + 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 + + ## https://github.com/vais-ral/CCPi-FISTA_Reconstruction/issues/4 + if (lambdaR_L1 > 0) : + sino_id2, sino_updt2 = astra.creators.create_sino3d_gpu( + X, proj_geom, vol_geom) + astra.matlab.data3d('delete', sino_id2) + + # subset loop + counterInd = 1 + for ss in range(fistaRecon.getParameter('subsets')): + print ("Subset {0}".format(ss)) + X_old = X.copy() + t_old = t + r_old = fistaRecon.r.copy() + + # the number of projections per subset + numProjSub = fistaRecon.getParameter('os_bins')[ss] + CurrSubIndices = fistaRecon.getParameter('os_indices')\ + [counterInd:counterInd+numProjSub-1] + proj_geomSUB['ProjectionAngles'] = angles[CurrSubIndeces] + +## 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 + residual2 = fistaRecon.residual2 + + 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") +## % the ring removal part (Group-Huber fidelity) +## % first 2 iterations do additional work reconstructing whole dataset to ensure +## % the stablility +## if (i < 3) +## [sino_id2, sino_updt2] = astra_create_sino3d_cuda(X_t, proj_geom, vol_geom); +## astra_mex_data3d('delete', sino_id2); +## else +## [sino_id, sino_updt] = astra_create_sino3d_cuda(X_t, proj_geomSUB, vol_geom); +## end + +## https://github.com/vais-ral/CCPi-FISTA_Reconstruction/issues/4 + if i < 3: + pass + else: + sino_id, sino_updt = astra.creators.create_sino3d_gpu( + X_t, proj_geomSUB, vol_geom) +## sino_id, sino_updt = astra.creators.create_sino3d_gpu( +## X, proj_geom, vol_geom) +## astra.matlab.data3d('delete', sino_id) + + for kkk in range(anglesNumb): + + residual2[:,kkk,:] = (weights[:,kkk,:]).squeeze() * \ + ((sino_updt2[:,kkk,:]).squeeze() - \ + (sino[:,kkk,:]).squeeze() -\ + (alpha_ring * r_x) + ) + shape = list(numpy.shape(fistaRecon.getParameter('input_sinogram'))) + shape[1] = numProjSub + fistaRecon.residual = numpy.zeros(shape) + if fistaRecon.residual.__hash__() != residual.__hash__(): + residual = fistaRecon.residual +## for kkk = 1:numProjSub +## indC = CurrSubIndeces(kkk); +## if (i < 3) +## residual(:,kkk,:) = squeeze(residual2(:,indC,:)); +## else +## residual(:,kkk,:) = squeeze(weights(:,indC,:)).*(squeeze(sino_updt(:,kkk,:)) - (squeeze(sino(:,indC,:)) - alpha_ring.*r_x)); +## end +## end + for kk in range(numProjSub): + indC = fistaRecon.getParameter('os_indices')[kkk] + if i < 3: + residual[:,kkk,:] = residual2[:,indC,:].squeeze() + else: + residual(:,kkk,:) = \ + weights[:,indC,:].squeeze() * sino_updt[:,kkk,:].squeeze() - \ + sino[:,indC,:].squeeze() - alpha_ring * fistaRecon.r_x + #squeeze(weights(:,indC,:)).* \ + # (squeeze(sino_updt(:,kkk,:)) - \ + #(squeeze(sino(:,indC,:)) - 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: + 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) + fistaRecon.prepareForIteration() + X = fistaRecon.iterate(numpy.load("X.npy")) diff --git a/src/Python/test_reconstructor.py b/src/Python/test_reconstructor.py index f8f6b3c..07668ba 100644 --- a/src/Python/test_reconstructor.py +++ b/src/Python/test_reconstructor.py @@ -11,10 +11,17 @@ import numpy from ccpi.fista.FISTAReconstructor import FISTAReconstructor import astra +import matplotlib.pyplot as plt -##def getEntry(nx, location): -## for item in nx[location].keys(): -## print (item) +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") @@ -68,7 +75,6 @@ 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) -#fistaRecon.setParameter(use_studentt_fidelity= True) ## Ordered subset if False: @@ -94,19 +100,34 @@ if False: counter = counter + binsDiscr[jj] - 1 -if True: - fistaRecon.prepareForIteration() +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) + #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)); @@ -114,30 +135,18 @@ if True: objective = numpy.zeros((iterFISTA)); - print ("line") t = 1 - print ("line") - 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 ("X_t copy") + + 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'] == 'parallel3d': + 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'); @@ -147,16 +156,15 @@ if True: vol_geomT['GridSliceCount'] = 1; sino_updt = numpy.zeros(numpy.shape(sino), dtype=numpy.float) for kkk in range(SlicesZ): - print (kkk) sino_id, sino_updt[kkk] = \ astra.creators.create_sino3d_gpu( - X_t[kkk:kkk+1], proj_geomT, vol_geomT) + 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.matlab.create_sino3d_gpu( + sino_id, sino_updt = astra.creators.create_sino3d_gpu( X_t, proj_geom, vol_geom) ## RING REMOVAL @@ -169,8 +177,9 @@ if True: SlicesZ, anglesNumb, Detectors = \ numpy.shape(fistaRecon.getParameter('input_sinogram')) if lambdaR_L1 > 0 : + print ("ring removal") for kkk in range(anglesNumb): - print ("angles {0}".format(kkk)) + residual[:,kkk,:] = (weights[:,kkk,:]).squeeze() * \ ((sino_updt[:,kkk,:]).squeeze() - \ (sino[:,kkk,:]).squeeze() -\ @@ -194,39 +203,16 @@ if True: ## r = r_x - (1./L_const).*vec; ## objective(i) = (0.5*sum(residual(:).^2)); % for the objective function output - else: - if fistaRecon.getParameter('use_studentt_fidelity'): - residual = weights * (sino_updt - sino) - for kkk in range(SlicesZ): - # reshape(residual(:,:,kkk), Detectors*anglesNumb, 1) - # 1D - res_vec = numpy.reshape(residual[kkk], (Detectors * anglesNumb,1)) - -## else -## if (studentt == 1) -## % artifacts removal with Students t penalty -## residual = weights.*(sino_updt - sino); -## for kkk = 1:SlicesZ -## res_vec = reshape(residual(:,:,kkk), Detectors*anglesNumb, 1); % 1D vectorized sinogram -## %s = 100; -## %gr = (2)*res_vec./(s*2 + conj(res_vec).*res_vec); -## [ff, gr] = studentst(res_vec, 1); -## residual(:,:,kkk) = reshape(gr, Detectors, anglesNumb); -## end -## objective(i) = ff; % for the objective function output -## else -## % no ring removal (LS model) -## residual = weights.*(sino_updt - sino); -## objective(i) = (0.5*sum(residual(:).^2)); % for the objective function output -## end -## end + # Projection/Backprojection Routine if fistaRecon.getParameter('projector_geometry')['type'] == 'parallel' or \ - fistaRecon.getParameter('projector_geometry')['type'] == 'parallel3d': + 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): - print ("Projection/Backprojection Routine {0}".format( kkk )) + x_id, x_temp[kkk] = \ astra.creators.create_backprojection3d_gpu( residual[kkk:kkk+1], @@ -248,9 +234,11 @@ if True: # 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( @@ -263,9 +251,16 @@ if True: fistaRecon.r_x = fistaRecon.r + \ (((t_old-1)/t) * (fistaRecon.r - r_old)) - if fistaRecon.getParameter('ideal_image') is None: + 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 @@ -291,3 +286,16 @@ if True: ## else ## fprintf('%s %i %s %s %f \n', 'Iteration Number:', i, '|', 'Objective:', objective(i)); ## end +else: + 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) + fistaRecon.prepareForIteration() + X = fistaRecon.iterate(numpy.load("X.npy")) |