From d24bad2c57b7470ae0e88faff9f122c2e070f180 Mon Sep 17 00:00:00 2001 From: Holger Kohr Date: Fri, 9 Jun 2017 00:24:44 +0200 Subject: [PATCH 01/24] ENH: add vector geometries, closes #1023 --- examples/tomo/ray_trafo_vec_geom_3d.py | 78 +++++ odl/tomo/backends/astra_setup.py | 126 ++++++- odl/tomo/geometry/conebeam.py | 174 +++++++++- odl/tomo/geometry/geometry.py | 450 ++++++++++++++++++++++++- odl/tomo/geometry/parallel.py | 110 +++++- odl/tomo/operators/ray_trafo.py | 26 +- 6 files changed, 922 insertions(+), 42 deletions(-) create mode 100644 examples/tomo/ray_trafo_vec_geom_3d.py diff --git a/examples/tomo/ray_trafo_vec_geom_3d.py b/examples/tomo/ray_trafo_vec_geom_3d.py new file mode 100644 index 00000000000..1e27bc29056 --- /dev/null +++ b/examples/tomo/ray_trafo_vec_geom_3d.py @@ -0,0 +1,78 @@ +"""Example using the ray transform a custom vector geometry. + +We manually build a "circle plus line trajectory" (CLT) geometry by +extracting the vectors from a circular geometry and extending it by +vertical shifts, starting at the initial position. +""" + +import numpy as np +import odl + +# Reconstruction space: discretized functions on the cube [-20, 20]^3 +# with 300 samples per dimension. +reco_space = odl.uniform_discr( + min_pt=[-20, -20, -20], max_pt=[20, 20, 20], shape=[300, 300, 300], + dtype='float32') + +# First part of the geometry: a 3D single-axis parallel beam geometry with +# flat detector +# Angles: uniformly spaced, n = 180, min = 0, max = 2 * pi +angle_partition = odl.uniform_partition(0, 2 * np.pi, 180) +# Detector: uniformly sampled, n = (512, 512), min = (-30, -30), max = (30, 30) +detector_partition = odl.uniform_partition([-30, -30], [30, 30], [512, 512]) +circle_geometry = odl.tomo.CircularConeFlatGeometry( + angle_partition, detector_partition, src_radius=1000, det_radius=100, + axis=[1, 0, 0]) + +circle_vecs = odl.tomo.astra_conebeam_3d_geom_to_vec(circle_geometry) + +# Cover the whole volume vertically, somewhat undersampled though +vert_shift_min = -22 +vert_shift_max = 22 +num_shifts = 180 +vert_shifts = np.linspace(vert_shift_min, vert_shift_max, num=num_shifts) +inital_vecs = circle_vecs[0] + +# Start from the initial position of the circle vectors and add the vertical +# shifts to the columns 2 and 5 (source and detector z positions) +line_vecs = np.repeat(circle_vecs[0][None, :], num_shifts, axis=0) +line_vecs[:, 2] += vert_shifts +line_vecs[:, 5] += vert_shifts + +# Build the composed geometry and the corresponding ray transform +# (= forward projection) +composed_vecs = np.vstack([circle_vecs, line_vecs]) +composed_geom = odl.tomo.ConeVecGeometry(detector_partition.shape, + composed_vecs) + +ray_trafo = odl.tomo.RayTransform(reco_space, composed_geom) + +# Create a Shepp-Logan phantom (modified version) and projection data +phantom = odl.phantom.shepp_logan(reco_space, True) +proj_data = ray_trafo(phantom) + +# Back-projection can be done by simply calling the adjoint operator on the +# projection data (or any element in the projection space). +backproj = ray_trafo.adjoint(proj_data) + +# Show the slice z=0 of phantom and backprojection, as well as a projection +# image at theta=0 and a sinogram at v=0 (middle detector row) +phantom.show(coords=[None, None, 0], title='Phantom, middle z slice') +backproj.show(coords=[None, None, 0], title='Back-projection, middle z slice') +proj_data.show(indices=[0, slice(None), slice(None)], + title='Projection 0 (circle start)') +proj_data.show(indices=[45, slice(None), slice(None)], + title='Projection 45 (circle 1/4)') +proj_data.show(indices=[90, slice(None), slice(None)], + title='Projection 90 (circle 1/2)') +proj_data.show(indices=[135, slice(None), slice(None)], + title='Projection 135 (circle 3/4)') +proj_data.show(indices=[179, slice(None), slice(None)], + title='Projection 179 (circle end)') +proj_data.show(indices=[180, slice(None), slice(None)], + title='Projection 180 (line start)') +proj_data.show(indices=[270, slice(None), slice(None)], + title='Projection 270 (line middle)') +proj_data.show(indices=[359, slice(None), slice(None)], + title='Projection 359 (line end)') +proj_data.show(coords=[None, None, 0], title='Sinogram, middle slice') diff --git a/odl/tomo/backends/astra_setup.py b/odl/tomo/backends/astra_setup.py index c847f1f2574..609fa387ba6 100644 --- a/odl/tomo/backends/astra_setup.py +++ b/odl/tomo/backends/astra_setup.py @@ -31,8 +31,8 @@ from odl.discr import DiscretizedSpace, DiscretizedSpaceElement from odl.tomo.geometry import ( - DivergentBeamGeometry, Flat1dDetector, Flat2dDetector, Geometry, - ParallelBeamGeometry) + ConeVecGeometry, DivergentBeamGeometry, Flat1dDetector, Flat2dDetector, + Geometry, ParallelBeamGeometry, ParallelVecGeometry) from odl.tomo.util.utility import euler_matrix try: @@ -43,6 +43,7 @@ else: ASTRA_AVAILABLE = True + # Make sure that ASTRA >= 1.7 is used if ASTRA_AVAILABLE: try: @@ -271,6 +272,75 @@ def astra_volume_geometry(vol_space): return vol_geom +def vecs_odl_axes_to_astra_axes(vecs): + """Convert geometry vectors from ODL axis convention to ASTRA. + + Parameters + ---------- + vecs : array-like, shape ``(N, 6)`` or ``(N, 12)`` + Vectors defining the geometric configuration in each + projection. The number ``N`` of rows determines the number + of projections, and the number of columns the spatial + dimension (6 for 2D, 12 for 3D). + + Returns + ------- + astra_vecs : `numpy.ndarray`, same shape as ``vecs`` + The converted geometry vectors. + """ + vecs = np.asarray(vecs, dtype=float) + + if vecs.shape[1] == 6: + # 2D geometry, nothing to do since the axes are the same + return vecs + elif vecs.shape[1] == 12: + # 3D geometry + # 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. + newind = [] + for i in range(4): + newind.extend([2 + 3 * i, 1 + 3 * i, 0 + 3 * i]) + return vecs[:, newind] + else: + raise ValueError('`vecs` must have shape (N, 6) or (N, 12), got ' + 'array with shape {}'.format(vecs.shape)) + + +def vecs_astra_axes_to_odl_axes(vecs): + """Convert geometry vectors from ASTRA axis convention to ODL. + + Parameters + ---------- + vecs : array-like, shape ``(N, 6)`` or ``(N, 12)`` + Vectors defining the geometric configuration in each + projection. The number ``N`` of rows determines the number + of projections, and the number of columns the spatial + dimension (6 for 2D, 12 for 3D). + + Returns + ------- + odl_vecs : `numpy.ndarray`, same shape as ``vecs`` + The converted geometry vectors. + """ + vecs = np.asarray(vecs, dtype=float) + + if vecs.shape[1] == 6: + # 2D geometry, nothing to do since the axes are the same + return vecs + elif vecs.shape[1] == 12: + # 3D geometry + # 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. + newind = [] + for i in range(4): + newind.extend([2 + 3 * i, 1 + 3 * i, 0 + 3 * i]) + newind = np.argsort(newind).tolist() + return vecs[:, newind] + else: + raise ValueError('`vecs` must have shape (N, 6) or (N, 12), got ' + 'array with shape {}'.format(vecs.shape)) + + def astra_conebeam_3d_geom_to_vec(geometry): """Create vectors for ASTRA projection geometries from ODL geometry. @@ -329,14 +399,7 @@ def astra_conebeam_3d_geom_to_vec(geometry): vectors[:, 9:12] = det_axes[0] * px_sizes[0] vectors[:, 6:9] = det_axes[1] * px_sizes[1] - # 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. - newind = [] - for i in range(4): - newind += [2 + 3 * i, 1 + 3 * i, 0 + 3 * i] - vectors = vectors[:, newind] - - return vectors + return vecs_odl_axes_to_astra_axes(vectors) def astra_conebeam_2d_geom_to_vec(geometry): @@ -395,7 +458,7 @@ def astra_conebeam_2d_geom_to_vec(geometry): px_size = geometry.det_partition.cell_sides[0] vectors[:, 4:6] = det_axis * px_size - return vectors + return vecs_odl_axes_to_astra_axes(vectors) def astra_parallel_3d_geom_to_vec(geometry): @@ -457,13 +520,7 @@ def astra_parallel_3d_geom_to_vec(geometry): vectors[:, 9:12] = det_axes[0] * px_sizes[0] vectors[:, 6:9] = det_axes[1] * px_sizes[1] - # 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. - new_ind = [] - for i in range(4): - new_ind += [2 + 3 * i, 1 + 3 * i, 0 + 3 * i] - vectors = vectors[:, new_ind] - return vectors + return vecs_odl_axes_to_astra_axes(vectors) def astra_projection_geometry(geometry): @@ -513,6 +570,22 @@ def astra_projection_geometry(geometry): vec = astra_conebeam_2d_geom_to_vec(geometry) proj_geom = astra.create_proj_geom('fanflat_vec', det_count, vec) + elif isinstance(geometry, ParallelVecGeometry) and geometry.ndim == 2: + det_count = geometry.detector.size + vec = geometry.vectors + # TODO: flip axes? + if not astra_supports('parallel2d_vec_geometry'): + raise NotImplementedError( + "'parallel_vec' geometry not supported by ASTRA " + 'v{}'.format(ASTRA_VERSION)) + proj_geom = astra.create_proj_geom('parallel_vec', det_count, vec) + + elif isinstance(geometry, ConeVecGeometry) and geometry.ndim == 2: + det_count = geometry.detector.size + vec = geometry.vectors + # TODO: flip axes? + proj_geom = astra.create_proj_geom('fanflat_vec', det_count, vec) + elif (isinstance(geometry, ParallelBeamGeometry) and isinstance(geometry.detector, (Flat1dDetector, Flat2dDetector)) and geometry.ndim == 3): @@ -532,6 +605,23 @@ def astra_projection_geometry(geometry): 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, ParallelVecGeometry) and geometry.ndim == 3: + det_row_count = geometry.det_partition.shape[1] + det_col_count = geometry.det_partition.shape[0] + vec = geometry.vectors + # TODO: flip axes? + proj_geom = astra.create_proj_geom('parallel3d_vec', det_row_count, + det_col_count, vec) + + elif isinstance(geometry, ConeVecGeometry) and geometry.ndim == 3: + det_row_count = geometry.det_partition.shape[1] + det_col_count = geometry.det_partition.shape[0] + vec = geometry.vectors + # TODO: flip axes? + proj_geom = astra.create_proj_geom('cone_vec', det_row_count, + det_col_count, vec) + else: raise NotImplementedError('unknown ASTRA geometry type {!r}' ''.format(geometry)) diff --git a/odl/tomo/geometry/conebeam.py b/odl/tomo/geometry/conebeam.py index 832d3f1a24e..15d05bbbb2d 100644 --- a/odl/tomo/geometry/conebeam.py +++ b/odl/tomo/geometry/conebeam.py @@ -14,16 +14,21 @@ from odl.discr import uniform_partition from odl.tomo.geometry.detector import ( - CircularDetector, CylindricalDetector, Flat1dDetector, Flat2dDetector, - SphericalDetector) + Flat1dDetector, Flat2dDetector, + CircularDetector, CylindricalDetector, SphericalDetector) from odl.tomo.geometry.geometry import ( - AxisOrientedGeometry, DivergentBeamGeometry) + AxisOrientedGeometry, DivergentBeamGeometry, VecGeometry) from odl.tomo.util.utility import ( euler_matrix, is_inside_bounds, transform_system) from odl.util import array_str, indent, signature_string -__all__ = ('FanBeamGeometry', 'ConeBeamGeometry', - 'cone_beam_geometry', 'helical_geometry') +__all__ = ( + 'FanBeamGeometry', + 'ConeBeamGeometry', + 'ConeVecGeometry', + 'cone_beam_geometry', + 'helical_geometry', +) class FanBeamGeometry(DivergentBeamGeometry): @@ -1919,6 +1924,165 @@ def helical_geometry(space, src_radius, det_radius, num_turns, pitch=pitch) +class ConeVecGeometry(VecGeometry): + + """Cone beam 2D or 3D geometry defined by a collection of vectors. + + This geometry gives maximal flexibility for representing locations + and orientations of source and detector for cone beam acquisition + schemes. It is defined by a set of vectors per projection, namely + + - source position (``src``), + - detector center (``d``), + - in 2D: vector (``u``) from the detector point with index ``0`` + to the point with index ``1`` + - in 3D: + + * vector (``u``) from detector point ``(0, 0)`` to ``(1, 0)`` and + * vector (``v``) from detector point ``(0, 0)`` to ``(0, 1)``. + + The vectors are given as a matrix with ``n_projs`` rows, where + ``n_projs`` is the number of projections. In 2D, 3 vectors have to + be specified, which makes for ``3 * 2 = 6`` columns:: + + vec2d = (srcX, srcY, dX, dY, uX, uY) + + 3D geometries require 4 vectors, resulting in ``4 * 3 = 12`` matrix + columns:: + + vec3d = (srcX, srcY, srcZ, dX, dY, dZ, uX, uY, uZ, vX, vY, vZ) + + This representation corresponds exactly to the ASTRA "vec" geometries, + see the `ASTRA documentation + `_. + + The geometry defined by these vectors is discrete in the motion + parameters ("angles") since they are merely indices for the individual + projections. For intermediate values, linear interpolation is used, + which corresponds to the assumption that the system moves on piecewise + linear paths. + """ + + @property + def _slice_src(self): + """Slice for the source position part of `vectors`.""" + if self.ndim == 2: + return slice(0, 2) + elif self.ndim == 3: + return slice(0, 3) + else: + raise RuntimeError('bad `ndim`') + + def src_position(self, index): + """Source position function. + + Parameters + ---------- + index : `motion_params` element + Index of the projection. Non-integer indices result in + interpolation, making the source point move on a piecewise + linear path. + + Returns + ------- + pos : `numpy.ndarray` (shape (`ndim`,)) + Source position, an `ndim`-dimensional vector. + + Examples + -------- + Initialize a 2D cone beam geometry with source positions ``(0, -1)`` + and ``(1, 0)`` (first two columns in the vectors): + + >>> vecs_2d = [[0, -1, 0, 1, 1, 0], + ... [1, 0, -1, 0, 0, 1]] + >>> det_shape_2d = 10 + >>> geom_2d = odl.tomo.ConeVecGeometry(det_shape_2d, vecs_2d) + >>> geom_2d.src_position(0) + array([ 0., -1.]) + >>> geom_2d.src_position(1) + array([ 1., 0.]) + >>> geom_2d.src_position(0.5) # mean value + array([ 0.5, -0.5]) + + In 3D, the first three columns determine the source position, + here ``(0, -1, 0)`` and ``(1, 0, 0)``: + + >>> vecs_3d = [[0, -1, 0, 0, 1, 0, 1, 0, 0, 0, 0, 1], + ... [1, 0, 0, -1, 0, 0, 0, 1, 0, 0, 0, 1]] + >>> det_shape_3d = (10, 20) + >>> geom_3d = odl.tomo.ConeVecGeometry(det_shape_3d, vecs_3d) + >>> geom_3d.src_position(0) + array([ 0., -1., 0.]) + >>> geom_3d.src_position(1) + array([ 1., 0., 0.]) + >>> geom_3d.src_position(0.5) # mean value + array([ 0.5, -0.5, 0. ]) + """ + if self.check_bounds and not self.motion_params.contains_all(index): + raise ValueError('`index` {} not in the valid range {}' + ''.format(index, self.motion_params)) + + index = np.array(index, dtype=float, copy=False, ndmin=1) + int_part = index.astype(int) + frac_part = index - int_part + + pos_vecs = np.empty((len(index), self.ndim)) + at_right_bdry = (int_part == self.motion_params.max_pt) + pos_vecs[at_right_bdry, :] = self.vectors[int_part[at_right_bdry], + self._slice_src] + + not_at_right_bdry = ~at_right_bdry + if np.any(not_at_right_bdry): + pt_left = self.vectors[int_part[not_at_right_bdry], + self._slice_src] + pt_right = self.vectors[int_part[not_at_right_bdry] + 1, + self._slice_src] + pos_vecs[not_at_right_bdry, :] = ( + pt_left + + frac_part[not_at_right_bdry, None] * (pt_right - pt_left)) + return pos_vecs.squeeze() + + def det_to_src(self, mparam, dparam, normalized=True): + """Vector pointing from a detector location to the source. + + A function of the motion and detector parameters. + + The default implementation uses the `det_point_position` and + `src_position` functions. Implementations can override this, for + example if no source position is given. + + Parameters + ---------- + mpar : `motion_params` element + Motion parameter at which to evaluate. + dpar : `det_params` element + Detector parameter at which to evaluate. + normalized : bool, optional + If ``True``, return a normalized (unit) vector. + + Returns + ------- + vec : `numpy.ndarray`, shape (`ndim`,) + (Unit) vector pointing from the detector to the source. + """ + if self.check_bounds: + if not self.motion_params.contains_all(mparam): + raise ValueError('`mparam` {} not in the valid range {}' + ''.format(mparam, self.motion_params)) + if not self.det_params.contains_all(dparam): + raise ValueError('`dparam` {} not in the valid range {}' + ''.format(dparam, self.det_params)) + + vec = (self.src_position(mparam) - + self.det_point_position(mparam, dparam)) + + if normalized: + # axis = -1 allows this to be vectorized + vec /= np.linalg.norm(vec, axis=-1) + + return vec + + if __name__ == '__main__': from odl.util.testutils import run_doctests run_doctests() diff --git a/odl/tomo/geometry/geometry.py b/odl/tomo/geometry/geometry.py index 65e4d45939f..05e656ff0ba 100644 --- a/odl/tomo/geometry/geometry.py +++ b/odl/tomo/geometry/geometry.py @@ -1,4 +1,4 @@ -# Copyright 2014-2017 The ODL contributors +# Copyright 2014-2019 The ODL contributors # # This file is part of ODL. # @@ -8,16 +8,24 @@ """Geometry base and mixin classes.""" -from __future__ import print_function, division, absolute_import +from __future__ import absolute_import, division, print_function + from builtins import object + import numpy as np -from odl.discr import RectPartition -from odl.tomo.geometry.detector import Detector +from odl.discr import RectPartition, uniform_partition +from odl.tomo.geometry.detector import Detector, Flat1dDetector, Flat2dDetector from odl.tomo.util import axis_rotation_matrix, is_inside_bounds +from odl.util import ( + indent_rows, normalized_scalar_param_list, safe_int_conv, signature_string) - -__all__ = ('Geometry', 'DivergentBeamGeometry', 'AxisOrientedGeometry') +__all__ = ( + 'Geometry', + 'DivergentBeamGeometry', + 'AxisOrientedGeometry', + 'VecGeometry', +) class Geometry(object): @@ -619,6 +627,436 @@ def rotation_matrix(self, angle): return matrix +class VecGeometry(Geometry): + + """Abstract 2D or 3D geometry defined by a collection of vectors. + + This geometry gives maximal flexibility for representing locations + and orientations of ray/source and detector for parallel or cone beam + acquisition schemes. It is defined by a set of vectors per projection, + namely + + - for parallel beam: ray direction (``ray``), + - for cone beam: source position (``src``), + - detector center (``d``), + - in 2D: vector (``u``) from the detector point with index ``0`` + to the point with index ``1`` + - in 3D: + + * vector (``u``) from detector point ``(0, 0)`` to ``(1, 0)`` and + * vector (``v``) from detector point ``(0, 0)`` to ``(0, 1)``. + + The vectors are given as a matrix with ``n_projs`` rows, where + ``n_projs`` is the number of projections. In 2D, 3 vectors have to + be specified, which makes for ``3 * 2 = 6`` columns:: + + vec_par2d = (rayX, rayY, dX, dY, uX, uY) + vec_cone2d = (srcX, srcY, dX, dY, uX, uY) + + 3D geometries require 4 vectors, resulting in ``4 * 3 = 12`` matrix + columns:: + + vec_par3d = (rayX, rayY, rayZ, dX, dY, dZ, uX, uY, uZ, vX, vY, vZ) + vec_cone3d = (srcX, srcY, srcZ, dX, dY, dZ, uX, uY, uZ, vX, vY, vZ) + + This representation corresponds exactly to the ASTRA "vec" geometries, + see the `ASTRA documentation + `_. + + The geometry defined by these vectors is discrete in the motion + parameters ("angles") since they are merely indices for the individual + projections. For intermediate values, linear interpolation is used, + which corresponds to the assumption that the system moves on piecewise + linear paths. + """ + + def __init__(self, det_shape, vectors): + """Initialize a new instance. + + Parameters + ---------- + det_shape : positive int or sequence of positive ints + Number of detector pixels per axis, can be given as single + integer for 2D geometries (with 1D detector). + vectors : array-like, shape ``(N, 6)`` or ``(N, 12)`` + Vectors defining the geometric configuration in each + projection. The number ``N`` of rows determines the number + of projections, and the number of columns the spatial + dimension (6 for 2D, 12 for 3D). + """ + self.__vectors = np.asarray(vectors, dtype=float) + if self.vectors.ndim != 2: + raise ValueError('`vectors` must be 2-dimensional, got array ' + 'with ndim = {}'.format(self.vectors.ndim)) + + num_projs = self.vectors.shape[0] + + if num_projs == 0: + raise ValueError('`vectors` is empty') + + if self.vectors.shape[1] == 6: + ndim = 2 + elif self.vectors.shape[1] == 12: + ndim = 3 + else: + raise ValueError('`vectors` must have 6 or 12 columns, got ' + 'array with {} columns' + ''.format(self.vectors.shape[1])) + + # Make angle partition with spacing 1 + mpart = uniform_partition(0, num_projs - 1, num_projs, + nodes_on_bdry=True) + + # Detector is initialized using the first set of vectors. This + # assumes that the detector pixels do not change their physical + # sizes during acquisition, but has no effect on the computations. + # For those, only the `vectors` are used. + if ndim == 2: + det_u = self.vectors[0, 4:6] + det_shape = normalized_scalar_param_list(det_shape, length=1, + param_conv=safe_int_conv) + det_cell_size = np.linalg.norm(det_u) + det_part = uniform_partition( + min_pt=-(det_cell_size * det_shape[0]) / 2, + max_pt=(det_cell_size * det_shape[0]) / 2, + shape=det_shape) + detector = Flat1dDetector(det_part, axis=det_u) + elif ndim == 3: + det_u = self.vectors[0, 6:9] + det_v = self.vectors[0, 9:12] + det_shape = normalized_scalar_param_list(det_shape, length=2, + param_conv=safe_int_conv) + det_cell_sides = np.array( + [np.linalg.norm(det_u), np.linalg.norm(det_v)]) + det_part = uniform_partition( + min_pt=-(det_cell_sides * det_shape) / 2, + max_pt=(det_cell_sides * det_shape) / 2, + shape=det_shape) + detector = Flat2dDetector(det_part, axes=[det_u, det_v]) + + Geometry.__init__(self, ndim, mpart, detector) + + @property + def vectors(self): + """Array of vectors defining this geometry.""" + return self.__vectors + + @property + def _slice_det_center(self): + """Slice for the detector center part of `vectors`.""" + if self.ndim == 2: + return slice(2, 4) + elif self.ndim == 3: + return slice(3, 6) + else: + raise RuntimeError('bad `ndim`') + + @property + def _slice_det_u(self): + """Slice for the detector u vector part of `vectors`.""" + if self.ndim == 2: + return slice(4, 6) + elif self.ndim == 3: + return slice(6, 9) + else: + raise RuntimeError('bad `ndim`') + + @property + def _slice_det_v(self): + """Slice for the detector v vector part of `vectors`.""" + if self.ndim == 3: + return slice(9, 12) + else: + raise RuntimeError('bad `ndim`') + + def det_refpoint(self, index): + """Detector reference point function. + + Parameters + ---------- + index : `motion_params` element + Index of the projection. Non-integer indices result in + interpolation, making the reference point move on a piecewise + linear path. + + Returns + ------- + point : `numpy.ndarray`, shape (`ndim`,) + The reference point, an `ndim`-dimensional vector. + + Examples + -------- + Initialize a 2D parallel beam geometry with detector positions + ``(0, 1)`` and ``(-1, 0)`` (columns 2 and 3 in the vectors, + counting from 0): + + >>> # d = refpoint + >>> # dx dy + >>> vecs_2d = [[0, -1, 0, 1, 1, 0], + ... [1, 0, -1, 0, 0, 1]] + >>> det_shape_2d = 10 + >>> geom_2d = odl.tomo.ParallelVecGeometry(det_shape_2d, vecs_2d) + >>> geom_2d.det_refpoint(0) + array([ 0., 1.]) + >>> geom_2d.det_refpoint(1) + array([-1., 0.]) + >>> geom_2d.det_refpoint(0.5) # mean value + array([-0.5, 0.5]) + + In 3D, columns 3 to 5 (starting at 0) determine the detector + reference point, here ``(0, 1, 0)`` and ``(-1, 0, 0)``: + + >>> # d = refpoint + >>> # dx dy dz + >>> vecs_3d = [[0, -1, 0, 0, 1, 0, 1, 0, 0, 0, 0, 1], + ... [1, 0, 0, -1, 0, 0, 0, 1, 0, 0, 0, 1]] + >>> det_shape_3d = (10, 20) + >>> geom_3d = odl.tomo.ParallelVecGeometry(det_shape_3d, vecs_3d) + >>> geom_3d.det_refpoint(0) + array([ 0., 1., 0.]) + >>> geom_3d.det_refpoint(1) + array([-1., 0., 0.]) + >>> geom_3d.det_refpoint(0.5) # mean value + array([-0.5, 0.5, 0. ]) + """ + if self.check_bounds and not self.motion_params.contains_all(index): + raise ValueError('`index` {} not in the valid range {}' + ''.format(index, self.motion_params)) + + index = np.array(index, dtype=float, copy=False, ndmin=1) + int_part = index.astype(int) + frac_part = index - int_part + + vecs = np.empty((len(index), self.ndim)) + at_right_bdry = (int_part == self.motion_params.max_pt) + vecs[at_right_bdry, :] = self.vectors[int_part[at_right_bdry], + self._slice_det_center] + + not_at_right_bdry = ~at_right_bdry + if np.any(not_at_right_bdry): + pt_left = self.vectors[int_part[not_at_right_bdry], + self._slice_det_center] + pt_right = self.vectors[int_part[not_at_right_bdry] + 1, + self._slice_det_center] + vecs[not_at_right_bdry, :] = ( + pt_left + + frac_part[not_at_right_bdry, None] * (pt_right - pt_left)) + return vecs.squeeze() + + def det_point_position(self, index, dparam): + """Return the detector point at ``(index, dparam)``. + + The projection ``index`` determines the location of the detector + reference point, and the detector parameter ``dparam`` defines + an intrinsic shift that is added to the reference point. + + Parameters + ---------- + index : `motion_params` element + Index of the projection. Non-integer indices result in + interpolated vectors. + dparam : `det_params` element + Detector parameter at which to evaluate. + + Returns + ------- + pos : `numpy.ndarray` (shape (`ndim`,)) + Source position, an `ndim`-dimensional vector. + + Examples + -------- + Initialize a 2D parallel beam geometry with detector positions + ``(0, 1)`` and ``(-1, 0)``, and detector axes ``(1, 0)`` and + ``(0, 1)``, respectively: + + >>> # d = refpoint, u = detector axis + >>> # dx dy ux uy + >>> vecs_2d = [[0, -1, 0, 1, 1, 0], + ... [1, 0, -1, 0, 0, 1]] + >>> det_shape_2d = 10 + >>> geom_2d = odl.tomo.ParallelVecGeometry(det_shape_2d, vecs_2d) + >>> # This is equal to d(0) = (0, 1) + >>> geom_2d.det_point_position(0, 0) + array([ 0., 1.]) + >>> # d(0) + 2 * u(0) = (0, 1) + 2 * (1, 0) + >>> geom_2d.det_point_position(0, 2) + array([ 2., 1.]) + >>> # d(1) + 2 * u(1) = (-1, 0) + 2 * (0, 1) + >>> geom_2d.det_point_position(1, 2) + array([-1., 2.]) + >>> # d(0.5) + 2 * u(0.5) = (-0.5, 0.5) + 2 * (0.5, 0.5) + >>> geom_2d.det_point_position(0.5, 2) + array([ 0.5, 1.5]) + + Do the same in 3D, with reference points ``(0, 1, 0)`` and + ``(-1, 0, 0)``, and horizontal ``u`` axis vectors ``(1, 0, 0)`` and + ``(0, 1, 0)``. The vertical axis ``v`` is always ``(0, 0, 1)``: + + >>> # d = refpoint, u = detector axis 0, v = detector axis 1 + >>> # dx dy dz ux uy uz vx vy vz + >>> vecs_3d = [[0, -1, 0, 0, 1, 0, 1, 0, 0, 0, 0, 1], + ... [1, 0, 0, -1, 0, 0, 0, 1, 0, 0, 0, 1]] + >>> det_shape_3d = (10, 20) + >>> geom_3d = odl.tomo.ParallelVecGeometry(det_shape_3d, vecs_3d) + >>> # This is equal to d(0) = (0, 1, 0) + >>> geom_3d.det_point_position(0, [0, 0]) + array([ 0., 1., 0.]) + >>> # d(0) + 2 * u(0) + 1 * v(0) = + >>> # (0, 1, 0) + 2 * (1, 0, 0) + 1 * (0, 0, 1) + >>> geom_3d.det_point_position(0, [2, 1]) + array([ 2., 1., 1.]) + >>> # d(1) + 2 * u(1) + 1 * v(1) = + >>> # (-1, 0, 0) + 2 * (0, 1, 0) + 1 * (0, 0, 1) + >>> geom_3d.det_point_position(1, [2, 1]) + array([-1., 2., 1.]) + >>> # d(0.5) + 2 * u(0.5) = (-0.5, 0.5) + 2 * (0.5, 0.5) + >>> geom_3d.det_point_position(0.5, [2, 1]) + array([ 0.5, 1.5, 1. ]) + """ + # TODO: vectorize! + + if index not in self.motion_params: + raise ValueError('`index` must be contained in `motion_params` ' + '{}, got {}'.format(self.motion_params, index)) + if dparam not in self.det_params: + raise ValueError('`dparam` must be contained in `det_params` ' + '{}, got {}'.format(self.det_params, dparam)) + + if self.ndim == 2: + det_shift = dparam * self.det_axis(index) + elif self.ndim == 3: + det_shift = sum(di * ax + for di, ax in zip(dparam, self.det_axes(index))) + + return self.det_refpoint(index) + det_shift + + def det_axis(self, index): + """Return the detector axis at ``index`` (for 2D geometries). + + Parameters + ---------- + index : `motion_params` element + Index of the projection. Non-integer indices result in + interpolated vectors. + + Returns + ------- + axis : `numpy.ndarray`, shape ``(2,)`` + The detector axis at ``index``. + + Examples + -------- + Initialize a 2D parallel beam geometry with detector axes + ``(1, 0)`` and ``(0, 1)`` (columns 4 and 5 in the vectors, + counting from 0): + + >>> # d = refpoint + >>> # ux uy + >>> vecs_2d = [[0, -1, 0, 1, 1, 0], + ... [1, 0, -1, 0, 0, 1]] + >>> det_shape_2d = 10 + >>> geom_2d = odl.tomo.ParallelVecGeometry(det_shape_2d, vecs_2d) + >>> geom_2d.det_axis(0) + array([ 1., 0.]) + >>> geom_2d.det_axis(1) + array([ 0., 1.]) + >>> geom_2d.det_axis(0.5) # mean value + array([ 0.5, 0.5]) + """ + # TODO: vectorize! + + if index not in self.motion_params: + raise ValueError('`index` must be contained in `motion_params` ' + '{}, got {}'.format(self.motion_params, index)) + + if self.ndim != 2: + raise NotImplementedError( + '`det_axis` only valid for 2D geometries, use `det_axes` ' + 'in 3D') + + index = float(index) + int_part = int(index) + frac_part = index - int_part + if int_part == self.motion_params.max_pt: + det_u = self.vectors[int_part, self._slice_det_u] + else: + det_u_left = self.vectors[int_part, self._slice_det_u] + det_u_right = self.vectors[int_part + 1, self._slice_det_u] + det_u = det_u_left + frac_part * (det_u_right - det_u_left) + + return det_u + + def det_axes(self, index): + """Return the detector axes at ``index`` (for 2D or 3D geometries). + + Parameters + ---------- + index : `motion_params` element + Index of the projection. Non-integer indices result in + interpolated vectors. + + Returns + ------- + axes : tuple of `numpy.ndarray`, shape ``(ndim,)`` + The detector axes at ``index``. + + Examples + -------- + Initialize a 3D parallel beam geometry with detector axis ``u`` + vectors ``(1, 0, 0)`` and ``(0, 1, 0)`` (columns 6 through 8 in the + vectors, counting from 0), and axis ``v`` equal to ``(0, 0, 1)`` + in both cases (columns 9 through 11): + + >>> # d = refpoint + >>> # ux uy uz vx vy vz + >>> vecs_3d = [[0, -1, 0, 0, 1, 0, 1, 0, 0, 0, 0, 1], + ... [1, 0, 0, -1, 0, 0, 0, 1, 0, 0, 0, 1]] + >>> det_shape_3d = (10, 20) + >>> geom_3d = odl.tomo.ParallelVecGeometry(det_shape_3d, vecs_3d) + >>> geom_3d.det_axes(0) + (array([ 1., 0., 0.]), array([ 0., 0., 1.])) + >>> geom_3d.det_axes(1) + (array([ 0., 1., 0.]), array([ 0., 0., 1.])) + >>> geom_3d.det_axes(0.5) # mean value + (array([ 0.5, 0.5, 0. ]), array([ 0., 0., 1.])) + """ + # TODO: vectorize! + + if index not in self.motion_params: + raise ValueError('`index` must be contained in `motion_params` ' + '{}, got {}'.format(self.motion_params, index)) + + index = float(index) + int_part = int(index) + frac_part = index - int_part + if int_part == self.motion_params.max_pt: + det_u = self.vectors[int_part, self._slice_det_u] + if self.ndim == 2: + return (det_u,) + elif self.ndim == 3: + det_v = self.vectors[int_part, self._slice_det_v] + return (det_u, det_v) + else: + det_u_left = self.vectors[int_part, self._slice_det_u] + det_u_right = self.vectors[int_part + 1, self._slice_det_u] + det_u = det_u_left + frac_part * (det_u_right - det_u_left) + + if self.ndim == 2: + return (det_u,) + elif self.ndim == 3: + det_v_left = self.vectors[int_part, self._slice_det_v] + det_v_right = self.vectors[int_part + 1, self._slice_det_v] + det_v = det_v_left + frac_part * (det_v_right - det_v_left) + return (det_u, det_v) + + def __repr__(self): + """Return ``repr(self)``.""" + posargs = [self.det_partition.shape, self.vectors] + inner_str = signature_string(posargs, [], sep=',\n') + return '{}(\n{}\n)'.format(self.__class__.__name__, + indent_rows(inner_str)) + + if __name__ == '__main__': from odl.util.testutils import run_doctests run_doctests() diff --git a/odl/tomo/geometry/parallel.py b/odl/tomo/geometry/parallel.py index 1d7d3dad6b8..1c57304da98 100644 --- a/odl/tomo/geometry/parallel.py +++ b/odl/tomo/geometry/parallel.py @@ -14,14 +14,19 @@ from odl.discr import uniform_partition from odl.tomo.geometry.detector import Flat1dDetector, Flat2dDetector -from odl.tomo.geometry.geometry import AxisOrientedGeometry, Geometry +from odl.tomo.geometry.geometry import ( + AxisOrientedGeometry, Geometry, VecGeometry) from odl.tomo.util import euler_matrix, is_inside_bounds, transform_system from odl.util import array_str, indent, signature_string -__all__ = ('ParallelBeamGeometry', - 'Parallel2dGeometry', - 'Parallel3dEulerGeometry', 'Parallel3dAxisGeometry', - 'parallel_beam_geometry') +__all__ = ( + 'ParallelBeamGeometry', + 'Parallel2dGeometry', + 'Parallel3dEulerGeometry', + 'Parallel3dAxisGeometry', + 'ParallelVecGeometry', + 'parallel_beam_geometry', +) class ParallelBeamGeometry(Geometry): @@ -1468,6 +1473,101 @@ def __getitem__(self, indices): rotation_matrix = AxisOrientedGeometry.rotation_matrix +class ParallelVecGeometry(VecGeometry): + + """Parallel beam 2D or 3D geometry defined by a collection of vectors. + + This geometry gives maximal flexibility for representing locations + and orientations of rays and detector for parallel beam acquisition + schemes. It is defined by a set of vectors per projection, namely + + - ray direction (``ray``), + - detector center (``d``), + - in 2D: vector (``u``) from the detector point with index ``0`` + to the point with index ``1`` + - in 3D: + + * vector (``u``) from detector point ``(0, 0)`` to ``(1, 0)`` and + * vector (``v``) from detector point ``(0, 0)`` to ``(0, 1)``. + + The vectors are given as a matrix with ``n_projs`` rows, where + ``n_projs`` is the number of projections. In 2D, 3 vectors have to + be specified, which makes for ``3 * 2 = 6`` columns:: + + vec2d = (rayX, rayY, dX, dY, uX, uY) + + 3D geometries require 4 vectors, resulting in ``4 * 3 = 12`` matrix + columns:: + + vec3d = (rayX, rayY, rayZ, dX, dY, dZ, uX, uY, uZ, vX, vY, vZ) + + This representation corresponds exactly to the ASTRA "vec" geometries, + see the `ASTRA documentation + `_. + + The geometry defined by these vectors is discrete in the motion + parameters ("angles") since they are merely indices for the individual + projections. For intermediate values, linear interpolation is used, + which corresponds to the assumption that the system moves on piecewise + linear paths. + """ + + @property + def _slice_ray(self): + """Slice for the ray direction part of `vectors`.""" + if self.ndim == 2: + return slice(0, 2) + elif self.ndim == 3: + return slice(0, 3) + else: + raise RuntimeError('bad `ndim`') + + def det_to_src(self, index, dparam, normalized=True): + """Vector pointing from a detector location to the source. + + A function of the motion and detector parameters. + + Parameters + ---------- + index : `motion_params` element + Index of the projection. Non-integer indices result in + interpolated vectors. + dparam : `det_params` element + Detector parameter at which to evaluate. + normalized : bool, optional + If ``True``, return a normalized (unit) vector. + The non-normalized variant is not available in parallel beam + geometries. + + Returns + ------- + vec : `numpy.ndarray`, shape (`ndim`,) + Unit vector pointing from the detector to the source + """ + if index not in self.motion_params: + raise ValueError('`index` must be contained in `motion_params` ' + '{}, got {}'.format(self.motion_params, index)) + if dparam not in self.det_params: + raise ValueError('`dparam` must be contained in `det_params` ' + '{}, got {}'.format(self.det_params, dparam)) + if not normalized: + raise NotImplementedError('non-normalized detector-to-source ' + 'vector not available for parallel ' + 'beam geometries') + + index = float(index) + int_part = int(index) + frac_part = index - int_part + if int_part == self.motion_params.max_pt: + ray_dir = self.vectors[int_part, self._slice_ray] + else: + ray_left = self.vectors[int_part, self._slice_ray] + ray_right = self.vectors[int_part + 1, self._slice_ray] + ray_dir = ray_left + frac_part * (ray_right - ray_left) + + return -ray_dir / np.linalg.norm(ray_dir) + + def parallel_beam_geometry(space, num_angles=None, det_shape=None): r"""Create default parallel beam geometry from ``space``. diff --git a/odl/tomo/operators/ray_trafo.py b/odl/tomo/operators/ray_trafo.py index 64419fd43b6..edc3e8c7bad 100644 --- a/odl/tomo/operators/ray_trafo.py +++ b/odl/tomo/operators/ray_trafo.py @@ -22,7 +22,7 @@ from odl.tomo.backends.astra_cpu import AstraCpuImpl from odl.tomo.backends.astra_cuda import AstraCudaImpl from odl.tomo.backends.skimage_radon import SkImageImpl -from odl.tomo.geometry import Geometry +from odl.tomo.geometry import Geometry, ParallelVecGeometry from odl.util import is_string # RAY_TRAFO_IMPLS are used by `RayTransform` when no `impl` is given. @@ -35,7 +35,9 @@ if ASTRA_CUDA_AVAILABLE: RAY_TRAFO_IMPLS['astra_cuda'] = AstraCudaImpl -__all__ = ('RayTransform',) +__all__ = ( + 'RayTransform', +) class RayTransform(Operator): @@ -101,8 +103,21 @@ def __init__(self, vol_space, geometry, **kwargs): if proj_space is None: dtype = vol_space.dtype + # TODO: use weighting that differentiates between angles and + # detector. We need an array of weights here, a different value + # for each angle, but constant in the detector axes. + # The required partition property is available since + # commit a551190d, but weighting is not adapted yet. + # See also issues #286 and #1001 + # The constant in a given angle has to be the average distance + # to previous and next angles. Probably there should also be + # an option to use the current implementation for faster + # runs (one weighting constant). if not vol_space.is_weighted: weighting = None + elif isinstance(self.geometry, ParallelVecGeometry): + # TODO: change to weighting constant per angle when available + weighting = 1.0 elif ( isinstance(vol_space.weighting, ConstWeighting) and np.isclose( @@ -110,12 +125,7 @@ def __init__(self, vol_space, geometry, **kwargs): ) ): # Approximate cell volume - # TODO: find a way to treat angles and detector differently - # regarding weighting. While the detector should be uniformly - # discretized, the angles do not have to and often are not. - # The needed partition property is available since - # commit a551190d, but weighting is not adapted yet. - # See also issue #286 + # TODO: change to weighting constant per angle when available extent = float(geometry.partition.extent.prod()) size = float(geometry.partition.size) weighting = extent / size From 452a9bf003a46159e137512fbd79b8ca1932b1c4 Mon Sep 17 00:00:00 2001 From: Holger Kohr Date: Fri, 9 Jun 2017 18:59:38 +0200 Subject: [PATCH 02/24] WIP: fix weighting for vec backprojection --- examples/tomo/ray_trafo_vec_geom_3d.py | 2 +- odl/tomo/backends/astra_cuda.py | 86 +++++++++++++++++++++++--- odl/tomo/backends/astra_setup.py | 4 -- odl/tomo/operators/ray_trafo.py | 15 +++-- 4 files changed, 84 insertions(+), 23 deletions(-) diff --git a/examples/tomo/ray_trafo_vec_geom_3d.py b/examples/tomo/ray_trafo_vec_geom_3d.py index 1e27bc29056..a6d4f076b36 100644 --- a/examples/tomo/ray_trafo_vec_geom_3d.py +++ b/examples/tomo/ray_trafo_vec_geom_3d.py @@ -1,4 +1,4 @@ -"""Example using the ray transform a custom vector geometry. +"""Example using the ray transform with a custom vector geometry. We manually build a "circle plus line trajectory" (CLT) geometry by extracting the vectors from a circular geometry and extending it by diff --git a/odl/tomo/backends/astra_cuda.py b/odl/tomo/backends/astra_cuda.py index 03b6eb877ae..6dd315e091c 100644 --- a/odl/tomo/backends/astra_cuda.py +++ b/odl/tomo/backends/astra_cuda.py @@ -23,8 +23,9 @@ astra_volume_geometry) from odl.tomo.backends.util import _add_default_complex_impl from odl.tomo.geometry import ( - ConeBeamGeometry, FanBeamGeometry, Geometry, Parallel2dGeometry, - Parallel3dAxisGeometry) + ConeBeamGeometry, ConeVecGeometry, FanBeamGeometry, Geometry, + Parallel2dGeometry, Parallel3dAxisGeometry, ParallelVecGeometry, + VecGeometry) try: import astra @@ -38,7 +39,7 @@ ) -class AstraCudaImpl: +class AstraCudaImpl(object): """`RayTransform` implementation for CUDA algorithms in ASTRA.""" algo_forward_id = None @@ -337,19 +338,22 @@ def astra_cuda_bp_scaling_factor(proj_space, vol_space, geometry): track of it and adapt the scaling accordingly. """ # Angular integration weighting factor - # angle interval weight by approximate cell volume - angle_extent = geometry.motion_partition.extent - num_angles = geometry.motion_partition.shape - # TODO: this gives the wrong factor for Parallel3dEulerGeometry with - # 2 angles - scaling_factor = (angle_extent / num_angles).prod() + if isinstance(geometry, VecGeometry): + angle_weighting = 1.0 + else: + # Average angle step, approximation in the case of non-uniform angles + angle_extent = geometry.motion_partition.extent + num_angles = geometry.motion_partition.shape + # TODO: this gives the wrong factor for Parallel3dEulerGeometry with + # 2 angles + angle_weighting = (angle_extent / num_angles).prod() # Correct in case of non-weighted spaces proj_extent = float(proj_space.partition.extent.prod()) proj_size = float(proj_space.partition.size) proj_weighting = proj_extent / proj_size - scaling_factor *= ( + scaling_factor = angle_weighting * ( proj_space.weighting.const / proj_weighting ) scaling_factor /= ( @@ -384,6 +388,27 @@ def astra_cuda_bp_scaling_factor(proj_space, vol_space, geometry): src_radius = geometry.src_radius det_radius = geometry.det_radius scaling_factor *= ((src_radius + det_radius) / src_radius) ** 2 + elif isinstance(geometry, VecGeometry): + scaling_factor = (geometry.det_partition.cell_volume / + vol_space.cell_volume) + elif isinstance(geometry, ParallelVecGeometry): + if geometry.ndim == 2: + # Scales with 1 / cell_volume + scaling_factor *= float(vol_space.cell_volume) + elif geometry.ndim == 3: + # Scales with cell volume + scaling_factor /= vol_space.cell_volume + elif isinstance(geometry, ConeVecGeometry): + # TODO: this should be based on the distances per angle, would + # probably be part of the weighting then + if geometry.ndim == 2: + # Scales with 1 / cell_volume + scaling_factor *= float(vol_space.cell_volume) + elif geometry.ndim == 3: + det_px_area = geometry.det_partition.cell_volume + scaling_factor *= (det_px_area ** 2 / + vol_space.cell_volume ** 2) + elif parse_version(ASTRA_VERSION) < parse_version('1.9.0dev'): # Scaling for the 1.8.x releases if isinstance(geometry, Parallel2dGeometry): @@ -409,6 +434,26 @@ def astra_cuda_bp_scaling_factor(proj_space, vol_space, geometry): src_radius = geometry.src_radius det_radius = geometry.det_radius scaling_factor *= ((src_radius + det_radius) / src_radius) ** 2 + elif isinstance(geometry, VecGeometry): + scaling_factor = (geometry.det_partition.cell_volume / + vol_space.cell_volume) + elif isinstance(geometry, ParallelVecGeometry): + if geometry.ndim == 2: + # Scales with 1 / cell_volume + scaling_factor *= float(vol_space.cell_volume) + elif geometry.ndim == 3: + # Scales with cell volume + scaling_factor /= vol_space.cell_volume + elif isinstance(geometry, ConeVecGeometry): + # TODO: this should be based on the distances per angle, would + # probably be part of the weighting then + if geometry.ndim == 2: + # Scales with 1 / cell_volume + scaling_factor *= float(vol_space.cell_volume) + elif geometry.ndim == 3: + det_px_area = geometry.det_partition.cell_volume + scaling_factor *= (det_px_area ** 2 / + vol_space.cell_volume ** 2) # Correction for scaled 1/r^2 factor in ASTRA's density weighting. # This compensates for scaled voxels and pixels, as well as a @@ -450,6 +495,27 @@ def astra_cuda_bp_scaling_factor(proj_space, vol_space, geometry): # density weighting. det_px_area = geometry.det_partition.cell_volume scaling_factor *= (src_radius ** 2 * det_px_area ** 2) + elif isinstance(geometry, VecGeometry): + scaling_factor = (geometry.det_partition.cell_volume / + vol_space.cell_volume) + elif isinstance(geometry, ParallelVecGeometry): + if geometry.ndim == 2: + # Scales with 1 / cell_volume + scaling_factor *= float(vol_space.cell_volume) + elif geometry.ndim == 3: + # Scales with cell volume + scaling_factor /= vol_space.cell_volume + elif isinstance(geometry, ConeVecGeometry): + # TODO: this should be based on the distances per angle, would + # probably be part of the weighting then + if geometry.ndim == 2: + # Scales with 1 / cell_volume + scaling_factor *= float(vol_space.cell_volume) + elif geometry.ndim == 3: + det_px_area = geometry.det_partition.cell_volume + scaling_factor *= (det_px_area ** 2 / + vol_space.cell_volume ** 2) + else: # Scaling for versions since 1.9.9.dev scaling_factor /= float(vol_space.cell_volume) diff --git a/odl/tomo/backends/astra_setup.py b/odl/tomo/backends/astra_setup.py index 609fa387ba6..5bed4200b31 100644 --- a/odl/tomo/backends/astra_setup.py +++ b/odl/tomo/backends/astra_setup.py @@ -573,7 +573,6 @@ def astra_projection_geometry(geometry): elif isinstance(geometry, ParallelVecGeometry) and geometry.ndim == 2: det_count = geometry.detector.size vec = geometry.vectors - # TODO: flip axes? if not astra_supports('parallel2d_vec_geometry'): raise NotImplementedError( "'parallel_vec' geometry not supported by ASTRA " @@ -583,7 +582,6 @@ def astra_projection_geometry(geometry): elif isinstance(geometry, ConeVecGeometry) and geometry.ndim == 2: det_count = geometry.detector.size vec = geometry.vectors - # TODO: flip axes? proj_geom = astra.create_proj_geom('fanflat_vec', det_count, vec) elif (isinstance(geometry, ParallelBeamGeometry) and @@ -610,7 +608,6 @@ def astra_projection_geometry(geometry): det_row_count = geometry.det_partition.shape[1] det_col_count = geometry.det_partition.shape[0] vec = geometry.vectors - # TODO: flip axes? proj_geom = astra.create_proj_geom('parallel3d_vec', det_row_count, det_col_count, vec) @@ -618,7 +615,6 @@ def astra_projection_geometry(geometry): det_row_count = geometry.det_partition.shape[1] det_col_count = geometry.det_partition.shape[0] vec = geometry.vectors - # TODO: flip axes? proj_geom = astra.create_proj_geom('cone_vec', det_row_count, det_col_count, vec) diff --git a/odl/tomo/operators/ray_trafo.py b/odl/tomo/operators/ray_trafo.py index edc3e8c7bad..4c1794fbe80 100644 --- a/odl/tomo/operators/ray_trafo.py +++ b/odl/tomo/operators/ray_trafo.py @@ -22,7 +22,7 @@ from odl.tomo.backends.astra_cpu import AstraCpuImpl from odl.tomo.backends.astra_cuda import AstraCudaImpl from odl.tomo.backends.skimage_radon import SkImageImpl -from odl.tomo.geometry import Geometry, ParallelVecGeometry +from odl.tomo.geometry import Geometry, VecGeometry from odl.util import is_string # RAY_TRAFO_IMPLS are used by `RayTransform` when no `impl` is given. @@ -115,9 +115,9 @@ def __init__(self, vol_space, geometry, **kwargs): # runs (one weighting constant). if not vol_space.is_weighted: weighting = None - elif isinstance(self.geometry, ParallelVecGeometry): - # TODO: change to weighting constant per angle when available - weighting = 1.0 + elif isinstance(self.geometry, VecGeometry): + # No angular weighting + weighting = geometry.det_partition.cell_volume elif ( isinstance(vol_space.weighting, ConstWeighting) and np.isclose( @@ -125,10 +125,9 @@ def __init__(self, vol_space, geometry, **kwargs): ) ): # Approximate cell volume - # TODO: change to weighting constant per angle when available - extent = float(geometry.partition.extent.prod()) - size = float(geometry.partition.size) - weighting = extent / size + angle_weight = float(geometry.motion_partition.extent.prod() / + geometry.motion_partition.size) + weighting = angle_weight * geometry.det_partition.cell_volume else: raise NotImplementedError('unknown weighting of domain') From 771b38d63a503773f999b5397594ebe06f0b677e Mon Sep 17 00:00:00 2001 From: Holger Kohr Date: Tue, 5 Sep 2017 18:12:30 +0200 Subject: [PATCH 03/24] WIP: vectorize VecGeometry methods --- odl/tomo/geometry/conebeam.py | 8 +- odl/tomo/geometry/geometry.py | 3 +- odl/tomo/geometry/parallel.py | 180 +++++++++++++++++++++++++++------- 3 files changed, 153 insertions(+), 38 deletions(-) diff --git a/odl/tomo/geometry/conebeam.py b/odl/tomo/geometry/conebeam.py index 15d05bbbb2d..ed991e9230e 100644 --- a/odl/tomo/geometry/conebeam.py +++ b/odl/tomo/geometry/conebeam.py @@ -2018,7 +2018,9 @@ def src_position(self, index): >>> geom_3d.src_position(0.5) # mean value array([ 0.5, -0.5, 0. ]) """ - if self.check_bounds and not self.motion_params.contains_all(index): + # TODO: vectorize + if (self.check_bounds and + not is_inside_bounds(index, self.motion_params)): raise ValueError('`index` {} not in the valid range {}' ''.format(index, self.motion_params)) @@ -2066,10 +2068,10 @@ def det_to_src(self, mparam, dparam, normalized=True): (Unit) vector pointing from the detector to the source. """ if self.check_bounds: - if not self.motion_params.contains_all(mparam): + if not is_inside_bounds(mparam, self.motion_params): raise ValueError('`mparam` {} not in the valid range {}' ''.format(mparam, self.motion_params)) - if not self.det_params.contains_all(dparam): + if not is_inside_bounds(dparam, self.det_params): raise ValueError('`dparam` {} not in the valid range {}' ''.format(dparam, self.det_params)) diff --git a/odl/tomo/geometry/geometry.py b/odl/tomo/geometry/geometry.py index 05e656ff0ba..9549e3f91cc 100644 --- a/odl/tomo/geometry/geometry.py +++ b/odl/tomo/geometry/geometry.py @@ -819,7 +819,8 @@ def det_refpoint(self, index): >>> geom_3d.det_refpoint(0.5) # mean value array([-0.5, 0.5, 0. ]) """ - if self.check_bounds and not self.motion_params.contains_all(index): + if (self.check_bounds and + not is_inside_bounds(index, self.motion_params)): raise ValueError('`index` {} not in the valid range {}' ''.format(index, self.motion_params)) diff --git a/odl/tomo/geometry/parallel.py b/odl/tomo/geometry/parallel.py index 1c57304da98..e8c1c2dd66e 100644 --- a/odl/tomo/geometry/parallel.py +++ b/odl/tomo/geometry/parallel.py @@ -1522,50 +1522,162 @@ def _slice_ray(self): else: raise RuntimeError('bad `ndim`') - def det_to_src(self, index, dparam, normalized=True): - """Vector pointing from a detector location to the source. + def det_to_src(self, index, dparam): + """Direction from a detector location to the source. - A function of the motion and detector parameters. + The direction vector is computed as follows:: + + dir = rotation_matrix(angle).dot(detector.surface_normal(dparam)) + + Note that for flat detectors, ``surface_normal`` does not depend + on the parameter ``dparam``, hence this function is constant in + that variable. Parameters ---------- - index : `motion_params` element - Index of the projection. Non-integer indices result in - interpolated vectors. - dparam : `det_params` element - Detector parameter at which to evaluate. - normalized : bool, optional - If ``True``, return a normalized (unit) vector. - The non-normalized variant is not available in parallel beam - geometries. + index : `array-like` or sequence + Index or indices of the projection. Non-integer indices + result in interpolated vectors. + dparam : `array-like` or sequence + Detector parameter(s) at which to evaluate. If + ``det_params.ndim == 2``, a sequence of that length must be + provided. Returns ------- - vec : `numpy.ndarray`, shape (`ndim`,) - Unit vector pointing from the detector to the source + det_to_src : `numpy.ndarray` + Vector(s) pointing from a detector point to the source (at + infinity). + The shape of the returned array is obtained from the + (broadcast) shapes of ``index`` and ``dparam``, and + broadcasting is supported within both parameters and between + them. The precise definition of the shape is + ``broadcast(index, bcast_dparam).shape + (ndim,)``, + where ``bcast_dparam`` is + + - ``dparam`` if `det_params` is 1D, + - ``broadcast(*dparam)`` otherwise. + + Examples + -------- + The method works with single parameter values, in which case + a single vector is returned: + + >>> # r = ray direction + >>> # rx ry + >>> vecs_2d = [[0, -1, 0, 1, 1, 0], + ... [1, 0, -1, 0, 0, 1]] + >>> det_shape_2d = 10 + >>> geom = odl.tomo.ParallelVecGeometry(det_shape_2d, vecs_2d) + >>> geom.det_to_src(0, 0) + array([-0., 1.]) + >>> geom.det_to_src(1, 0) + array([-1., -0.]) + >>> geom.det_to_src(1, 1) # dparam has no effect + array([-1., -0.]) + >>> dir = geom.det_to_src(0.5, 0) # interpolate + >>> np.allclose(dir, [-1, 1] / np.linalg.norm([-1, 1])) + True + + Both variables support vectorized calls, i.e., stacks of + parameters can be provided. The order of axes in the output (left + of the ``ndim`` axis for the vector dimension) corresponds to the + order of arguments: + TODO: adapt this + + >>> dirs = geom.det_to_src(0, [-1, 0, 0.5, 1]) + >>> dirs + array([[ 0., -1.], + [ 0., -1.], + [ 0., -1.], + [ 0., -1.]]) + >>> dirs.shape # (num_dparams, ndim) + (4, 2) + >>> dirs = geom.det_to_src([0, np.pi / 2, np.pi], 0) + >>> np.allclose(dirs, [[0, -1], + ... [1, 0], + ... [0, 1]]) + True + >>> dirs.shape # (num_angles, ndim) + (3, 2) + >>> # Providing 3 pairs of parameters, resulting in 3 vectors + >>> dirs = geom.det_to_src([0, np.pi / 2, np.pi], [-1, 0, 1]) + >>> dirs[0] # Corresponds to angle = 0, dparam = -1 + array([ 0., -1.]) + >>> dirs.shape + (3, 2) + >>> # Pairs of parameters arranged in arrays of same size + >>> geom.det_to_src(np.zeros((4, 5)), np.zeros((4, 5))).shape + (4, 5, 2) + >>> # "Outer product" type evaluation using broadcasting + >>> geom.det_to_src(np.zeros((4, 1)), np.zeros((1, 5))).shape + (4, 5, 2) """ - if index not in self.motion_params: - raise ValueError('`index` must be contained in `motion_params` ' - '{}, got {}'.format(self.motion_params, index)) - if dparam not in self.det_params: - raise ValueError('`dparam` must be contained in `det_params` ' - '{}, got {}'.format(self.det_params, dparam)) - if not normalized: - raise NotImplementedError('non-normalized detector-to-source ' - 'vector not available for parallel ' - 'beam geometries') - - index = float(index) - int_part = int(index) - frac_part = index - int_part - if int_part == self.motion_params.max_pt: - ray_dir = self.vectors[int_part, self._slice_ray] + squeeze_index = (np.shape(index) == ()) + index_in = index + index = np.array(index, dtype=float, copy=False, ndmin=1) + + dparam_in = dparam + if self.det_params.ndim == 1: + squeeze_dparam = (np.shape(dparam) == ()) + dparam = np.array(dparam, dtype=float, copy=False, ndmin=1) else: - ray_left = self.vectors[int_part, self._slice_ray] - ray_right = self.vectors[int_part + 1, self._slice_ray] - ray_dir = ray_left + frac_part * (ray_right - ray_left) + squeeze_dparam = (np.broadcast(*dparam).shape == ()) + dparam = tuple(np.array(p, dtype=float, copy=False, ndmin=1) + for p in dparam) + + if self.check_bounds: + if not is_inside_bounds(index, self.motion_params): + raise ValueError('`index` {} not in the valid range {}' + ''.format(index_in, self.motion_params)) - return -ray_dir / np.linalg.norm(ray_dir) + if not is_inside_bounds(dparam, self.det_params): + raise ValueError('`dparam` {} not in the valid range {}' + ''.format(dparam_in, self.det_params)) + + at_max_flat = (index == self.motion_params.max_pt).ravel() + + if self.det_params.ndim == 1: + bcast_shape = np.broadcast(index, dparam).shape + else: + bcast_shape = np.broadcast(index, *dparam).shape + + # TODO: shapes are not quite right here. Scalar works, but + # vectorized goes wrong + ray_dir = np.empty(bcast_shape + (self.ndim,), dtype=float) + print('result shape:', ray_dir.shape) + + # Values at rightmost boundary, cannot interpolate due to IndexError + if np.any(at_max_flat): + print('at max bool arr:', at_max_flat) + vec_at_max = self.vectors[-1] + ray_dir[at_max_flat, ...] = vec_at_max[self._slice_ray] + + # Inner part, use linear interpolation + if np.any(~at_max_flat): + print('inner bool arr:', ~at_max_flat) + index_int_part = index[~at_max_flat, ...].astype(int) + index_int_part_flat = index_int_part.ravel() + # Keep shape for broadcasting + index_frac_part = index[~at_max_flat, ...] - index_int_part + # Get ray direction vectors left and right of the given points + vecs_left = np.take(self.vectors, index_int_part_flat, axis=0) + ray_left = vecs_left[..., self._slice_ray] + vecs_right = np.take(self.vectors, index_int_part_flat + 1, axis=0) + ray_right = vecs_right[..., self._slice_ray] + print('ray_left shape:', ray_left.shape) + print('ray_right shape:', ray_right.shape) + print('index_frac_part shape:', index_frac_part.shape) + ray_dir[~at_max_flat, ...] = ( + ray_left + + index_frac_part * (ray_right - ray_left)) + + ray_dir *= -1 / np.linalg.norm(ray_dir, axis=-1, keepdims=True) + + if squeeze_index and squeeze_dparam: + ray_dir = ray_dir.squeeze() + + return ray_dir def parallel_beam_geometry(space, num_angles=None, det_shape=None): From 2d99692aedd15966b99d08870d3642973b4050b6 Mon Sep 17 00:00:00 2001 From: Holger Kohr Date: Mon, 20 Aug 2018 10:30:45 +0200 Subject: [PATCH 04/24] MAINT: fix old-style repr of vec geometry --- odl/tomo/geometry/geometry.py | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/odl/tomo/geometry/geometry.py b/odl/tomo/geometry/geometry.py index 9549e3f91cc..44af6b8ff3a 100644 --- a/odl/tomo/geometry/geometry.py +++ b/odl/tomo/geometry/geometry.py @@ -18,7 +18,8 @@ from odl.tomo.geometry.detector import Detector, Flat1dDetector, Flat2dDetector from odl.tomo.util import axis_rotation_matrix, is_inside_bounds from odl.util import ( - indent_rows, normalized_scalar_param_list, safe_int_conv, signature_string) + REPR_PRECISION, normalized_scalar_param_list, npy_printoptions, + repr_string, safe_int_conv, signature_string_parts) __all__ = ( 'Geometry', @@ -1053,9 +1054,10 @@ def det_axes(self, index): def __repr__(self): """Return ``repr(self)``.""" posargs = [self.det_partition.shape, self.vectors] - inner_str = signature_string(posargs, [], sep=',\n') - return '{}(\n{}\n)'.format(self.__class__.__name__, - indent_rows(inner_str)) + with npy_printoptions(precision=REPR_PRECISION): + inner_parts = signature_string_parts(posargs, []) + return repr_string(self.__class__.__name__, inner_parts, + allow_mixed_seps=False) if __name__ == '__main__': From 195f2dc076d7956b80c1953bbc941205d525df44 Mon Sep 17 00:00:00 2001 From: Holger Kohr Date: Mon, 20 Aug 2018 15:33:26 +0200 Subject: [PATCH 05/24] ENH: add parallel 2d version of geom->vectors converter --- odl/tomo/backends/astra_setup.py | 59 ++++++++++++++++++++++++++++++++ 1 file changed, 59 insertions(+) diff --git a/odl/tomo/backends/astra_setup.py b/odl/tomo/backends/astra_setup.py index 5bed4200b31..5615068d505 100644 --- a/odl/tomo/backends/astra_setup.py +++ b/odl/tomo/backends/astra_setup.py @@ -523,6 +523,65 @@ def astra_parallel_3d_geom_to_vec(geometry): return vecs_odl_axes_to_astra_axes(vectors) +def astra_parallel_2d_geom_to_vec(geometry): + """Create vectors for ASTRA projection geometries from ODL geometry. + + The 2D vectors are used to create an ASTRA projection geometry for + fan beam geometries, see ``'parallel2d_vec'`` in the + `ASTRA projection geometry documentation`_. + + Each row of the returned vectors corresponds to a single projection + and consists of :: + + (rayX, rayY, dX, dY, uX, uY) + + with + + - ``ray``: the ray direction + - ``d`` : the center of the detector + - ``u`` : the vector from detector pixel 0 to 1 + + Parameters + ---------- + geometry : `Geometry` + ODL projection geometry from which to create the ASTRA geometry. + + Returns + ------- + vectors : `numpy.ndarray` + Array of shape ``(num_angles, 6)`` containing the vectors. + + References + ---------- + .. _ASTRA projection geometry documentation: + http://www.astra-toolbox.com/docs/geom2d.html#projection-geometries + """ + # Instead of rotating the data by 90 degrees counter-clockwise, + # we subtract pi/2 from the geometry angles, thereby rotating the + # geometry by 90 degrees clockwise + rot_minus_90 = euler_matrix(-np.pi / 2) + angles = geometry.angles + mid_pt = geometry.det_params.mid_pt + vectors = np.zeros((angles.size, 6)) + + # Ray direction = -(detector-to-source normal vector) + ray = -geometry.det_to_src(angles, mid_pt) + vectors[:, 0:2] = rot_minus_90.dot(ray.T).T # dot along 2nd axis + + # Center of detector + # Need to cast `mid_pt` to float since otherwise the empty axis is + # not removed + centers = geometry.det_point_position(angles, float(mid_pt)) + vectors[:, 2:4] = rot_minus_90.dot(centers.T).T + + # Vector from detector pixel 0 to 1 + det_axis = rot_minus_90.dot(geometry.det_axis(angles).T).T + px_size = geometry.det_partition.cell_sides[0] + vectors[:, 4:6] = det_axis * px_size + + return vecs_odl_axes_to_astra_axes(vectors) + + def astra_projection_geometry(geometry): """Create an ASTRA projection geometry from an ODL geometry object. From 6e8b102ed166320a1eaf677b7e947260018d0767 Mon Sep 17 00:00:00 2001 From: Holger Kohr Date: Tue, 21 Aug 2018 18:40:12 +0200 Subject: [PATCH 06/24] ENH: allow different discr in irrelevant axes in resizing op --- odl/discr/discr_ops.py | 8 +++++--- 1 file changed, 5 insertions(+), 3 deletions(-) diff --git a/odl/discr/discr_ops.py b/odl/discr/discr_ops.py index bebaa195973..0ae01f9b648 100644 --- a/odl/discr/discr_ops.py +++ b/odl/discr/discr_ops.py @@ -320,10 +320,12 @@ def __init__(self, domain, range=None, ran_shp=None, **kwargs): '`ran_shp`') for i in range(domain.ndim): - if (ran.is_uniform_byaxis[i] and + if ( + domain.shape[i] != ran.shape[i] and # axis i relevant + ran.is_uniform_byaxis[i] and domain.is_uniform_byaxis[i] and - not np.isclose(ran.cell_sides[i], - domain.cell_sides[i])): + not np.isclose(ran.cell_sides[i], domain.cell_sides[i]) + ): raise ValueError( 'in axis {}: cell sides of domain and range differ ' 'significantly: (difference {})' From 36f09112407de98ca7fce49f6784adc541493156 Mon Sep 17 00:00:00 2001 From: Holger Kohr Date: Tue, 21 Aug 2018 18:44:28 +0200 Subject: [PATCH 07/24] WIP: fix adjoint weighting in par2d vec geometry --- odl/tomo/backends/astra_cuda.py | 14 ++++++++++---- 1 file changed, 10 insertions(+), 4 deletions(-) diff --git a/odl/tomo/backends/astra_cuda.py b/odl/tomo/backends/astra_cuda.py index 6dd315e091c..52eb8fda150 100644 --- a/odl/tomo/backends/astra_cuda.py +++ b/odl/tomo/backends/astra_cuda.py @@ -239,8 +239,14 @@ def _call_forward_real(self, vol_data, out=None, **kwargs): # Fix scaling to weight by pixel size if ( - isinstance(self.geometry, Parallel2dGeometry) - and parse_version(ASTRA_VERSION) < parse_version('1.9.9.dev') + parse_version(ASTRA_VERSION) < parse_version('1.9.9.dev') + and ( + isinstance(self.geometry, Parallel2dGeometry) + or ( + isinstance(self.geometry, ParallelVecGeometry) + and self.geometry.ndim == 2 + ) + ) ): # parallel2d scales with pixel stride out *= 1 / float(self.geometry.det_partition.cell_volume) @@ -496,8 +502,8 @@ def astra_cuda_bp_scaling_factor(proj_space, vol_space, geometry): det_px_area = geometry.det_partition.cell_volume scaling_factor *= (src_radius ** 2 * det_px_area ** 2) elif isinstance(geometry, VecGeometry): - scaling_factor = (geometry.det_partition.cell_volume / - vol_space.cell_volume) + # TODO: correct in 1.8? With differently weighted spaces as well? + scaling_factor = vol_space.cell_volume elif isinstance(geometry, ParallelVecGeometry): if geometry.ndim == 2: # Scales with 1 / cell_volume From f891d9d126feb5e1fb733b92f33ec287216d3e59 Mon Sep 17 00:00:00 2001 From: Holger Kohr Date: Tue, 21 Aug 2018 18:45:29 +0200 Subject: [PATCH 08/24] ENH: add support for ASTRA parallel_vec --- odl/tomo/backends/astra_setup.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/odl/tomo/backends/astra_setup.py b/odl/tomo/backends/astra_setup.py index 5615068d505..78ca25b9ef9 100644 --- a/odl/tomo/backends/astra_setup.py +++ b/odl/tomo/backends/astra_setup.py @@ -632,7 +632,7 @@ def astra_projection_geometry(geometry): elif isinstance(geometry, ParallelVecGeometry) and geometry.ndim == 2: det_count = geometry.detector.size vec = geometry.vectors - if not astra_supports('parallel2d_vec_geometry'): + if not astra_supports('par2d_vec_geometry'): raise NotImplementedError( "'parallel_vec' geometry not supported by ASTRA " 'v{}'.format(ASTRA_VERSION)) From 567bef667c91385d607e435d8627b3a00b8fc01f Mon Sep 17 00:00:00 2001 From: Holger Kohr Date: Tue, 21 Aug 2018 18:45:56 +0200 Subject: [PATCH 09/24] WIP: add axis check examples for vec geometries --- .../checks/check_axes_parallel2d_vec_fp.py | 95 +++++++++++++++++++ 1 file changed, 95 insertions(+) create mode 100644 examples/tomo/checks/check_axes_parallel2d_vec_fp.py diff --git a/examples/tomo/checks/check_axes_parallel2d_vec_fp.py b/examples/tomo/checks/check_axes_parallel2d_vec_fp.py new file mode 100644 index 00000000000..c69acd62b8c --- /dev/null +++ b/examples/tomo/checks/check_axes_parallel2d_vec_fp.py @@ -0,0 +1,95 @@ +"""Parallel2D_vec example for checking that orientations are handled correctly. + +Due to differing axis conventions between ODL and the ray transform +back-ends, a check is needed to confirm that the translation steps are +done correctly. + +Both pairs of plots of ODL projections and NumPy axis sums should look +the same in the sense that they should show the same features in the +right arrangement (not flipped, rotated, etc.). + +This example is best run in Spyder section-by-section (CTRL-Enter). +""" + +# %% Set up the things that never change + +import matplotlib.pyplot as plt +import numpy as np +import odl + +# Set back-end here (for `None` the fastest available is chosen) +impl = None +# Set a volume shift. This should move the projections in the same direction. +shift = np.array([0.0, 25.0]) + +img_shape = (100, 150) +img_max_pt = np.array(img_shape, dtype=float) / 2 +img_min_pt = -img_max_pt +reco_space = odl.uniform_discr(img_min_pt + shift, img_max_pt + shift, + img_shape, dtype='float32') +phantom = odl.phantom.indicate_proj_axis(reco_space) + +assert np.allclose(reco_space.cell_sides, 1) + +# Generate ASTRA vectors, but in the ODL geometry convention. +# We use 0, 90, 180 and 270 degrees as angles, with the detector starting +# with reference point (0, 1) and axis (1, 0), rotating counter-clockwise. +# The vector to the detector reference point coincides with the ray direction. +# The detector `u` vector from pixel 0 to pixel 1 is equal to the detector +# axis, since we choose the pixel size to be equal to 1. +det_refpoints = np.array([(0, 1), (-1, 0), (0, -1), (1, 0)]) +ray_dirs = det_refpoints +det_axes = np.array([(1, 0), (0, 1), (-1, 0), (0, -1)]) +vectors = np.empty((4, 6)) +vectors[:, 0:2] = ray_dirs +vectors[:, 2:4] = det_refpoints +vectors[:, 4:6] = det_axes + +# Choose enough pixels to cover the object projections +det_size = np.floor(1.1 * np.sqrt(np.sum(np.square(img_shape)))) +det_shape = int(det_size) + +# Sum manually using Numpy +sum_along_x = np.sum(phantom, axis=0) +sum_along_y = np.sum(phantom, axis=1) + + +# %% Test forward projection along y axis + + +geometry = odl.tomo.ParallelVecGeometry(det_shape, vectors) +# Check initial configuration +assert np.allclose(geometry.det_axes(0), [[1, 0]]) +assert np.allclose(geometry.det_refpoint(0), [[0, 1]]) + +# Create projections +ray_trafo = odl.tomo.RayTransform(reco_space, geometry, impl=impl) +proj_data = ray_trafo(phantom) + +# Axis in this image is x. This corresponds to 0 degrees. +proj_data.show(indices=[0, None], + title='Projection at 0 degrees ~ Sum along y axis') +fig, ax = plt.subplots() +ax.plot(sum_along_y) +ax.set_xlabel('x') +plt.title('Sum along y axis') +plt.show() +# Check axes in geometry +axis_sum_y = geometry.det_axis(np.deg2rad(0)) +assert np.allclose(axis_sum_y, [1, 0]) + + +# %% Test forward projection along x axis + + +# Axis in this image is y. This corresponds to 90 degrees. +proj_data.show(indices=[1, None], + title='Projection at 90 degrees ~ Sum along x axis') +fig, ax = plt.subplots() +ax.plot(sum_along_x) +ax.set_xlabel('y') +plt.title('Sum along x axis') +plt.show() +# Check axes in geometry +axis_sum_x = geometry.det_axis(np.deg2rad(90)) +assert np.allclose(axis_sum_x, [0, 1]) From 87ae3f1a8f3ffb50a6eb1660e7eaee13f1474422 Mon Sep 17 00:00:00 2001 From: Holger Kohr Date: Wed, 22 Aug 2018 17:13:04 +0200 Subject: [PATCH 10/24] ENH: fix coord sys handling in vec geometries --- examples/tomo/checks/README.md | 19 +- .../tomo/checks/check_axes_cone2d_vec_fp.py | 98 ++++ .../checks/check_axes_parallel2d_vec_bp.py | 42 ++ .../checks/check_axes_parallel2d_vec_fp.py | 8 +- odl/test/tomo/backends/astra_setup_test.py | 35 +- odl/tomo/backends/astra_setup.py | 472 ++++++++++-------- odl/tomo/geometry/geometry.py | 30 +- 7 files changed, 448 insertions(+), 256 deletions(-) create mode 100644 examples/tomo/checks/check_axes_cone2d_vec_fp.py create mode 100644 examples/tomo/checks/check_axes_parallel2d_vec_bp.py diff --git a/examples/tomo/checks/README.md b/examples/tomo/checks/README.md index 480849c2f4b..6adfa8a1b89 100644 --- a/examples/tomo/checks/README.md +++ b/examples/tomo/checks/README.md @@ -7,11 +7,14 @@ They are primarily intended for developers who change back-end code or introduce Example | Purpose | Complexity ------- | ------- | ---------- -[`check_axes_cone2d_bp.py`](https://github.com/odlgroup/odl/blob/master/examples/tomo/checks/check_axes_cone2d_bp.py) | Check axis conventions and shifts for back-projection in 2D cone beam (= fan beam) geometry | middle -[`check_axes_cone2d_fp.py`](https://github.com/odlgroup/odl/blob/master/examples/tomo/checks/check_axes_cone2d_fp.py) | Check axis conventions and shifts for forward projection in 2D cone beam (= fan beam) geometry | middle -[`check_axes_cone3d_bp.py`](https://github.com/odlgroup/odl/blob/master/examples/tomo/checks/check_axes_cone3d_bp.py) | Check axis conventions and shifts for back-projection in 3D circular cone beam geometry | high -[`check_axes_cone3d_fp.py`](https://github.com/odlgroup/odl/blob/master/examples/tomo/checks/check_axes_cone3d_fp.py) | Check axis conventions and shifts for forward projection in 3D circular cone beam geometry | high -[`check_axes_parallel2d_bp.py`](https://github.com/odlgroup/odl/blob/master/examples/tomo/checks/check_axes_parallel2d_bp.py) | Check axis conventions and shifts for back-projection in 2D parallel beam geometry | middle -[`check_axes_parallel2d_fp.py`](https://github.com/odlgroup/odl/blob/master/examples/tomo/checks/check_axes_parallel2d_fp.py) | Check axis conventions and shifts for forward projection in 2D parallel beam geometry | middle -[`check_axes_parallel3d_bp.py`](https://github.com/odlgroup/odl/blob/master/examples/tomo/checks/check_axes_parallel3d_bp.py) | Check axis conventions and shifts for back-projection in 3D parallel beam geometry | high -[`check_axes_parallel3d_fp.py`](https://github.com/odlgroup/odl/blob/master/examples/tomo/checks/check_axes_parallel3d_fp.py) | Check axis conventions and shifts for forward projection in 3D parallel beam geometry | high +[`check_axes_cone2d_bp.py`](check_axes_cone2d_bp.py) | Check axis conventions and shifts for back-projection in 2D cone beam (= fan beam) geometry | middle +[`check_axes_cone2d_fp.py`](check_axes_cone2d_fp.py) | Check axis conventions and shifts for forward projection in 2D cone beam (= fan beam) geometry | middle +[`check_axes_cone2d_vec_fp.py`](check_axes_cone2d_vec_fp.py) | Check axis conventions and shifts for forward projection in 2D cone beam "vec" geometry | middle +[`check_axes_cone3d_bp.py`](check_axes_cone3d_bp.py) | Check axis conventions and shifts for back-projection in 3D circular cone beam geometry | high +[`check_axes_cone3d_fp.py`](check_axes_cone3d_fp.py) | Check axis conventions and shifts for forward projection in 3D circular cone beam geometry | high +[`check_axes_parallel2d_bp.py`](check_axes_parallel2d_bp.py) | Check axis conventions and shifts for back-projection in 2D parallel beam geometry | middle +[`check_axes_parallel2d_fp.py`](check_axes_parallel2d_fp.py) | Check axis conventions and shifts for forward projection in 2D parallel beam geometry | middle +[`check_axes_parallel2d_vec_bp.py`](check_axes_parallel2d_vec_bp.py) | Check axis conventions and shifts for back-projection in 2D parallel beam "vec" geometry | middle +[`check_axes_parallel2d_vec_fp.py`](check_axes_parallel2d_vec_fp.py) | Check axis conventions and shifts for forward projection in 2D parallel beam "vec" geometry | middle +[`check_axes_parallel3d_bp.py`](check_axes_parallel3d_bp.py) | Check axis conventions and shifts for back-projection in 3D parallel beam geometry | high +[`check_axes_parallel3d_fp.py`](check_axes_parallel3d_fp.py) | Check axis conventions and shifts for forward projection in 3D parallel beam geometry | high diff --git a/examples/tomo/checks/check_axes_cone2d_vec_fp.py b/examples/tomo/checks/check_axes_cone2d_vec_fp.py new file mode 100644 index 00000000000..947b74c174c --- /dev/null +++ b/examples/tomo/checks/check_axes_cone2d_vec_fp.py @@ -0,0 +1,98 @@ +"""Cone2D_vec example for checking that orientations are handled correctly. + +Due to differing axis conventions between ODL and the ray transform +back-ends, a check is needed to confirm that the translation steps are +done correctly. + +Both pairs of plots of ODL projections and NumPy axis sums should look +the same in the sense that they should show the same features in the +right arrangement (not flipped, rotated, etc.). + +This example is best run in Spyder section-by-section (CTRL-Enter). +""" + +# %% Set up the things that never change + +import matplotlib.pyplot as plt +import numpy as np +import odl + +# Set back-end here (for `None` the fastest available is chosen) +impl = None +# Set a volume shift. This should move the projections in the same direction. +shift = np.array([0.0, 25.0]) + +img_shape = (100, 150) +img_max_pt = np.array(img_shape, dtype=float) / 2 +img_min_pt = -img_max_pt +reco_space = odl.uniform_discr(img_min_pt + shift, img_max_pt + shift, + img_shape, dtype='float32') +phantom = odl.phantom.indicate_proj_axis(reco_space) + +assert np.allclose(reco_space.cell_sides, 1) + +# Generate ASTRA vectors, but in the ODL geometry convention. +# We use 0, 90, 180 and 270 degrees as angles, with the detector starting +# with reference point (0, 1000) and axis (1, 0), rotating counter-clockwise. +# The source points are chosen opposite to the detector at radius 500. +# The detector `u` vector from pixel 0 to pixel 1 is equal to the detector +# axis, since we choose the pixel size to be equal to 1. +src_radius = 500 +det_radius = 1000 +det_refpoints = np.array([(0, 1), (-1, 0), (0, -1), (1, 0)]) * det_radius +src_points = -det_refpoints / det_radius * src_radius +det_axes = np.array([(1, 0), (0, 1), (-1, 0), (0, -1)]) +vectors = np.empty((4, 6)) +vectors[:, 0:2] = src_points +vectors[:, 2:4] = det_refpoints +vectors[:, 4:6] = det_axes + +# Choose enough pixels to cover the object projections +fan_angle = np.arctan(img_max_pt[1] / src_radius) +det_size = np.floor(2 * (src_radius + det_radius) * np.sin(fan_angle)) +det_shape = int(det_size) + +# Sum manually using Numpy +sum_along_x = np.sum(phantom, axis=0) +sum_along_y = np.sum(phantom, axis=1) + + +# %% Test forward projection along y axis + + +geometry = odl.tomo.ConeVecGeometry(det_shape, vectors) +# Check initial configuration +assert np.allclose(geometry.det_axes(0), [[1, 0]]) +assert np.allclose(geometry.det_refpoint(0), [[0, det_radius]]) + +# Create projections +ray_trafo = odl.tomo.RayTransform(reco_space, geometry, impl=impl) +proj_data = ray_trafo(phantom) + +# Axis in this image is x. This corresponds to 0 degrees or index 0. +proj_data.show(indices=[0, None], + title='Projection at 0 degrees ~ Sum along y axis') +fig, ax = plt.subplots() +ax.plot(sum_along_y) +ax.set_xlabel('x') +plt.title('Sum along y axis') +plt.show() +# Check axes in geometry +axis_sum_y = geometry.det_axis(0) +assert np.allclose(axis_sum_y, [1, 0]) + + +# %% Test forward projection along x axis + + +# Axis in this image is y. This corresponds to 90 degrees or index 1. +proj_data.show(indices=[1, None], + title='Projection at 90 degrees ~ Sum along x axis') +fig, ax = plt.subplots() +ax.plot(sum_along_x) +ax.set_xlabel('y') +plt.title('Sum along x axis') +plt.show() +# Check axes in geometry +axis_sum_x = geometry.det_axis(1) +assert np.allclose(axis_sum_x, [0, 1]) diff --git a/examples/tomo/checks/check_axes_parallel2d_vec_bp.py b/examples/tomo/checks/check_axes_parallel2d_vec_bp.py new file mode 100644 index 00000000000..31904413ec6 --- /dev/null +++ b/examples/tomo/checks/check_axes_parallel2d_vec_bp.py @@ -0,0 +1,42 @@ +"""Parallel2D_vec example for checking that orientations are handled correctly. + +Due to differing axis conventions between ODL and the ray transform +back-ends, a check is needed to confirm that the translation steps are +done correctly. + +The back-projected data should be a blurry version of the phantom, with +all features in the correct positions, not flipped or rotated. +""" + +import numpy as np +import odl + +# Set back-end here (for `None` the fastest available is chosen) +impl = None +# Set a volume shift. This should not have any influence on the back-projected +# data. +shift = np.array([0.0, 25.0]) + +img_shape = (100, 150) +img_max_pt = np.array(img_shape, dtype=float) / 2 +img_min_pt = -img_max_pt +reco_space = odl.uniform_discr(img_min_pt + shift, img_max_pt + shift, + img_shape, dtype='float32') +phantom = odl.phantom.indicate_proj_axis(reco_space) + +assert np.allclose(reco_space.cell_sides, 1) + +# Make standard parallel beam geometry with 360 angles and cast it to a +# vec geometry +geometry = odl.tomo.parallel_beam_geometry(reco_space, num_angles=360) +geometry = odl.tomo.ParallelVecGeometry( + det_shape=geometry.det_partition.shape, + vectors=odl.tomo.parallel_2d_geom_to_astra_vecs(geometry) +) + +# Test back-projection +ray_trafo = odl.tomo.RayTransform(reco_space, geometry, impl=impl) +proj_data = ray_trafo(phantom) +backproj = ray_trafo.adjoint(proj_data) +backproj.show('Back-projection') +phantom.show('Phantom') diff --git a/examples/tomo/checks/check_axes_parallel2d_vec_fp.py b/examples/tomo/checks/check_axes_parallel2d_vec_fp.py index c69acd62b8c..655270bf365 100644 --- a/examples/tomo/checks/check_axes_parallel2d_vec_fp.py +++ b/examples/tomo/checks/check_axes_parallel2d_vec_fp.py @@ -66,7 +66,7 @@ ray_trafo = odl.tomo.RayTransform(reco_space, geometry, impl=impl) proj_data = ray_trafo(phantom) -# Axis in this image is x. This corresponds to 0 degrees. +# Axis in this image is x. This corresponds to 0 degrees or index 0. proj_data.show(indices=[0, None], title='Projection at 0 degrees ~ Sum along y axis') fig, ax = plt.subplots() @@ -75,14 +75,14 @@ plt.title('Sum along y axis') plt.show() # Check axes in geometry -axis_sum_y = geometry.det_axis(np.deg2rad(0)) +axis_sum_y = geometry.det_axis(0) assert np.allclose(axis_sum_y, [1, 0]) # %% Test forward projection along x axis -# Axis in this image is y. This corresponds to 90 degrees. +# Axis in this image is y. This corresponds to 90 degrees or index 1. proj_data.show(indices=[1, None], title='Projection at 90 degrees ~ Sum along x axis') fig, ax = plt.subplots() @@ -91,5 +91,5 @@ plt.title('Sum along x axis') plt.show() # Check axes in geometry -axis_sum_x = geometry.det_axis(np.deg2rad(90)) +axis_sum_x = geometry.det_axis(1) assert np.allclose(axis_sum_x, [0, 1]) diff --git a/odl/test/tomo/backends/astra_setup_test.py b/odl/test/tomo/backends/astra_setup_test.py index d3502cf5bc2..ac11c428718 100644 --- a/odl/test/tomo/backends/astra_setup_test.py +++ b/odl/test/tomo/backends/astra_setup_test.py @@ -17,7 +17,6 @@ from odl.tomo.backends.astra_setup import ( astra_algorithm, astra_data, astra_projection_geometry, astra_projector, astra_supports, astra_volume_geometry) -from odl.util.testutils import is_subdict try: import astra @@ -151,23 +150,6 @@ def test_vol_geom_3d(): astra_volume_geometry(discr_dom) -def test_proj_geom_parallel_2d(): - """Create ASTRA 2D projection geometry.""" - - apart = odl.uniform_partition(0, 2, 5) - dpart = odl.uniform_partition(-1, 1, 10) - geom = odl.tomo.Parallel2dGeometry(apart, dpart) - - proj_geom = astra_projection_geometry(geom) - - correct_subdict = { - 'type': 'parallel', - 'DetectorCount': 10, 'DetectorWidth': 0.2} - - assert is_subdict(correct_subdict, proj_geom) - assert 'ProjectionAngles' in proj_geom - - def test_astra_projection_geometry(): """Create ASTRA projection geometry from geometry objects.""" @@ -191,7 +173,10 @@ def test_astra_projection_geometry(): # Parallel 2D geometry geom_p2d = odl.tomo.Parallel2dGeometry(apart, dpart) astra_geom = astra_projection_geometry(geom_p2d) - assert astra_geom['type'] == 'parallel' + if odl.tomo.astra_supports('par2d_vec_geometry'): + assert astra_geom['type'] == 'parallel_vec' + else: + assert astra_geom['type'] == 'parallel' # Fan flat src_rad = 10 @@ -386,22 +371,22 @@ def test_geom_to_vec(): src_rad = 10 det_rad = 5 geom_ff = odl.tomo.FanBeamGeometry(apart, dpart, src_rad, det_rad) - vec = odl.tomo.astra_conebeam_2d_geom_to_vec(geom_ff) + vecs = odl.tomo.cone_2d_geom_to_astra_vecs(geom_ff) - assert vec.shape == (apart.size, 6) + assert vecs.shape == (apart.size, 6) # Circular cone flat dpart = odl.uniform_partition([-40, -3], [40, 3], (10, 5)) geom_ccf = odl.tomo.ConeBeamGeometry(apart, dpart, src_rad, det_rad) - vec = odl.tomo.astra_conebeam_3d_geom_to_vec(geom_ccf) - assert vec.shape == (apart.size, 12) + vecs = odl.tomo.cone_3d_geom_to_astra_vecs(geom_ccf) + assert vecs.shape == (apart.size, 12) # Helical cone flat pitch = 1 geom_hcf = odl.tomo.ConeBeamGeometry(apart, dpart, src_rad, det_rad, pitch=pitch) - vec = odl.tomo.astra_conebeam_3d_geom_to_vec(geom_hcf) - assert vec.shape == (apart.size, 12) + vecs = odl.tomo.cone_3d_geom_to_astra_vecs(geom_hcf) + assert vecs.shape == (apart.size, 12) if __name__ == '__main__': diff --git a/odl/tomo/backends/astra_setup.py b/odl/tomo/backends/astra_setup.py index 78ca25b9ef9..cb8ec0b25c6 100644 --- a/odl/tomo/backends/astra_setup.py +++ b/odl/tomo/backends/astra_setup.py @@ -65,10 +65,13 @@ 'astra_supports', 'astra_versions_supporting', 'astra_volume_geometry', - 'astra_conebeam_3d_geom_to_vec', - 'astra_conebeam_2d_geom_to_vec', - 'astra_parallel_3d_geom_to_vec', 'astra_projection_geometry', + 'vecs_odl_to_astra_coords', + 'vecs_astra_to_odl_coords', + 'parallel_2d_geom_to_astra_vecs', + 'parallel_3d_geom_to_astra_vecs', + 'cone_2d_geom_to_astra_vecs', + 'cone_3d_geom_to_astra_vecs', 'astra_data', 'astra_projector', 'astra_algorithm', @@ -272,8 +275,139 @@ def astra_volume_geometry(vol_space): return vol_geom -def vecs_odl_axes_to_astra_axes(vecs): - """Convert geometry vectors from ODL axis convention to ASTRA. +def astra_projection_geometry(geometry): + """Create an ASTRA projection geometry from an ODL geometry object. + + Parameters + ---------- + geometry : `Geometry` + ODL projection geometry from which to create the ASTRA geometry. + + Returns + ------- + proj_geom : dict + Dictionary defining the ASTRA projection geometry. + """ + if not isinstance(geometry, Geometry): + raise TypeError('`geometry` {!r} is not a `Geometry` instance' + ''.format(geometry)) + + if 'astra' in geometry.implementation_cache: + # Shortcut, reuse already computed value. + return geometry.implementation_cache['astra'] + + if not geometry.det_partition.is_uniform: + raise ValueError('non-uniform detector sampling is not supported') + + # Parallel 2D + if (isinstance(geometry, ParallelBeamGeometry) and + isinstance(geometry.detector, (Flat1dDetector, Flat2dDetector)) and + geometry.ndim == 2): + det_count = geometry.detector.size + if astra_supports('par2d_vec_geometry'): + vecs = parallel_2d_geom_to_astra_vecs(geometry) + proj_geom = astra.create_proj_geom('parallel_vec', det_count, vecs) + else: + det_width = geometry.det_partition.cell_sides[0] + # Instead of rotating the data by 90 degrees counter-clockwise, + # we subtract pi/2 from the geometry angles, thereby rotating the + # geometry by 90 degrees clockwise + angles = geometry.angles - np.pi / 2 + proj_geom = astra.create_proj_geom( + 'parallel', det_width, det_count, angles) + + # Parallel 2D vec + elif isinstance(geometry, ParallelVecGeometry) and geometry.ndim == 2: + if not astra_supports('par2d_vec_geometry'): + raise NotImplementedError( + "'parallel_vec' geometry not supported by ASTRA {}" + ''.format(ASTRA_VERSION)) + det_count = geometry.detector.size + vecs = vecs_odl_to_astra_coords(geometry.vectors) + proj_geom = astra.create_proj_geom('parallel_vec', det_count, vecs) + + # Cone 2D (aka fan beam) + elif (isinstance(geometry, DivergentBeamGeometry) and + isinstance(geometry.detector, (Flat1dDetector, Flat2dDetector)) and + geometry.ndim == 2): + det_count = geometry.detector.size + vecs = cone_2d_geom_to_astra_vecs(geometry) + proj_geom = astra.create_proj_geom('fanflat_vec', det_count, vecs) + + # Cone 2D vec + elif isinstance(geometry, ConeVecGeometry) and geometry.ndim == 2: + det_count = geometry.detector.size + vecs = vecs_odl_to_astra_coords(geometry.vectors) + proj_geom = astra.create_proj_geom('fanflat_vec', det_count, vecs) + + # Parallel 3D + elif (isinstance(geometry, ParallelBeamGeometry) and + isinstance(geometry.detector, (Flat1dDetector, Flat2dDetector)) and + geometry.ndim == 3): + # Swap detector axes (see `parallel_3d_geom_to_astra_vecs`) + det_row_count = geometry.det_partition.shape[0] + det_col_count = geometry.det_partition.shape[1] + vecs = parallel_3d_geom_to_astra_vecs(geometry) + proj_geom = astra.create_proj_geom( + 'parallel3d_vec', det_row_count, det_col_count, vecs) + + # Parallel 3D vec + elif isinstance(geometry, ParallelVecGeometry) and geometry.ndim == 3: + det_row_count = geometry.det_partition.shape[1] + det_col_count = geometry.det_partition.shape[0] + vecs = vecs_odl_to_astra_coords(geometry.vectors) + proj_geom = astra.create_proj_geom('parallel3d_vec', det_row_count, + det_col_count, vecs) + + # Cone 3D + elif (isinstance(geometry, DivergentBeamGeometry) and + isinstance(geometry.detector, (Flat1dDetector, Flat2dDetector)) and + geometry.ndim == 3): + # Swap detector axes (see `conebeam_3d_geom_to_astra_vecs`) + det_row_count = geometry.det_partition.shape[0] + det_col_count = geometry.det_partition.shape[1] + vecs = cone_3d_geom_to_astra_vecs(geometry) + proj_geom = astra.create_proj_geom( + 'cone_vec', det_row_count, det_col_count, vecs) + + # Cone 3D vec + elif isinstance(geometry, ConeVecGeometry) and geometry.ndim == 3: + det_row_count = geometry.det_partition.shape[1] + det_col_count = geometry.det_partition.shape[0] + vecs = vecs_odl_to_astra_coords(geometry.vectors) + proj_geom = astra.create_proj_geom('cone_vec', det_row_count, + det_col_count, vecs) + + else: + raise NotImplementedError('unknown ASTRA geometry type {!r}' + ''.format(geometry)) + + if 'astra' not in geometry.implementation_cache: + # Save computed value for later + geometry.implementation_cache['astra'] = proj_geom + + return proj_geom + + +def vecs_odl_to_astra_coords(vecs): + """Convert geometry vectors from ODL coordinates to ASTRA. + + We have the following relation between ODL and ASTRA coordinates: + + +----------------+----------------+ + | 2D | 3D | + +-------+--------+-------+--------+ + | ODL | ASTRA | ODL | ASTRA | + +=======+========+=======+========+ + | ``X`` | ``-Y`` | ``X`` | ``Z`` | + +-------+--------+-------+--------+ + | ``Y`` | ``X`` | ``Y`` | ``Y`` | + +-------+--------+-------+--------+ + | | | ``Z`` | ``X`` | + +-------+--------+-------+--------+ + + Thus, the conversion ODL→ASTRA is a rotation by +90 degrees in 2D, and + an XZ axis swap in 3D. Parameters ---------- @@ -288,11 +422,19 @@ def vecs_odl_axes_to_astra_axes(vecs): astra_vecs : `numpy.ndarray`, same shape as ``vecs`` The converted geometry vectors. """ - vecs = np.asarray(vecs, dtype=float) + vecs = np.array(vecs, dtype=float, copy=True) if vecs.shape[1] == 6: - # 2D geometry, nothing to do since the axes are the same + # 2D geometry + # ODL x = ASTRA -y, ODL y = ASTRA x + # Thus: ODL->ASTRA conversion is a rotation by +90 degrees + rot_90 = euler_matrix(np.pi / 2) + # Bulk matrix-vector product + vecs[:, 0:2] = np.einsum('kj,ik->ij', rot_90, vecs[:, 0:2]) + vecs[:, 2:4] = np.einsum('kj,ik->ij', rot_90, vecs[:, 2:4]) + vecs[:, 4:6] = np.einsum('kj,ik->ij', rot_90, vecs[:, 4:6]) return vecs + elif vecs.shape[1] == 12: # 3D geometry # ASTRA has (z, y, x) axis convention, in contrast to (x, y, z) in ODL, @@ -301,13 +443,31 @@ def vecs_odl_axes_to_astra_axes(vecs): for i in range(4): newind.extend([2 + 3 * i, 1 + 3 * i, 0 + 3 * i]) return vecs[:, newind] + else: raise ValueError('`vecs` must have shape (N, 6) or (N, 12), got ' 'array with shape {}'.format(vecs.shape)) -def vecs_astra_axes_to_odl_axes(vecs): - """Convert geometry vectors from ASTRA axis convention to ODL. +def vecs_astra_to_odl_coords(vecs): + """Convert geometry vectors from ASTRA coordinates to ODL. + + We have the following relation between ODL and ASTRA coordinates: + + +----------------+----------------+ + | 2D | 3D | + +-------+--------+-------+--------+ + | ODL | ASTRA | ODL | ASTRA | + +=======+========+=======+========+ + | ``X`` | ``-Y`` | ``X`` | ``Z`` | + +-------+--------+-------+--------+ + | ``Y`` | ``X`` | ``Y`` | ``Y`` | + +-------+--------+-------+--------+ + | | | ``Z`` | ``X`` | + +-------+--------+-------+--------+ + + Thus, the conversion ASTRA→ODL is a rotation by -90 degrees in 2D, and + an XZ axis swap in 3D. Parameters ---------- @@ -319,14 +479,22 @@ def vecs_astra_axes_to_odl_axes(vecs): Returns ------- - odl_vecs : `numpy.ndarray`, same shape as ``vecs`` + astra_vecs : `numpy.ndarray`, same shape as ``vecs`` The converted geometry vectors. """ vecs = np.asarray(vecs, dtype=float) if vecs.shape[1] == 6: - # 2D geometry, nothing to do since the axes are the same + # 2D geometry + # ODL x = ASTRA -y, ODL y = ASTRA x + # Thus: ODL->ASTRA conversion is a rotation by -90 degrees + rot_minus_90 = euler_matrix(-np.pi / 2) + # Bulk matrix-vector product + vecs[:, 0:2] = np.einsum('kj,ik->ij', rot_minus_90, vecs[:, 0:2]) + vecs[:, 2:4] = np.einsum('kj,ik->ij', rot_minus_90, vecs[:, 2:4]) + vecs[:, 4:6] = np.einsum('kj,ik->ij', rot_minus_90, vecs[:, 4:6]) return vecs + elif vecs.shape[1] == 12: # 3D geometry # ASTRA has (z, y, x) axis convention, in contrast to (x, y, z) in ODL, @@ -341,85 +509,28 @@ def vecs_astra_axes_to_odl_axes(vecs): 'array with shape {}'.format(vecs.shape)) -def astra_conebeam_3d_geom_to_vec(geometry): - """Create vectors for ASTRA projection geometries from ODL geometry. - - The 3D vectors are used to create an ASTRA projection geometry for - cone beam geometries, see ``'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) - - with - - - ``src``: the ray source position - - ``d`` : the center of the detector - - ``u`` : the vector from detector pixel ``(0,0)`` to ``(0,1)`` - - ``v`` : the vector from detector pixel ``(0,0)`` to ``(1,0)`` - - Parameters - ---------- - geometry : `Geometry` - ODL projection geometry from which to create the ASTRA geometry. - - Returns - ------- - vectors : `numpy.ndarray` - Array of shape ``(num_angles, 12)`` 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, 12)) - - # Source position - vectors[:, 0:3] = geometry.src_position(angles) - - # Center of detector in 3D space - mid_pt = geometry.det_params.mid_pt - vectors[:, 3:6] = geometry.det_point_position(angles, mid_pt) - - # Vectors from detector pixel (0, 0) to (1, 0) and (0, 0) to (0, 1) - # `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 - # Swap detector axes to have better memory layout in projection data. - # ASTRA produces `(v, theta, u)` layout, and to map to ODL layout - # `(theta, u, v)` a complete roll must be performed, which is the - # worst case (compeltely discontiguous). - # Instead we swap `u` and `v`, resulting in the effective ASTRA result - # `(u, theta, v)`. Here we only need to swap axes 0 and 1, which - # keeps at least contiguous blocks in `v`. - vectors[:, 9:12] = det_axes[0] * px_sizes[0] - vectors[:, 6:9] = det_axes[1] * px_sizes[1] - - return vecs_odl_axes_to_astra_axes(vectors) - - -def astra_conebeam_2d_geom_to_vec(geometry): - """Create vectors for ASTRA projection geometries from ODL geometry. +def parallel_2d_geom_to_astra_vecs(geometry): + """Create vectors for ASTRA projection geometries from an ODL geometry. The 2D vectors are used to create an ASTRA projection geometry for - fan beam geometries, see ``'fanflat_vec'`` in the + parallel beam geometries, see ``'parallel_vec'`` in the `ASTRA projection geometry documentation`_. Each row of the returned vectors corresponds to a single projection and consists of :: - (srcX, srcY, dX, dY, uX, uY) + (rayX, rayY, dX, dY, uX, uY) with - - ``src``: the ray source position + - ``ray``: the ray direction - ``d`` : the center of the detector - ``u`` : the vector from detector pixel 0 to 1 + .. note:: + The returned vectors are the result of converting from ODL coordinates + to ASTRA coordinates, see `vecs_odl_to_astra_coords` for details. + Parameters ---------- geometry : `Geometry` @@ -435,34 +546,27 @@ def astra_conebeam_2d_geom_to_vec(geometry): .. _ASTRA projection geometry documentation: http://www.astra-toolbox.com/docs/geom2d.html#projection-geometries """ - # Instead of rotating the data by 90 degrees counter-clockwise, - # we subtract pi/2 from the geometry angles, thereby rotating the - # geometry by 90 degrees clockwise - rot_minus_90 = euler_matrix(-np.pi / 2) angles = geometry.angles + mid_pt = geometry.det_params.mid_pt vectors = np.zeros((angles.size, 6)) - # Source position - src_pos = geometry.src_position(angles) - vectors[:, 0:2] = rot_minus_90.dot(src_pos.T).T # dot along 2nd axis + # Ray direction = -(detector-to-source normal vector) + vectors[:, 0:2] = -geometry.det_to_src(angles, mid_pt) # Center of detector - mid_pt = geometry.det_params.mid_pt - # Need to cast `mid_pt` to float since otherwise the empty axis is - # not removed - centers = geometry.det_point_position(angles, float(mid_pt)) - vectors[:, 2:4] = rot_minus_90.dot(centers.T).T + # Cast `mid_pt` to float to avoid broadcasting + vectors[:, 2:4] = geometry.det_point_position(angles, float(mid_pt)) # Vector from detector pixel 0 to 1 - det_axis = rot_minus_90.dot(geometry.det_axis(angles).T).T + det_axis = geometry.det_axis(angles) px_size = geometry.det_partition.cell_sides[0] vectors[:, 4:6] = det_axis * px_size - return vecs_odl_axes_to_astra_axes(vectors) + return vecs_odl_to_astra_coords(vectors) -def astra_parallel_3d_geom_to_vec(geometry): - """Create vectors for ASTRA projection geometries from ODL geometry. +def parallel_3d_geom_to_astra_vecs(geometry): + """Create vectors for ASTRA projection geometries from an ODL geometry. The 3D vectors are used to create an ASTRA projection geometry for parallel beam geometries, see ``'parallel3d_vec'`` in the @@ -480,6 +584,10 @@ def astra_parallel_3d_geom_to_vec(geometry): - ``u`` : the vector from detector pixel ``(0,0)`` to ``(0,1)`` - ``v`` : the vector from detector pixel ``(0,0)`` to ``(1,0)`` + .. note:: + The returned vectors are the result of converting from ODL coordinates + to ASTRA coordinates, see `vecs_odl_to_astra_coords` for details. + Parameters ---------- geometry : `Geometry` @@ -520,31 +628,35 @@ def astra_parallel_3d_geom_to_vec(geometry): vectors[:, 9:12] = det_axes[0] * px_sizes[0] vectors[:, 6:9] = det_axes[1] * px_sizes[1] - return vecs_odl_axes_to_astra_axes(vectors) + return vecs_odl_to_astra_coords(vectors) -def astra_parallel_2d_geom_to_vec(geometry): - """Create vectors for ASTRA projection geometries from ODL geometry. +def cone_2d_geom_to_astra_vecs(geometry): + """Create vectors for ASTRA projection geometries from an ODL geometry. The 2D vectors are used to create an ASTRA projection geometry for - fan beam geometries, see ``'parallel2d_vec'`` in the + fan beam geometries, see ``'fanflat_vec'`` in the `ASTRA projection geometry documentation`_. Each row of the returned vectors corresponds to a single projection and consists of :: - (rayX, rayY, dX, dY, uX, uY) + (srcX, srcY, dX, dY, uX, uY) with - - ``ray``: the ray direction + - ``src``: the ray source position - ``d`` : the center of the detector - ``u`` : the vector from detector pixel 0 to 1 + .. note:: + The returned vectors are the result of converting from ODL coordinates + to ASTRA coordinates, see `vecs_odl_to_astra_coords` for details. + Parameters ---------- geometry : `Geometry` - ODL projection geometry from which to create the ASTRA geometry. + ODL projection geometry from which to create the ASTRA vectors. Returns ------- @@ -556,136 +668,88 @@ def astra_parallel_2d_geom_to_vec(geometry): .. _ASTRA projection geometry documentation: http://www.astra-toolbox.com/docs/geom2d.html#projection-geometries """ - # Instead of rotating the data by 90 degrees counter-clockwise, - # we subtract pi/2 from the geometry angles, thereby rotating the - # geometry by 90 degrees clockwise - rot_minus_90 = euler_matrix(-np.pi / 2) angles = geometry.angles mid_pt = geometry.det_params.mid_pt vectors = np.zeros((angles.size, 6)) - # Ray direction = -(detector-to-source normal vector) - ray = -geometry.det_to_src(angles, mid_pt) - vectors[:, 0:2] = rot_minus_90.dot(ray.T).T # dot along 2nd axis + # Source position + vectors[:, 0:2] = geometry.src_position(angles) # Center of detector - # Need to cast `mid_pt` to float since otherwise the empty axis is - # not removed - centers = geometry.det_point_position(angles, float(mid_pt)) - vectors[:, 2:4] = rot_minus_90.dot(centers.T).T + # Cast `mid_pt` to float to avoid broadcasting + vectors[:, 2:4] = geometry.det_point_position(angles, float(mid_pt)) # Vector from detector pixel 0 to 1 - det_axis = rot_minus_90.dot(geometry.det_axis(angles).T).T + det_axis = geometry.det_axis(angles) px_size = geometry.det_partition.cell_sides[0] vectors[:, 4:6] = det_axis * px_size - return vecs_odl_axes_to_astra_axes(vectors) - - -def astra_projection_geometry(geometry): - """Create an ASTRA projection geometry from an ODL geometry object. - - As of ASTRA version 1.7, the length values are not required any more to be - rescaled for 3D geometries and non-unit (but isotropic) voxel sizes. + return vecs_odl_to_astra_coords(vectors) - Parameters - ---------- - geometry : `Geometry` - ODL projection geometry from which to create the ASTRA geometry. - Returns - ------- - proj_geom : dict - Dictionary defining the ASTRA projection geometry. - """ - if not isinstance(geometry, Geometry): - raise TypeError('`geometry` {!r} is not a `Geometry` instance' - ''.format(geometry)) +def cone_3d_geom_to_astra_vecs(geometry): + """Create vectors for ASTRA projection geometries from an ODL geometry. - if 'astra' in geometry.implementation_cache: - # Shortcut, reuse already computed value. - return geometry.implementation_cache['astra'] + The 3D vectors are used to create an ASTRA projection geometry for + cone beam geometries, see ``'cone_vec'`` in the + `ASTRA projection geometry documentation`_. - if not geometry.det_partition.is_uniform: - raise ValueError('non-uniform detector sampling is not supported') + Each row of the returned vectors corresponds to a single projection + and consists of :: - if (isinstance(geometry, ParallelBeamGeometry) and - isinstance(geometry.detector, (Flat1dDetector, Flat2dDetector)) and - geometry.ndim == 2): - # TODO: change to parallel_vec when available - det_width = geometry.det_partition.cell_sides[0] - det_count = geometry.detector.size - # Instead of rotating the data by 90 degrees counter-clockwise, - # we subtract pi/2 from the geometry angles, thereby rotating the - # geometry by 90 degrees clockwise - angles = geometry.angles - np.pi / 2 - proj_geom = astra.create_proj_geom('parallel', det_width, det_count, - angles) + (srcX, srcY, srcZ, dX, dY, dZ, uX, uY, uZ, vX, vY, vZ) - elif (isinstance(geometry, DivergentBeamGeometry) and - isinstance(geometry.detector, (Flat1dDetector, Flat2dDetector)) and - geometry.ndim == 2): - det_count = geometry.detector.size - vec = astra_conebeam_2d_geom_to_vec(geometry) - proj_geom = astra.create_proj_geom('fanflat_vec', det_count, vec) + with - elif isinstance(geometry, ParallelVecGeometry) and geometry.ndim == 2: - det_count = geometry.detector.size - vec = geometry.vectors - if not astra_supports('par2d_vec_geometry'): - raise NotImplementedError( - "'parallel_vec' geometry not supported by ASTRA " - 'v{}'.format(ASTRA_VERSION)) - proj_geom = astra.create_proj_geom('parallel_vec', det_count, vec) + - ``src``: the ray source position + - ``d`` : the center of the detector + - ``u`` : the vector from detector pixel ``(0,0)`` to ``(0,1)`` + - ``v`` : the vector from detector pixel ``(0,0)`` to ``(1,0)`` - elif isinstance(geometry, ConeVecGeometry) and geometry.ndim == 2: - det_count = geometry.detector.size - vec = geometry.vectors - proj_geom = astra.create_proj_geom('fanflat_vec', det_count, vec) + .. note:: + The returned vectors are the result of converting from ODL coordinates + to ASTRA coordinates, see `vecs_odl_to_astra_coords` for details. - elif (isinstance(geometry, ParallelBeamGeometry) and - isinstance(geometry.detector, (Flat1dDetector, Flat2dDetector)) and - geometry.ndim == 3): - # Swap detector axes (see astra_*_3d_to_vec) - det_row_count = geometry.det_partition.shape[0] - det_col_count = geometry.det_partition.shape[1] - vec = astra_parallel_3d_geom_to_vec(geometry) - proj_geom = astra.create_proj_geom('parallel3d_vec', det_row_count, - det_col_count, vec) + Parameters + ---------- + geometry : `Geometry` + ODL projection geometry from which to create the ASTRA vectors. - elif (isinstance(geometry, DivergentBeamGeometry) and - isinstance(geometry.detector, (Flat1dDetector, Flat2dDetector)) and - geometry.ndim == 3): - # Swap detector axes (see astra_*_3d_to_vec) - det_row_count = geometry.det_partition.shape[0] - det_col_count = geometry.det_partition.shape[1] - vec = astra_conebeam_3d_geom_to_vec(geometry) - proj_geom = astra.create_proj_geom('cone_vec', det_row_count, - det_col_count, vec) + Returns + ------- + vectors : `numpy.ndarray` + Array of shape ``(num_angles, 12)`` containing the vectors. - elif isinstance(geometry, ParallelVecGeometry) and geometry.ndim == 3: - det_row_count = geometry.det_partition.shape[1] - det_col_count = geometry.det_partition.shape[0] - vec = geometry.vectors - proj_geom = astra.create_proj_geom('parallel3d_vec', det_row_count, - det_col_count, vec) + References + ---------- + .. _ASTRA projection geometry documentation: + http://www.astra-toolbox.com/docs/geom3d.html#projection-geometries + """ + angles = geometry.angles + mid_pt = geometry.det_params.mid_pt + vectors = np.zeros((angles.size, 12)) - elif isinstance(geometry, ConeVecGeometry) and geometry.ndim == 3: - det_row_count = geometry.det_partition.shape[1] - det_col_count = geometry.det_partition.shape[0] - vec = geometry.vectors - proj_geom = astra.create_proj_geom('cone_vec', det_row_count, - det_col_count, vec) + # Source position + vectors[:, 0:3] = geometry.src_position(angles) - else: - raise NotImplementedError('unknown ASTRA geometry type {!r}' - ''.format(geometry)) + # Center of detector in 3D space + vectors[:, 3:6] = geometry.det_point_position(angles, mid_pt) - if 'astra' not in geometry.implementation_cache: - # Save computed value for later - geometry.implementation_cache['astra'] = proj_geom + # Vectors from detector pixel (0, 0) to (1, 0) and (0, 0) to (0, 1) + # `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 + # Swap detector axes to have better memory layout in projection data. + # ASTRA produces `(v, theta, u)` layout, and to map to ODL layout + # `(theta, u, v)` a complete roll must be performed, which is the + # worst case (compeltely discontiguous). + # Instead we swap `u` and `v`, resulting in the effective ASTRA result + # `(u, theta, v)`. Here we only need to swap axes 0 and 1, which + # keeps at least contiguous blocks in `v`. + vectors[:, 9:12] = det_axes[0] * px_sizes[0] + vectors[:, 6:9] = det_axes[1] * px_sizes[1] - return proj_geom + return vecs_odl_to_astra_coords(vectors) def astra_data(astra_geom, datatype, data=None, ndim=2, allow_copy=False): diff --git a/odl/tomo/geometry/geometry.py b/odl/tomo/geometry/geometry.py index 44af6b8ff3a..dec2ffb3846 100644 --- a/odl/tomo/geometry/geometry.py +++ b/odl/tomo/geometry/geometry.py @@ -680,29 +680,28 @@ def __init__(self, det_shape, vectors): Number of detector pixels per axis, can be given as single integer for 2D geometries (with 1D detector). vectors : array-like, shape ``(N, 6)`` or ``(N, 12)`` - Vectors defining the geometric configuration in each - projection. The number ``N`` of rows determines the number - of projections, and the number of columns the spatial - dimension (6 for 2D, 12 for 3D). + Vectors defining the geometric configuration (in ASTRA's + coordinate system) in each projection. The number ``N`` of rows + determines the number of projections, and the number of columns + the spatial dimension (6 for 2D, 12 for 3D). """ self.__vectors = np.asarray(vectors, dtype=float) - if self.vectors.ndim != 2: + if self.__vectors.ndim != 2: raise ValueError('`vectors` must be 2-dimensional, got array ' - 'with ndim = {}'.format(self.vectors.ndim)) - - num_projs = self.vectors.shape[0] + 'with ndim = {}'.format(self.__vectors.ndim)) + num_projs = self.__vectors.shape[0] if num_projs == 0: raise ValueError('`vectors` is empty') - if self.vectors.shape[1] == 6: + if self.__vectors.shape[1] == 6: ndim = 2 - elif self.vectors.shape[1] == 12: + elif self.__vectors.shape[1] == 12: ndim = 3 else: raise ValueError('`vectors` must have 6 or 12 columns, got ' 'array with {} columns' - ''.format(self.vectors.shape[1])) + ''.format(self.__vectors.shape[1])) # Make angle partition with spacing 1 mpart = uniform_partition(0, num_projs - 1, num_projs, @@ -713,7 +712,7 @@ def __init__(self, det_shape, vectors): # sizes during acquisition, but has no effect on the computations. # For those, only the `vectors` are used. if ndim == 2: - det_u = self.vectors[0, 4:6] + det_u = self.__vectors[0, 4:6] det_shape = normalized_scalar_param_list(det_shape, length=1, param_conv=safe_int_conv) det_cell_size = np.linalg.norm(det_u) @@ -723,8 +722,8 @@ def __init__(self, det_shape, vectors): shape=det_shape) detector = Flat1dDetector(det_part, axis=det_u) elif ndim == 3: - det_u = self.vectors[0, 6:9] - det_v = self.vectors[0, 9:12] + det_u = self.__vectors[0, 6:9] + det_v = self.__vectors[0, 9:12] det_shape = normalized_scalar_param_list(det_shape, length=2, param_conv=safe_int_conv) det_cell_sides = np.array( @@ -842,7 +841,8 @@ def det_refpoint(self, index): self._slice_det_center] vecs[not_at_right_bdry, :] = ( pt_left + - frac_part[not_at_right_bdry, None] * (pt_right - pt_left)) + frac_part[not_at_right_bdry, None] * (pt_right - pt_left) + ) return vecs.squeeze() def det_point_position(self, index, dparam): From a05ebd3c9cafb4b8783ec4cc4da3479741d9851d Mon Sep 17 00:00:00 2001 From: Holger Kohr Date: Thu, 23 Aug 2018 16:51:06 +0200 Subject: [PATCH 11/24] ENH: ditch pkg_resources for `pkg_supports` function --- odl/util/utility.py | 15 +++++---------- 1 file changed, 5 insertions(+), 10 deletions(-) diff --git a/odl/util/utility.py b/odl/util/utility.py index 9b5732d91e8..ca41e88dd67 100644 --- a/odl/util/utility.py +++ b/odl/util/utility.py @@ -1402,7 +1402,7 @@ def pkg_supports(feature, pkg_version, pkg_feat_dict): >>> pkg_supports('feat5', '1.0', feat_dict) False """ - from pkg_resources import parse_requirements + from packaging.requirements import Requirement feature = str(feature) pkg_version = str(pkg_version) @@ -1418,16 +1418,11 @@ def pkg_supports(feature, pkg_version, pkg_feat_dict): ver_specs = ['pkg' + supp_ver for supp_ver in supp_versions] # Each parse_requirements list contains only one entry since we specify # only one package - ver_reqs = [list(parse_requirements(ver_spec))[0] - for ver_spec in ver_specs] + ver_reqs = [Requirement(ver_spec) for ver_spec in ver_specs] - # If one of the requirements in the list is met, return True - for req in ver_reqs: - if req.specifier.contains(pkg_version, prereleases=True): - return True - - # No match - return False + # If one of the requirements in the list is met, return True, else False + return any(req.specifier.contains(pkg_version, prereleases=True) + for req in ver_reqs) @contextmanager From bd0a2a547d06ece3cc8d64b8017abd70f9f1c04e Mon Sep 17 00:00:00 2001 From: Holger Kohr Date: Thu, 23 Aug 2018 18:28:01 +0200 Subject: [PATCH 12/24] WIP: add `coords` arg to geom->vec converters --- odl/tomo/backends/astra_setup.py | 70 ++++++++++++++++++++++++++------ 1 file changed, 58 insertions(+), 12 deletions(-) diff --git a/odl/tomo/backends/astra_setup.py b/odl/tomo/backends/astra_setup.py index cb8ec0b25c6..350caa1b961 100644 --- a/odl/tomo/backends/astra_setup.py +++ b/odl/tomo/backends/astra_setup.py @@ -305,7 +305,7 @@ def astra_projection_geometry(geometry): geometry.ndim == 2): det_count = geometry.detector.size if astra_supports('par2d_vec_geometry'): - vecs = parallel_2d_geom_to_astra_vecs(geometry) + vecs = parallel_2d_geom_to_astra_vecs(geometry, coords='ASTRA') proj_geom = astra.create_proj_geom('parallel_vec', det_count, vecs) else: det_width = geometry.det_partition.cell_sides[0] @@ -331,7 +331,7 @@ def astra_projection_geometry(geometry): isinstance(geometry.detector, (Flat1dDetector, Flat2dDetector)) and geometry.ndim == 2): det_count = geometry.detector.size - vecs = cone_2d_geom_to_astra_vecs(geometry) + vecs = cone_2d_geom_to_astra_vecs(geometry, coords='ASTRA') proj_geom = astra.create_proj_geom('fanflat_vec', det_count, vecs) # Cone 2D vec @@ -347,7 +347,7 @@ def astra_projection_geometry(geometry): # Swap detector axes (see `parallel_3d_geom_to_astra_vecs`) det_row_count = geometry.det_partition.shape[0] det_col_count = geometry.det_partition.shape[1] - vecs = parallel_3d_geom_to_astra_vecs(geometry) + vecs = parallel_3d_geom_to_astra_vecs(geometry, coords='ASTRA') proj_geom = astra.create_proj_geom( 'parallel3d_vec', det_row_count, det_col_count, vecs) @@ -366,7 +366,7 @@ def astra_projection_geometry(geometry): # Swap detector axes (see `conebeam_3d_geom_to_astra_vecs`) det_row_count = geometry.det_partition.shape[0] det_col_count = geometry.det_partition.shape[1] - vecs = cone_3d_geom_to_astra_vecs(geometry) + vecs = cone_3d_geom_to_astra_vecs(geometry, coords='ASTRA') proj_geom = astra.create_proj_geom( 'cone_vec', det_row_count, det_col_count, vecs) @@ -509,7 +509,7 @@ def vecs_astra_to_odl_coords(vecs): 'array with shape {}'.format(vecs.shape)) -def parallel_2d_geom_to_astra_vecs(geometry): +def parallel_2d_geom_to_astra_vecs(geometry, coords='ODL'): """Create vectors for ASTRA projection geometries from an ODL geometry. The 2D vectors are used to create an ASTRA projection geometry for @@ -534,18 +534,29 @@ def parallel_2d_geom_to_astra_vecs(geometry): Parameters ---------- geometry : `Geometry` - ODL projection geometry from which to create the ASTRA geometry. + ODL projection geometry from which to create the ASTRA vectors. + coords : {'ODL', 'ASTRA'}, optional + Coordinate system in which the resulting vectors should be + represented. Returns ------- vectors : `numpy.ndarray` Array of shape ``(num_angles, 6)`` containing the vectors. + See Also + -------- + vecs_odl_to_astra_coords : Conversion to ASTRA coordinates + References ---------- .. _ASTRA projection geometry documentation: http://www.astra-toolbox.com/docs/geom2d.html#projection-geometries """ + coords, coords_in = str(coords).upper(), coords + if coords not in ('ODL', 'ASTRA'): + raise ValueError('`coords` {!r} not understood'.format(coords_in)) + angles = geometry.angles mid_pt = geometry.det_params.mid_pt vectors = np.zeros((angles.size, 6)) @@ -562,10 +573,13 @@ def parallel_2d_geom_to_astra_vecs(geometry): px_size = geometry.det_partition.cell_sides[0] vectors[:, 4:6] = det_axis * px_size - return vecs_odl_to_astra_coords(vectors) + if coords == 'ASTRA': + vectors = vecs_odl_to_astra_coords(vectors) + + return vectors -def parallel_3d_geom_to_astra_vecs(geometry): +def parallel_3d_geom_to_astra_vecs(geometry, coords='ODL'): """Create vectors for ASTRA projection geometries from an ODL geometry. The 3D vectors are used to create an ASTRA projection geometry for @@ -591,7 +605,10 @@ def parallel_3d_geom_to_astra_vecs(geometry): Parameters ---------- geometry : `Geometry` - ODL projection geometry from which to create the ASTRA geometry. + ODL projection geometry from which to create the ASTRA vectors. + coords : {'ODL', 'ASTRA'}, optional + Coordinate system in which the resulting vectors should be + represented. Returns ------- @@ -603,6 +620,12 @@ def parallel_3d_geom_to_astra_vecs(geometry): .. _ASTRA projection geometry documentation: http://www.astra-toolbox.com/docs/geom3d.html#projection-geometries """ + coords, coords_in = str(coords).upper(), coords + if coords not in ('ODL', 'ASTRA'): + raise ValueError('`coords` {!r} not understood'.format(coords_in)) + + # TODO: use `coords` + angles = geometry.angles mid_pt = geometry.det_params.mid_pt @@ -625,13 +648,15 @@ def parallel_3d_geom_to_astra_vecs(geometry): # Instead we swap `u` and `v`, resulting in the effective ASTRA result # `(u, theta, v)`. Here we only need to swap axes 0 and 1, which # keeps at least contiguous blocks in `v`. + # + # TODO: this swapping should not happen in this function! vectors[:, 9:12] = det_axes[0] * px_sizes[0] vectors[:, 6:9] = det_axes[1] * px_sizes[1] return vecs_odl_to_astra_coords(vectors) -def cone_2d_geom_to_astra_vecs(geometry): +def cone_2d_geom_to_astra_vecs(geometry, coords='ODL'): """Create vectors for ASTRA projection geometries from an ODL geometry. The 2D vectors are used to create an ASTRA projection geometry for @@ -657,6 +682,9 @@ def cone_2d_geom_to_astra_vecs(geometry): ---------- geometry : `Geometry` ODL projection geometry from which to create the ASTRA vectors. + coords : {'ODL', 'ASTRA'}, optional + Coordinate system in which the resulting vectors should be + represented. Returns ------- @@ -668,6 +696,10 @@ def cone_2d_geom_to_astra_vecs(geometry): .. _ASTRA projection geometry documentation: http://www.astra-toolbox.com/docs/geom2d.html#projection-geometries """ + coords, coords_in = str(coords).upper(), coords + if coords not in ('ODL', 'ASTRA'): + raise ValueError('`coords` {!r} not understood'.format(coords_in)) + angles = geometry.angles mid_pt = geometry.det_params.mid_pt vectors = np.zeros((angles.size, 6)) @@ -684,10 +716,13 @@ def cone_2d_geom_to_astra_vecs(geometry): px_size = geometry.det_partition.cell_sides[0] vectors[:, 4:6] = det_axis * px_size - return vecs_odl_to_astra_coords(vectors) + if coords == 'ASTRA': + vectors = vecs_odl_to_astra_coords(vectors) + return vectors -def cone_3d_geom_to_astra_vecs(geometry): + +def cone_3d_geom_to_astra_vecs(geometry, coords='ODL'): """Create vectors for ASTRA projection geometries from an ODL geometry. The 3D vectors are used to create an ASTRA projection geometry for @@ -714,6 +749,9 @@ def cone_3d_geom_to_astra_vecs(geometry): ---------- geometry : `Geometry` ODL projection geometry from which to create the ASTRA vectors. + coords : {'ODL', 'ASTRA'}, optional + Coordinate system in which the resulting vectors should be + represented. Returns ------- @@ -725,6 +763,12 @@ def cone_3d_geom_to_astra_vecs(geometry): .. _ASTRA projection geometry documentation: http://www.astra-toolbox.com/docs/geom3d.html#projection-geometries """ + coords, coords_in = str(coords).upper(), coords + if coords not in ('ODL', 'ASTRA'): + raise ValueError('`coords` {!r} not understood'.format(coords_in)) + + # TODO: use `coords` + angles = geometry.angles mid_pt = geometry.det_params.mid_pt vectors = np.zeros((angles.size, 12)) @@ -746,6 +790,8 @@ def cone_3d_geom_to_astra_vecs(geometry): # Instead we swap `u` and `v`, resulting in the effective ASTRA result # `(u, theta, v)`. Here we only need to swap axes 0 and 1, which # keeps at least contiguous blocks in `v`. + # + # TODO: this swapping should not happen in this function! vectors[:, 9:12] = det_axes[0] * px_sizes[0] vectors[:, 6:9] = det_axes[1] * px_sizes[1] From d487ce3710f09dca550b9be9dc206fc2929dff3f Mon Sep 17 00:00:00 2001 From: Holger Kohr Date: Mon, 1 Oct 2018 18:20:31 +0200 Subject: [PATCH 13/24] ENH: add VecGeometry.__getitem__ --- odl/tomo/geometry/geometry.py | 68 +++++++++++++++++++++++++++++++++++ 1 file changed, 68 insertions(+) diff --git a/odl/tomo/geometry/geometry.py b/odl/tomo/geometry/geometry.py index dec2ffb3846..e87cb0ec9f5 100644 --- a/odl/tomo/geometry/geometry.py +++ b/odl/tomo/geometry/geometry.py @@ -1051,6 +1051,74 @@ def det_axes(self, index): det_v = det_v_left + frac_part * (det_v_right - det_v_left) return (det_u, det_v) + def __getitem__(self, indices): + """Return ``self[indices]``. + + The returned geometry is the vec geometry gained by slicing + `vectors` with the given ``indices`` along the first axis and + using the original `det_shape`. + + Parameters + ---------- + indices : index expression + Object used to slice `vectors` along the first axis. + + + Raises + ------ + ValueError + If ``indices`` slices along other axes than the first. + + Examples + -------- + >>> vecs_2d = [[0, 1, 0, 1, 1, 0], + ... [0, 1, 0, 1, -1, 0], + ... [-1, 0, -1, 0, 0, 1], + ... [-1, 0, -1, 0, 0, -1]] + >>> det_shape = (10,) + >>> geom = odl.tomo.ParallelVecGeometry(det_shape, vecs_2d) + >>> geom[0] # first angle only + ParallelVecGeometry((10,), array([[ 0., 1., 0., 1., 1., 0.]])) + >>> geom[:3] # first 3 angles + ParallelVecGeometry( + (10,), + array([[ 0., 1., 0., 1., 1., 0.], + [ 0., 1., 0., 1., -1., 0.], + [-1., 0., -1., 0., 0., 1.]]) + ) + >>> geom[::2] # every second angle, starting at the first + ParallelVecGeometry( + (10,), + array([[ 0., 1., 0., 1., 1., 0.], + [-1., 0., -1., 0., 0., 1.]]) + ) + >>> geom[[0, -1]] # first and last + ParallelVecGeometry( + (10,), + array([[ 0., 1., 0., 1., 1., 0.], + [-1., 0., -1., 0., 0., -1.]]) + ) + """ + # Index vectors and make sure that we still have a valid shape + # for ASTRA vecs of the same type + sub_vecs = self.vectors[indices] + if np.isscalar(sub_vecs) or not hasattr(sub_vecs, 'shape'): + raise ValueError( + 'indexing of `vectors` results in scalar or non-array ' + '(type {})'.format(type(sub_vecs))) + + if sub_vecs.ndim == 1 and sub_vecs.shape[0] == self.vectors.shape[1]: + # Reduced to single vector, add missing dimension + sub_vecs = sub_vecs[None, :] + + if sub_vecs.ndim != 2 or sub_vecs.shape[1] != self.vectors.shape[1]: + raise ValueError( + 'indexing of `vectors` results in array of shape {0}, ' + 'expected (N, {1}) or ({1},)' + ''.format(sub_vecs.shape, self.vectors.shape[1])) + + return type(self)(self.detector.shape, sub_vecs) + def __repr__(self): """Return ``repr(self)``.""" posargs = [self.det_partition.shape, self.vectors] From d3a75a84a33e6b7fdf2a1b31faecd377f2794068 Mon Sep 17 00:00:00 2001 From: Holger Kohr Date: Wed, 13 Feb 2019 14:10:05 +0100 Subject: [PATCH 14/24] MAINT: update 3d vec geom example --- examples/tomo/ray_trafo_vec_geom_3d.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/examples/tomo/ray_trafo_vec_geom_3d.py b/examples/tomo/ray_trafo_vec_geom_3d.py index a6d4f076b36..d30316a6753 100644 --- a/examples/tomo/ray_trafo_vec_geom_3d.py +++ b/examples/tomo/ray_trafo_vec_geom_3d.py @@ -20,11 +20,11 @@ angle_partition = odl.uniform_partition(0, 2 * np.pi, 180) # Detector: uniformly sampled, n = (512, 512), min = (-30, -30), max = (30, 30) detector_partition = odl.uniform_partition([-30, -30], [30, 30], [512, 512]) -circle_geometry = odl.tomo.CircularConeFlatGeometry( +circle_geometry = odl.tomo.ConeFlatGeometry( angle_partition, detector_partition, src_radius=1000, det_radius=100, axis=[1, 0, 0]) -circle_vecs = odl.tomo.astra_conebeam_3d_geom_to_vec(circle_geometry) +circle_vecs = odl.tomo.cone_3d_geom_to_astra_vecs(circle_geometry) # Cover the whole volume vertically, somewhat undersampled though vert_shift_min = -22 From e084e957f358d748cdb4b77df5ff5fd541932eeb Mon Sep 17 00:00:00 2001 From: Holger Date: Tue, 5 Nov 2019 14:24:58 +0100 Subject: [PATCH 15/24] Fix imports and names in tomo subpackage --- odl/tomo/__init__.py | 12 +++++++++--- 1 file changed, 9 insertions(+), 3 deletions(-) diff --git a/odl/tomo/__init__.py b/odl/tomo/__init__.py index 94ac8f4913a..58845a61c2f 100644 --- a/odl/tomo/__init__.py +++ b/odl/tomo/__init__.py @@ -13,7 +13,9 @@ from .analytic import * from .backends import ( ASTRA_AVAILABLE, ASTRA_CUDA_AVAILABLE, SKIMAGE_AVAILABLE, - astra_conebeam_2d_geom_to_vec, astra_conebeam_3d_geom_to_vec) + cone_2d_geom_to_astra_vecs, cone_3d_geom_to_astra_vecs, + parallel_2d_geom_to_astra_vecs, parallel_3d_geom_to_astra_vecs, + vecs_astra_to_odl_coords, vecs_odl_to_astra_coords) from .geometry import * from .operators import * from .util import * @@ -27,6 +29,10 @@ 'ASTRA_AVAILABLE', 'ASTRA_CUDA_AVAILABLE', 'SKIMAGE_AVAILABLE', - 'astra_conebeam_2d_geom_to_vec', - 'astra_conebeam_3d_geom_to_vec' + 'vecs_astra_to_odl_coords', + 'vecs_odl_to_astra_coords', + 'parallel_2d_geom_to_astra_vecs', + 'parallel_3d_geom_to_astra_vecs', + 'cone_2d_geom_to_astra_vecs', + 'cone_3d_geom_to_astra_vecs' ) From f21af6b2dc2d7e3c88efddb9c2ec9043602b561d Mon Sep 17 00:00:00 2001 From: Holger Date: Tue, 5 Nov 2019 14:25:14 +0100 Subject: [PATCH 16/24] Fix astra_setup tests --- odl/test/tomo/backends/astra_setup_test.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/odl/test/tomo/backends/astra_setup_test.py b/odl/test/tomo/backends/astra_setup_test.py index ac11c428718..5d5112c518a 100644 --- a/odl/test/tomo/backends/astra_setup_test.py +++ b/odl/test/tomo/backends/astra_setup_test.py @@ -173,7 +173,7 @@ def test_astra_projection_geometry(): # Parallel 2D geometry geom_p2d = odl.tomo.Parallel2dGeometry(apart, dpart) astra_geom = astra_projection_geometry(geom_p2d) - if odl.tomo.astra_supports('par2d_vec_geometry'): + if astra_supports('par2d_vec_geometry'): assert astra_geom['type'] == 'parallel_vec' else: assert astra_geom['type'] == 'parallel' From 6b8eb8bc8651ef3ba93fe816bbc66531c7e43ed4 Mon Sep 17 00:00:00 2001 From: Holger Kohr Date: Sun, 19 Jan 2020 23:03:01 +0100 Subject: [PATCH 17/24] WIP: simplify index interpolation to rounding --- .../tomo/checks/check_axes_cone2d_vec_fp.py | 9 +- odl/tomo/backends/astra_cpu.py | 11 +- odl/tomo/backends/astra_setup.py | 53 ++-- odl/tomo/geometry/conebeam.py | 27 +- odl/tomo/geometry/geometry.py | 257 +++++++++++------- odl/tomo/geometry/parallel.py | 52 ++-- 6 files changed, 249 insertions(+), 160 deletions(-) diff --git a/examples/tomo/checks/check_axes_cone2d_vec_fp.py b/examples/tomo/checks/check_axes_cone2d_vec_fp.py index 947b74c174c..0be2dc41755 100644 --- a/examples/tomo/checks/check_axes_cone2d_vec_fp.py +++ b/examples/tomo/checks/check_axes_cone2d_vec_fp.py @@ -70,12 +70,13 @@ proj_data = ray_trafo(phantom) # Axis in this image is x. This corresponds to 0 degrees or index 0. -proj_data.show(indices=[0, None], - title='Projection at 0 degrees ~ Sum along y axis') +proj_data.show( + indices=[0, None], + title='Projection at 0 degrees ~ Sum along y axis' +) fig, ax = plt.subplots() ax.plot(sum_along_y) -ax.set_xlabel('x') -plt.title('Sum along y axis') +ax.set(xlabel="x", title='Sum along y axis') plt.show() # Check axes in geometry axis_sum_y = geometry.det_axis(0) diff --git a/odl/tomo/backends/astra_cpu.py b/odl/tomo/backends/astra_cpu.py index d844be51da8..299c1e011e9 100644 --- a/odl/tomo/backends/astra_cpu.py +++ b/odl/tomo/backends/astra_cpu.py @@ -20,7 +20,8 @@ astra_volume_geometry) from odl.tomo.backends.util import _add_default_complex_impl from odl.tomo.geometry import ( - DivergentBeamGeometry, Geometry, ParallelBeamGeometry) + ConeVecGeometry, DivergentBeamGeometry, Geometry, ParallelBeamGeometry, + ParallelVecGeometry) from odl.util import writable_array try: @@ -52,16 +53,24 @@ def default_astra_proj_type(geom): - `ParallelBeamGeometry`: ``'linear'`` - `DivergentBeamGeometry`: ``'line_fanflat'`` + - `ParallelVecGeometry`: ``'linear'`` + - `ConeVecGeometry`: ``'line_fanflat'`` In 3D: - `ParallelBeamGeometry`: ``'linear3d'`` - `DivergentBeamGeometry`: ``'linearcone'`` + - `ParallelVecGeometry`: ``'linear3d'`` + - `ConeVecGeometry`: ``'linearcone'`` """ if isinstance(geom, ParallelBeamGeometry): return 'linear' if geom.ndim == 2 else 'linear3d' elif isinstance(geom, DivergentBeamGeometry): return 'line_fanflat' if geom.ndim == 2 else 'linearcone' + elif isinstance(geom, ParallelVecGeometry): + return 'linear' if geom.ndim == 2 else 'linear3d' + elif isinstance(geom, ConeVecGeometry): + return 'line_fanflat' if geom.ndim == 2 else 'linearcone' else: raise TypeError( 'no default exists for {}, `astra_proj_type` must be given ' diff --git a/odl/tomo/backends/astra_setup.py b/odl/tomo/backends/astra_setup.py index 350caa1b961..9078dd8b385 100644 --- a/odl/tomo/backends/astra_setup.py +++ b/odl/tomo/backends/astra_setup.py @@ -300,9 +300,11 @@ def astra_projection_geometry(geometry): raise ValueError('non-uniform detector sampling is not supported') # Parallel 2D - if (isinstance(geometry, ParallelBeamGeometry) and - isinstance(geometry.detector, (Flat1dDetector, Flat2dDetector)) and - geometry.ndim == 2): + if ( + isinstance(geometry, ParallelBeamGeometry) + and isinstance(geometry.detector, (Flat1dDetector, Flat2dDetector)) + and geometry.ndim == 2 + ): det_count = geometry.detector.size if astra_supports('par2d_vec_geometry'): vecs = parallel_2d_geom_to_astra_vecs(geometry, coords='ASTRA') @@ -327,9 +329,11 @@ def astra_projection_geometry(geometry): proj_geom = astra.create_proj_geom('parallel_vec', det_count, vecs) # Cone 2D (aka fan beam) - elif (isinstance(geometry, DivergentBeamGeometry) and - isinstance(geometry.detector, (Flat1dDetector, Flat2dDetector)) and - geometry.ndim == 2): + elif ( + isinstance(geometry, DivergentBeamGeometry) + and isinstance(geometry.detector, (Flat1dDetector, Flat2dDetector)) + and geometry.ndim == 2 + ): det_count = geometry.detector.size vecs = cone_2d_geom_to_astra_vecs(geometry, coords='ASTRA') proj_geom = astra.create_proj_geom('fanflat_vec', det_count, vecs) @@ -341,46 +345,55 @@ def astra_projection_geometry(geometry): proj_geom = astra.create_proj_geom('fanflat_vec', det_count, vecs) # Parallel 3D - elif (isinstance(geometry, ParallelBeamGeometry) and - isinstance(geometry.detector, (Flat1dDetector, Flat2dDetector)) and - geometry.ndim == 3): + elif ( + isinstance(geometry, ParallelBeamGeometry) + and isinstance(geometry.detector, (Flat1dDetector, Flat2dDetector)) + and geometry.ndim == 3 + ): # Swap detector axes (see `parallel_3d_geom_to_astra_vecs`) det_row_count = geometry.det_partition.shape[0] det_col_count = geometry.det_partition.shape[1] vecs = parallel_3d_geom_to_astra_vecs(geometry, coords='ASTRA') proj_geom = astra.create_proj_geom( - 'parallel3d_vec', det_row_count, det_col_count, vecs) + 'parallel3d_vec', det_row_count, det_col_count, vecs + ) # Parallel 3D vec elif isinstance(geometry, ParallelVecGeometry) and geometry.ndim == 3: det_row_count = geometry.det_partition.shape[1] det_col_count = geometry.det_partition.shape[0] vecs = vecs_odl_to_astra_coords(geometry.vectors) - proj_geom = astra.create_proj_geom('parallel3d_vec', det_row_count, - det_col_count, vecs) + proj_geom = astra.create_proj_geom( + 'parallel3d_vec', det_row_count, det_col_count, vecs + ) # Cone 3D - elif (isinstance(geometry, DivergentBeamGeometry) and - isinstance(geometry.detector, (Flat1dDetector, Flat2dDetector)) and - geometry.ndim == 3): + elif ( + isinstance(geometry, DivergentBeamGeometry) + and isinstance(geometry.detector, (Flat1dDetector, Flat2dDetector)) + and geometry.ndim == 3 + ): # Swap detector axes (see `conebeam_3d_geom_to_astra_vecs`) det_row_count = geometry.det_partition.shape[0] det_col_count = geometry.det_partition.shape[1] vecs = cone_3d_geom_to_astra_vecs(geometry, coords='ASTRA') proj_geom = astra.create_proj_geom( - 'cone_vec', det_row_count, det_col_count, vecs) + 'cone_vec', det_row_count, det_col_count, vecs + ) # Cone 3D vec elif isinstance(geometry, ConeVecGeometry) and geometry.ndim == 3: det_row_count = geometry.det_partition.shape[1] det_col_count = geometry.det_partition.shape[0] vecs = vecs_odl_to_astra_coords(geometry.vectors) - proj_geom = astra.create_proj_geom('cone_vec', det_row_count, - det_col_count, vecs) + proj_geom = astra.create_proj_geom( + 'cone_vec', det_row_count, det_col_count, vecs + ) else: - raise NotImplementedError('unknown ASTRA geometry type {!r}' - ''.format(geometry)) + raise NotImplementedError( + 'unknown ASTRA geometry type {!r}'.format(geometry) + ) if 'astra' not in geometry.implementation_cache: # Save computed value for later diff --git a/odl/tomo/geometry/conebeam.py b/odl/tomo/geometry/conebeam.py index ed991e9230e..9faee63fd57 100644 --- a/odl/tomo/geometry/conebeam.py +++ b/odl/tomo/geometry/conebeam.py @@ -32,7 +32,6 @@ class FanBeamGeometry(DivergentBeamGeometry): - """Fan beam (2d cone beam) geometry. The source moves on a circle with radius ``src_radius``, and the @@ -41,7 +40,7 @@ class FanBeamGeometry(DivergentBeamGeometry): radii can be chosen as 0, which corresponds to a stationary source or detector, respectively. - The motion parameter is the 1d rotation angle parameterizing source + The motion parameter is the 1d rotation angle that parametrizes source and detector positions simultaneously. In the standard configuration, the detector is perpendicular to the @@ -713,7 +712,6 @@ def __getitem__(self, indices): class ConeBeamGeometry(DivergentBeamGeometry, AxisOrientedGeometry): - """Cone beam geometry with circular/helical source curve. The source moves along a spiral oriented along a fixed ``axis``, with @@ -1212,7 +1210,7 @@ def det_axes(self, angle): Parameters ---------- - angles : float or `array-like` + angle : float or `array-like` Angle(s) in radians describing the counter-clockwise rotation of the detector around `axis`. @@ -1707,7 +1705,7 @@ def cone_beam_geometry(space, src_radius, det_radius, num_angles=None, # used here is (w/2)/(rs+rd) = rho/rs since both are equal to tan(alpha), # where alpha is the half fan angle. rs = float(src_radius) - if (rs <= rho): + if rs <= rho: raise ValueError('source too close to the object, resulting in ' 'infinite detector for full coverage') rd = float(det_radius) @@ -1749,6 +1747,10 @@ def cone_beam_geometry(space, src_radius, det_radius, num_angles=None, det_max_pt = [w / 2, h / 2] if det_shape is None: det_shape = [num_px_horiz, num_px_vert] + else: + raise ValueError( + '`space.ndim` must be 2 or 3, got {}'.format(space.ndim) + ) fan_angle = 2 * np.arctan(rho / rs) if short_scan: @@ -1877,7 +1879,7 @@ def helical_geometry(space, src_radius, det_radius, num_turns, # used here is (w/2)/(rs+rd) = rho/rs since both are equal to tan(alpha), # where alpha is the half fan angle. rs = float(src_radius) - if (rs <= rho): + if rs <= rho: raise ValueError('source too close to the object, resulting in ' 'infinite detector for full coverage') rd = float(det_radius) @@ -1925,7 +1927,6 @@ def helical_geometry(space, src_radius, det_radius, num_turns, class ConeVecGeometry(VecGeometry): - """Cone beam 2D or 3D geometry defined by a collection of vectors. This geometry gives maximal flexibility for representing locations @@ -1963,6 +1964,8 @@ class ConeVecGeometry(VecGeometry): linear paths. """ + # `rotation_matrix` not implemented; reason: missing + @property def _slice_src(self): """Slice for the source position part of `vectors`.""" @@ -2055,9 +2058,9 @@ def det_to_src(self, mparam, dparam, normalized=True): Parameters ---------- - mpar : `motion_params` element + mparam : `motion_params` element Motion parameter at which to evaluate. - dpar : `det_params` element + dparam : `det_params` element Detector parameter at which to evaluate. normalized : bool, optional If ``True``, return a normalized (unit) vector. @@ -2075,8 +2078,10 @@ def det_to_src(self, mparam, dparam, normalized=True): raise ValueError('`dparam` {} not in the valid range {}' ''.format(dparam, self.det_params)) - vec = (self.src_position(mparam) - - self.det_point_position(mparam, dparam)) + vec = ( + self.src_position(mparam) + - self.det_point_position(mparam, dparam) + ) if normalized: # axis = -1 allows this to be vectorized diff --git a/odl/tomo/geometry/geometry.py b/odl/tomo/geometry/geometry.py index e87cb0ec9f5..b3f5ffbc96c 100644 --- a/odl/tomo/geometry/geometry.py +++ b/odl/tomo/geometry/geometry.py @@ -1,4 +1,4 @@ -# Copyright 2014-2019 The ODL contributors +# Copyright 2014-2020 The ODL contributors # # This file is part of ODL. # @@ -30,7 +30,6 @@ class Geometry(object): - """Abstract geometry class. A geometry is described by @@ -419,7 +418,6 @@ def implementation_cache(self): class DivergentBeamGeometry(Geometry): - """Abstract divergent beam geometry class. A geometry characterized by the presence of a point-like ray source. @@ -462,6 +460,8 @@ def det_to_src(self, angle, dparam, normalized=True): Detector parameter(s) at which to evaluate. If ``det_params.ndim >= 2``, a sequence of that length must be provided. + normalized : bool, optional + If ``True``, return a normalized (unit) vector. Returns ------- @@ -552,8 +552,10 @@ def det_to_src(self, angle, dparam, normalized=True): dparam = tuple(np.array(p, dtype=float, copy=False, ndmin=1) for p in dparam) - det_to_src = (self.src_position(angle) - - self.det_point_position(angle, dparam)) + det_to_src = ( + self.src_position(angle) + - self.det_point_position(angle, dparam) + ) if normalized: det_to_src /= np.linalg.norm(det_to_src, axis=-1, keepdims=True) @@ -565,8 +567,11 @@ def det_to_src(self, angle, dparam, normalized=True): class AxisOrientedGeometry(object): + """Mixin class for 3d geometries oriented along an axis. - """Mixin class for 3d geometries oriented along an axis.""" + Makes use of `Geometry` attributes and should thus only be used for + `Geometry` subclasses. + """ def __init__(self, axis): """Initialize a new instance. @@ -616,8 +621,10 @@ def rotation_matrix(self, angle): """ squeeze_out = (np.shape(angle) == ()) angle = np.array(angle, dtype=float, copy=False, ndmin=1) - if (self.check_bounds and - not is_inside_bounds(angle, self.motion_params)): + if ( + self.check_bounds + and not is_inside_bounds(angle, self.motion_params) + ): raise ValueError('`angle` {} not in the valid range {}' ''.format(angle, self.motion_params)) @@ -629,7 +636,6 @@ def rotation_matrix(self, angle): class VecGeometry(Geometry): - """Abstract 2D or 3D geometry defined by a collection of vectors. This geometry gives maximal flexibility for representing locations @@ -671,6 +677,9 @@ class VecGeometry(Geometry): linear paths. """ + # `rotation_matrix` not implemented; reason: missing + # `det_to_src` not implemented; reason: depends on subclass + def __init__(self, det_shape, vectors): """Initialize a new instance. @@ -733,6 +742,8 @@ def __init__(self, det_shape, vectors): max_pt=(det_cell_sides * det_shape) / 2, shape=det_shape) detector = Flat2dDetector(det_part, axes=[det_u, det_v]) + else: + raise RuntimeError('invalid `ndim`') Geometry.__init__(self, ndim, mpart, detector) @@ -800,8 +811,9 @@ def det_refpoint(self, index): array([ 0., 1.]) >>> geom_2d.det_refpoint(1) array([-1., 0.]) - >>> geom_2d.det_refpoint(0.5) # mean value - array([-0.5, 0.5]) + >>> geom_2d.det_refpoint([0.4, 0.6]) # values at closest indices + array([[ 0., 1.], + [-1., 0.]]) In 3D, columns 3 to 5 (starting at 0) determine the detector reference point, here ``(0, 1, 0)`` and ``(-1, 0, 0)``: @@ -816,35 +828,33 @@ def det_refpoint(self, index): array([ 0., 1., 0.]) >>> geom_3d.det_refpoint(1) array([-1., 0., 0.]) - >>> geom_3d.det_refpoint(0.5) # mean value - array([-0.5, 0.5, 0. ]) + >>> geom_3d.det_refpoint([0.4, 0.6]) # values at closest indices + array([[ 0., 1., 0.], + [-1., 0., 0.]]) """ - if (self.check_bounds and - not is_inside_bounds(index, self.motion_params)): - raise ValueError('`index` {} not in the valid range {}' - ''.format(index, self.motion_params)) - - index = np.array(index, dtype=float, copy=False, ndmin=1) - int_part = index.astype(int) - frac_part = index - int_part - - vecs = np.empty((len(index), self.ndim)) - at_right_bdry = (int_part == self.motion_params.max_pt) - vecs[at_right_bdry, :] = self.vectors[int_part[at_right_bdry], - self._slice_det_center] - - not_at_right_bdry = ~at_right_bdry - if np.any(not_at_right_bdry): - pt_left = self.vectors[int_part[not_at_right_bdry], - self._slice_det_center] - pt_right = self.vectors[int_part[not_at_right_bdry] + 1, - self._slice_det_center] - vecs[not_at_right_bdry, :] = ( - pt_left + - frac_part[not_at_right_bdry, None] * (pt_right - pt_left) + if ( + self.check_bounds + and not is_inside_bounds(index, self.motion_params) + ): + raise ValueError( + '`index` {} not in the valid range {}' + ''.format(index, self.motion_params) ) - return vecs.squeeze() + index_int = np.round(index).astype(int) + if index_int.shape == (): + index_int = index_int[None] + squeeze_index = True + else: + squeeze_index = False + + det_refpts = self.vectors[index_int, self._slice_det_center] + if squeeze_index: + return det_refpts[0] + else: + return det_refpts + + # Overrides implementation in `Geometry` def det_point_position(self, index, dparam): """Return the detector point at ``(index, dparam)``. @@ -886,9 +896,16 @@ def det_point_position(self, index, dparam): >>> # d(1) + 2 * u(1) = (-1, 0) + 2 * (0, 1) >>> geom_2d.det_point_position(1, 2) array([-1., 2.]) - >>> # d(0.5) + 2 * u(0.5) = (-0.5, 0.5) + 2 * (0.5, 0.5) - >>> geom_2d.det_point_position(0.5, 2) - array([ 0.5, 1.5]) + >>> # d(0.4) + 2 * u(0.4) = d(0) + 2 * u(0) + >>> geom_2d.det_point_position(0.4, 2) + array([ 2., 1.]) + + Broadcasting of arguments: + + >>> idcs = np.array([0.4, 0.6])[:, None] + >>> dpar = np.array([2.0])[None, :] + >>> geom_2d.det_point_position(idcs, dpar) + Do the same in 3D, with reference points ``(0, 1, 0)`` and ``(-1, 0, 0)``, and horizontal ``u`` axis vectors ``(1, 0, 0)`` and @@ -911,24 +928,39 @@ def det_point_position(self, index, dparam): >>> # (-1, 0, 0) + 2 * (0, 1, 0) + 1 * (0, 0, 1) >>> geom_3d.det_point_position(1, [2, 1]) array([-1., 2., 1.]) - >>> # d(0.5) + 2 * u(0.5) = (-0.5, 0.5) + 2 * (0.5, 0.5) - >>> geom_3d.det_point_position(0.5, [2, 1]) - array([ 0.5, 1.5, 1. ]) - """ - # TODO: vectorize! + >>> # d(0.4) + 2 * u(0.4) + 1 * u(0.4) = d(0) + 2 * u(0) + 1 * v(0) + >>> geom_3d.det_point_position(0.4, [2, 1]) + array([ 2., 1., 1.]) + + Broadcasting of arguments: - if index not in self.motion_params: - raise ValueError('`index` must be contained in `motion_params` ' - '{}, got {}'.format(self.motion_params, index)) - if dparam not in self.det_params: - raise ValueError('`dparam` must be contained in `det_params` ' - '{}, got {}'.format(self.det_params, dparam)) + >>> idcs = np.array([0.4, 0.6])[:, None] + >>> dpar = np.array([2.0, 1.0])[None, :] + >>> geom_3d.det_point_position(idcs, dpar) + """ + if self.check_bounds: + if not is_inside_bounds(index, self.motion_params): + raise ValueError( + '`index` {} not in the valid range {}' + ''.format(index, self.motion_params) + ) + + if not is_inside_bounds(dparam, self.det_params): + raise ValueError( + '`dparam` {} not in the valid range {}' + ''.format(dparam, self.det_params) + ) + + # TODO: broadcast correctly if self.ndim == 2: det_shift = dparam * self.det_axis(index) elif self.ndim == 3: - det_shift = sum(di * ax - for di, ax in zip(dparam, self.det_axes(index))) + det_shift = sum( + di * ax for di, ax in zip(dparam, self.det_axes(index)) + ) + else: + raise RuntimeError('invalid `ndim`') return self.det_refpoint(index) + det_shift @@ -962,31 +994,38 @@ def det_axis(self, index): array([ 1., 0.]) >>> geom_2d.det_axis(1) array([ 0., 1.]) - >>> geom_2d.det_axis(0.5) # mean value - array([ 0.5, 0.5]) + >>> geom_2d.det_axis([0.4, 0.6]) # values at closest indices + array([[ 1., 0.], + [ 0., 1.]]) """ - # TODO: vectorize! - - if index not in self.motion_params: - raise ValueError('`index` must be contained in `motion_params` ' - '{}, got {}'.format(self.motion_params, index)) + if ( + self.check_bounds + and not is_inside_bounds(index, self.motion_params) + ): + raise ValueError( + '`index` {} not in the valid range {}' + ''.format(index, self.motion_params) + ) if self.ndim != 2: raise NotImplementedError( '`det_axis` only valid for 2D geometries, use `det_axes` ' - 'in 3D') + 'in 3D' + ) - index = float(index) - int_part = int(index) - frac_part = index - int_part - if int_part == self.motion_params.max_pt: - det_u = self.vectors[int_part, self._slice_det_u] + index_int = np.round(index).astype(int) + if index_int.shape == (): + index_int = index_int[None] + squeeze_index = True else: - det_u_left = self.vectors[int_part, self._slice_det_u] - det_u_right = self.vectors[int_part + 1, self._slice_det_u] - det_u = det_u_left + frac_part * (det_u_right - det_u_left) + squeeze_index = False - return det_u + vectors = self.vectors[index_int] + det_us = vectors[:, self._slice_det_u] + if squeeze_index: + return det_us[0] + else: + return det_us def det_axes(self, index): """Return the detector axes at ``index`` (for 2D or 3D geometries). @@ -994,13 +1033,14 @@ def det_axes(self, index): Parameters ---------- index : `motion_params` element - Index of the projection. Non-integer indices result in - interpolated vectors. + Index of the projection. Non-integer indices are rounded to + closest integer. Returns ------- axes : tuple of `numpy.ndarray`, shape ``(ndim,)`` - The detector axes at ``index``. + The detector axes at ``index``, 1 for ``ndim == 2`` and + 2 for ``ndim == 3``. Examples -------- @@ -1019,37 +1059,45 @@ def det_axes(self, index): (array([ 1., 0., 0.]), array([ 0., 0., 1.])) >>> geom_3d.det_axes(1) (array([ 0., 1., 0.]), array([ 0., 0., 1.])) - >>> geom_3d.det_axes(0.5) # mean value - (array([ 0.5, 0.5, 0. ]), array([ 0., 0., 1.])) + >>> axs = geom_3d.det_axes([0.4, 0.6]) # values at closest indices + >>> axs[0] # first axis + array([[ 1., 0., 0.], + [ 0., 1., 0.]]) + >>> axs[1] # second axis + array([[ 0., 0., 1.], + [ 0., 0., 1.]]) """ - # TODO: vectorize! - - if index not in self.motion_params: - raise ValueError('`index` must be contained in `motion_params` ' - '{}, got {}'.format(self.motion_params, index)) - - index = float(index) - int_part = int(index) - frac_part = index - int_part - if int_part == self.motion_params.max_pt: - det_u = self.vectors[int_part, self._slice_det_u] - if self.ndim == 2: - return (det_u,) - elif self.ndim == 3: - det_v = self.vectors[int_part, self._slice_det_v] - return (det_u, det_v) + if ( + self.check_bounds + and not is_inside_bounds(index, self.motion_params) + ): + raise ValueError( + '`index` {} not in the valid range {}' + ''.format(index, self.motion_params) + ) + + index_int = np.round(index).astype(int) + if index_int.shape == (): + index_int = index_int[None] + squeeze_index = True else: - det_u_left = self.vectors[int_part, self._slice_det_u] - det_u_right = self.vectors[int_part + 1, self._slice_det_u] - det_u = det_u_left + frac_part * (det_u_right - det_u_left) - - if self.ndim == 2: - return (det_u,) - elif self.ndim == 3: - det_v_left = self.vectors[int_part, self._slice_det_v] - det_v_right = self.vectors[int_part + 1, self._slice_det_v] - det_v = det_v_left + frac_part * (det_v_right - det_v_left) - return (det_u, det_v) + squeeze_index = False + + vectors = self.vectors[index_int] + if self.ndim == 2: + det_us = vectors[:, self._slice_det_u] + retval_lst = [det_us[0]] if squeeze_index else [det_us] + elif self.ndim == 3: + det_us = vectors[:, self._slice_det_u] + det_vs = vectors[:, self._slice_det_v] + if squeeze_index: + retval_lst = [det_us[0], det_vs[0]] + else: + retval_lst = [det_us, det_vs] + else: + raise RuntimeError('invalid `ndim`') + + return tuple(retval_lst) def __getitem__(self, indices): """Return ``self[indices]``. @@ -1124,8 +1172,9 @@ def __repr__(self): posargs = [self.det_partition.shape, self.vectors] with npy_printoptions(precision=REPR_PRECISION): inner_parts = signature_string_parts(posargs, []) - return repr_string(self.__class__.__name__, inner_parts, - allow_mixed_seps=False) + return repr_string( + self.__class__.__name__, inner_parts, allow_mixed_seps=False + ) if __name__ == '__main__': diff --git a/odl/tomo/geometry/parallel.py b/odl/tomo/geometry/parallel.py index e8c1c2dd66e..11e20123129 100644 --- a/odl/tomo/geometry/parallel.py +++ b/odl/tomo/geometry/parallel.py @@ -30,7 +30,6 @@ class ParallelBeamGeometry(Geometry): - """Abstract parallel beam geometry in 2 or 3 dimensions. Parallel geometries are characterized by a virtual source at @@ -331,7 +330,6 @@ def det_to_src(self, angle, dparam): class Parallel2dGeometry(ParallelBeamGeometry): - """Parallel beam geometry in 2d. The motion parameter is the counter-clockwise rotation angle around @@ -635,8 +633,10 @@ def rotation_matrix(self, angle): """ squeeze_out = (np.shape(angle) == ()) angle = np.array(angle, dtype=float, copy=False, ndmin=1) - if (self.check_bounds and - not is_inside_bounds(angle, self.motion_params)): + if ( + self.check_bounds + and not is_inside_bounds(angle, self.motion_params) + ): raise ValueError('`angle` {} not in the valid range {}' ''.format(angle, self.motion_params)) @@ -669,7 +669,7 @@ def __repr__(self): return '{}(\n{}\n)'.format(self.__class__.__name__, indent(sig_str)) def __getitem__(self, indices): - """Return self[slc] + """Return self[slc]. This is defined by:: @@ -705,7 +705,6 @@ def __getitem__(self, indices): class Parallel3dEulerGeometry(ParallelBeamGeometry): - """Parallel beam geometry in 3d. The motion parameters are two or three Euler angles, and the detector @@ -1044,8 +1043,10 @@ def rotation_matrix(self, angles): angles_in = angles angles = tuple(np.array(angle, dtype=float, copy=False, ndmin=1) for angle in angles) - if (self.check_bounds and - not is_inside_bounds(angles, self.motion_params)): + if ( + self.check_bounds + and not is_inside_bounds(angles, self.motion_params) + ): raise ValueError('`angles` {} not in the valid range ' '{}'.format(angles_in, self.motion_params)) @@ -1079,7 +1080,6 @@ def __repr__(self): class Parallel3dAxisGeometry(ParallelBeamGeometry, AxisOrientedGeometry): - """Parallel beam geometry in 3d with single rotation axis. The motion parameter is the rotation angle around the specified @@ -1474,7 +1474,6 @@ def __getitem__(self, indices): class ParallelVecGeometry(VecGeometry): - """Parallel beam 2D or 3D geometry defined by a collection of vectors. This geometry gives maximal flexibility for representing locations @@ -1512,6 +1511,8 @@ class ParallelVecGeometry(VecGeometry): linear paths. """ + # `rotation_matrix` not implemented; reason: missing + @property def _slice_ray(self): """Slice for the ray direction part of `vectors`.""" @@ -1614,6 +1615,25 @@ def det_to_src(self, index, dparam): (4, 5, 2) """ squeeze_index = (np.shape(index) == ()) + + if self.check_bounds: + if not is_inside_bounds(index, self.motion_params): + raise ValueError( + '`index` {} not in the valid range {}' + ''.format(index, self.motion_params) + ) + + if not is_inside_bounds(dparam, self.det_params): + raise ValueError( + '`dparam` {} not in the valid range {}' + ''.format(dparam, self.det_params) + ) + + index_int = np.round(index).astype(int) + vectors = self.vectors[index_int] + ray_dirs = vectors[self._slice_ray] + det_centers = vectors[self._slice_det_center] + index_in = index index = np.array(index, dtype=float, copy=False, ndmin=1) @@ -1626,14 +1646,6 @@ def det_to_src(self, index, dparam): dparam = tuple(np.array(p, dtype=float, copy=False, ndmin=1) for p in dparam) - if self.check_bounds: - if not is_inside_bounds(index, self.motion_params): - raise ValueError('`index` {} not in the valid range {}' - ''.format(index_in, self.motion_params)) - - if not is_inside_bounds(dparam, self.det_params): - raise ValueError('`dparam` {} not in the valid range {}' - ''.format(dparam_in, self.det_params)) at_max_flat = (index == self.motion_params.max_pt).ravel() @@ -1669,8 +1681,8 @@ def det_to_src(self, index, dparam): print('ray_right shape:', ray_right.shape) print('index_frac_part shape:', index_frac_part.shape) ray_dir[~at_max_flat, ...] = ( - ray_left + - index_frac_part * (ray_right - ray_left)) + ray_left + index_frac_part * (ray_right - ray_left) + ) ray_dir *= -1 / np.linalg.norm(ray_dir, axis=-1, keepdims=True) From 252b7d0d8b61501407043b008e74a55549bb5658 Mon Sep 17 00:00:00 2001 From: Holger Kohr Date: Tue, 21 Jan 2020 00:06:29 +0100 Subject: [PATCH 18/24] WIP: Fix broadcasting in vec geometry methods --- odl/tomo/geometry/geometry.py | 60 +++++++++++++++++++---------------- 1 file changed, 33 insertions(+), 27 deletions(-) diff --git a/odl/tomo/geometry/geometry.py b/odl/tomo/geometry/geometry.py index b3f5ffbc96c..a222ba8c26e 100644 --- a/odl/tomo/geometry/geometry.py +++ b/odl/tomo/geometry/geometry.py @@ -889,16 +889,16 @@ def det_point_position(self, index, dparam): >>> geom_2d = odl.tomo.ParallelVecGeometry(det_shape_2d, vecs_2d) >>> # This is equal to d(0) = (0, 1) >>> geom_2d.det_point_position(0, 0) - array([ 0., 1.]) + array([ 0., 1.]) >>> # d(0) + 2 * u(0) = (0, 1) + 2 * (1, 0) >>> geom_2d.det_point_position(0, 2) - array([ 2., 1.]) + array([ 2., 1.]) >>> # d(1) + 2 * u(1) = (-1, 0) + 2 * (0, 1) >>> geom_2d.det_point_position(1, 2) - array([-1., 2.]) + array([-1., 2.]) >>> # d(0.4) + 2 * u(0.4) = d(0) + 2 * u(0) >>> geom_2d.det_point_position(0.4, 2) - array([ 2., 1.]) + array([ 2., 1.]) Broadcasting of arguments: @@ -935,7 +935,7 @@ def det_point_position(self, index, dparam): Broadcasting of arguments: >>> idcs = np.array([0.4, 0.6])[:, None] - >>> dpar = np.array([2.0, 1.0])[None, :] + >>> dpar = np.array([2.0, 1.0])[:, None] >>> geom_3d.det_point_position(idcs, dpar) """ @@ -952,10 +952,10 @@ def det_point_position(self, index, dparam): ''.format(dparam, self.det_params) ) - # TODO: broadcast correctly if self.ndim == 2: det_shift = dparam * self.det_axis(index) elif self.ndim == 3: + axes = self.det_axes(index) det_shift = sum( di * ax for di, ax in zip(dparam, self.det_axes(index)) ) @@ -1038,9 +1038,15 @@ def det_axes(self, index): Returns ------- - axes : tuple of `numpy.ndarray`, shape ``(ndim,)`` - The detector axes at ``index``, 1 for ``ndim == 2`` and - 2 for ``ndim == 3``. + axes : `numpy.ndarray` + Unit vectors along which the detector is aligned, an array + with following shape: + + - In 2D: If ``index`` is a single parameter, the shape is + ``(2,)``, otherwise ``index.shape + (2,)``. + + - In 3D: If ``index`` is a single parameter, the shape is + ``(2, 3)``, otherwise ``index.shape + (2, 3)``. Examples -------- @@ -1056,16 +1062,17 @@ def det_axes(self, index): >>> det_shape_3d = (10, 20) >>> geom_3d = odl.tomo.ParallelVecGeometry(det_shape_3d, vecs_3d) >>> geom_3d.det_axes(0) - (array([ 1., 0., 0.]), array([ 0., 0., 1.])) + array([[ 1., 0., 0.], + [ 0., 0., 1.]]) >>> geom_3d.det_axes(1) - (array([ 0., 1., 0.]), array([ 0., 0., 1.])) - >>> axs = geom_3d.det_axes([0.4, 0.6]) # values at closest indices - >>> axs[0] # first axis - array([[ 1., 0., 0.], - [ 0., 1., 0.]]) - >>> axs[1] # second axis - array([[ 0., 0., 1.], - [ 0., 0., 1.]]) + array([[ 0., 1., 0.], + [ 0., 0., 1.]]) + >>> geom_3d.det_axes([0.4, 0.6]) # values at closest indices + array([[[ 1., 0., 0.], + [ 0., 0., 1.]], + + [[ 0., 1., 0.], + [ 0., 0., 1.]]]) """ if ( self.check_bounds @@ -1085,19 +1092,18 @@ def det_axes(self, index): vectors = self.vectors[index_int] if self.ndim == 2: - det_us = vectors[:, self._slice_det_u] - retval_lst = [det_us[0]] if squeeze_index else [det_us] + axes = np.empty(index_int.shape + (2,)) + axes[:] = vectors[:, self._slice_det_u] elif self.ndim == 3: - det_us = vectors[:, self._slice_det_u] - det_vs = vectors[:, self._slice_det_v] - if squeeze_index: - retval_lst = [det_us[0], det_vs[0]] - else: - retval_lst = [det_us, det_vs] + axes = np.empty(index_int.shape + (2, 3)) + axes[..., 0, :] = vectors[..., self._slice_det_u] + axes[..., 1, :] = vectors[..., self._slice_det_v] else: raise RuntimeError('invalid `ndim`') - return tuple(retval_lst) + if squeeze_index: + axes = axes[0] + return axes def __getitem__(self, indices): """Return ``self[indices]``. From 047a79dc1074ff4e664e3f47ab446a997c60a683 Mon Sep 17 00:00:00 2001 From: Holger Kohr Date: Sun, 6 Sep 2020 13:16:12 +0200 Subject: [PATCH 19/24] Check dimension in parallel_beam_geometry --- odl/tomo/geometry/conebeam.py | 1 + odl/tomo/geometry/parallel.py | 5 +++++ 2 files changed, 6 insertions(+) diff --git a/odl/tomo/geometry/conebeam.py b/odl/tomo/geometry/conebeam.py index 9faee63fd57..25eca19a2af 100644 --- a/odl/tomo/geometry/conebeam.py +++ b/odl/tomo/geometry/conebeam.py @@ -2070,6 +2070,7 @@ def det_to_src(self, mparam, dparam, normalized=True): vec : `numpy.ndarray`, shape (`ndim`,) (Unit) vector pointing from the detector to the source. """ + # TODO: vectorize if self.check_bounds: if not is_inside_bounds(mparam, self.motion_params): raise ValueError('`mparam` {} not in the valid range {}' diff --git a/odl/tomo/geometry/parallel.py b/odl/tomo/geometry/parallel.py index 11e20123129..5b632220679 100644 --- a/odl/tomo/geometry/parallel.py +++ b/odl/tomo/geometry/parallel.py @@ -1796,6 +1796,11 @@ def parallel_beam_geometry(space, num_angles=None, det_shape=None): det_max_pt = [rho, max_h] if det_shape is None: det_shape = [num_px_horiz, num_px_vert] + else: + raise ValueError( + "`space` must be 2- or 3-dimensional, got space.ndim={}" + "".format(space.ndim) + ) if num_angles is None: num_angles = int(np.ceil(omega * rho)) From 1445d133909df8731ca37d5cd3a2a26710937e41 Mon Sep 17 00:00:00 2001 From: Holger Kohr Date: Sun, 6 Sep 2020 15:54:24 +0200 Subject: [PATCH 20/24] Minor fixes in docs and comments of astra_setup.py --- odl/tomo/backends/astra_setup.py | 43 ++++++++++++++++---------------- 1 file changed, 21 insertions(+), 22 deletions(-) diff --git a/odl/tomo/backends/astra_setup.py b/odl/tomo/backends/astra_setup.py index 9078dd8b385..e19894abf86 100644 --- a/odl/tomo/backends/astra_setup.py +++ b/odl/tomo/backends/astra_setup.py @@ -419,8 +419,8 @@ def vecs_odl_to_astra_coords(vecs): | | | ``Z`` | ``X`` | +-------+--------+-------+--------+ - Thus, the conversion ODL→ASTRA is a rotation by +90 degrees in 2D, and - an XZ axis swap in 3D. + Thus, the coordinate conversion ODL→ASTRA is a transposed rotation + by +90 degrees in 2D, and an XZ axis swap in 3D. Parameters ---------- @@ -440,9 +440,9 @@ def vecs_odl_to_astra_coords(vecs): if vecs.shape[1] == 6: # 2D geometry # ODL x = ASTRA -y, ODL y = ASTRA x - # Thus: ODL->ASTRA conversion is a rotation by +90 degrees + # Thus: ODL->ASTRA conversion is a transposed rotation by +90 degrees rot_90 = euler_matrix(np.pi / 2) - # Bulk matrix-vector product + # Bulk transposed-matrix-vector product vecs[:, 0:2] = np.einsum('kj,ik->ij', rot_90, vecs[:, 0:2]) vecs[:, 2:4] = np.einsum('kj,ik->ij', rot_90, vecs[:, 2:4]) vecs[:, 4:6] = np.einsum('kj,ik->ij', rot_90, vecs[:, 4:6]) @@ -451,7 +451,7 @@ def vecs_odl_to_astra_coords(vecs): elif vecs.shape[1] == 12: # 3D geometry # 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. + # so we need to adapt to this by swapping x and z newind = [] for i in range(4): newind.extend([2 + 3 * i, 1 + 3 * i, 0 + 3 * i]) @@ -479,8 +479,8 @@ def vecs_astra_to_odl_coords(vecs): | | | ``Z`` | ``X`` | +-------+--------+-------+--------+ - Thus, the conversion ASTRA→ODL is a rotation by -90 degrees in 2D, and - an XZ axis swap in 3D. + Thus, the coordinate conversion ASTRA->ODL is a transposed rotation + by -90 degrees in 2D, and an XZ axis swap in 3D. Parameters ---------- @@ -500,9 +500,9 @@ def vecs_astra_to_odl_coords(vecs): if vecs.shape[1] == 6: # 2D geometry # ODL x = ASTRA -y, ODL y = ASTRA x - # Thus: ODL->ASTRA conversion is a rotation by -90 degrees + # Thus: ODL->ASTRA conversion is a transposed rotation by -90 degrees rot_minus_90 = euler_matrix(-np.pi / 2) - # Bulk matrix-vector product + # Bulk transposed-matrix-vector product vecs[:, 0:2] = np.einsum('kj,ik->ij', rot_minus_90, vecs[:, 0:2]) vecs[:, 2:4] = np.einsum('kj,ik->ij', rot_minus_90, vecs[:, 2:4]) vecs[:, 4:6] = np.einsum('kj,ik->ij', rot_minus_90, vecs[:, 4:6]) @@ -511,11 +511,10 @@ def vecs_astra_to_odl_coords(vecs): elif vecs.shape[1] == 12: # 3D geometry # 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. + # so we need to adapt to this by swapping x and z newind = [] for i in range(4): newind.extend([2 + 3 * i, 1 + 3 * i, 0 + 3 * i]) - newind = np.argsort(newind).tolist() return vecs[:, newind] else: raise ValueError('`vecs` must have shape (N, 6) or (N, 12), got ' @@ -529,16 +528,16 @@ def parallel_2d_geom_to_astra_vecs(geometry, coords='ODL'): parallel beam geometries, see ``'parallel_vec'`` in the `ASTRA projection geometry documentation`_. - Each row of the returned vectors corresponds to a single projection - and consists of :: + Each row of the returned array is a vector corresponding to a single + projection and consists of :: (rayX, rayY, dX, dY, uX, uY) with - - ``ray``: the ray direction - - ``d`` : the center of the detector - - ``u`` : the vector from detector pixel 0 to 1 + - ``ray``: the ray direction + - ``d`` : the center of the detector + - ``u`` : the vector from detector pixel 0 to 1 .. note:: The returned vectors are the result of converting from ODL coordinates @@ -599,8 +598,8 @@ def parallel_3d_geom_to_astra_vecs(geometry, coords='ODL'): parallel beam geometries, see ``'parallel3d_vec'`` in the `ASTRA projection geometry documentation`_. - Each row of the returned vectors corresponds to a single projection - and consists of :: + Each row of the returned array is a vector corresponding to a single + projection and consists of :: (rayX, rayY, rayZ, dX, dY, dZ, uX, uY, uZ, vX, vY, vZ) @@ -676,8 +675,8 @@ def cone_2d_geom_to_astra_vecs(geometry, coords='ODL'): fan beam geometries, see ``'fanflat_vec'`` in the `ASTRA projection geometry documentation`_. - Each row of the returned vectors corresponds to a single projection - and consists of :: + Each row of the returned array is a vector corresponding to a single + projection and consists of :: (srcX, srcY, dX, dY, uX, uY) @@ -742,8 +741,8 @@ def cone_3d_geom_to_astra_vecs(geometry, coords='ODL'): cone beam geometries, see ``'cone_vec'`` in the `ASTRA projection geometry documentation`_. - Each row of the returned vectors corresponds to a single projection - and consists of :: + Each row of the returned array is a vector corresponding to a single + projection and consists of :: (srcX, srcY, srcZ, dX, dY, dZ, uX, uY, uZ, vX, vY, vZ) From 5562656bb60c59743f1d819bcade795b63aff07b Mon Sep 17 00:00:00 2001 From: Holger Kohr Date: Sun, 8 Nov 2020 01:16:35 +0100 Subject: [PATCH 21/24] Replace `is` by `==` in pspace ops comparison --- odl/operator/pspace_ops.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/odl/operator/pspace_ops.py b/odl/operator/pspace_ops.py index e1ddcd7802a..219cea714c1 100644 --- a/odl/operator/pspace_ops.py +++ b/odl/operator/pspace_ops.py @@ -263,7 +263,7 @@ def _convert_to_spmatrix(operators): '{}'.format(len(row), i, ncols)) for j, col in enumerate(row): - if col is None or col is 0: + if col is None or col == 0: pass elif isinstance(col, Operator): irow.append(i) From 6ca631df9d0019e03f2a183124f70f701a1a1109 Mon Sep 17 00:00:00 2001 From: Holger Kohr Date: Sun, 8 Nov 2020 01:17:09 +0100 Subject: [PATCH 22/24] Improve formatting in `pkg_supports` --- odl/util/utility.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/odl/util/utility.py b/odl/util/utility.py index ca41e88dd67..14c3f1d9ff2 100644 --- a/odl/util/utility.py +++ b/odl/util/utility.py @@ -1421,8 +1421,10 @@ def pkg_supports(feature, pkg_version, pkg_feat_dict): ver_reqs = [Requirement(ver_spec) for ver_spec in ver_specs] # If one of the requirements in the list is met, return True, else False - return any(req.specifier.contains(pkg_version, prereleases=True) - for req in ver_reqs) + return any( + req.specifier.contains(pkg_version, prereleases=True) + for req in ver_reqs + ) @contextmanager From accaf92bab3a6a3e5b5b3f41ad55dfa73f9065cf Mon Sep 17 00:00:00 2001 From: Holger Kohr Date: Sun, 8 Nov 2020 01:18:17 +0100 Subject: [PATCH 23/24] Temporary weight changes and minor fixes --- odl/tomo/backends/astra_cuda.py | 45 ++++--------- odl/tomo/backends/astra_setup.py | 109 +++++++++++++++++++------------ odl/tomo/geometry/conebeam.py | 3 +- odl/tomo/operators/ray_trafo.py | 6 +- 4 files changed, 88 insertions(+), 75 deletions(-) diff --git a/odl/tomo/backends/astra_cuda.py b/odl/tomo/backends/astra_cuda.py index 52eb8fda150..4388dd7e47a 100644 --- a/odl/tomo/backends/astra_cuda.py +++ b/odl/tomo/backends/astra_cuda.py @@ -351,11 +351,11 @@ def astra_cuda_bp_scaling_factor(proj_space, vol_space, geometry): angle_extent = geometry.motion_partition.extent num_angles = geometry.motion_partition.shape # TODO: this gives the wrong factor for Parallel3dEulerGeometry with - # 2 angles + # 2 angles angle_weighting = (angle_extent / num_angles).prod() # Correct in case of non-weighted spaces - proj_extent = float(proj_space.partition.extent.prod()) + proj_extent = float(np.prod(proj_space.partition.extent)) proj_size = float(proj_space.partition.size) proj_weighting = proj_extent / proj_size @@ -394,9 +394,6 @@ def astra_cuda_bp_scaling_factor(proj_space, vol_space, geometry): src_radius = geometry.src_radius det_radius = geometry.det_radius scaling_factor *= ((src_radius + det_radius) / src_radius) ** 2 - elif isinstance(geometry, VecGeometry): - scaling_factor = (geometry.det_partition.cell_volume / - vol_space.cell_volume) elif isinstance(geometry, ParallelVecGeometry): if geometry.ndim == 2: # Scales with 1 / cell_volume @@ -405,15 +402,14 @@ def astra_cuda_bp_scaling_factor(proj_space, vol_space, geometry): # Scales with cell volume scaling_factor /= vol_space.cell_volume elif isinstance(geometry, ConeVecGeometry): - # TODO: this should be based on the distances per angle, would - # probably be part of the weighting then if geometry.ndim == 2: # Scales with 1 / cell_volume scaling_factor *= float(vol_space.cell_volume) elif geometry.ndim == 3: det_px_area = geometry.det_partition.cell_volume - scaling_factor *= (det_px_area ** 2 / - vol_space.cell_volume ** 2) + scaling_factor *= ( + det_px_area ** 2 / vol_space.cell_volume ** 2 + ) elif parse_version(ASTRA_VERSION) < parse_version('1.9.0dev'): # Scaling for the 1.8.x releases @@ -440,9 +436,6 @@ def astra_cuda_bp_scaling_factor(proj_space, vol_space, geometry): src_radius = geometry.src_radius det_radius = geometry.det_radius scaling_factor *= ((src_radius + det_radius) / src_radius) ** 2 - elif isinstance(geometry, VecGeometry): - scaling_factor = (geometry.det_partition.cell_volume / - vol_space.cell_volume) elif isinstance(geometry, ParallelVecGeometry): if geometry.ndim == 2: # Scales with 1 / cell_volume @@ -451,24 +444,15 @@ def astra_cuda_bp_scaling_factor(proj_space, vol_space, geometry): # Scales with cell volume scaling_factor /= vol_space.cell_volume elif isinstance(geometry, ConeVecGeometry): - # TODO: this should be based on the distances per angle, would - # probably be part of the weighting then if geometry.ndim == 2: # Scales with 1 / cell_volume scaling_factor *= float(vol_space.cell_volume) elif geometry.ndim == 3: det_px_area = geometry.det_partition.cell_volume - scaling_factor *= (det_px_area ** 2 / - vol_space.cell_volume ** 2) + scaling_factor *= ( + det_px_area ** 2 / vol_space.cell_volume ** 2 + ) - # Correction for scaled 1/r^2 factor in ASTRA's density weighting. - # This compensates for scaled voxels and pixels, as well as a - # missing factor src_radius ** 2 in the ASTRA BP with - # density weighting. - det_px_area = geometry.det_partition.cell_volume - scaling_factor *= ( - src_radius ** 2 * det_px_area ** 2 / vol_space.cell_volume ** 2 - ) elif parse_version(ASTRA_VERSION) < parse_version('1.9.9.dev'): # Scaling for intermediate dev releases between 1.8.3 and 1.9.9.dev if isinstance(geometry, Parallel2dGeometry): @@ -501,9 +485,6 @@ def astra_cuda_bp_scaling_factor(proj_space, vol_space, geometry): # density weighting. det_px_area = geometry.det_partition.cell_volume scaling_factor *= (src_radius ** 2 * det_px_area ** 2) - elif isinstance(geometry, VecGeometry): - # TODO: correct in 1.8? With differently weighted spaces as well? - scaling_factor = vol_space.cell_volume elif isinstance(geometry, ParallelVecGeometry): if geometry.ndim == 2: # Scales with 1 / cell_volume @@ -512,15 +493,17 @@ def astra_cuda_bp_scaling_factor(proj_space, vol_space, geometry): # Scales with cell volume scaling_factor /= vol_space.cell_volume elif isinstance(geometry, ConeVecGeometry): - # TODO: this should be based on the distances per angle, would - # probably be part of the weighting then + # TODO: the distance weighting correction should be based on the + # distances per angle; part of the weighting scheme probably if geometry.ndim == 2: # Scales with 1 / cell_volume scaling_factor *= float(vol_space.cell_volume) elif geometry.ndim == 3: det_px_area = geometry.det_partition.cell_volume - scaling_factor *= (det_px_area ** 2 / - vol_space.cell_volume ** 2) + scaling_factor *= ( + det_px_area ** 2 / vol_space.cell_volume ** 2 + ) + else: # Scaling for versions since 1.9.9.dev diff --git a/odl/tomo/backends/astra_setup.py b/odl/tomo/backends/astra_setup.py index e19894abf86..a64c287b70b 100644 --- a/odl/tomo/backends/astra_setup.py +++ b/odl/tomo/backends/astra_setup.py @@ -238,9 +238,11 @@ def astra_volume_geometry(vol_space): # NOTE: We need to flip the sign of the (ODL) x component since # ASTRA seems to move it in the other direction. Not quite clear # why. - vol_geom = astra.create_vol_geom(vol_shp[0], vol_shp[1], - vol_min[1], vol_max[1], - -vol_max[0], -vol_min[0]) + vol_geom = astra.create_vol_geom( + vol_shp[0], vol_shp[1], + vol_min[1], vol_max[1], + -vol_max[0], -vol_min[0], + ) elif vol_space.ndim == 3: # Not supported in all versions of ASTRA if ( @@ -265,10 +267,12 @@ def astra_volume_geometry(vol_space): # 'WindowMaxY': y_min, # 'WindowMinZ': x_min, # 'WindowMaxZ': x_min} - vol_geom = astra.create_vol_geom(vol_shp[1], vol_shp[2], vol_shp[0], - vol_min[2], vol_max[2], - vol_min[1], vol_max[1], - vol_min[0], vol_max[0]) + vol_geom = astra.create_vol_geom( + vol_shp[1], vol_shp[2], vol_shp[0], + vol_min[2], vol_max[2], + vol_min[1], vol_max[1], + vol_min[0], vol_max[0], + ) else: raise ValueError('{}-dimensional volume geometries not supported ' 'by ASTRA'.format(vol_space.ndim)) @@ -316,7 +320,8 @@ def astra_projection_geometry(geometry): # geometry by 90 degrees clockwise angles = geometry.angles - np.pi / 2 proj_geom = astra.create_proj_geom( - 'parallel', det_width, det_count, angles) + 'parallel', det_width, det_count, angles + ) # Parallel 2D vec elif isinstance(geometry, ParallelVecGeometry) and geometry.ndim == 2: @@ -458,8 +463,10 @@ def vecs_odl_to_astra_coords(vecs): return vecs[:, newind] else: - raise ValueError('`vecs` must have shape (N, 6) or (N, 12), got ' - 'array with shape {}'.format(vecs.shape)) + raise ValueError( + '`vecs` must have shape (N, 6) or (N, 12), got array with shape {}' + ''.format(vecs.shape) + ) def vecs_astra_to_odl_coords(vecs): @@ -516,9 +523,12 @@ def vecs_astra_to_odl_coords(vecs): for i in range(4): newind.extend([2 + 3 * i, 1 + 3 * i, 0 + 3 * i]) return vecs[:, newind] + else: - raise ValueError('`vecs` must have shape (N, 6) or (N, 12), got ' - 'array with shape {}'.format(vecs.shape)) + raise ValueError( + '`vecs` must have shape (N, 6) or (N, 12), got array with shape {}' + ''.format(vecs.shape) + ) def parallel_2d_geom_to_astra_vecs(geometry, coords='ODL'): @@ -653,10 +663,10 @@ def parallel_3d_geom_to_astra_vecs(geometry, coords='ODL'): # `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 - # Swap detector axes to have better memory layout in projection data. + # Swap detector axes to have better memory layout in projection data. # ASTRA produces `(v, theta, u)` layout, and to map to ODL layout # `(theta, u, v)` a complete roll must be performed, which is the - # worst case (compeltely discontiguous). + # worst case (completely non-contiguous). # Instead we swap `u` and `v`, resulting in the effective ASTRA result # `(u, theta, v)`. Here we only need to swap axes 0 and 1, which # keeps at least contiguous blocks in `v`. @@ -840,8 +850,10 @@ def astra_data(astra_geom, datatype, data=None, ndim=2, allow_copy=False): if isinstance(data, (DiscretizedSpaceElement, np.ndarray)): ndim = data.ndim else: - raise TypeError('`data` {!r} is neither DiscretizedSpaceElement ' - 'instance nor a `numpy.ndarray`'.format(data)) + raise TypeError( + '`data` must be `DiscretizedSpaceElement` or `numpy.ndarray`, ' + 'got {!r}'.format(data) + ) else: ndim = int(ndim) @@ -860,8 +872,9 @@ def astra_data(astra_geom, datatype, data=None, ndim=2, allow_copy=False): link = astra.data3d.link create = astra.data3d.create else: - raise ValueError('{}-dimensional data not supported' - ''.format(ndim)) + raise ValueError( + '{}-dimensional data not supported'.format(ndim) + ) # ASTRA checks if data is c-contiguous and aligned if data is not None: @@ -875,9 +888,10 @@ def astra_data(astra_geom, datatype, data=None, ndim=2, allow_copy=False): return link(astra_dtype_str, astra_geom, data.asarray()) else: # Something else than NumPy data representation - raise NotImplementedError('ASTRA supports data wrapping only ' - 'for `numpy.ndarray` instances, got ' - '{!r}'.format(data)) + raise NotImplementedError( + 'ASTRA supports data wrapping only for `numpy.ndarray`, ' + 'got {!r}'.format(data) + ) else: return create(astra_dtype_str, astra_geom) @@ -904,8 +918,9 @@ def astra_projector(astra_proj_type, astra_vol_geom, astra_proj_geom, ndim): Handle for the created ASTRA internal projector object. """ if 'type' not in astra_proj_geom: - raise ValueError('invalid projection geometry dict {}' - ''.format(astra_proj_geom)) + raise ValueError( + 'invalid projection geometry dict {}'.format(astra_proj_geom) + ) ndim = int(ndim) @@ -954,11 +969,12 @@ def astra_projector(astra_proj_type, astra_vol_geom, astra_proj_geom, ndim): ) # Create config dict - proj_cfg = {} - proj_cfg['type'] = astra_proj_type - proj_cfg['VolumeGeometry'] = astra_vol_geom - proj_cfg['ProjectionGeometry'] = astra_proj_geom - proj_cfg['options'] = {} + proj_cfg = { + 'type': astra_proj_type, + 'VolumeGeometry': astra_vol_geom, + 'ProjectionGeometry': astra_proj_geom, + 'options': {} + } # Add the approximate 1/r^2 weighting exposed in intermediate versions of # ASTRA @@ -1001,25 +1017,36 @@ def astra_algorithm(direction, ndim, vol_id, sino_id, proj_id, impl): if direction not in ('forward', 'backward'): raise ValueError("`direction` '{}' not understood".format(direction)) if ndim not in (2, 3): - raise ValueError('{}-dimensional projectors not supported' - ''.format(ndim)) + raise ValueError( + '{}-dimensional projectors not supported'.format(ndim) + ) if impl not in ('cpu', 'cuda'): - raise ValueError("`impl` type '{}' not understood" - ''.format(impl)) + raise ValueError( + '`impl` type {!r} not understood'.format(impl) + ) if ndim == 3 and impl == 'cpu': raise NotImplementedError( - '3d algorithms for CPU not supported by ASTRA') + '3d algorithms for CPU not supported by ASTRA' + ) if proj_id is None and impl == 'cpu': raise ValueError("'cpu' implementation requires projector ID") - algo_map = {'forward': {2: {'cpu': 'FP', 'cuda': 'FP_CUDA'}, - 3: {'cpu': None, 'cuda': 'FP3D_CUDA'}}, - 'backward': {2: {'cpu': 'BP', 'cuda': 'BP_CUDA'}, - 3: {'cpu': None, 'cuda': 'BP3D_CUDA'}}} - - algo_cfg = {'type': algo_map[direction][ndim][impl], - 'ProjectorId': proj_id, - 'ProjectionDataId': sino_id} + algo_map = { + 'forward': { + 2: {'cpu': 'FP', 'cuda': 'FP_CUDA'}, + 3: {'cpu': None, 'cuda': 'FP3D_CUDA'} + }, + 'backward': { + 2: {'cpu': 'BP', 'cuda': 'BP_CUDA'}, + 3: {'cpu': None, 'cuda': 'BP3D_CUDA'} + } + } + + algo_cfg = { + 'type': algo_map[direction][ndim][impl], + 'ProjectorId': proj_id, + 'ProjectionDataId': sino_id + } if direction == 'forward': algo_cfg['VolumeDataId'] = vol_id else: diff --git a/odl/tomo/geometry/conebeam.py b/odl/tomo/geometry/conebeam.py index 25eca19a2af..a68ebe3dc6a 100644 --- a/odl/tomo/geometry/conebeam.py +++ b/odl/tomo/geometry/conebeam.py @@ -1964,7 +1964,7 @@ class ConeVecGeometry(VecGeometry): linear paths. """ - # `rotation_matrix` not implemented; reason: missing + # `rotation_matrix` not implemented; reason: does not apply @property def _slice_src(self): @@ -2021,6 +2021,7 @@ def src_position(self, index): >>> geom_3d.src_position(0.5) # mean value array([ 0.5, -0.5, 0. ]) """ + # TODO: Use NN interpolation instead of linear # TODO: vectorize if (self.check_bounds and not is_inside_bounds(index, self.motion_params)): diff --git a/odl/tomo/operators/ray_trafo.py b/odl/tomo/operators/ray_trafo.py index 4c1794fbe80..060cee6fd9e 100644 --- a/odl/tomo/operators/ray_trafo.py +++ b/odl/tomo/operators/ray_trafo.py @@ -125,8 +125,10 @@ def __init__(self, vol_space, geometry, **kwargs): ) ): # Approximate cell volume - angle_weight = float(geometry.motion_partition.extent.prod() / - geometry.motion_partition.size) + angle_weight = float( + np.prod(geometry.motion_partition.extent) + / geometry.motion_partition.size + ) weighting = angle_weight * geometry.det_partition.cell_volume else: raise NotImplementedError('unknown weighting of domain') From 4bf2d7a7a31d36ac66128bc478221028eaf69eb8 Mon Sep 17 00:00:00 2001 From: Holger Kohr Date: Sun, 8 Nov 2020 01:19:20 +0100 Subject: [PATCH 24/24] Fix broadcasting in base VecGeometry methods --- odl/tomo/geometry/geometry.py | 79 ++++++++++++++++++++++------------- 1 file changed, 50 insertions(+), 29 deletions(-) diff --git a/odl/tomo/geometry/geometry.py b/odl/tomo/geometry/geometry.py index a222ba8c26e..ba15552fa08 100644 --- a/odl/tomo/geometry/geometry.py +++ b/odl/tomo/geometry/geometry.py @@ -638,7 +638,7 @@ def rotation_matrix(self, angle): class VecGeometry(Geometry): """Abstract 2D or 3D geometry defined by a collection of vectors. - This geometry gives maximal flexibility for representing locations + This geometry gives maximum flexibility for representing locations and orientations of ray/source and detector for parallel or cone beam acquisition schemes. It is defined by a set of vectors per projection, namely @@ -677,7 +677,7 @@ class VecGeometry(Geometry): linear paths. """ - # `rotation_matrix` not implemented; reason: missing + # `rotation_matrix` not implemented; reason: does not apply # `det_to_src` not implemented; reason: depends on subclass def __init__(self, det_shape, vectors): @@ -694,12 +694,15 @@ def __init__(self, det_shape, vectors): determines the number of projections, and the number of columns the spatial dimension (6 for 2D, 12 for 3D). """ + vectors = self.__vectors = np.asarray(vectors, dtype=float) self.__vectors = np.asarray(vectors, dtype=float) - if self.__vectors.ndim != 2: - raise ValueError('`vectors` must be 2-dimensional, got array ' - 'with ndim = {}'.format(self.__vectors.ndim)) + if vectors.ndim != 2: + raise ValueError( + '`vectors` must be 2-dimensional, got array with ndim = {}' + ''.format(vectors.ndim) + ) - num_projs = self.__vectors.shape[0] + num_projs = vectors.shape[0] if num_projs == 0: raise ValueError('`vectors` is empty') @@ -708,35 +711,41 @@ def __init__(self, det_shape, vectors): elif self.__vectors.shape[1] == 12: ndim = 3 else: - raise ValueError('`vectors` must have 6 or 12 columns, got ' - 'array with {} columns' - ''.format(self.__vectors.shape[1])) + raise ValueError( + '`vectors` must have 6 or 12 columns, got array with ' + '{} columns'.format(vectors.shape[1]) + ) # Make angle partition with spacing 1 - mpart = uniform_partition(0, num_projs - 1, num_projs, - nodes_on_bdry=True) + mpart = uniform_partition( + 0, num_projs - 1, num_projs, nodes_on_bdry=True + ) # Detector is initialized using the first set of vectors. This # assumes that the detector pixels do not change their physical # sizes during acquisition, but has no effect on the computations. # For those, only the `vectors` are used. if ndim == 2: - det_u = self.__vectors[0, 4:6] - det_shape = normalized_scalar_param_list(det_shape, length=1, - param_conv=safe_int_conv) + det_u = vectors[0, 4:6] + det_shape = normalized_scalar_param_list( + det_shape, length=1, param_conv=safe_int_conv + ) det_cell_size = np.linalg.norm(det_u) det_part = uniform_partition( min_pt=-(det_cell_size * det_shape[0]) / 2, max_pt=(det_cell_size * det_shape[0]) / 2, - shape=det_shape) + shape=det_shape + ) detector = Flat1dDetector(det_part, axis=det_u) elif ndim == 3: - det_u = self.__vectors[0, 6:9] - det_v = self.__vectors[0, 9:12] - det_shape = normalized_scalar_param_list(det_shape, length=2, - param_conv=safe_int_conv) + det_u = vectors[0, 6:9] + det_v = vectors[0, 9:12] + det_shape = normalized_scalar_param_list( + det_shape, length=2, param_conv=safe_int_conv + ) det_cell_sides = np.array( - [np.linalg.norm(det_u), np.linalg.norm(det_v)]) + [np.linalg.norm(det_u), np.linalg.norm(det_v)] + ) det_part = uniform_partition( min_pt=-(det_cell_sides * det_shape) / 2, max_pt=(det_cell_sides * det_shape) / 2, @@ -902,10 +911,14 @@ def det_point_position(self, index, dparam): Broadcasting of arguments: + >>> # d(0.4) + 2 * u(0.4) = d(0) + 2 * u(0) + >>> # d(0.6) + 2 * u(0.6) = d(1) + 2 * u(1) >>> idcs = np.array([0.4, 0.6])[:, None] >>> dpar = np.array([2.0])[None, :] >>> geom_2d.det_point_position(idcs, dpar) - + array([[[ 2., 1.]], + + [[-1., 2.]]]) Do the same in 3D, with reference points ``(0, 1, 0)`` and ``(-1, 0, 0)``, and horizontal ``u`` axis vectors ``(1, 0, 0)`` and @@ -934,11 +947,16 @@ def det_point_position(self, index, dparam): Broadcasting of arguments: + >>> # d(0.4) + 2 * u(0.4) + 1 * u(0.4) = d(0) + 2 * u(0) + 1 * v(0) + >>> # d(0.6) + 2 * u(0.6) + 1 * u(0.6) = d(1) + 2 * u(1) + 1 * v(1) >>> idcs = np.array([0.4, 0.6])[:, None] - >>> dpar = np.array([2.0, 1.0])[:, None] + >>> dpar = np.array([2.0, 1.0])[None, :] >>> geom_3d.det_point_position(idcs, dpar) - + array([[[ 2., 1., 1.]], + + [[-1., 2., 1.]]]) """ + dparam = np.asarray(dparam) if self.check_bounds: if not is_inside_bounds(index, self.motion_params): raise ValueError( @@ -955,10 +973,12 @@ def det_point_position(self, index, dparam): if self.ndim == 2: det_shift = dparam * self.det_axis(index) elif self.ndim == 3: + dpar0 = dparam[..., 0] + dpar1 = dparam[..., 1] axes = self.det_axes(index) - det_shift = sum( - di * ax for di, ax in zip(dparam, self.det_axes(index)) - ) + ax0 = axes[..., 0, :] + ax1 = axes[..., 1, :] + det_shift = dpar0 * ax0 + dpar1 * ax1 else: raise RuntimeError('invalid `ndim`') @@ -1021,7 +1041,7 @@ def det_axis(self, index): squeeze_index = False vectors = self.vectors[index_int] - det_us = vectors[:, self._slice_det_u] + det_us = vectors[..., self._slice_det_u] if squeeze_index: return det_us[0] else: @@ -1093,7 +1113,7 @@ def det_axes(self, index): vectors = self.vectors[index_int] if self.ndim == 2: axes = np.empty(index_int.shape + (2,)) - axes[:] = vectors[:, self._slice_det_u] + axes[:] = vectors[..., self._slice_det_u] elif self.ndim == 3: axes = np.empty(index_int.shape + (2, 3)) axes[..., 0, :] = vectors[..., self._slice_det_u] @@ -1159,7 +1179,8 @@ def __getitem__(self, indices): if np.isscalar(sub_vecs) or not hasattr(sub_vecs, 'shape'): raise ValueError( 'indexing of `vectors` results in scalar or non-array ' - '(type {})'.format(type(sub_vecs))) + '(type {})'.format(type(sub_vecs)) + ) if sub_vecs.ndim == 1 and sub_vecs.shape[0] == self.vectors.shape[1]: # Reduced to single vector, add missing dimension