From c563c660764915865e369d8e1ddb7798feaaca78 Mon Sep 17 00:00:00 2001
From: Edoardo Pasca <edo.paskino@gmail.com>
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(-)
(limited to 'src/Python')
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