Merge pull request #1157 from adler-j/issue-1054__helical_helper

Add helical_geometry

Author: adler-j

examples/tomo/filtered_backprojection_helical_3d.py

@@ -17,29 +17,21 @@
 
 
 # Reconstruction space: discretized functions on the cube
-# [-20, 20]^2 x [0, 40] with 300 samples per dimension.
-reco_space = odl.uniform_discr(
-    min_pt=[-20, -20, 0], max_pt=[20, 20, 40], shape=[300, 300, 300],
+# [-20, 20]^3  with 200 samples per dimension.
+space = odl.uniform_discr(
+    min_pt=[-20, -20, -20], max_pt=[20, 20, 20], shape=[200, 200, 200],
     dtype='float32')
 
-# Make a helical cone beam geometry with flat detector
-# Angles: uniformly spaced, n = 2000, min = 0, max = 8 * 2 * pi
-# This gives 8 full turns of the helix.
-angle_partition = odl.uniform_partition(0, 8 * 2 * np.pi, 2000)
-# Detector: uniformly sampled with a small height,
-# n = (512, 64), min = (-30, -4), max = (30, 4)
-detector_partition = odl.uniform_partition([-40, -4], [40, 4], [512, 64])
-# Create geometry
-geometry = odl.tomo.ConeFlatGeometry(
-    angle_partition, detector_partition, src_radius=100, det_radius=100,
-    pitch=5.0)
-
+# Create helical geometry
+geometry = odl.tomo.helical_geometry(space,
+                                     src_radius=100, det_radius=100,
+                                     num_turns=7.5, num_angles=1000)
 
 # --- Create Filtered Back-Projection (FBP) operator --- #
 
 
 # Ray transform (= forward projection).
-ray_trafo = odl.tomo.RayTransform(reco_space, geometry)
+ray_trafo = odl.tomo.RayTransform(space, geometry)
 
 # Unwindowed fbp
 # We select a Hamming filter, and only use the lowest 80% of frequencies to
@@ -54,7 +46,7 @@
 
 
 # Create a discrete Shepp-Logan phantom (modified version)
-phantom = odl.phantom.shepp_logan(reco_space, modified=True)
+phantom = odl.phantom.shepp_logan(space, modified=True)
 
 # Create projection data by calling the ray transform on the phantom
 proj_data = ray_trafo(phantom)

odl/test/solvers/functional/default_functionals_test.py

@@ -56,8 +56,8 @@ def test_L1_norm(space, sigma):
     x = noise_element(space)
 
     # Test functional evaluation
-    expected_result = func.domain.element(np.abs(x)).inner(space.one())
-    assert pytest.approx(func(x), expected_result)
+    expected_result = np.abs(x).inner(space.one())
+    assert func(x) == pytest.approx(expected_result)
 
     # Test gradient - expecting sign function
     expected_result = func.domain.element(np.sign(x))
@@ -125,7 +125,7 @@ def test_L2_norm(space, sigma):
 
     # Test functional evaluation
     expected_result = np.sqrt((x ** 2).inner(space.one()))
-    assert pytest.approx(func(x), expected_result)
+    assert func(x) == pytest.approx(expected_result)
 
     # Test gradient
     if x_norm > 0:
@@ -172,7 +172,7 @@ def test_L2_norm(space, sigma):
 
     # Verify that the biconjugate is the functional itself
     func_cc_cc = func_cc.convex_conj
-    assert pytest.approx(func_cc_cc(x), func(x))
+    assert func_cc_cc(x) == pytest.approx(func(x))
 
 
 def test_L2_norm_squared(space, sigma):
@@ -183,7 +183,7 @@ def test_L2_norm_squared(space, sigma):
 
     # Test functional evaluation
     expected_result = x_norm ** 2
-    assert pytest.approx(func(x), expected_result)
+    assert func(x) == pytest.approx(expected_result)
 
     # Test gradient
     expected_result = 2.0 * x
@@ -198,7 +198,7 @@ def test_L2_norm_squared(space, sigma):
 
     # Test evaluation of the convex conjugate
     expected_result = x_norm ** 2 / 4.0
-    assert pytest.approx(func_cc(x), expected_result)
+    assert func_cc(x) == pytest.approx(expected_result)
 
     # Test gradient of the convex conjugate
     expected_result = x / 2.0
@@ -212,7 +212,7 @@ def test_L2_norm_squared(space, sigma):
     func_cc_cc = func_cc.convex_conj
 
     # Check that they evaluate to the same value
-    assert pytest.approx(func_cc_cc(x), func(x))
+    assert func_cc_cc(x) == pytest.approx(func(x))
 
     # Check that their gradients evaluate to the same value
     assert all_almost_equal(func_cc_cc.gradient(x), func.gradient(x))
@@ -281,7 +281,7 @@ def test_kullback_leibler(space):
         # Evaluation of the functional
         expected_result = ((x - prior + prior * np.log(prior / x))
                            .inner(one_elem))
-        assert pytest.approx(func(x), expected_result)
+        assert func(x) == pytest.approx(expected_result)
 
     # Check property for prior
     assert all_almost_equal(func.prior, prior)
@@ -317,7 +317,7 @@ def test_kullback_leibler(space):
     # Evaluation of convex conjugate
     with np.errstate(all='ignore'):
         expected_result = - (prior * np.log(1 - x)).inner(one_elem)
-        assert pytest.approx(cc_func(x), expected_result)
+        assert cc_func(x) == pytest.approx(expected_result)
 
     x_wrong = noise_element(space)
     x_wrong = x_wrong - x_wrong.ufuncs.max() + 1.01
@@ -340,7 +340,7 @@ def test_kullback_leibler(space):
 
     # Check that they evaluate the same
     with np.errstate(all='ignore'):
-        assert pytest.approx(cc_cc_func(x), func(x))
+        assert cc_cc_func(x) == pytest.approx(func(x))
 
 
 def test_kullback_leibler_cross_entorpy(space):
@@ -359,7 +359,7 @@ def test_kullback_leibler_cross_entorpy(space):
     # Evaluation of the functional
     expected_result = ((prior - x + x * np.log(x / prior))
                        .inner(one_elem))
-    assert pytest.approx(func(x), expected_result)
+    assert func(x) == pytest.approx(expected_result)
 
     # Check property for prior
     assert all_almost_equal(func.prior, prior)
@@ -390,7 +390,7 @@ def test_kullback_leibler_cross_entorpy(space):
 
     # Evaluation of convex conjugate
     expected_result = (prior * (np.exp(x) - 1)).inner(one_elem)
-    assert pytest.approx(cc_func(x), expected_result)
+    assert cc_func(x) == pytest.approx(expected_result)
 
     # The gradient of the convex conjugate
     expected_result = prior * np.exp(x)
@@ -406,7 +406,7 @@ def test_kullback_leibler_cross_entorpy(space):
     cc_cc_func = cc_func.convex_conj
 
     # Check that they evaluate the same
-    assert pytest.approx(cc_cc_func(x), func(x))
+    assert cc_cc_func(x) == pytest.approx(func(x))
 
 
 def test_quadratic_form(space):
@@ -425,7 +425,7 @@ def test_quadratic_form(space):
 
     # Evaluation of the functional
     expected_result = x.inner(operator(x)) + vector.inner(x) + constant
-    assert pytest.approx(func(x), expected_result)
+    assert func(x) == pytest.approx(expected_result)
 
     # The gradient
     expected_gradient = 2 * operator(x) + vector
@@ -438,7 +438,7 @@ def test_quadratic_form(space):
     func_no_operator = odl.solvers.QuadraticForm(vector=vector,
                                                  constant=constant)
     expected_result = vector.inner(x) + constant
-    assert pytest.approx(func_no_operator(x), expected_result)
+    assert func_no_operator(x) == pytest.approx(expected_result)
 
     expected_gradient = vector
     assert all_almost_equal(func_no_operator.gradient(x), expected_gradient)
@@ -455,7 +455,7 @@ def test_quadratic_form(space):
     # Test with no offset
     func_no_offset = odl.solvers.QuadraticForm(operator, constant=constant)
     expected_result = x.inner(operator(x)) + constant
-    assert pytest.approx(func_no_offset(x), expected_result)
+    assert func_no_offset(x) == pytest.approx(expected_result)
 
 
 def test_separable_sum(space):

odl/test/tomo/geometry/geometry_test.py

@@ -376,12 +376,12 @@ def test_parallel_beam_geometry_helper():
 
     # Validate angles
     assert geometry.motion_partition.is_uniform
-    assert pytest.approx(geometry.motion_partition.extent, np.pi)
+    assert geometry.motion_partition.extent == pytest.approx(np.pi)
     assert geometry.motion_partition.cell_sides <= np.pi / (rho * omega)
 
     # Validate detector
     assert geometry.det_partition.cell_sides <= np.pi / omega
-    assert pytest.approx(geometry.det_partition.extent, 2 * rho)
+    assert geometry.det_partition.extent == pytest.approx(2 * rho)
 
     # --- 3d case ---
 
@@ -390,18 +390,18 @@ def test_parallel_beam_geometry_helper():
 
     # Validate angles
     assert geometry.motion_partition.is_uniform
-    assert pytest.approx(geometry.motion_partition.extent, np.pi)
+    assert geometry.motion_partition.extent == pytest.approx(np.pi)
     assert geometry.motion_partition.cell_sides <= np.pi / (rho * omega)
 
     # Validate detector
     assert geometry.det_partition.cell_sides[0] <= np.pi / omega
-    assert pytest.approx(geometry.det_partition.cell_sides[0], 0.05)
-    assert pytest.approx(geometry.det_partition.extent[0], 2 * rho)
+    assert geometry.det_partition.cell_sides[1] == pytest.approx(0.05)
+    assert geometry.det_partition.extent[0] == pytest.approx(2 * rho)
 
     # Validate that new detector axis is correctly aligned with the rotation
     # axis and that there is one detector row per slice
-    assert pytest.approx(geometry.det_partition.min_pt[1], 0.0)
-    assert pytest.approx(geometry.det_partition.max_pt[1], 2.0)
+    assert geometry.det_partition.min_pt[1] == pytest.approx(0.0)
+    assert geometry.det_partition.max_pt[1] == pytest.approx(2.0)
     assert geometry.det_partition.shape[1] == space.shape[2]
     assert all_equal(geometry.det_axes_init[1], geometry.axis)
 
@@ -414,12 +414,12 @@ def test_parallel_beam_geometry_helper():
 
     # Validate angles
     assert geometry.motion_partition.is_uniform
-    assert pytest.approx(geometry.motion_partition.extent, np.pi)
+    assert geometry.motion_partition.extent == pytest.approx(np.pi)
     assert geometry.motion_partition.cell_sides <= np.pi / (rho * omega)
 
     # Validate detector
     assert geometry.det_partition.cell_sides <= np.pi / omega
-    assert pytest.approx(geometry.det_partition.extent, 2 * rho)
+    assert geometry.det_partition.extent == pytest.approx(2 * rho)
 
 
 def test_fanflat_props(shift):
@@ -653,21 +653,21 @@ def test_cone_beam_geometry_helper():
 
     # Validate angles
     assert geometry.motion_partition.is_uniform
-    assert pytest.approx(geometry.motion_partition.extent, 2 * np.pi)
+    assert geometry.motion_partition.extent == pytest.approx(2 * np.pi)
     assert (geometry.motion_partition.cell_sides <=
             (r + rho) / r * np.pi / (rho * omega))
 
     # Validate detector
     det_width = 2 * magnification * rho
     R = np.hypot(r, det_width / 2)
     assert geometry.det_partition.cell_sides <= np.pi * R / (r * omega)
-    assert pytest.approx(geometry.det_partition.extent, det_width)
+    assert geometry.det_partition.extent == pytest.approx(det_width)
 
     # Short scan option
     fan_angle = 2 * np.arctan(det_width / (2 * r))
     geometry = odl.tomo.cone_beam_geometry(space, src_radius, det_radius,
                                            short_scan=True)
-    assert pytest.approx(geometry.motion_params.extent, np.pi + fan_angle)
+    assert geometry.motion_params.extent == pytest.approx(np.pi + fan_angle)
 
     # --- 3d case ---
 
@@ -676,7 +676,7 @@ def test_cone_beam_geometry_helper():
 
     # Validate angles
     assert geometry.motion_partition.is_uniform
-    assert pytest.approx(geometry.motion_partition.extent, 2 * np.pi)
+    assert geometry.motion_partition.extent == pytest.approx(2 * np.pi)
     assert (geometry.motion_partition.cell_sides <=
             (r + rho) / r * np.pi / (rho * omega))
 
@@ -699,15 +699,51 @@ def test_cone_beam_geometry_helper():
 
     # Validate angles
     assert geometry.motion_partition.is_uniform
-    assert pytest.approx(geometry.motion_partition.extent, 2 * np.pi)
+    assert geometry.motion_partition.extent == pytest.approx(2 * np.pi)
     assert (geometry.motion_partition.cell_sides <=
             (r + rho) / r * np.pi / (rho * omega))
 
     # Validate detector
     det_width = 2 * magnification * rho
     R = np.hypot(r, det_width / 2)
     assert geometry.det_partition.cell_sides <= np.pi * R / (r * omega)
-    assert pytest.approx(geometry.det_partition.extent, det_width)
+    assert geometry.det_partition.extent == pytest.approx(det_width)
+
+
+def test_helical_geometry_helper():
+    """Test that helical_geometry satisfies the sampling conditions.
+
+    See the `helical_geometry` documentation for the exact conditions.
+    """
+    # Parameters
+    src_radius = 3
+    det_radius = 9
+    num_turns = 4
+
+    magnification = (src_radius + det_radius) / src_radius
+    rho = np.sqrt(2)
+    omega = np.pi * 10.0
+    r = src_radius + det_radius
+
+    # Create object
+    space = odl.uniform_discr([-1, -1, -2], [1, 1, 2], [20, 20, 40])
+    geometry = odl.tomo.helical_geometry(space, src_radius, det_radius,
+                                         num_turns=num_turns)
+
+    # Validate angles
+    assert geometry.motion_partition.is_uniform
+    assert (geometry.motion_partition.extent ==
+            pytest.approx(num_turns * 2 * np.pi))
+    assert (geometry.motion_partition.cell_sides <=
+            (r + rho) / r * np.pi / (rho * omega))
+
+    # Validate detector
+    det_width = 2 * magnification * rho
+    R = np.hypot(r, det_width / 2)
+    assert geometry.det_partition.cell_sides[0] <= np.pi * R / (r * omega)
+    mag = (3 + 9) / (3 + rho)
+    delta_h = space.cell_sides[2] * mag
+    assert geometry.det_partition.cell_sides[1] <= delta_h
 
 
 if __name__ == '__main__':