Flying Focal Spot for FanBeam and ConeBeam geometries

odl/test/tomo/geometry/geometry_test.py

@@ -10,6 +10,7 @@
 
 from __future__ import division
 from itertools import permutations, product
+from functools import partial
 import pytest
 import numpy as np
 
@@ -549,10 +550,61 @@ def test_fanbeam_frommatrix():
                                                    det_rad, sing_mat)
 
 
+def test_fanbeam_flying_focal_spot(init1=None):
+    """Test the flying focal spot in 2d fan beam geometry."""
+    full_angle = np.pi
+    n_angles = 2 * 6
+    apart = odl.uniform_partition(0, full_angle, n_angles)
+    dpart = odl.uniform_partition(-1, 1, 11)
+    src_rad = 10
+    det_rad = 5
+    # Source positions with flying focal spot should correspond to
+    # source positions of 2 geometries with different starting positions
+    shift1 = np.array([2.0, -3.0])
+    shift2 = np.array([-2.0, 3.0])
+    init = np.array([1, 0], dtype=np.float32)
+
+    ffs = partial(odl.tomo.flying_focal_spot,
+                  apart=apart,
+                  shifts=[shift1, shift2])
+    geom_ffs = odl.tomo.FanBeamGeometry(
+        apart, dpart,
+        src_rad, det_rad,
+        src_to_det_init=init,
+        src_shift_func=ffs)
+    # angles must be shifted to match discretization of apart
+    ang1 = -full_angle / (n_angles * 2)
+    apart1 = odl.uniform_partition(ang1, full_angle + ang1, n_angles // 2)
+    ang2 = full_angle / (n_angles * 2)
+    apart2 = odl.uniform_partition(ang2, full_angle + ang2, n_angles // 2)
+
+    init1 = init + np.array([0, shift1[1]]) / (src_rad + shift1[0])
+    init2 = init + np.array([0, shift2[1]]) / (src_rad + shift2[0])
+    # radius also changes when a shift is applied
+    src_rad1 = np.linalg.norm(np.array([src_rad, 0]) + shift1)
+    src_rad2 = np.linalg.norm(np.array([src_rad, 0]) + shift2)
+    geom1 = odl.tomo.FanBeamGeometry(apart1, dpart, src_rad1, det_rad,
+                                     src_to_det_init=init1)
+    geom2 = odl.tomo.FanBeamGeometry(apart2, dpart, src_rad2, det_rad,
+                                     src_to_det_init=init2)
+
+    sp1 = geom1.src_position(geom1.angles)
+    sp2 = geom2.src_position(geom2.angles)
+    sp = geom_ffs.src_position(geom_ffs.angles)
+    assert all_almost_equal(sp[0::2], sp1)
+    assert all_almost_equal(sp[1::2], sp2)
+
+    # detector positions are not affected by flying focal spot
+    geom = odl.tomo.FanBeamGeometry(apart, dpart,
+                                    src_rad, det_rad,
+                                    src_to_det_init=init)
+    assert all_almost_equal(geom.det_refpoint(geom.angles),
+                            geom_ffs.det_refpoint(geom_ffs.angles))
+
 def test_helical_cone_beam_props(detector_type, shift):
     """Test basic properties of 3D helical cone beam geometries."""
     full_angle = 2 * np.pi
-    apart = odl.uniform_partition(0, full_angle, 10)
+    apart = odl.uniform_partition(0, full_angle, 13)
     dpart = odl.uniform_partition([0, 0], [1, 1], (10, 10))
     src_rad = 10
     det_rad = 5
@@ -686,6 +738,59 @@ def test_helical_cone_beam_props(detector_type, shift):
     assert repr(geom)
 
 
+def test_conebeam_flying_focal_spot():
+    """Test the flying focal spot in 3d cone beam geometry."""
+    full_angle = np.pi
+    n_angles = 2 * 7
+    apart = odl.uniform_partition(0, full_angle, n_angles)
+    dpart = odl.uniform_partition([-1, -1], [1, 1], (10, 10))
+    src_rad = 10
+    det_rad = 5
+    # Source positions with flying focal spot should correspond to
+    # source positions of 2 geometries with different starting positions
+    shift1 = np.array([2.0, -3.0, 1.0])
+    shift2 = np.array([-2.0, 3.0, -1.0])
+    init = np.array([1, 0, 0], dtype=np.float32)
+    ffs = partial(odl.tomo.flying_focal_spot,
+                  apart=apart,
+                  shifts=[shift1, shift2])
+    geom_ffs = odl.tomo.ConeBeamGeometry(
+        apart, dpart,
+        src_rad, det_rad,
+        src_to_det_init=init,
+        src_shift_func=ffs)
+    # angles must be shifted to match discretization of apart
+    ang1 = -full_angle / (n_angles * 2)
+    apart1 = odl.uniform_partition(ang1, full_angle + ang1, n_angles // 2)
+    ang2 = full_angle / (n_angles * 2)
+    apart2 = odl.uniform_partition(ang2, full_angle + ang2, n_angles // 2)
+
+    init1 = init + np.array([0, shift1[1], 0]) / (src_rad + shift1[0])
+    init2 = init + np.array([0, shift2[1], 0]) / (src_rad + shift2[0])
+    # radius also changes when a shift is applied
+    src_rad1 = np.linalg.norm(np.array([src_rad + shift1[0], shift1[1], 0]))
+    src_rad2 = np.linalg.norm(np.array([src_rad + shift2[0], shift2[1], 0]))
+    geom1 = odl.tomo.ConeBeamGeometry(apart1, dpart, src_rad1, det_rad,
+                                      src_to_det_init=init1,
+                                      offset_along_axis=shift1[2])
+    geom2 = odl.tomo.ConeBeamGeometry(apart2, dpart, src_rad2, det_rad,
+                                      src_to_det_init=init2,
+                                      offset_along_axis=shift2[2])
+
+    sp1 = geom1.src_position(geom1.angles)
+    sp2 = geom2.src_position(geom2.angles)
+    sp = geom_ffs.src_position(geom_ffs.angles)
+    assert all_almost_equal(sp[0::2], sp1)
+    assert all_almost_equal(sp[1::2], sp2)
+
+    # detector positions are not affected by flying focal spot
+    geom = odl.tomo.ConeBeamGeometry(apart, dpart,
+                                     src_rad, det_rad,
+                                     src_to_det_init=init)
+    assert all_almost_equal(geom.det_refpoint(geom.angles),
+                            geom_ffs.det_refpoint(geom_ffs.angles))
+
+
 def test_cone_beam_slanted_detector():
     """Check if non-standard detector axes are handled correctly."""
     full_angle = np.pi
@@ -846,5 +951,41 @@ def test_helical_geometry_helper():
     assert geometry.det_partition.cell_sides[1] <= delta_h
 
 
+def test_source_detector_shifts():
+    """Test source-detector shift functions, e.g. flying focal spot.
+
+    See the `flying_focal_spot` documentation for the exact conditions.
+    """
+    n_angles = np.random.randint(1, 100)
+    apart = odl.uniform_partition(0, np.pi, n_angles)
+    part_angles = apart.meshgrid[0]
+
+    # shifts are periodic
+    def check_shifts(ffs, shifts):
+        i = 0
+        while i < part_angles.size:
+            j = min(len(ffs), i + len(shifts))
+            assert all_almost_equal(ffs[i:j], shifts[:(j - i)])
+            i = j
+
+    # shifts define ffs at partition points
+    n_shifts = np.random.randint(1, n_angles)
+    shift_dim = 3
+    shifts = np.random.uniform(size=(n_shifts, shift_dim))
+    ffs = odl.tomo.flying_focal_spot(part_angles, apart, shifts)
+    check_shifts(ffs, shifts)
+
+    shift_dim = 2
+    shifts = np.random.uniform(size=(n_shifts, shift_dim))
+    ffs = odl.tomo.flying_focal_spot(part_angles, apart, shifts)
+    check_shifts(ffs, shifts)
+
+    # shifts at other angles ar defined by nearest neighbor interpolation
+    d = np.random.uniform(-0.49, 0.49) * apart.cell_volume
+    shifts = np.random.uniform(size=(n_shifts, shift_dim))
+    ffs = odl.tomo.flying_focal_spot(part_angles + d, apart, shifts)
+    check_shifts(ffs, shifts)
+
+
 if __name__ == '__main__':
     odl.util.test_file(__file__)

odl/test/tomo/operators/ray_trafo_test.py

@@ -13,6 +13,7 @@
 import numpy as np
 import pytest
 from packaging.version import parse as parse_version
+from functools import partial
 
 import odl
 from odl.tomo.backends import ASTRA_VERSION
@@ -475,5 +476,138 @@ def test_shifted_volume(geometry_type):
     assert np.max(proj[3, 15:]) > 5
 
 
+def test_source_detector_shifts_2d():
+    """Check that source/detector shifts are handled correctly.
+
+    We forward project a Shepp-Logan phantom and check that reconstruction
+    with flying focal spot is equal to a sum of reconstructions with two
+    geometries which mimic ffs by using initial offsets
+    (the detector must be large enough, not to be influenced by shifts)
+    """
+
+    # If no implementation is available, skip
+    if not odl.tomo.ASTRA_AVAILABLE:
+        pytest.skip(msg='ASTRA not available, skipping 2d test')
+
+    space = odl.uniform_discr([-1] * 2, [1] * 2, [10] * 2)
+    phantom = odl.phantom.shepp_logan(space)
+
+    full_angle = np.pi
+    n_angles = 2 * 6
+    apart = odl.uniform_partition(0, full_angle, n_angles)
+    dpart = odl.uniform_partition(-2, 2, 11)
+    src_rad = 10
+    det_rad = 10
+    # Source positions with flying focal spot should correspond to
+    # source positions of 2 geometries with different starting positions
+    shift1 = np.array([2.0, -3.0])
+    shift2 = np.array([-2.0, 3.0])
+    init = np.array([1, 0], dtype=np.float32)
+
+    ffs = partial(odl.tomo.flying_focal_spot,
+                  apart=apart,
+                  shifts=[shift1, shift2])
+    geom_ffs = odl.tomo.FanBeamGeometry(
+        apart, dpart,
+        src_rad, det_rad,
+        src_to_det_init=init,
+        src_shift_func=ffs)
+    # angles must be shifted to match discretization of apart
+    ang1 = -full_angle / (n_angles * 2)
+    apart1 = odl.uniform_partition(ang1, full_angle + ang1, n_angles // 2)
+    ang2 = full_angle / (n_angles * 2)
+    apart2 = odl.uniform_partition(ang2, full_angle + ang2, n_angles // 2)
+
+    init1 = init + np.array([0, shift1[1]]) / (src_rad + shift1[0])
+    init2 = init + np.array([0, shift2[1]]) / (src_rad + shift2[0])
+    # radius also changes when a shift is applied
+    src_rad1 = np.linalg.norm(np.array([src_rad, 0]) + shift1)
+    src_rad2 = np.linalg.norm(np.array([src_rad, 0]) + shift2)
+    geom1 = odl.tomo.FanBeamGeometry(apart1, dpart, src_rad1, det_rad,
+                                     src_to_det_init=init1)
+    geom2 = odl.tomo.FanBeamGeometry(apart2, dpart, src_rad2, det_rad,
+                                     src_to_det_init=init2)
+
+    sp1 = geom1.src_position(geom1.angles)
+    sp2 = geom2.src_position(geom2.angles)
+    sp = geom_ffs.src_position(geom_ffs.angles)
+    assert all_almost_equal(sp[0::2], sp1)
+    assert all_almost_equal(sp[1::2], sp2)
+
+    op_ffs = odl.tomo.RayTransform(space, geom_ffs)
+    op1 = odl.tomo.RayTransform(space, geom1)
+    op2 = odl.tomo.RayTransform(space, geom2)
+    im = op_ffs.adjoint(op_ffs(phantom))
+    im_combined = op1.adjoint(op1(phantom)) + op2.adjoint(op2(phantom))
+    assert all_almost_equal(im.asarray(), im_combined.asarray())
+
+
+def test_source_detector_shifts_3d():
+    """Check that source/detector shifts are handled correctly.
+
+    We forward project a Shepp-Logan phantom and check that reconstruction
+    with flying focal spot is equal to a sum of reconstructions with two
+    geometries which mimic ffs by using initial offsets
+    (the detector must be large enough, not to be influenced by shifts)
+    """
+
+    # If no implementation is available, skip
+    if not odl.tomo.ASTRA_CUDA_AVAILABLE:
+        pytest.skip(msg='ASTRA_CUDA not available, skipping 3d test')
+
+    space = odl.uniform_discr([-1] * 3, [1] * 3, [10] * 3)
+    phantom = odl.phantom.shepp_logan(space)
+
+    full_angle = np.pi
+    n_angles = 2 * 7
+    apart = odl.uniform_partition(0, full_angle, n_angles)
+    dpart = odl.uniform_partition([-2, -2], [2, 2], (11, 11))
+    src_rad = 10
+    det_rad = 10
+    # Source positions with flying focal spot should correspond to
+    # source positions of 2 geometries with different starting positions
+    shift1 = np.array([2.0, -3.0, 1.0])
+    shift2 = np.array([-2.0, 3.0, -1.0])
+    init = np.array([1, 0, 0], dtype=np.float32)
+    ffs = partial(odl.tomo.flying_focal_spot,
+                  apart=apart,
+                  shifts=[shift1, shift2])
+    geom_ffs = odl.tomo.ConeBeamGeometry(
+        apart, dpart,
+        src_rad, det_rad,
+        src_to_det_init=init,
+        src_shift_func=ffs)
+    # angles must be shifted to match discretization of apart
+    ang1 = -full_angle / (n_angles * 2)
+    apart1 = odl.uniform_partition(ang1, full_angle + ang1, n_angles // 2)
+    ang2 = full_angle / (n_angles * 2)
+    apart2 = odl.uniform_partition(ang2, full_angle + ang2, n_angles // 2)
+
+    init1 = init + np.array([0, shift1[1], 0]) / (src_rad + shift1[0])
+    init2 = init + np.array([0, shift2[1], 0]) / (src_rad + shift2[0])
+    # radius also changes when a shift is applied
+    src_rad1 = np.linalg.norm(np.array([src_rad + shift1[0], shift1[1], 0]))
+    src_rad2 = np.linalg.norm(np.array([src_rad + shift2[0], shift2[1], 0]))
+    geom1 = odl.tomo.ConeBeamGeometry(apart1, dpart, src_rad1, det_rad,
+                                      src_to_det_init=init1,
+                                      offset_along_axis=shift1[2])
+    geom2 = odl.tomo.ConeBeamGeometry(apart2, dpart, src_rad2, det_rad,
+                                      src_to_det_init=init2,
+                                      offset_along_axis=shift2[2])
+
+    sp1 = geom1.src_position(geom1.angles)
+    sp2 = geom2.src_position(geom2.angles)
+    sp = geom_ffs.src_position(geom_ffs.angles)
+    assert all_almost_equal(sp[0::2], sp1)
+    assert all_almost_equal(sp[1::2], sp2)
+
+    op_ffs = odl.tomo.RayTransform(space, geom_ffs)
+    op1 = odl.tomo.RayTransform(space, geom1)
+    op2 = odl.tomo.RayTransform(space, geom2)
+    im = op_ffs.adjoint(op_ffs(phantom))
+    im_combined = op1.adjoint(op1(phantom)) + op2.adjoint(op2(phantom))
+    assert all_almost_equal(im.asarray(), im_combined.asarray())
+
+
 if __name__ == '__main__':
     odl.util.test_file(__file__)

odl/tomo/__init__.py

@@ -16,15 +16,17 @@
     astra_conebeam_2d_geom_to_vec, astra_conebeam_3d_geom_to_vec)
 from .geometry import *
 from .operators import *
+from .util import *
 
 __all__ = ()
 __all__ += analytic.__all__
 __all__ += geometry.__all__
 __all__ += operators.__all__
+__all__ += util.source_detector_shifts.__all__
 __all__ += (
     'ASTRA_AVAILABLE',
     'ASTRA_CUDA_AVAILABLE',
     'SKIMAGE_AVAILABLE',
     'astra_conebeam_2d_geom_to_vec',
-    'astra_conebeam_3d_geom_to_vec',
+    'astra_conebeam_3d_geom_to_vec'
 )