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..0be2dc41755 --- /dev/null +++ b/examples/tomo/checks/check_axes_cone2d_vec_fp.py @@ -0,0 +1,99 @@ +"""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", 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 new file mode 100644 index 00000000000..655270bf365 --- /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 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/ray_trafo_vec_geom_3d.py b/examples/tomo/ray_trafo_vec_geom_3d.py new file mode 100644 index 00000000000..d30316a6753 --- /dev/null +++ b/examples/tomo/ray_trafo_vec_geom_3d.py @@ -0,0 +1,78 @@ +"""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 +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.ConeFlatGeometry( + angle_partition, detector_partition, src_radius=1000, det_radius=100, + axis=[1, 0, 0]) + +circle_vecs = odl.tomo.cone_3d_geom_to_astra_vecs(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/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 {})' 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) diff --git a/odl/test/tomo/backends/astra_setup_test.py b/odl/test/tomo/backends/astra_setup_test.py index d3502cf5bc2..5d5112c518a 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 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/__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' ) 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_cuda.py b/odl/tomo/backends/astra_cuda.py index 03b6eb877ae..4388dd7e47a 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 @@ -238,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) @@ -337,19 +344,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_extent = float(np.prod(proj_space.partition.extent)) 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 +394,23 @@ 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, 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): + 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,15 +436,23 @@ 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, 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): + 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 - # 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): @@ -450,6 +485,26 @@ 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, 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: 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 + ) + + 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 c847f1f2574..a64c287b70b 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: @@ -64,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', @@ -234,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 ( @@ -261,39 +267,370 @@ 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)) return vol_geom -def astra_conebeam_3d_geom_to_vec(geometry): - """Create vectors for ASTRA projection geometries from ODL geometry. +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, coords='ASTRA') + 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, coords='ASTRA') + 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, coords='ASTRA') + 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, coords='ASTRA') + 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 coordinate conversion ODL→ASTRA is a transposed rotation + by +90 degrees in 2D, and an XZ axis swap in 3D. + + 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.array(vecs, dtype=float, copy=True) + + if vecs.shape[1] == 6: + # 2D geometry + # ODL x = ASTRA -y, ODL y = ASTRA x + # Thus: ODL->ASTRA conversion is a transposed rotation by +90 degrees + rot_90 = euler_matrix(np.pi / 2) + # 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]) + 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 swapping x and z + 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_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 coordinate conversion ASTRA->ODL is a transposed rotation + by -90 degrees in 2D, and an XZ axis swap in 3D. + + 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 + # ODL x = ASTRA -y, ODL y = ASTRA x + # Thus: ODL->ASTRA conversion is a transposed rotation by -90 degrees + rot_minus_90 = euler_matrix(-np.pi / 2) + # 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]) + 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 swapping x and z + 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 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 + parallel beam geometries, see ``'parallel_vec'`` in the + `ASTRA projection geometry documentation`_. + + 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 + + .. 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 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)) + + # Ray direction = -(detector-to-source normal vector) + vectors[:, 0:2] = -geometry.det_to_src(angles, mid_pt) + + # Center of detector + # 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 = geometry.det_axis(angles) + px_size = geometry.det_partition.cell_sides[0] + vectors[:, 4:6] = det_axis * px_size + + if coords == 'ASTRA': + vectors = vecs_odl_to_astra_coords(vectors) + + return vectors + + +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 - cone beam geometries, see ``'cone_vec'`` in the + 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 :: - (srcX, srcY, srcZ, dX, dY, dZ, uX, uY, uZ, vX, vY, vZ) + (rayX, rayY, rayZ, dX, dY, dZ, uX, uY, uZ, vX, vY, vZ) with - - ``src``: the ray source position + - ``ray``: the ray direction - ``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)`` + .. 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. + coords : {'ODL', 'ASTRA'}, optional + Coordinate system in which the resulting vectors should be + represented. Returns ------- @@ -305,49 +642,51 @@ def astra_conebeam_3d_geom_to_vec(geometry): .. _ASTRA projection geometry documentation: http://www.astra-toolbox.com/docs/geom3d.html#projection-geometries """ - angles = geometry.angles - vectors = np.zeros((angles.size, 12)) + coords, coords_in = str(coords).upper(), coords + if coords not in ('ODL', 'ASTRA'): + raise ValueError('`coords` {!r} not understood'.format(coords_in)) - # Source position - vectors[:, 0:3] = geometry.src_position(angles) + # TODO: use `coords` - # Center of detector in 3D space + angles = geometry.angles mid_pt = geometry.det_params.mid_pt + + vectors = np.zeros((angles.shape[-1], 12)) + + # Ray direction = -(detector-to-source normal vector) + vectors[:, 0:3] = -geometry.det_to_src(angles, mid_pt) + + # Center of the detector in 3D space 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. + # 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`. + # + # 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] - # 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_to_astra_coords(vectors) -def astra_conebeam_2d_geom_to_vec(geometry): - """Create vectors for ASTRA projection geometries from ODL 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 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) @@ -357,10 +696,17 @@ def astra_conebeam_2d_geom_to_vec(geometry): - ``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. + coords : {'ODL', 'ASTRA'}, optional + Coordinate system in which the resulting vectors should be + represented. Returns ------- @@ -372,55 +718,62 @@ 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) + 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)) # Source position - src_pos = geometry.src_position(angles) - vectors[:, 0:2] = rot_minus_90.dot(src_pos.T).T # dot along 2nd axis + vectors[:, 0:2] = geometry.src_position(angles) # 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 + if coords == 'ASTRA': + vectors = vecs_odl_to_astra_coords(vectors) + return vectors -def astra_parallel_3d_geom_to_vec(geometry): - """Create vectors for ASTRA projection geometries from ODL 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 - parallel beam geometries, see ``'parallel3d_vec'`` in the + 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 :: - (rayX, rayY, rayZ, dX, dY, dZ, uX, uY, uZ, vX, vY, vZ) + (srcX, srcY, srcZ, dX, dY, dZ, uX, uY, uZ, vX, vY, vZ) with - - ``ray``: the ray direction + - ``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)`` + .. 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. + coords : {'ODL', 'ASTRA'}, optional + Coordinate system in which the resulting vectors should be + represented. Returns ------- @@ -432,15 +785,20 @@ def astra_parallel_3d_geom_to_vec(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)) - vectors = np.zeros((angles.shape[-1], 12)) - - # Ray direction = -(detector-to-source normal vector) - vectors[:, 0:3] = -geometry.det_to_src(angles, mid_pt) + # Source position + vectors[:, 0:3] = geometry.src_position(angles) - # Center of the detector in 3D space + # Center of detector in 3D space 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) @@ -454,93 +812,12 @@ def astra_parallel_3d_geom_to_vec(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] - # 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 - - -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. - - 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') - - 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) - - 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) - - 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) - - 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) - 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 + return vecs_odl_to_astra_coords(vectors) def astra_data(astra_geom, datatype, data=None, ndim=2, allow_copy=False): @@ -573,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) @@ -593,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: @@ -608,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) @@ -637,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) @@ -687,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 @@ -734,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 832d3f1a24e..a68ebe3dc6a 100644 --- a/odl/tomo/geometry/conebeam.py +++ b/odl/tomo/geometry/conebeam.py @@ -14,20 +14,24 @@ 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): - """Fan beam (2d cone beam) geometry. The source moves on a circle with radius ``src_radius``, and the @@ -36,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 @@ -708,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 @@ -1207,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`. @@ -1702,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) @@ -1744,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: @@ -1872,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) @@ -1919,6 +1926,172 @@ 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. + """ + + # `rotation_matrix` not implemented; reason: does not apply + + @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. ]) + """ + # TODO: Use NN interpolation instead of linear + # 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)) + + 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 + ---------- + mparam : `motion_params` element + Motion parameter at which to evaluate. + dparam : `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. + """ + # TODO: vectorize + if self.check_bounds: + if not is_inside_bounds(mparam, self.motion_params): + raise ValueError('`mparam` {} not in the valid range {}' + ''.format(mparam, 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)) + + 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..ba15552fa08 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-2020 The ODL contributors # # This file is part of ODL. # @@ -8,20 +8,28 @@ """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 ( + REPR_PRECISION, normalized_scalar_param_list, npy_printoptions, + repr_string, safe_int_conv, signature_string_parts) - -__all__ = ('Geometry', 'DivergentBeamGeometry', 'AxisOrientedGeometry') +__all__ = ( + 'Geometry', + 'DivergentBeamGeometry', + 'AxisOrientedGeometry', + 'VecGeometry', +) class Geometry(object): - """Abstract geometry class. A geometry is described by @@ -410,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. @@ -453,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 ------- @@ -543,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) @@ -556,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. @@ -607,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)) @@ -619,6 +635,575 @@ def rotation_matrix(self, angle): return matrix +class VecGeometry(Geometry): + """Abstract 2D or 3D geometry defined by a collection of vectors. + + 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 + + - 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. + """ + + # `rotation_matrix` not implemented; reason: does not apply + # `det_to_src` not implemented; reason: depends on subclass + + 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 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). + """ + vectors = self.__vectors = np.asarray(vectors, dtype=float) + self.__vectors = np.asarray(vectors, dtype=float) + if vectors.ndim != 2: + raise ValueError( + '`vectors` must be 2-dimensional, got array with ndim = {}' + ''.format(vectors.ndim) + ) + + num_projs = 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(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 = 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 = 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)] + ) + 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]) + else: + raise RuntimeError('invalid `ndim`') + + 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.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)``: + + >>> # 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.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_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)``. + + 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.4) + 2 * u(0.4) = d(0) + 2 * u(0) + >>> geom_2d.det_point_position(0.4, 2) + array([ 2., 1.]) + + 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 + ``(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.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: + + >>> # 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, :] + >>> 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( + '`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) + ) + + 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) + ax0 = axes[..., 0, :] + ax1 = axes[..., 1, :] + det_shift = dpar0 * ax0 + dpar1 * ax1 + else: + raise RuntimeError('invalid `ndim`') + + 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.4, 0.6]) # values at closest indices + array([[ 1., 0.], + [ 0., 1.]]) + """ + 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' + ) + + index_int = np.round(index).astype(int) + if index_int.shape == (): + index_int = index_int[None] + squeeze_index = True + else: + squeeze_index = False + + 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). + + Parameters + ---------- + index : `motion_params` element + Index of the projection. Non-integer indices are rounded to + closest integer. + + Returns + ------- + 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 + -------- + 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.], + [ 0., 0., 1.]]) + >>> geom_3d.det_axes(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 + 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: + squeeze_index = False + + vectors = self.vectors[index_int] + if self.ndim == 2: + axes = np.empty(index_int.shape + (2,)) + 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] + axes[..., 1, :] = vectors[..., self._slice_det_v] + else: + raise RuntimeError('invalid `ndim`') + + if squeeze_index: + axes = axes[0] + return axes + + 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] + 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 odl.util.testutils import run_doctests run_doctests() diff --git a/odl/tomo/geometry/parallel.py b/odl/tomo/geometry/parallel.py index 1d7d3dad6b8..5b632220679 100644 --- a/odl/tomo/geometry/parallel.py +++ b/odl/tomo/geometry/parallel.py @@ -14,18 +14,22 @@ 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): - """Abstract parallel beam geometry in 2 or 3 dimensions. Parallel geometries are characterized by a virtual source at @@ -326,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 @@ -630,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)) @@ -664,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:: @@ -700,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 @@ -1039,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)) @@ -1074,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 @@ -1468,6 +1473,225 @@ 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. + """ + + # `rotation_matrix` not implemented; reason: missing + + @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): + """Direction from a detector location to the source. + + 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 : `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 + ------- + 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) + """ + 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) + + 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: + squeeze_dparam = (np.broadcast(*dparam).shape == ()) + dparam = tuple(np.array(p, dtype=float, copy=False, ndmin=1) + for p in dparam) + + + 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): r"""Create default parallel beam geometry from ``space``. @@ -1572,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)) diff --git a/odl/tomo/operators/ray_trafo.py b/odl/tomo/operators/ray_trafo.py index 64419fd43b6..060cee6fd9e 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, VecGeometry 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, VecGeometry): + # No angular weighting + weighting = geometry.det_partition.cell_volume elif ( isinstance(vol_space.weighting, ConstWeighting) and np.isclose( @@ -110,15 +125,11 @@ 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 - extent = float(geometry.partition.extent.prod()) - size = float(geometry.partition.size) - weighting = extent / 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') diff --git a/odl/util/utility.py b/odl/util/utility.py index 9b5732d91e8..14c3f1d9ff2 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,13 @@ 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