From abaa1f0be4c3f7042260dba0ac23dbc53e7d9c46 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Julius=20H=C3=A4ger?= Date: Tue, 5 May 2026 13:10:37 +0200 Subject: [PATCH 1/7] Added support for ASTRA cylindrical detector geometry. --- .../applications/tomo/backends/astra_cuda.py | 6 +- .../applications/tomo/backends/astra_setup.py | 94 ++++++++++++++++++- .../applications/tomo/geometry/conebeam.py | 1 + 3 files changed, 99 insertions(+), 2 deletions(-) diff --git a/src/odl/applications/tomo/backends/astra_cuda.py b/src/odl/applications/tomo/backends/astra_cuda.py index 8477906adbd..b713f4013a2 100644 --- a/src/odl/applications/tomo/backends/astra_cuda.py +++ b/src/odl/applications/tomo/backends/astra_cuda.py @@ -113,8 +113,12 @@ def __init__(self, geometry, vol_space, proj_space): if self.geometry.ndim == 3: if vol_space.impl == 'numpy': - self.transpose_tuple = (1,0,2) + self.transpose_tuple = (1,0,2) if self.geometry.det_curvature_radius is None else (2, 0, 1) elif vol_space.impl == 'pytorch': + # FIXME: if self.geometry.det_curvature_radius is None + # We can't use a single PyTorch transpose... + if self.geometry.det_curvature_radius is not None: + raise NotImplementedError("Curved detectors currently do not support pytorch") self.transpose_tuple = (1,0) else: raise NotImplementedError("Not implemented for another backend") diff --git a/src/odl/applications/tomo/backends/astra_setup.py b/src/odl/applications/tomo/backends/astra_setup.py index c5d3c194dc6..fb86f8c9765 100644 --- a/src/odl/applications/tomo/backends/astra_setup.py +++ b/src/odl/applications/tomo/backends/astra_setup.py @@ -30,7 +30,7 @@ from odl.core.discr import DiscretizedSpace, DiscretizedSpaceElement from odl.applications.tomo.geometry import ( - DivergentBeamGeometry, Flat1dDetector, Flat2dDetector, Geometry, + ConeBeamGeometry, CylindricalDetector, DivergentBeamGeometry, Flat1dDetector, Flat2dDetector, Geometry, ParallelBeamGeometry) from odl.applications.tomo.util.utility import euler_matrix from odl.core.array_API_support import get_array_and_backend @@ -124,6 +124,10 @@ # next release after 1.8.3, see # https://github.com/astra-toolbox/astra-toolbox/pull/183 'par2d_distance_driven_proj': '>1.8.3', + + # Cylidrical detector geometry, see + # https://github.com/astra-toolbox/astra-toolbox/pull/444 + 'cyl_cone_vec': ">= 2.4.0" } ODL_TO_ASTRA_INDEX_PERMUTATIONS = [ @@ -354,6 +358,75 @@ def astra_conebeam_3d_geom_to_vec(geometry:Geometry): return vectors +def astra_cyl_conebeam_3d_geom_to_vec(geometry:DivergentBeamGeometry): + """Create vectors for ASTRA projection geometries from ODL geometry. + + The 3D vectors are used to create an ASTRA projection geometry for + cone beam geometries with a cylindrical detector, see ``'cyl_cone_vec'`` + in the `ASTRA projection geometry documentation`_. + + Each row of the returned vectors corresponds to a single projection + and consists of :: + + (srcX, srcY, srcZ, dX, dY, dZ, uX, uY, uZ, vX, vY, vZ, R) + + with + + - ``src``: the ray source position + - ``d`` : the center of the detector + - ``u`` : tangential direction at center of detector; + the length of u is the arc length of a detector pixel + - ``v`` : the vector from detector pixel ``(0,0)`` to ``(1,0)`` + - ``R`` : the radius of the detector cylinder + + Parameters + ---------- + geometry : `Geometry` + ODL projection geometry from which to create the ASTRA geometry. + + Returns + ------- + vectors : `numpy.ndarray` + Array of shape ``(num_angles, 13)`` containing the vectors. + + References + ---------- + .. _ASTRA projection geometry documentation: + http://www.astra-toolbox.com/docs/geom3d.html#projection-geometries + """ + angles = geometry.angles + vectors = np.zeros((angles.size, 13)) + + # Source position + vectors[:, 0:3] = geometry.src_position(angles) + + # Center of detector in 3D space + # FIXME: This is not correct: det_point_position returns the zero-point of + # the detector, and not the center of the detector. For quarter-pixel-shifted + # detector these two do not coincide. + mid_pt = geometry.det_params.mid_pt + vectors[:, 3:6] = geometry.det_point_position(angles, mid_pt) + + # `det_axes` gives shape (N, 2, 3), swap to get (2, N, 3) + det_axes = np.moveaxis(geometry.det_axes(angles), -2, 0) + px_sizes = geometry.det_partition.cell_sides + + # `px_sizes[0]` is angular partition; scale by radius to get arc length + # NB: For flat panel detector we swap the u and v axes to get a better + # memory layout. For cylindrical detectors this is (currently) not possible + # since both ODL and Astra have the v direction along the axial direction. + vectors[:, 6:9] = det_axes[0] * px_sizes[0] * geometry.det_curvature_radius + vectors[:, 9:12] = det_axes[1] * px_sizes[1] + + # detector curvature radius + vectors[:, 12] = geometry.det_curvature_radius + + # ASTRA has (z, y, x) axis convention, in contrast to (x, y, z) in ODL, + # so we need to adapt to this by changing the order. + vectors = vectors[:, ODL_TO_ASTRA_INDEX_PERMUTATIONS] + + return vectors + def astra_fanflat_2d_geom_to_conebeam_vec(geometry:Geometry): """ Create vectors for ASTRA projection geometry. This is required for the CUDA implementation of fanflat geometry. @@ -609,6 +682,23 @@ def astra_projection_geometry(geometry: Geometry, astra_impl: str): vec = astra_conebeam_3d_geom_to_vec(geometry) proj_geom = astra.create_proj_geom('cone_vec', det_row_count, det_col_count, vec) + + elif (isinstance(geometry, DivergentBeamGeometry) and + isinstance(geometry.detector, CylindricalDetector) and + geometry.ndim == 3): + + if not astra_supports('cyl_cone_vec'): + req_ver = astra_versions_supporting('cyl_cone_vec') + raise NotImplementedError( + f"support for cylindrical detector geometry requires ASTRA {req_ver}" + ) + # Do NOT swap detector axes (see astra_cyl_conebeam_3d_geom_to_vec) + det_row_count = geometry.det_partition.shape[1] + det_col_count = geometry.det_partition.shape[0] + vec = astra_cyl_conebeam_3d_geom_to_vec(geometry) + proj_geom = astra.create_proj_geom('cyl_cone_vec', det_row_count, + det_col_count, vec) + else: raise NotImplementedError(f"unknown ASTRA geometry type {geometry}") @@ -750,6 +840,8 @@ def astra_projector( valid_proj_types = ['linear3d', 'cuda3d'] elif astra_geom in {'cone', 'cone_vec'}: valid_proj_types = ['linearcone', 'cuda3d'] + elif astra_geom in {'cyl_cone_vec'}: + valid_proj_types = ['cuda3d'] else: raise ValueError(f"invalid geometry type {astra_geom}") diff --git a/src/odl/applications/tomo/geometry/conebeam.py b/src/odl/applications/tomo/geometry/conebeam.py index 7ecebb858c7..1865cad2fb3 100644 --- a/src/odl/applications/tomo/geometry/conebeam.py +++ b/src/odl/applications/tomo/geometry/conebeam.py @@ -1471,6 +1471,7 @@ def __repr__(self): posargs = [self.motion_partition, self.det_partition] optargs = [('src_radius', self.src_radius, -1), ('det_radius', self.det_radius, -1), + ('det_curvature_radius', self.det_curvature_radius, None), ('pitch', self.pitch, 0) ] From c5f83c0fc59a6b446da5f92f19adea9170e60251 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Julius=20H=C3=A4ger?= Date: Tue, 5 May 2026 13:36:16 +0200 Subject: [PATCH 2/7] Fix permutation issue. --- src/odl/applications/tomo/backends/astra_setup.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/odl/applications/tomo/backends/astra_setup.py b/src/odl/applications/tomo/backends/astra_setup.py index fb86f8c9765..0bb0093d9a8 100644 --- a/src/odl/applications/tomo/backends/astra_setup.py +++ b/src/odl/applications/tomo/backends/astra_setup.py @@ -423,7 +423,7 @@ def astra_cyl_conebeam_3d_geom_to_vec(geometry:DivergentBeamGeometry): # ASTRA has (z, y, x) axis convention, in contrast to (x, y, z) in ODL, # so we need to adapt to this by changing the order. - vectors = vectors[:, ODL_TO_ASTRA_INDEX_PERMUTATIONS] + vectors = vectors[:, [*ODL_TO_ASTRA_INDEX_PERMUTATIONS, 12]] return vectors From 1e76c2111d84403a4ae168ff7c380c69f04e1b69 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Julius=20H=C3=A4ger?= Date: Tue, 5 May 2026 15:30:10 +0200 Subject: [PATCH 3/7] Switch to permute_dims. --- .../applications/tomo/backends/astra_cuda.py | 17 ++++------------- 1 file changed, 4 insertions(+), 13 deletions(-) diff --git a/src/odl/applications/tomo/backends/astra_cuda.py b/src/odl/applications/tomo/backends/astra_cuda.py index b713f4013a2..fb595f272fc 100644 --- a/src/odl/applications/tomo/backends/astra_cuda.py +++ b/src/odl/applications/tomo/backends/astra_cuda.py @@ -112,16 +112,7 @@ def __init__(self, geometry, vol_space, proj_space): ), f"Volume space ({vol_space.impl}) != Projection space ({proj_space.impl})" if self.geometry.ndim == 3: - if vol_space.impl == 'numpy': - self.transpose_tuple = (1,0,2) if self.geometry.det_curvature_radius is None else (2, 0, 1) - elif vol_space.impl == 'pytorch': - # FIXME: if self.geometry.det_curvature_radius is None - # We can't use a single PyTorch transpose... - if self.geometry.det_curvature_radius is not None: - raise NotImplementedError("Curved detectors currently do not support pytorch") - self.transpose_tuple = (1,0) - else: - raise NotImplementedError("Not implemented for another backend") + self.transpose_tuple = (1,0,2) if self.geometry.det_curvature_radius is None else (2, 0, 1) self.fp_scaling_factor = astra_cuda_fp_scaling_factor(self.geometry) self.bp_scaling_factor = astra_cuda_bp_scaling_factor( @@ -217,7 +208,7 @@ def _call_forward_real(self, vol_data:DiscretizedSpaceElement, out=None, **kwarg ) proj_data = out.data[None] if self.proj_ndim == 2 else out.data if self.geometry.ndim == 3: - proj_data = proj_data.transpose(*self.transpose_tuple) + proj_data = proj_data.__array_namespace__().permute_dims(proj_data, self.transpose_tuple) else: proj_data = empty( @@ -253,7 +244,7 @@ def _call_forward_real(self, vol_data:DiscretizedSpaceElement, out=None, **kwarg proj_data = ( proj_data[0] if self.geometry.ndim == 2 - else proj_data.transpose(*self.transpose_tuple) + else proj_data.__array_namespace__().permute_dims(proj_data, self.transpose_tuple) ) if out is not None: @@ -327,7 +318,7 @@ def _call_backward_real(self, proj_data:DiscretizedSpaceElement, out=None, **kwa if self.proj_ndim == 2: proj_data = proj_data.data[None] elif self.proj_ndim == 3: - proj_data = proj_data.data.transpose(*self.transpose_tuple) + proj_data = proj_data.data.__array_namespace__().permute_dims(proj_data, self.transpose_tuple) else: raise NotImplementedError From da505f5f4c51f578aecb47604a0bd18bac77688a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Julius=20H=C3=A4ger?= Date: Tue, 5 May 2026 15:36:05 +0200 Subject: [PATCH 4/7] Fix typo. --- src/odl/applications/tomo/backends/astra_cuda.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/odl/applications/tomo/backends/astra_cuda.py b/src/odl/applications/tomo/backends/astra_cuda.py index fb595f272fc..99a609ba9f9 100644 --- a/src/odl/applications/tomo/backends/astra_cuda.py +++ b/src/odl/applications/tomo/backends/astra_cuda.py @@ -318,7 +318,7 @@ def _call_backward_real(self, proj_data:DiscretizedSpaceElement, out=None, **kwa if self.proj_ndim == 2: proj_data = proj_data.data[None] elif self.proj_ndim == 3: - proj_data = proj_data.data.__array_namespace__().permute_dims(proj_data, self.transpose_tuple) + proj_data = proj_data.data.__array_namespace__().permute_dims(proj_data.data, self.transpose_tuple) else: raise NotImplementedError From 4a29b2deb53ca07bdba8ac03caae8cc00fa662ed Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Julius=20H=C3=A4ger?= Date: Tue, 5 May 2026 16:33:15 +0200 Subject: [PATCH 5/7] Use proj_space.array_namespace. --- src/odl/applications/tomo/backends/astra_cuda.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/odl/applications/tomo/backends/astra_cuda.py b/src/odl/applications/tomo/backends/astra_cuda.py index 99a609ba9f9..37c6eb3081a 100644 --- a/src/odl/applications/tomo/backends/astra_cuda.py +++ b/src/odl/applications/tomo/backends/astra_cuda.py @@ -208,7 +208,7 @@ def _call_forward_real(self, vol_data:DiscretizedSpaceElement, out=None, **kwarg ) proj_data = out.data[None] if self.proj_ndim == 2 else out.data if self.geometry.ndim == 3: - proj_data = proj_data.__array_namespace__().permute_dims(proj_data, self.transpose_tuple) + proj_data = self._proj_space.array_namespace.permute_dims(proj_data, self.transpose_tuple) else: proj_data = empty( @@ -244,7 +244,7 @@ def _call_forward_real(self, vol_data:DiscretizedSpaceElement, out=None, **kwarg proj_data = ( proj_data[0] if self.geometry.ndim == 2 - else proj_data.__array_namespace__().permute_dims(proj_data, self.transpose_tuple) + else self._proj_space.array_namespace.permute_dims(proj_data, self.transpose_tuple) ) if out is not None: @@ -318,7 +318,7 @@ def _call_backward_real(self, proj_data:DiscretizedSpaceElement, out=None, **kwa if self.proj_ndim == 2: proj_data = proj_data.data[None] elif self.proj_ndim == 3: - proj_data = proj_data.data.__array_namespace__().permute_dims(proj_data.data, self.transpose_tuple) + proj_data = self._proj_space.array_namespace.permute_dims(proj_data.data, self.transpose_tuple) else: raise NotImplementedError From d5fefa57a25b4af1639437b310bb73eb2f6ae7a7 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Julius=20H=C3=A4ger?= Date: Thu, 7 May 2026 11:02:11 +0200 Subject: [PATCH 6/7] Try and fix the dimention missmatch problem. --- src/odl/applications/tomo/backends/astra_cuda.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/odl/applications/tomo/backends/astra_cuda.py b/src/odl/applications/tomo/backends/astra_cuda.py index 37c6eb3081a..9155e99c324 100644 --- a/src/odl/applications/tomo/backends/astra_cuda.py +++ b/src/odl/applications/tomo/backends/astra_cuda.py @@ -113,6 +113,7 @@ def __init__(self, geometry, vol_space, proj_space): if self.geometry.ndim == 3: self.transpose_tuple = (1,0,2) if self.geometry.det_curvature_radius is None else (2, 0, 1) + self.inverse_transpose_tuple = (1,0,2) if self.geometry.det_curvature_radius is None else (1, 2, 0) self.fp_scaling_factor = astra_cuda_fp_scaling_factor(self.geometry) self.bp_scaling_factor = astra_cuda_bp_scaling_factor( @@ -244,7 +245,7 @@ def _call_forward_real(self, vol_data:DiscretizedSpaceElement, out=None, **kwarg proj_data = ( proj_data[0] if self.geometry.ndim == 2 - else self._proj_space.array_namespace.permute_dims(proj_data, self.transpose_tuple) + else self._proj_space.array_namespace.permute_dims(proj_data, self.invserse_transpose_tuple) ) if out is not None: From 900090a9e4633678a592940457e682fd15496e99 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Julius=20H=C3=A4ger?= Date: Thu, 7 May 2026 11:07:05 +0200 Subject: [PATCH 7/7] Fix typo --- src/odl/applications/tomo/backends/astra_cuda.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/odl/applications/tomo/backends/astra_cuda.py b/src/odl/applications/tomo/backends/astra_cuda.py index 9155e99c324..5475bd7d5f1 100644 --- a/src/odl/applications/tomo/backends/astra_cuda.py +++ b/src/odl/applications/tomo/backends/astra_cuda.py @@ -245,7 +245,7 @@ def _call_forward_real(self, vol_data:DiscretizedSpaceElement, out=None, **kwarg proj_data = ( proj_data[0] if self.geometry.ndim == 2 - else self._proj_space.array_namespace.permute_dims(proj_data, self.invserse_transpose_tuple) + else self._proj_space.array_namespace.permute_dims(proj_data, self.inverse_transpose_tuple) ) if out is not None: