From ec373635a4b3e095cfcc87ae03bd52b05389e5d1 Mon Sep 17 00:00:00 2001 From: Edoardo Pasca Date: Mon, 30 Oct 2017 11:09:53 +0000 Subject: bugfixes for AstraDevice use --- src/Python/ccpi/reconstruction/AstraDevice.py | 36 +++++++++++--------- src/Python/ccpi/reconstruction/DeviceModel.py | 3 +- .../ccpi/reconstruction/FISTAReconstructor.py | 28 ++++++++-------- src/Python/test/test_reconstructor.py | 39 +++++++++++++++++++--- 4 files changed, 71 insertions(+), 35 deletions(-) (limited to 'src/Python') diff --git a/src/Python/ccpi/reconstruction/AstraDevice.py b/src/Python/ccpi/reconstruction/AstraDevice.py index 999533b..ac9206e 100644 --- a/src/Python/ccpi/reconstruction/AstraDevice.py +++ b/src/Python/ccpi/reconstruction/AstraDevice.py @@ -1,5 +1,6 @@ import astra from ccpi.reconstruction.DeviceModel import DeviceModel +import numpy class AstraDevice(DeviceModel): '''Concrete class for Astra Device''' @@ -17,7 +18,7 @@ class AstraDevice(DeviceModel): self.proj_geom = astra.creators.create_proj_geom( device_type, self.acquisition_data_geometry['detectorSpacingX'], - 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'], @@ -34,10 +35,15 @@ class AstraDevice(DeviceModel): Uses Astra-toolbox ''' - sino_id, y = astra.creators.create_sino3d_gpu( - volume, self.proj_geom, self.vol_geom) - astra.matlab.data3d('delete', sino_id) - return y + + 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 @@ -49,27 +55,25 @@ Uses Astra-toolbox projections, self.proj_geom, self.vol_geom) + astra.matlab.data3d('delete', idx) return volume def createReducedDevice(self): - proj_geom = astra.creators.create_proj_geom( - device_type, - self.acquisition_data_geometry['detectorSpacingX'], - self.acquisition_data_geometry['detectorSpacingX'], + proj_geom = [ self.acquisition_data_geometry['cameraX'], 1, - self.acquisition_data_geometry['angles'], - ) + self.acquisition_data_geometry['detectorSpacingX'], + self.acquisition_data_geometry['detectorSpacingY'], + self.acquisition_data_geometry['angles'] + ] - vol_geom = astra.creators.create_vol_geom( + vol_geom = [ self.reconstructed_volume_geometry['X'], self.reconstructed_volume_geometry['Y'], 1 - ) - return AstraDevice(proj_geom['type'] , - proj_geom, - vol_geom) + ] + return AstraDevice(self.type, proj_geom, vol_geom) diff --git a/src/Python/ccpi/reconstruction/DeviceModel.py b/src/Python/ccpi/reconstruction/DeviceModel.py index 6214437..eeb9a34 100644 --- a/src/Python/ccpi/reconstruction/DeviceModel.py +++ b/src/Python/ccpi/reconstruction/DeviceModel.py @@ -28,7 +28,8 @@ HELICAL''' Mandatory parameters are: device_type from DeviceType Enum -data_acquisition_geometry: tuple (camera_X, camera_Y) +data_acquisition_geometry: tuple (camera_X, camera_Y, detectorSpacingX, + detectorSpacingY, angles) reconstructed_volume_geometry: tuple (dimX,dimY,dimZ) ''' self.device_geometry = device_type diff --git a/src/Python/ccpi/reconstruction/FISTAReconstructor.py b/src/Python/ccpi/reconstruction/FISTAReconstructor.py index d2eefc0..0607003 100644 --- a/src/Python/ccpi/reconstruction/FISTAReconstructor.py +++ b/src/Python/ccpi/reconstruction/FISTAReconstructor.py @@ -27,7 +27,7 @@ import numpy from enum import Enum import astra -#from ccpi.reconstruction.AstraDevice import AstraDevice +from ccpi.reconstruction.AstraDevice import AstraDevice @@ -90,11 +90,7 @@ class FISTAReconstructor(): self.pars['number_of_angles'] = nangles self.pars['SlicesZ'] = sliceZ self.pars['output_volume'] = None - self.pars['device'] = device - - - reduced_device = device.createReducedDevice() - self.setParameter(reduced_device_model=reduced_device) + self.pars['device_model'] = device self.use_device = True @@ -183,6 +179,9 @@ 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): @@ -455,22 +454,24 @@ class FISTAReconstructor(): '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'); - if self.use_device: + if self.use_device : self.sino_updt = numpy.zeros(numpy.shape(sino), dtype=numpy.float) for kkk in range(SlicesZ): @@ -491,6 +492,7 @@ 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); + if self.use_device: self.sino_updt = device.doForwardProject(X_t) else: @@ -563,7 +565,7 @@ class FISTAReconstructor(): #sino_updt = zeros(size(sino),'single'); x_temp = numpy.zeros(numpy.shape(X),dtype=numpy.float32) - if use_device: + if self.use_device: proj_geomT = proj_geom.copy() proj_geomT['DetectorRowCount'] = 1 vol_geomT = vol_geom.copy() @@ -581,7 +583,7 @@ class FISTAReconstructor(): x_temp[kkk] = \ reduced_device.doBackwardProject(residual[kkk:kkk+1]) else: - if use_device: + if self.use_device: x_id, x_temp = \ astra.creators.create_backprojection3d_gpu( residual, proj_geom, vol_geom) diff --git a/src/Python/test/test_reconstructor.py b/src/Python/test/test_reconstructor.py index 343b9bb..f4627d7 100644 --- a/src/Python/test/test_reconstructor.py +++ b/src/Python/test/test_reconstructor.py @@ -30,10 +30,10 @@ def createAstraDevice(projector_geometry, output_geometry): '''TODO remove''' device = AstraDevice(DeviceModel.DeviceType.PARALLEL3D.value, - [projector_geometry['DetectorSpacingX'] , - projector_geometry['DetectorSpacingY'] , + [projector_geometry['DetectorRowCount'] , projector_geometry['DetectorColCount'] , - projector_geometry['DetectorRowCount'] , + projector_geometry['DetectorSpacingX'] , + projector_geometry['DetectorSpacingY'] , projector_geometry['ProjectionAngles'] ], [ @@ -332,5 +332,34 @@ else: 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) -- cgit v1.2.3