Merge pull request #825 from odlgroup/issue-824__parker_weighting

Issue 824  parker weighting

examples/tomo/filtered_backprojection_cone_2d_partial_scan.py

@@ -0,0 +1,86 @@
+# Copyright 2014-2016 The ODL development group
+#
+# This file is part of ODL.
+#
+# ODL is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# ODL is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with ODL.  If not, see <http://www.gnu.org/licenses/>.
+
+"""
+Example using a filtered back-projection (FBP) in fan beam using `fbp_op`.
+
+Note that the FBP is only approximate in this geometry, but still gives a
+decent reconstruction that can be used as an initial guess in more complicated
+methods.
+
+Here we look at a partial scan, where the angular interval is not 2 * pi.
+This caues issues for the regular FBP reconstruction, but can be improved
+via a Parker weighting.
+"""
+
+import numpy as np
+import odl
+
+
+# --- Set-up geometry of the problem --- #
+
+
+# Discrete reconstruction space: discretized functions on the cube
+# [-20, 20]^2 with 300 samples per dimension.
+reco_space = odl.uniform_discr(
+    min_pt=[-20, -20], max_pt=[20, 20], shape=[300, 300],
+    dtype='float32')
+
+# Make a circular cone beam geometry with flat detector
+# Angles: uniformly spaced, n = 360, min = 0, max = pi + fan angle
+angle_partition = odl.uniform_partition(0, np.pi + 0.7, 360)
+# Detector: uniformly sampled, n = 558, min = -40, max = 40
+detector_partition = odl.uniform_partition(-40, 40, 558)
+# Geometry with large fan angle
+geometry = odl.tomo.FanFlatGeometry(
+    angle_partition, detector_partition, src_radius=80, det_radius=40)
+
+
+# --- Create Filtered Back-Projection (FBP) operator --- #
+
+
+# Ray transform (= forward projection). We use the ASTRA CUDA backend.
+ray_trafo = odl.tomo.RayTransform(reco_space, geometry, impl='astra_cuda')
+
+# Create FBP operator using utility function
+# We select a Hann filter, and only use the lowest 80% of frequencies to avoid
+# high frequency noise.
+fbp = odl.tomo.fbp_op(ray_trafo, filter_type='Hann', frequency_scaling=0.8)
+
+# Apply parker weighting in order to improve reconstruction
+parker_weighting = odl.tomo.parker_weighting(ray_trafo)
+parker_weighting.show()
+parker_weighted_fbp = fbp * parker_weighting
+
+# --- Show some examples --- #
+
+
+# Create a discrete Shepp-Logan phantom (modified version)
+phantom = odl.phantom.shepp_logan(reco_space, modified=True)
+
+# Create projection data by calling the ray transform on the phantom
+proj_data = ray_trafo(phantom)
+
+# Calculate filtered back-projection of data
+fbp_reconstruction = fbp(proj_data)
+pw_fbp_reconstruction = parker_weighted_fbp(proj_data)
+
+# Shows a slice of the phantom, projections, and reconstruction
+phantom.show(title='Phantom')
+proj_data.show(title='Projection data (sinogram)')
+fbp_reconstruction.show(title='Filtered back-projection')
+pw_fbp_reconstruction.show(title='Parker weighted filtered back-projection')

examples/tomo/filtered_backprojection_cone_3d_partial_scan.py

@@ -0,0 +1,89 @@
+# Copyright 2014-2016 The ODL development group
+#
+# This file is part of ODL.
+#
+# ODL is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# ODL is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with ODL.  If not, see <http://www.gnu.org/licenses/>.
+
+"""
+Example using a filtered back-projection (FBP) in cone-beam 3d using `fbp_op`.
+
+Note that the FBP is only approximate in this geometry, but still gives a
+decent reconstruction that can be used as an initial guess in more complicated
+methods.
+
+Here we look at a partial scan, where the angular interval is not 2 * pi.
+This caues issues for the regular FBP reconstruction, but can be improved
+via a Parker weighting.
+
+Note that since this is a fully 3d example, it may take some time to run,
+about ~20s.
+"""
+
+import numpy as np
+import odl
+
+
+# --- Set-up geometry of the problem --- #
+
+
+# Discrete 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')
+
+# Make a circular cone beam geometry with flat detector and with a short scan
+# Angles: uniformly spaced, n = 360, min = 0, max = 1.3 * pi
+angle_partition = odl.uniform_partition(0, 1.3 * np.pi, 360)
+# Detector: uniformly sampled, n = (558, 558), min = (-60, -60), max = (60, 60)
+detector_partition = odl.uniform_partition([-60, -60], [60, 60], [558, 558])
+# Geometry with large cone and fan angle and tilted axis.
+geometry = odl.tomo.CircularConeFlatGeometry(
+    angle_partition, detector_partition, src_radius=80, det_radius=40)
+
+
+# --- Create Filteredback-projection (FBP) operator --- #
+
+
+# Ray transform (= forward projection). We use the ASTRA CUDA backend.
+ray_trafo = odl.tomo.RayTransform(reco_space, geometry, impl='astra_cuda')
+
+# Create FBP operator using utility function
+# We select a Shepp-Logan filter, and only use the lowest 80% of frequencies to
+# avoid high frequency noise.
+fbp = odl.tomo.fbp_op(ray_trafo,
+                      filter_type='Shepp-Logan', frequency_scaling=0.8)
+
+# Apply parker weighting in order to improve reconstruction
+parker_weighting = odl.tomo.parker_weighting(ray_trafo)
+parker_weighted_fbp = fbp * parker_weighting
+
+# --- Show some examples --- #
+
+
+# Create a discrete Shepp-Logan phantom (modified version)
+phantom = odl.phantom.shepp_logan(reco_space, modified=True)
+
+# Create projection data by calling the ray transform on the phantom
+proj_data = ray_trafo(phantom)
+
+# Calculate filtered back-projection of data
+fbp_reconstruction = fbp(proj_data)
+pw_fbp_reconstruction = parker_weighted_fbp(proj_data)
+
+# Shows a slice of the phantom, projections, and reconstruction
+phantom.show(title='Phantom')
+proj_data.show(title='Simulated data (sinogram)')
+fbp_reconstruction.show(title='Filtered back-projection')
+pw_fbp_reconstruction.show(title='Parker weighted filtered back-projection')

odl/tomo/analytic/filtered_back_projection.py

@@ -26,7 +26,8 @@
 from odl.trafos import FourierTransform, PYFFTW_AVAILABLE
 
 
-__all__ = ('fbp_op', 'tam_danielson_window')
+__all__ = ('fbp_op', 'fbp_filter_op', 'tam_danielson_window',
+           'parker_weighting')
 
 
 def _axis_in_detector(geometry):
@@ -134,6 +135,7 @@ def tam_danielson_window(ray_trafo, smoothing_width=0.05, n_half_rot=1):
     See Also
     --------
     fbp_op : Filtered back-projection operator from `RayTransform`
+    tam_danielson_window : Weighting for short scan data
     HelicalConeFlatGeometry : The primary use case for this window function.
 
     References
@@ -199,6 +201,97 @@ def window_fcn(x):
     return ray_trafo.range.element(window_fcn) / n_half_rot
 
 
+def parker_weighting(ray_trafo, q=0.25):
+    """Create parker weighting for a `RayTransform`.
+
+    Parker weighting is a weighting function that ensures that oversampled
+    fan/cone beam data are weighted such that each line has unit weight. It is
+    useful in analytic reconstruction methods such as FBP to give a more
+    accurate result and can improve convergence rates for iterative methods.
+
+    See the article `Parker weights revisited`_ for more information.
+
+    Parameters
+    ----------
+    ray_trafo : `RayTransform`
+        The ray transform for which to compute the weights.
+    q : float
+        Parameter controlling the speed of the roll-off at the edges of the
+        weighting. 1.0 gives the classical Parker weighting, while smaller
+        values in general lead to lower noise but stronger discretization
+        artifacts.
+
+    Returns
+    -------
+    parker_weighting : ``ray_trafo.range`` element
+
+    See Also
+    --------
+    fbp_op : Filtered back-projection operator from `RayTransform`
+    tam_danielson_window : Indicator function for helical data
+    FanFlatGeometry : Use case in 2d
+    CircularConeFlatGeometry : Use case in 3d
+
+    References
+    ----------
+    .. _Parker weights revisited: https://www.ncbi.nlm.nih.gov/pubmed/11929021
+    """
+    # Note: Parameter names taken from WES2002
+
+    # Extract parameters
+    src_radius = ray_trafo.geometry.src_radius
+    det_radius = ray_trafo.geometry.det_radius
+    ndim = ray_trafo.geometry.ndim
+    angles = ray_trafo.range.meshgrid[0]
+    min_rot_angle = ray_trafo.geometry.motion_partition.min_pt
+    alen = ray_trafo.geometry.motion_params.length
+
+    # Parker weightings are not defined for helical geometries
+    if ray_trafo.geometry.ndim != 2:
+        pitch = ray_trafo.geometry.pitch
+        if pitch != 0:
+            raise ValueError('Parker weighting window is only defined with '
+                             '`pitch==0`')
+
+    # Find distance from projection of rotation axis for each pixel
+    if ndim == 2:
+        dx = ray_trafo.range.meshgrid[1]
+    elif ndim == 3:
+        # Find projection of axis on detector
+        rot_dir = _rotation_direction_in_detector(ray_trafo.geometry)
+        dx = (rot_dir[0] * ray_trafo.range.meshgrid[1] +
+              rot_dir[1] * ray_trafo.range.meshgrid[2])
+
+    # Compute parameters
+    dx_abs_max = np.max(np.abs(dx))
+    max_fan_angle = 2 * np.arctan2(dx_abs_max, src_radius + det_radius)
+    delta = max_fan_angle / 2
+    epsilon = alen - np.pi - max_fan_angle
+
+    if epsilon < 0:
+        raise Exception('data not sufficiently sampled for parker weighting')
+
+    # Define utility functions
+    def S(betap):
+        return (0.5 * (1.0 + np.sin(np.pi * betap)) * (np.abs(betap) < 0.5) +
+                (betap >= 0.5))
+
+    def b(alpha):
+        return q * (2 * delta - 2 * alpha + epsilon)
+
+    # Create weighting function
+    beta = angles - min_rot_angle  # rotation angle
+    alpha = np.arctan2(dx, src_radius + det_radius)
+
+    S1 = S(beta / b(alpha) - 0.5)
+    S2 = S((beta - 2 * delta + 2 * alpha - epsilon) / b(alpha) + 0.5)
+    S3 = S((beta - np.pi + 2 * alpha) / b(-alpha) - 0.5)
+    S4 = S((beta - np.pi - 2 * delta - epsilon) / b(-alpha) + 0.5)
+
+    scale = 0.5 * alen / np.pi
+    return ray_trafo.range.element((S1 + S2 - S3 - S4) * scale)
+
+
 def fbp_filter_op(ray_trafo, padding=True, filter_type='Ram-Lak',
                   frequency_scaling=1.0):
     """Create a filter operator for FBP from a `RayTransform`.
@@ -396,12 +489,14 @@ def fbp_op(ray_trafo, padding=True, filter_type='Ram-Lak',
 
 
 if __name__ == '__main__':
-    # Display the various filters
     import odl
     import matplotlib.pyplot as plt
+
+    # Display the various filters
     x = np.linspace(0, 1, 100)
     cutoff = 0.7
 
+    plt.figure('fbp filter')
     for filter_name in ['Ram-Lak', 'Shepp-Logan', 'Cosine', 'Hamming', 'Hann']:
         plt.plot(x, x * _fbp_filter(x, filter_name, cutoff), label=filter_name)
 
@@ -424,6 +519,19 @@ def fbp_op(ray_trafo, padding=True, filter_type='Ram-Lak',
     td_window = tam_danielson_window(ray_trafo, smoothing_width=0)
     td_window.show('Tam-Danielson window', coords=[0, None, None])
 
+    # Show the Parker weighting
+
+    # Create Ray Transform in fan beam geometry
+    angle_partition = odl.uniform_partition(0, np.pi + 0.8, 360)
+    detector_partition = odl.uniform_partition(-40, 40, 558)
+    geometry = odl.tomo.FanFlatGeometry(
+        angle_partition, detector_partition, src_radius=80, det_radius=40)
+    ray_trafo = odl.tomo.RayTransform(reco_space, geometry, impl='astra_cuda')
+
+    # Crete and show parker weighting
+    parker_weighting = parker_weighting(ray_trafo)
+    parker_weighting.show('Parker weighting')
+
     # pylint: disable=wrong-import-position
     from odl.util.testutils import run_doctests
     run_doctests()