From 99685257792688f27502c73f51fc39ad4de84728 Mon Sep 17 00:00:00 2001 From: Willem Jan Palenstijn Date: Tue, 8 Apr 2025 15:39:40 +0200 Subject: [PATCH 1/3] Switch to ASTRA direct_FP3D/BP3D interface --- .../astra/processors/AstraBackProjector3D.py | 85 +++--------------- .../processors/AstraForwardProjector3D.py | 86 +++---------------- 2 files changed, 24 insertions(+), 147 deletions(-) diff --git a/Wrappers/Python/cil/plugins/astra/processors/AstraBackProjector3D.py b/Wrappers/Python/cil/plugins/astra/processors/AstraBackProjector3D.py index a9d61a23a6..45f18ed811 100644 --- a/Wrappers/Python/cil/plugins/astra/processors/AstraBackProjector3D.py +++ b/Wrappers/Python/cil/plugins/astra/processors/AstraBackProjector3D.py @@ -19,7 +19,7 @@ from cil.framework.labels import AcquisitionDimension, ImageDimension from cil.plugins.astra.utilities import convert_geometry_to_astra_vec_3D import astra -from astra import astra_dict, algorithm, data3d +import astra.experimental import numpy as np @@ -44,8 +44,7 @@ def __init__(self, kwargs = { 'volume_geometry' : volume_geometry, 'sinogram_geometry' : sinogram_geometry, - 'proj_geom' : None, - 'vol_geom' : None} + 'projector_id' : None} #DataProcessor.__init__(self, **kwargs) super(AstraBackProjector3D, self).__init__(**kwargs) @@ -53,7 +52,11 @@ def __init__(self, self.set_ImageGeometry(volume_geometry) self.set_AcquisitionGeometry(sinogram_geometry) - self.vol_geom, self.proj_geom = convert_geometry_to_astra_vec_3D(self.volume_geometry, self.sinogram_geometry) + vol_geom, proj_geom = convert_geometry_to_astra_vec_3D(self.volume_geometry, self.sinogram_geometry) + self.projector_id = astra.create_projector('cuda3d', proj_geom, vol_geom) + + def __del__(self): + astra.projector3d.delete(self.projector_id) def check_input(self, dataset): @@ -95,20 +98,15 @@ def process(self, out=None): new_shape_ag = [self.sinogram_geometry.pixel_num_v,self.sinogram_geometry.num_projections,self.sinogram_geometry.pixel_num_h] data_temp = DATA.as_array().reshape(new_shape_ag) + new_shape_ig = [self.volume_geometry.voxel_num_z,self.volume_geometry.voxel_num_y,self.volume_geometry.voxel_num_x] + new_shape_ig = [x if x>0 else 1 for x in new_shape_ig] + if out is None: - rec_id, arr_out = astra.create_backprojection3d_gpu(data_temp, - self.proj_geom, - self.vol_geom) + arr_out = np.zeros(new_shape_ig, dtype=np.float32) else: - new_shape_ig = [self.volume_geometry.voxel_num_z,self.volume_geometry.voxel_num_y,self.volume_geometry.voxel_num_x] - new_shape_ig = [x if x>0 else 1 for x in new_shape_ig] arr_out = out.as_array().reshape(new_shape_ig) - rec_id = astra.data3d.link('-vol', self.vol_geom, arr_out) - self.create_backprojection3d_gpu( data_temp, self.proj_geom, self.vol_geom, False, rec_id) - - # delete the GPU copy - astra.data3d.delete(rec_id) + astra.experimental.direct_BP3D(self.projector_id, arr_out, data_temp) arr_out = np.squeeze(arr_out) @@ -117,62 +115,3 @@ def process(self, out=None): else: out.fill(arr_out) return out - - def create_backprojection3d_gpu(self, data, proj_geom, vol_geom, returnData=True, vol_id=None): - - """ - Call to ASTRA to create a backward projection of an image (3D) - - Parameters - ---------- - - data : numpy.ndarray or int - Image data or ID. - - proj_geom : dict - Projection geometry. - - vol_geom : dict - Volume geometry. - - returnData : bool - If False, only return the ID of the forward projection. - - vol_id : int, default None - ID of the np array linked with astra. - - Returns - ------- - - proj_geom : int or (int, numpy.ndarray) - If ``returnData=False``, returns the ID of the back projection. Otherwise returns a tuple containing the ID of the back projection and the back projection itself. - """ - - if isinstance(data, np.ndarray): - sino_id = data3d.create('-sino', proj_geom, data) - else: - sino_id = data - - if vol_id is None: - vol_id = data3d.create('-vol', vol_geom, 0) - - cfg = astra_dict('BP3D_CUDA') - cfg['ProjectionDataId'] = sino_id - cfg['ReconstructionDataId'] = vol_id - alg_id = algorithm.create(cfg) - algorithm.run(alg_id) - algorithm.delete(alg_id) - - if isinstance(data, np.ndarray): - data3d.delete(sino_id) - - if vol_id is not None: - if returnData: - return vol_id, data3d.get_shared(vol_id) - else: - return vol_id - else: - if returnData: - return vol_id, data3d.get(vol_id) - else: - return vol_id diff --git a/Wrappers/Python/cil/plugins/astra/processors/AstraForwardProjector3D.py b/Wrappers/Python/cil/plugins/astra/processors/AstraForwardProjector3D.py index e03959d7db..9f85472395 100644 --- a/Wrappers/Python/cil/plugins/astra/processors/AstraForwardProjector3D.py +++ b/Wrappers/Python/cil/plugins/astra/processors/AstraForwardProjector3D.py @@ -19,7 +19,7 @@ from cil.framework.labels import ImageDimension, AcquisitionDimension from cil.plugins.astra.utilities import convert_geometry_to_astra_vec_3D import astra -from astra import astra_dict, algorithm, data3d +import astra.experimental import numpy as np class AstraForwardProjector3D(DataProcessor): @@ -43,8 +43,7 @@ def __init__(self, kwargs = { 'volume_geometry' : volume_geometry, 'sinogram_geometry' : sinogram_geometry, - 'proj_geom' : None, - 'vol_geom' : None, + 'projector_id' : None } super(AstraForwardProjector3D, self).__init__(**kwargs) @@ -52,7 +51,12 @@ def __init__(self, self.set_ImageGeometry(volume_geometry) self.set_AcquisitionGeometry(sinogram_geometry) - self.vol_geom, self.proj_geom = convert_geometry_to_astra_vec_3D(self.volume_geometry, self.sinogram_geometry) + vol_geom, proj_geom = convert_geometry_to_astra_vec_3D(self.volume_geometry, self.sinogram_geometry) + self.projector_id = astra.create_projector('cuda3d', proj_geom, vol_geom) + + def __del__(self): + astra.projector3d.delete(self.projector_id) + def check_input(self, dataset): @@ -97,19 +101,14 @@ def process(self, out=None): data_temp = IM.as_array().reshape(new_shape_ig) + new_shape_ag = [self.sinogram_geometry.pixel_num_v,self.sinogram_geometry.num_projections,self.sinogram_geometry.pixel_num_h] + if out is None: - sinogram_id, arr_out = astra.create_sino3d_gpu(data_temp, - self.proj_geom, - self.vol_geom) + arr_out = np.zeros(new_shape_ag, dtype=np.float32) else: - new_shape_ag = [self.sinogram_geometry.pixel_num_v,self.sinogram_geometry.num_projections,self.sinogram_geometry.pixel_num_h] arr_out = out.as_array().reshape(new_shape_ag) - sinogram_id = astra.data3d.link('-sino', self.proj_geom, arr_out) - self.create_sino3d_gpu(data_temp, self.proj_geom, self.vol_geom, False, sino_id=sinogram_id) - - #clear the memory on GPU - astra.data3d.delete(sinogram_id) + astra.experimental.direct_FP3D(self.projector_id, data_temp, arr_out) arr_out = np.squeeze(arr_out) @@ -118,64 +117,3 @@ def process(self, out=None): else: out.fill(arr_out) return out - - def create_sino3d_gpu(self, data, proj_geom, vol_geom, returnData=True, gpuIndex=None, sino_id=None): - """ - Call to ASTRA to create a forward projection of an image (3D) - - Parameters - ---------- - - data : numpy.ndarray or int - Image data or ID. - - proj_geom : dict - Projection geometry. - - vol_geom : dict - Volume geometry. - - returnData : bool - If False, only return the ID of the forward projection. - - gpuIndex : int, optional - Optional GPU index. - - Returns - ------- - - proj_geom : int or (int, numpy.ndarray) - If ``returnData=False``, returns the ID of the forward projection. Otherwise returns a tuple containing the ID of the forward projection and the forward projection itself. - - """ - - - if isinstance(data, np.ndarray): - volume_id = data3d.create('-vol', vol_geom, data) - else: - volume_id = data - if sino_id is None: - sino_id = data3d.create('-sino', proj_geom, 0) - - algString = 'FP3D_CUDA' - cfg = astra_dict(algString) - if not gpuIndex==None: - cfg['option']={'GPUindex':gpuIndex} - cfg['ProjectionDataId'] = sino_id - cfg['VolumeDataId'] = volume_id - alg_id = algorithm.create(cfg) - algorithm.run(alg_id) - algorithm.delete(alg_id) - - if isinstance(data, np.ndarray): - data3d.delete(volume_id) - if sino_id is None: - if returnData: - return sino_id, data3d.get(sino_id) - else: - return sino_id - else: - if returnData: - return sino_id, data3d.get_shared(sino_id) - else: - return sino_id From 51cc36d5b512efa1d790672359b03ed2e667ab0f Mon Sep 17 00:00:00 2001 From: Willem Jan Palenstijn Date: Tue, 8 Apr 2025 17:50:48 +0200 Subject: [PATCH 2/3] Add self to contributors --- NOTICE.txt | 2 ++ 1 file changed, 2 insertions(+) diff --git a/NOTICE.txt b/NOTICE.txt index f8cad6aeda..e4d5f9b819 100644 --- a/NOTICE.txt +++ b/NOTICE.txt @@ -35,6 +35,7 @@ Institutions Key: 12 - Australian e-Health Research, Australia 13 - KU Leuven 14 - Independent Contributor +15 - Universiteit Leiden CIL Developers in date order: Edoardo Pasca (2017 – present) - 1 @@ -72,6 +73,7 @@ Nicholas Whyatt (2024) - 1 Rasmia Kulan (2024) - 1 Emmanuel Ferdman (2025) - 14 Mariam Demir (2025) - 1 +Willem Jan Palenstijn (2025) - 15 CIL Advisory Board: Llion Evans - 9 From ec987cf41316bee4fd4941774a7aa860328a2112 Mon Sep 17 00:00:00 2001 From: Willem Jan Palenstijn Date: Tue, 8 Apr 2025 17:42:58 +0200 Subject: [PATCH 3/3] Use 2d_weighting for ASTRA adjoint when in 2D mode This option requires ASTRA version 2.3. --- .../plugins/astra/processors/AstraBackProjector3D.py | 10 ++++++++-- 1 file changed, 8 insertions(+), 2 deletions(-) diff --git a/Wrappers/Python/cil/plugins/astra/processors/AstraBackProjector3D.py b/Wrappers/Python/cil/plugins/astra/processors/AstraBackProjector3D.py index 45f18ed811..12b322ba98 100644 --- a/Wrappers/Python/cil/plugins/astra/processors/AstraBackProjector3D.py +++ b/Wrappers/Python/cil/plugins/astra/processors/AstraBackProjector3D.py @@ -16,7 +16,7 @@ # Authors: # CIL Developers, listed at: https://github.com/TomographicImaging/CIL/blob/master/NOTICE.txt from cil.framework import DataProcessor, ImageData -from cil.framework.labels import AcquisitionDimension, ImageDimension +from cil.framework.labels import AcquisitionDimension, AcquisitionType, ImageDimension from cil.plugins.astra.utilities import convert_geometry_to_astra_vec_3D import astra import astra.experimental @@ -53,7 +53,13 @@ def __init__(self, self.set_AcquisitionGeometry(sinogram_geometry) vol_geom, proj_geom = convert_geometry_to_astra_vec_3D(self.volume_geometry, self.sinogram_geometry) - self.projector_id = astra.create_projector('cuda3d', proj_geom, vol_geom) + + proj_cfg = astra.astra_dict('cuda3d') + proj_cfg['ProjectionGeometry'] = proj_geom + proj_cfg['VolumeGeometry'] = vol_geom + if AcquisitionType.DIM2 & self.sinogram_geometry.dimension: + proj_cfg['ProjectionKernel'] = '2d_weighting' + self.projector_id = astra.projector3d.create(proj_cfg) def __del__(self): astra.projector3d.delete(self.projector_id)