diff options
| -rw-r--r-- | src/Python/CMakeLists.txt | 6 | ||||
| -rw-r--r-- | src/Python/FindAnacondaEnvironment.cmake | 7 | ||||
| -rw-r--r-- | src/Python/ccpi/imaging/Regularizer.py | 6 | ||||
| -rw-r--r-- | src/Python/ccpi/reconstruction/AstraDevice.py | 83 | ||||
| -rw-r--r-- | src/Python/ccpi/reconstruction/DeviceModel.py | 63 | ||||
| -rw-r--r-- | src/Python/ccpi/reconstruction/FISTAReconstructor.py | 119 | ||||
| -rw-r--r-- | src/Python/fista_module.cpp | 5 | ||||
| -rw-r--r-- | src/Python/test/test_reconstructor.py | 100 |
8 files changed, 320 insertions, 69 deletions
diff --git a/src/Python/CMakeLists.txt b/src/Python/CMakeLists.txt index 1b73380..506159a 100644 --- a/src/Python/CMakeLists.txt +++ b/src/Python/CMakeLists.txt @@ -114,6 +114,8 @@ endif() # fista reconstructor file(COPY ${CMAKE_CURRENT_SOURCE_DIR}/ccpi/reconstruction/FISTAReconstructor.py DESTINATION ${CMAKE_CURRENT_BINARY_DIR}/ccpi/reconstruction) file(COPY ${CMAKE_CURRENT_SOURCE_DIR}/ccpi/reconstruction/__init__.py DESTINATION ${CMAKE_CURRENT_BINARY_DIR}/ccpi/reconstruction) +file(COPY ${CMAKE_CURRENT_SOURCE_DIR}/ccpi/reconstruction/DeviceModel.py DESTINATION ${CMAKE_CURRENT_BINARY_DIR}/ccpi/reconstruction) +file(COPY ${CMAKE_CURRENT_SOURCE_DIR}/ccpi/reconstruction/AstraDevice.py DESTINATION ${CMAKE_CURRENT_BINARY_DIR}/ccpi/reconstruction) configure_file(${CMAKE_CURRENT_SOURCE_DIR}/setup-fista.py.in ${CMAKE_CURRENT_BINARY_DIR}/setup-fista.py) file(MAKE_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}/fista-recipe) @@ -158,13 +160,13 @@ add_custom_target(regularizers ) add_custom_target(install-fista - COMMAND conda + COMMAND ${CONDA_EXECUTABLE} install --force --use-local ccpi-fista=${CIL_VERSION} -c ccpi -c conda-forge WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR} ) add_custom_target(install-regularizers - COMMAND conda + COMMAND ${CONDA_EXECUTABLE} install --force --use-local ccpi-regularizers=${CIL_VERSION} -c ccpi -c conda-forge WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR} ) diff --git a/src/Python/FindAnacondaEnvironment.cmake b/src/Python/FindAnacondaEnvironment.cmake index fa4637a..6475128 100644 --- a/src/Python/FindAnacondaEnvironment.cmake +++ b/src/Python/FindAnacondaEnvironment.cmake @@ -97,7 +97,12 @@ function (findPythonForAnacondaEnvironment env) set (PYTHON_VERSION_MINOR ${_PYTHON_VERSION_MINOR} PARENT_SCOPE) set (PYTHON_VERSION_PATCH ${_PYTHON_VERSION_PATCH} PARENT_SCOPE) message("My version found " ${PYTHON_VERSION_STRING}) - + ## find conda executable + if (WIN32) + set (CONDA_EXECUTABLE ${env}/Script/conda PARENT_SCOPE) + elseif(UNIX) + set (CONDA_EXECUTABLE ${env}/bin/conda PARENT_SCOPE) + endif() endfunction() diff --git a/src/Python/ccpi/imaging/Regularizer.py b/src/Python/ccpi/imaging/Regularizer.py index fb9ae08..8ab6c6a 100644 --- a/src/Python/ccpi/imaging/Regularizer.py +++ b/src/Python/ccpi/imaging/Regularizer.py @@ -216,13 +216,13 @@ class Regularizer(): #assuming it's 3D # run independent calls on each slice out3d = input.copy() - for i in range(np.shape(input)[2]): - out = self.algorithm(input, regularization_parameter, + for i in range(np.shape(input)[0]): + out = self.algorithm(input[i], regularization_parameter, self.pars['first_order_term'] , self.pars['second_order_term'] , self.pars['number_of_iterations']) # copy the result in the 3D image - out3d.T[i] = out[0].copy() + out3d[i] = out[0].copy() # append the rest of the info that the algorithm returns output = [out3d] for i in range(1,len(out)): diff --git a/src/Python/ccpi/reconstruction/AstraDevice.py b/src/Python/ccpi/reconstruction/AstraDevice.py new file mode 100644 index 0000000..ac9206e --- /dev/null +++ b/src/Python/ccpi/reconstruction/AstraDevice.py @@ -0,0 +1,83 @@ +import astra +from ccpi.reconstruction.DeviceModel import DeviceModel +import numpy + +class AstraDevice(DeviceModel): + '''Concrete class for Astra Device''' + + def __init__(self, + device_type, + data_aquisition_geometry, + reconstructed_volume_geometry): + + super(AstraDevice, self).__init__(device_type, + data_aquisition_geometry, + reconstructed_volume_geometry) + + self.type = device_type + self.proj_geom = astra.creators.create_proj_geom( + device_type, + self.acquisition_data_geometry['detectorSpacingX'], + self.acquisition_data_geometry['detectorSpacingY'], + self.acquisition_data_geometry['cameraX'], + self.acquisition_data_geometry['cameraY'], + self.acquisition_data_geometry['angles'], + ) + + self.vol_geom = astra.creators.create_vol_geom( + self.reconstructed_volume_geometry['X'], + self.reconstructed_volume_geometry['Y'], + self.reconstructed_volume_geometry['Z'] + ) + + def doForwardProject(self, volume): + '''Forward projects the volume according to the device geometry + +Uses Astra-toolbox +''' + + try: + sino_id, y = astra.creators.create_sino3d_gpu( + volume, self.proj_geom, self.vol_geom) + astra.matlab.data3d('delete', sino_id) + return y + except Exception(e): + print(e) + print("Value Error: ", self.proj_geom, self.vol_geom) + + def doBackwardProject(self, projections): + '''Backward projects the projections according to the device geometry + +Uses Astra-toolbox +''' + idx, volume = \ + astra.creators.create_backprojection3d_gpu( + projections, + self.proj_geom, + self.vol_geom) + + astra.matlab.data3d('delete', idx) + return volume + + def createReducedDevice(self): + proj_geom = [ + self.acquisition_data_geometry['cameraX'], + 1, + self.acquisition_data_geometry['detectorSpacingX'], + self.acquisition_data_geometry['detectorSpacingY'], + self.acquisition_data_geometry['angles'] + ] + + vol_geom = [ + self.reconstructed_volume_geometry['X'], + self.reconstructed_volume_geometry['Y'], + 1 + ] + return AstraDevice(self.type, proj_geom, vol_geom) + + + +if __name__=="main": + a = AstraDevice() + + diff --git a/src/Python/ccpi/reconstruction/DeviceModel.py b/src/Python/ccpi/reconstruction/DeviceModel.py new file mode 100644 index 0000000..eeb9a34 --- /dev/null +++ b/src/Python/ccpi/reconstruction/DeviceModel.py @@ -0,0 +1,63 @@ +from abc import ABCMeta, abstractmethod +from enum import Enum + +class DeviceModel(metaclass=ABCMeta): + '''Abstract class that defines the device for projection and backprojection + +This class defines the methods that must be implemented by concrete classes. + + ''' + + class DeviceType(Enum): + '''Type of device +PARALLEL BEAM +PARALLEL BEAM 3D +CONE BEAM +HELICAL''' + + PARALLEL = 'parallel' + PARALLEL3D = 'parallel3d' + CONE_BEAM = 'cone-beam' + HELICAL = 'helical' + + def __init__(self, + device_type, + data_aquisition_geometry, + reconstructed_volume_geometry): + '''Initializes the class + +Mandatory parameters are: +device_type from DeviceType Enum +data_acquisition_geometry: tuple (camera_X, camera_Y, detectorSpacingX, + detectorSpacingY, angles) +reconstructed_volume_geometry: tuple (dimX,dimY,dimZ) +''' + self.device_geometry = device_type + self.acquisition_data_geometry = { + 'cameraX': data_aquisition_geometry[0], + 'cameraY': data_aquisition_geometry[1], + 'detectorSpacingX' : data_aquisition_geometry[2], + 'detectorSpacingY' : data_aquisition_geometry[3], + 'angles' : data_aquisition_geometry[4],} + self.reconstructed_volume_geometry = { + 'X': reconstructed_volume_geometry[0] , + 'Y': reconstructed_volume_geometry[1] , + 'Z': reconstructed_volume_geometry[2] } + + @abstractmethod + def doForwardProject(self, volume): + '''Forward projects the volume according to the device geometry''' + return NotImplemented + + + @abstractmethod + def doBackwardProject(self, projections): + '''Backward projects the projections according to the device geometry''' + return NotImplemented + + @abstractmethod + def createReducedDevice(self): + '''Create a Device to do forward/backward projections on 2D slices''' + return NotImplemented + + diff --git a/src/Python/ccpi/reconstruction/FISTAReconstructor.py b/src/Python/ccpi/reconstruction/FISTAReconstructor.py index c903712..0607003 100644 --- a/src/Python/ccpi/reconstruction/FISTAReconstructor.py +++ b/src/Python/ccpi/reconstruction/FISTAReconstructor.py @@ -27,6 +27,7 @@ import numpy from enum import Enum import astra +from ccpi.reconstruction.AstraDevice import AstraDevice @@ -73,7 +74,10 @@ class FISTAReconstructor(): # 3. "A novel tomographic reconstruction method based on the robust # Student's t function for suppressing data outliers" D. Kazantsev et.al. # D. Kazantsev, 2016-17 - def __init__(self, projector_geometry, output_geometry, input_sinogram, + def __init__(self, projector_geometry, + output_geometry, + input_sinogram, + device, **kwargs): # handle parmeters: # obligatory parameters @@ -86,7 +90,10 @@ class FISTAReconstructor(): self.pars['number_of_angles'] = nangles self.pars['SlicesZ'] = sliceZ self.pars['output_volume'] = None + self.pars['device_model'] = device + self.use_device = True + print (self.pars) # handle optional input parameters (at instantiation) @@ -113,7 +120,9 @@ class FISTAReconstructor(): 'output_volume', 'os_subsets', 'os_indices', - 'os_bins') + 'os_bins', + 'device_model', + 'reduced_device_model') self.acceptedInputKeywords = list(kw) # handle keyworded parameters @@ -170,9 +179,11 @@ class FISTAReconstructor(): if not 'initialize' in kwargs.keys(): self.pars['initialize'] = False + reduced_device = device.createReducedDevice() + self.setParameter(reduced_device_model=reduced_device) + - - + def setParameter(self, **kwargs): '''set named parameter for the reconstructor engine @@ -436,46 +447,66 @@ class FISTAReconstructor(): X_t = X.copy() # convenience variable storage proj_geom , vol_geom, sino , \ - SlicesZ = self.getParameter([ 'projector_geometry' , + SlicesZ , ring_lambda_R_L1 = self.getParameter([ 'projector_geometry' , 'output_geometry', 'input_sinogram', - 'SlicesZ' ]) + 'SlicesZ' , + 'ring_lambda_R_L1']) t = 1 + + device = self.getParameter('device_model') + reduced_device = self.getParameter('reduced_device_model') for i in range(self.getParameter('number_of_iterations')): + print("iteration", i) 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'] == 'fanflat' or \ - self.getParameter('projector_geometry')['type'] == 'fanflat_vec': + pg = self.getParameter('projector_geometry')['type'] + if pg == 'parallel' or \ + pg == 'fanflat' or \ + pg == '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; - self.sino_updt = numpy.zeros(numpy.shape(sino), dtype=numpy.float) - for kkk in range(SlicesZ): - 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) + + if self.use_device : + self.sino_updt = numpy.zeros(numpy.shape(sino), dtype=numpy.float) + + for kkk in range(SlicesZ): + self.sino_updt[kkk] = \ + reduced_device.doForwardProject( X_t[kkk:kkk+1] ) + else: + proj_geomT = proj_geom.copy() + proj_geomT['DetectorRowCount'] = 1 + vol_geomT = vol_geom.copy() + vol_geomT['GridSliceCount'] = 1; + self.sino_updt = numpy.zeros(numpy.shape(sino), dtype=numpy.float) + for kkk in range(SlicesZ): + 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) 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, self.sino_updt = astra.creators.create_sino3d_gpu( - X_t, proj_geom, vol_geom) + + if self.use_device: + self.sino_updt = device.doForwardProject(X_t) + else: + sino_id, self.sino_updt = astra.creators.create_sino3d_gpu( + X_t, proj_geom, vol_geom) + astra.matlab.data3d('delete', sino_id) ## RING REMOVAL - self.ringRemoval(i) + if ring_lambda_R_L1 != 0: + self.ringRemoval(i) ## Projection/Backprojection Routine self.projectionBackprojection(X, X_t) - astra.matlab.data3d('delete', sino_id) + ## REGULARIZATION X = self.regularize(X) ## Update Loop @@ -523,6 +554,8 @@ class FISTAReconstructor(): 'output_geometry', 'Lipschitz_constant']) + device, reduced_device = self.getParameter(['device_model', + 'reduced_device_model']) if self.getParameter('projector_geometry')['type'] == 'parallel' or \ self.getParameter('projector_geometry')['type'] == 'fanflat' or \ @@ -530,27 +563,39 @@ class FISTAReconstructor(): # 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) + if self.use_device: + proj_geomT = proj_geom.copy() + proj_geomT['DetectorRowCount'] = 1 + vol_geomT = vol_geom.copy() + vol_geomT['GridSliceCount'] = 1; + + 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: + for kkk in range(SliceZ): + x_temp[kkk] = \ + reduced_device.doBackwardProject(residual[kkk:kkk+1]) else: - x_id, x_temp = \ + if self.use_device: + x_id, x_temp = \ astra.creators.create_backprojection3d_gpu( - residual, proj_geom, vol_geom) + residual, proj_geom, vol_geom) + astra.matlab.data3d('delete', x_id) + else: + x_temp = \ + device.doBackwardProject(residual) + 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") diff --git a/src/Python/fista_module.cpp b/src/Python/fista_module.cpp index aca3be0..f3add76 100644 --- a/src/Python/fista_module.cpp +++ b/src/Python/fista_module.cpp @@ -328,7 +328,6 @@ bp::list FGP_TV(np::ndarray input, double d_mu, int iter, double d_epsil, int me else { dim_array[2] = input.shape(2); } - // Parameter handling is be done in Python ///*Handling Matlab input data*/ //if ((nrhs < 2) || (nrhs > 5)) mexErrMsgTxt("At least 2 parameters is required: Image(2D/3D), Regularization parameter. The full list of parameters: Image(2D/3D), Regularization parameter, iterations number, tolerance, penalty type ('iso' or 'l1')"); @@ -513,13 +512,11 @@ bp::list FGP_TV(np::ndarray input, double d_mu, int iter, double d_epsil, int me P3_old = reinterpret_cast<float *>(npP3_old.get_data()); R1 = reinterpret_cast<float *>(npR1.get_data()); R2 = reinterpret_cast<float *>(npR2.get_data()); - R2 = reinterpret_cast<float *>(npR3.get_data()); + R3 = reinterpret_cast<float *>(npR3.get_data()); /* begin iterations */ for (ll = 0; ll<iter; ll++) { - /* computing the gradient of the objective function */ Obj_func3D(A, D, R1, R2, R3, lambda, dimX, dimY, dimZ); - /*Taking a step towards minus of the gradient*/ Grad_func3D(P1, P2, P3, D, R1, R2, R3, lambda, dimX, dimY, dimZ); diff --git a/src/Python/test/test_reconstructor.py b/src/Python/test/test_reconstructor.py index 3342301..f4627d7 100644 --- a/src/Python/test/test_reconstructor.py +++ b/src/Python/test/test_reconstructor.py @@ -12,6 +12,9 @@ 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''' @@ -22,7 +25,23 @@ def RMSE(signal1, signal2): 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, '/') @@ -65,23 +84,23 @@ vol_geom = astra.creators.create_vol_geom( image_size_x, ## First pass the arguments to the FISTAReconstructor and test the ## Lipschitz constant -fistaRecon = FISTAReconstructor(proj_geom, - vol_geom, - Sino3D , - weights=Weights3D) - -print ("Lipschitz Constant {0}".format(fistaRecon.pars['Lipschitz_constant'])) -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) +##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: @@ -294,16 +313,53 @@ if False: ## 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) fistaRecon = FISTAReconstructor(proj_geom, vol_geom, Sino3D , - weights=Weights3D) + device = astradevice, + weights=Weights3D + ) print ("Lipschitz Constant {0}".format(fistaRecon.pars['Lipschitz_constant'])) - fistaRecon.setParameter(number_of_iterations = 12) + fistaRecon.setParameter(number_of_iterations = 3) 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")) - numpy.save("X_out.npy", X) + X = numpy.load("X.npy") +## rd = astradevice.createReducedDevice() +## print ("rd proj_geom" , rd.proj_geom) +## +## +## rd.doForwardProject(X[0:1]) +## proj_geomT = proj_geom.copy() +## for ekey in rd.proj_geom.keys(): +## if ekey == 'ProjectionAngles': +## valrd = numpy.shape(rd.proj_geom[ekey]) +## valg = numpy.shape(proj_geomT[ekey]) +## else: +## valrd = rd.proj_geom[ekey] +## valg = proj_geomT[ekey] +## +## print ("key {0}: RD {1} geomT {2}".format(ekey, valrd, valg)) +## +## +## proj_geomT['DetectorRowCount'] = 1 +## vol_geomT = vol_geom.copy() +## vol_geomT['GridSliceCount'] = 1; +## rd.proj_geom = proj_geomT.copy() +## rd.vol_geom = vol_geomT.copy() +## +## +## +## sino_id, y = astra.creators.create_sino3d_gpu( +## X[0:1], rd.proj_geom, rd.vol_geom) + + X = fistaRecon.iterate(X) + #numpy.save("X_out.npy", X) |