From 2014650ab9fbf5a7d1c7334fa54ac0b1c5908915 Mon Sep 17 00:00:00 2001 From: Edoardo Pasca Date: Tue, 17 Oct 2017 17:01:12 +0100 Subject: Progress in pythonization --- src/Python/ccpi/fista/FISTAReconstructor.py | 104 +++++++++---- src/Python/test_reconstructor.py | 229 +++++++++++++++++++++++++++- 2 files changed, 298 insertions(+), 35 deletions(-) diff --git a/src/Python/ccpi/fista/FISTAReconstructor.py b/src/Python/ccpi/fista/FISTAReconstructor.py index 8318ea6..87dd2c0 100644 --- a/src/Python/ccpi/fista/FISTAReconstructor.py +++ b/src/Python/ccpi/fista/FISTAReconstructor.py @@ -81,7 +81,7 @@ class FISTAReconstructor(): self.pars['projector_geometry'] = projector_geometry # proj_geom self.pars['output_geometry'] = output_geometry # vol_geom self.pars['input_sinogram'] = input_sinogram # sino - detectors, nangles, sliceZ = numpy.shape(input_sinogram) + sliceZ, nangles, detectors = numpy.shape(input_sinogram) self.pars['detectors'] = detectors self.pars['number_of_angles'] = nangles self.pars['SlicesZ'] = sliceZ @@ -108,7 +108,9 @@ class FISTAReconstructor(): 'regularizer' , 'ring_lambda_R_L1', 'ring_alpha', - 'subsets') + 'subsets', + 'use_studentt_fidelity', + 'studentt') self.acceptedInputKeywords = list(kw) # handle keyworded parameters @@ -143,16 +145,18 @@ class FISTAReconstructor(): else: self.pars['region_of_interest'] = numpy.nonzero( self.pars['ideal_image']>0.0) - + + # the regularizer must be a correctly instantiated object if not 'regularizer' in kwargs.keys() : self.pars['regularizer'] = None - else: - # the regularizer must be a correctly instantiated object - if not 'ring_lambda_R_L1' in kwargs.keys(): - self.pars['ring_lambda_R_L1'] = 0 - if not 'ring_alpha' in kwargs.keys(): - self.pars['ring_alpha'] = 1 + #RING REMOVAL + if not 'ring_lambda_R_L1' in kwargs.keys(): + self.pars['ring_lambda_R_L1'] = 0 + if not 'ring_alpha' in kwargs.keys(): + self.pars['ring_alpha'] = 1 + + # ORDERED SUBSET if not 'subsets' in kwargs.keys(): self.pars['subsets'] = 0 else: @@ -160,6 +164,15 @@ 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): @@ -170,6 +183,8 @@ 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) + @@ -354,10 +369,28 @@ class FISTAReconstructor(): #return subsets angles = self.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 + + + return IndicesReorg def prepareForIteration(self): @@ -368,23 +401,24 @@ class FISTAReconstructor(): detectors, nangles, sliceZ = numpy.shape(self.pars['input_sinogram']) self.r = numpy.zeros((detectors, sliceZ), dtype=numpy.float) # another ring variable - self.rx = self.r.copy() + self.r_x = self.r.copy() self.residual = numpy.zeros(numpy.shape(self.pars['input_sinogram'])) if self.getParameter('Lipschitz_constant') is None: self.pars['Lipschitz_constant'] = \ self.calculateLipschitzConstantWithPowerMethod() + # prepareForIteration def iterate(self, Xin=None): # convenience variable storage proj_geom , vol_geom, sino , \ - SlicesZ = self.getParameter(['projector_geometry' , - 'output_geometry', - 'input_sinogram', - 'SlicesZ']) + SlicesZ = self.getParameter([ 'projector_geometry' , + 'output_geometry', + 'input_sinogram', + 'SlicesZ']) t = 1 if Xin is None: @@ -394,7 +428,8 @@ class FISTAReconstructor(): N = vol_geom['GridColCount'] X = numpy.zeros((N,N,SlicesZ), dtype=numpy.float) else: - X = Xin.copy() + # copy by reference + X = Xin X_t = X.copy() @@ -402,28 +437,31 @@ class FISTAReconstructor(): X_old = X.copy() t_old = t r_old = self.r.copy() - if self.pars['projector_geometry']['type'] == 'parallel' or \ - self.pars['projector_geometry']['type'] == 'parallel3d': + if self.getParameter('projector_geometry')['type'] == 'parallel' or \ + self.getParameter('projector_geometry')['type'] == 'parallel3d': # 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 = 1:SlicesZ - # [sino_id, sino_updt(:,:,kkk)] = - # astra_create_sino3d_cuda(X_t(:,:,kkk), proj_geomT, vol_geomT); - # astra_mex_data3d('delete', sino_id); for kkk in range(SlicesZ): + print (kkk) sino_id, sino_updt[kkk] = \ astra.creators.create_sino3d_gpu( - X_t[kkk], proj_geomT, vol_geomT) - + X_t[kkk:kkk+1], proj_geomT, vol_geomT) + astra.matlab.data3d('delete', sino_id) else: - # for divergent 3D geometry (watch GPU memory overflow in - # Astra < 1.8 - sino_id, y = astra.creators.create_sino3d_gpu(X_t, - proj_geom, - vol_geom) - + # 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( + X_t, proj_geom, vol_geom) + + + ## RING REMOVAL + ## REGULARIZATION diff --git a/src/Python/test_reconstructor.py b/src/Python/test_reconstructor.py index a338d34..f8f6b3c 100644 --- a/src/Python/test_reconstructor.py +++ b/src/Python/test_reconstructor.py @@ -31,7 +31,7 @@ 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 = 3 +Z_slices = 20 det_row_count = Z_slices # next definition is just for consistency of naming det_col_count = size_det @@ -64,5 +64,230 @@ fistaRecon = FISTAReconstructor(proj_geom, 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.setParameter(use_studentt_fidelity= True) -## the calculation of the lipschitz constant should not start by itself +## 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 True: + fistaRecon.prepareForIteration() + print ("Lipschitz Constant {0}".format(fistaRecon.pars['Lipschitz_constant'])) + + + + 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)); + + + 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") +## % 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': + # 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): + print (kkk) + sino_id, sino_updt[kkk] = \ + astra.creators.create_sino3d_gpu( + X_t[kkk:kkk+1], proj_geomT, vol_geomT) + 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( + 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 : + for kkk in range(anglesNumb): + print ("angles {0}".format(kkk)) + 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 + + 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': + x_temp = numpy.zeros(numpy.shape(X),dtype=numpy.float32) + 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], + 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) + + + ## 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('ideal_image') is None: + string = 'Iteration Number {0} | Objective {1} \n' + print (string.format( 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 -- cgit v1.2.3