From 6fa45d5eecf270de9b51e70d4b64b4b6ae1c8fd0 Mon Sep 17 00:00:00 2001 From: Holger Kohr Date: Fri, 9 Jun 2017 00:24:44 +0200 Subject: [PATCH] ENH: add vector geometries, closes #1023 --- examples/tomo/ray_trafo_vec_geom_3d.py | 78 +++++ odl/tomo/backends/astra_setup.py | 124 ++++++- odl/tomo/geometry/conebeam.py | 170 +++++++++- odl/tomo/geometry/geometry.py | 442 ++++++++++++++++++++++++- odl/tomo/geometry/parallel.py | 99 +++++- odl/tomo/operators/ray_trafo.py | 23 +- 6 files changed, 904 insertions(+), 32 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 174c36a8bb6..dc4c06a4dca 100644 --- a/odl/tomo/backends/astra_setup.py +++ b/odl/tomo/backends/astra_setup.py @@ -1,4 +1,4 @@ -# Copyright 2014-2017 The ODL contributors +# Copyright 2014-2019 The ODL contributors # # This file is part of ODL. # @@ -37,6 +37,7 @@ from odl.discr import DiscreteLp, DiscreteLpElement from odl.tomo.geometry import ( Geometry, DivergentBeamGeometry, ParallelBeamGeometry, + ParallelVecGeometry, ConeVecGeometry, Flat1dDetector, Flat2dDetector) from odl.tomo.util.utility import euler_matrix from odl.util.npy_compat import moveaxis @@ -230,6 +231,75 @@ def astra_volume_geometry(reco_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. @@ -288,14 +358,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): @@ -354,7 +417,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): @@ -416,13 +479,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): @@ -472,6 +529,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): @@ -491,6 +564,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 e4da61f4e23..63df47242a8 100644 --- a/odl/tomo/geometry/conebeam.py +++ b/odl/tomo/geometry/conebeam.py @@ -16,14 +16,19 @@ Flat1dDetector, Flat2dDetector, CircularDetector, CylindricalDetector, SphericalDetector) from odl.tomo.geometry.geometry import ( - DivergentBeamGeometry, AxisOrientedGeometry) + DivergentBeamGeometry, AxisOrientedGeometry, VecGeometry) from odl.tomo.util.utility import ( euler_matrix, transform_system, is_inside_bounds) from odl.util import signature_string, indent, array_str -__all__ = ('FanBeamGeometry', 'ConeBeamGeometry', - 'cone_beam_geometry', 'helical_geometry') +__all__ = ( + 'FanBeamGeometry', + 'ConeBeamGeometry', + 'ConeVecGeometry', + 'cone_beam_geometry', + 'helical_geometry', +) class FanBeamGeometry(DivergentBeamGeometry): @@ -1678,6 +1683,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..1101c2f1ce6 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. # @@ -12,12 +12,16 @@ 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 ( + normalized_scalar_param_list, safe_int_conv, signature_string, + indent_rows) -__all__ = ('Geometry', 'DivergentBeamGeometry', 'AxisOrientedGeometry') +__all__ = ('Geometry', 'DivergentBeamGeometry', 'AxisOrientedGeometry', + 'VecGeometry') class Geometry(object): @@ -619,6 +623,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 f39e9c1db28..9f52a910fda 100644 --- a/odl/tomo/geometry/parallel.py +++ b/odl/tomo/geometry/parallel.py @@ -13,7 +13,8 @@ from odl.discr import uniform_partition from odl.tomo.geometry.detector import Flat1dDetector, Flat2dDetector -from odl.tomo.geometry.geometry import Geometry, AxisOrientedGeometry +from odl.tomo.geometry.geometry import ( + Geometry, AxisOrientedGeometry, VecGeometry) from odl.tomo.util import euler_matrix, transform_system, is_inside_bounds from odl.util import signature_string, indent, array_str @@ -21,6 +22,7 @@ __all__ = ('ParallelBeamGeometry', 'Parallel2dGeometry', 'Parallel3dEulerGeometry', 'Parallel3dAxisGeometry', + 'ParallelVecGeometry', 'parallel_beam_geometry') @@ -1468,6 +1470,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 8175f78e1e3..382b60cfe43 100644 --- a/odl/tomo/operators/ray_trafo.py +++ b/odl/tomo/operators/ray_trafo.py @@ -16,7 +16,8 @@ from odl.operator import Operator from odl.space import FunctionSpace from odl.tomo.geometry import ( - Geometry, Parallel2dGeometry, Parallel3dAxisGeometry) + Geometry, Parallel2dGeometry, Parallel3dAxisGeometry, + ParallelVecGeometry) from odl.space.weighting import ConstWeighting from odl.tomo.backends import ( ASTRA_AVAILABLE, ASTRA_CUDA_AVAILABLE, SKIMAGE_AVAILABLE, @@ -214,18 +215,26 @@ def __init__(self, reco_space, geometry, variant, **kwargs): dtype = reco_space.dtype proj_fspace = FunctionSpace(geometry.params, out_dtype=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 reco_space.is_weighted: weighting = None + if isinstance(self.geometry, ParallelVecGeometry): + # TODO: change to weighting constant per angle when available + weighting = 1.0 elif (isinstance(reco_space.weighting, ConstWeighting) and np.isclose(reco_space.weighting.const, reco_space.cell_volume)): # 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