ENH: add vector geometries, closes odlgroup#1023

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.CircularConeFlatGeometry(
+    angle_partition, detector_partition, src_radius=1000, det_radius=100,
+    axis=[1, 0, 0])
+
+circle_vecs = odl.tomo.astra_conebeam_3d_geom_to_vec(circle_geometry)
+
+# Cover the whole volume vertically, somewhat undersampled though
+vert_shift_min = -22
+vert_shift_max = 22
+num_shifts = 180
+vert_shifts = np.linspace(vert_shift_min, vert_shift_max, num=num_shifts)
+inital_vecs = circle_vecs[0]
+
+# Start from the initial position of the circle vectors and add the vertical
+# shifts to the columns 2 and 5 (source and detector z positions)
+line_vecs = np.repeat(circle_vecs[0][None, :], num_shifts, axis=0)
+line_vecs[:, 2] += vert_shifts
+line_vecs[:, 5] += vert_shifts
+
+# Build the composed geometry and the corresponding ray transform
+# (= forward projection)
+composed_vecs = np.vstack([circle_vecs, line_vecs])
+composed_geom = odl.tomo.ConeVecGeometry(detector_partition.shape,
+                                         composed_vecs)
+
+ray_trafo = odl.tomo.RayTransform(reco_space, composed_geom)
+
+# Create a Shepp-Logan phantom (modified version) and projection data
+phantom = odl.phantom.shepp_logan(reco_space, True)
+proj_data = ray_trafo(phantom)
+
+# Back-projection can be done by simply calling the adjoint operator on the
+# projection data (or any element in the projection space).
+backproj = ray_trafo.adjoint(proj_data)
+
+# Show the slice z=0 of phantom and backprojection, as well as a projection
+# image at theta=0 and a sinogram at v=0 (middle detector row)
+phantom.show(coords=[None, None, 0], title='Phantom, middle z slice')
+backproj.show(coords=[None, None, 0], title='Back-projection, middle z slice')
+proj_data.show(indices=[0, slice(None), slice(None)],
+               title='Projection 0 (circle start)')
+proj_data.show(indices=[45, slice(None), slice(None)],
+               title='Projection 45 (circle 1/4)')
+proj_data.show(indices=[90, slice(None), slice(None)],
+               title='Projection 90 (circle 1/2)')
+proj_data.show(indices=[135, slice(None), slice(None)],
+               title='Projection 135 (circle 3/4)')
+proj_data.show(indices=[179, slice(None), slice(None)],
+               title='Projection 179 (circle end)')
+proj_data.show(indices=[180, slice(None), slice(None)],
+               title='Projection 180 (line start)')
+proj_data.show(indices=[270, slice(None), slice(None)],
+               title='Projection 270 (line middle)')
+proj_data.show(indices=[359, slice(None), slice(None)],
+               title='Projection 359 (line end)')
+proj_data.show(coords=[None, None, 0], title='Sinogram, middle slice')

odl/tomo/backends/astra_cuda.py

@@ -27,8 +27,9 @@
     astra_projection_geometry, astra_volume_geometry, astra_projector,
     astra_data, astra_algorithm)
 from odl.tomo.geometry import (
-    Geometry, Parallel2dGeometry, FanFlatGeometry, Parallel3dAxisGeometry,
-    HelicalConeFlatGeometry)
+    Geometry, VecGeometry, ParallelVecGeometry, ConeVecGeometry,
+    Parallel2dGeometry, Parallel3dAxisGeometry,
+    FanFlatGeometry, HelicalConeFlatGeometry)
 
 
 __all__ = ('ASTRA_CUDA_AVAILABLE',
@@ -370,18 +371,21 @@ def astra_cuda_bp_scaling_factor(proj_space, reco_space, geometry):
     track of it and adapt the scaling accordingly.
     """
     # Angular integration weighting factor
-    # angle interval weight by approximate cell volume
-    angle_extent = float(geometry.motion_partition.extent)
-    num_angles = float(geometry.motion_partition.size)
-    scaling_factor = angle_extent / num_angles
+    if isinstance(geometry, VecGeometry):
+        angle_weighting = 1.0
+    else:
+        # Average angle step, approximation in the case of non-uniform angles
+        angle_extent = float(geometry.motion_partition.extent)
+        num_angles = float(geometry.motion_partition.size)
+        angle_weighting = angle_extent / num_angles
 
     # Correct in case of non-weighted spaces
     proj_extent = float(proj_space.partition.extent.prod())
     proj_size = float(proj_space.partition.size)
     proj_weighting = proj_extent / proj_size
 
-    scaling_factor *= (proj_space.weighting.const /
-                       proj_weighting)
+    scaling_factor = angle_weighting * (proj_space.weighting.const /
+                                        proj_weighting)
     scaling_factor /= (reco_space.weighting.const /
                        reco_space.cell_volume)
 
@@ -410,6 +414,26 @@ def astra_cuda_bp_scaling_factor(proj_space, reco_space, geometry):
             src_radius = geometry.src_radius
             det_radius = geometry.det_radius
             scaling_factor *= ((src_radius + det_radius) / src_radius) ** 2
+        elif isinstance(geometry, VecGeometry):
+            scaling_factor = (geometry.det_partition.cell_volume /
+                              reco_space.cell_volume)
+        elif isinstance(geometry, ParallelVecGeometry):
+            if geometry.ndim == 2:
+                # Scales with 1 / cell_volume
+                scaling_factor *= float(reco_space.cell_volume)
+            elif geometry.ndim == 3:
+                # Scales with cell volume
+                scaling_factor /= reco_space.cell_volume
+        elif isinstance(geometry, ConeVecGeometry):
+            # TODO: this should be based on the distances per angle, would
+            # probably be part of the weighting then
+            if geometry.ndim == 2:
+                # Scales with 1 / cell_volume
+                scaling_factor *= float(reco_space.cell_volume)
+            elif geometry.ndim == 3:
+                det_px_area = geometry.det_partition.cell_volume
+                scaling_factor *= (det_px_area ** 2 /
+                                   reco_space.cell_volume ** 2)
 
     else:
         if isinstance(geometry, Parallel2dGeometry):
@@ -424,7 +448,6 @@ def astra_cuda_bp_scaling_factor(proj_space, reco_space, geometry):
             scaling_factor *= ((src_radius + det_radius) / src_radius)
         elif isinstance(geometry, Parallel3dAxisGeometry):
             # Scales with cell volume
-            # currently only square voxels are supported
             scaling_factor /= reco_space.cell_volume
         elif isinstance(geometry, HelicalConeFlatGeometry):
             # Scales with cell volume
@@ -441,6 +464,26 @@ def astra_cuda_bp_scaling_factor(proj_space, reco_space, geometry):
             det_px_area = geometry.det_partition.cell_volume
             scaling_factor *= (src_radius ** 2 * det_px_area ** 2 /
                                reco_space.cell_volume ** 2)
+        elif isinstance(geometry, VecGeometry):
+            scaling_factor = (geometry.det_partition.cell_volume /
+                              reco_space.cell_volume)
+        elif isinstance(geometry, ParallelVecGeometry):
+            if geometry.ndim == 2:
+                # Scales with 1 / cell_volume
+                scaling_factor *= float(reco_space.cell_volume)
+            elif geometry.ndim == 3:
+                # Scales with cell volume
+                scaling_factor /= reco_space.cell_volume
+        elif isinstance(geometry, ConeVecGeometry):
+            # TODO: this should be based on the distances per angle, would
+            # probably be part of the weighting then
+            if geometry.ndim == 2:
+                # Scales with 1 / cell_volume
+                scaling_factor *= float(reco_space.cell_volume)
+            elif geometry.ndim == 3:
+                det_px_area = geometry.det_partition.cell_volume
+                scaling_factor *= (det_px_area ** 2 /
+                                   reco_space.cell_volume ** 2)
 
         # TODO: add case with new ASTRA release