Merge pull request #584 from josteinbf/iradon-algebraic

SART algorithm for tomography reconstruction

CONTRIBUTORS.txt

@@ -144,3 +144,4 @@
 
 - Jostein Bø Fløystad
   Reconstruction circle mode for Radon transform
+  Simultaneous Algebraic Reconstruction Technique for inverse Radon transform

doc/examples/plot_radon_transform.py

@@ -1,57 +1,200 @@
+# -*- coding: utf-8 -*-
 """
 ===============
 Radon transform
 ===============
 
-The radon transform is a technique widely used in tomography to
-reconstruct an object from different projections. A projection is, for
-example, the scattering data obtained as the output of a tomographic
-scan.
+In computed tomography, the tomography reconstruction problem is to obtain
+a tomographic slice image from a set of projections [1]_. A projection is formed
+by drawing a set of parallel rays through the 2D object of interest, assigning
+the integral of the object's contrast along each ray to a single pixel in the
+projection. A single projection of a 2D object is one dimensional. To
+enable computed tomography reconstruction of the object, several projections
+must be acquired, each of them corresponding to a different angle between the
+rays with respect to the object. A collection of projections at several angles
+is called a sinogram, which is a linear transform of the original image.
 
-For more information see:
+The inverse Radon transform is used in computed tomography to reconstruct
+a 2D image from the measured projections (the sinogram). A practical, exact
+implementation of the inverse Radon transform does not exist, but there are
+several good approximate algorithms available.
 
-  - http://en.wikipedia.org/wiki/Radon_transform
-  - http://www.clear.rice.edu/elec431/projects96/DSP/bpanalysis.html
+As the inverse Radon transform reconstructs the object from a set of
+projections, the (forward) Radon transform can be used to simulate a
+tomography experiment.
 
-This script performs the radon transform, and reconstructs the
-input image based on the resulting sinogram.
+This script performs the Radon transform to simulate a tomography experiment
+and reconstructs the input image based on the resulting sinogram formed by
+the simulation. Two methods for performing the inverse Radon transform
+and reconstructing the original image are compared: The Filtered Back
+Projection (FBP) and the Simultaneous Algebraic Reconstruction
+Technique (SART).
 
+.. seealso::
+
+    - AC Kak, M Slaney, "Principles of Computerized Tomographic Imaging",
+      http://www.slaney.org/pct/pct-toc.html
+    - http://en.wikipedia.org/wiki/Radon_transform
+
+The forward transform
+=====================
+
+As our original image, we will use the Shepp-Logan phantom. When calculating
+the Radon transform, we need to decide how many projection angles we wish
+to use. As a rule of thumb, the number of projections should be about the
+same as the number of pixels there are across the object (to see why this
+is so, consider how many unknown pixel values must be determined in the
+reconstruction process and compare this to the number of measurements
+provided by the projections), and we follow that rule here. Below is the
+original image and its Radon transform, often known as its _sinogram_:
 """
 
+from __future__ import print_function, division
+
+import numpy as np
 import matplotlib.pyplot as plt
 
 from skimage.io import imread
 from skimage import data_dir
-from skimage.transform import radon, iradon, rescale
-
+from skimage.transform import radon, rescale
 
 image = imread(data_dir + "/phantom.png", as_grey=True)
 image = rescale(image, scale=0.4)
 
-plt.figure(figsize=(8, 8.5))
+plt.figure(figsize=(8, 4.5))
 
-plt.subplot(221)
+plt.subplot(121)
 plt.title("Original")
 plt.imshow(image, cmap=plt.cm.Greys_r)
 
-plt.subplot(222)
-projections = radon(image, theta=[0, 45, 90])
-plt.plot(projections)
-plt.title("Projections at\n0, 45 and 90 degrees")
-plt.xlabel("Projection axis")
-plt.ylabel("Intensity")
-
-projections = radon(image)
-plt.subplot(223)
+theta = np.linspace(0., 180., max(image.shape), endpoint=True)
+sinogram = radon(image, theta=theta, circle=True)
+plt.subplot(122)
 plt.title("Radon transform\n(Sinogram)")
-plt.xlabel("Projection angle (degrees)")
-plt.ylabel("Projection axis")
-plt.imshow(projections, aspect='auto')
-
-reconstruction = iradon(projections)
-plt.subplot(224)
-plt.title("Reconstruction\nfrom sinogram")
-plt.imshow(reconstruction, cmap=plt.cm.Greys_r)
+plt.xlabel("Projection angle (deg)")
+plt.ylabel("Projection position (pixels)")
+plt.imshow(sinogram, cmap=plt.cm.Greys_r,
+           extent=(0, 180, 0, sinogram.shape[0]), aspect='auto')
 
 plt.subplots_adjust(hspace=0.4, wspace=0.5)
 plt.show()
+
+"""
+.. image:: PLOT2RST.current_figure
+
+Reconstruction with the Filtered Back Projection (FBP)
+======================================================
+
+The mathematical foundation of the filtered back projection is the Fourier
+slice theorem [2]_. It uses Fourier transform of the projection and
+interpolation in Fourier space to obtain the 2D Fourier transform of the image,
+which is then inverted to form the reconstructed image. The filtered back
+projection is among the fastest methods of performing the inverse Radon
+transform. The only tunable parameter for the FBP is the filter, which is
+applied to the Fourier transformed projections. It may be used to suppress
+high frequency noise in the reconstruction. ``skimage`` provides a few
+different options for the filter.
+
+"""
+
+from skimage.transform import iradon
+
+reconstruction_fbp = iradon(sinogram, theta=theta, circle=True)
+error = reconstruction_fbp - image
+print('FBP rms reconstruction error: %.3g' % np.sqrt(np.mean(error**2)))
+
+imkwargs = dict(vmin=-0.2, vmax=0.2)
+plt.figure(figsize=(8, 4.5))
+plt.subplot(121)
+plt.title("Reconstruction\nFiltered back projection")
+plt.imshow(reconstruction_fbp, cmap=plt.cm.Greys_r)
+plt.subplot(122)
+plt.title("Reconstruction error\nFiltered back projection")
+plt.imshow(reconstruction_fbp - image, cmap=plt.cm.Greys_r, **imkwargs)
+plt.show()
+
+"""
+.. image:: PLOT2RST.current_figure
+
+Reconstruction with the Simultaneous Algebraic Reconstruction Technique
+=======================================================================
+
+Algebraic reconstruction techniques for tomography are based on a
+straightforward idea: for a pixelated image the value of a single ray in a
+particular projection is simply a sum of all the pixels the ray passes through
+on its way through the object. This is a way of expressing the forward Radon
+transform. The inverse Radon transform can then be formulated as a (large) set
+of linear equations. As each ray passes through a small fraction of the pixels
+in the image, this set of equations is sparse, allowing iterative solvers for
+sparse linear systems to tackle the system of equations. One iterative method
+has been particularly popular, namely Kaczmarz' method [3]_, which has the
+property that the solution will approach a least-squares solution of the
+equation set.
+
+The combination of the formulation of the reconstruction problem as a set
+of linear equations and an iterative solver makes algebraic techniques
+relatively flexible, hence some forms of prior knowledge can be incorporated
+with relative ease.
+
+``skimage`` provides one of the more popular variations of the algebraic
+reconstruction techniques: the Simultaneous Algebraic Reconstruction Technique
+(SART) [1]_ [4]_. It uses Kaczmarz' method [3]_ as the iterative solver. A good
+reconstruction is normally obtained in a single iteration, making the method
+computationally effective. Running one or more extra iterations will normally
+improve the reconstruction of sharp, high frequency features and reduce the
+mean squared error at the expense of increased high frequency noise (the user
+will need to decide on what number of iterations is best suited to the problem
+at hand. The implementation in ``skimage`` allows prior information of the
+form of a lower and upper threshold on the reconstructed values to be supplied
+to the reconstruction.
+
+"""
+
+from skimage.transform import iradon_sart
+
+reconstruction_sart = iradon_sart(sinogram, theta=theta)
+error = reconstruction_sart - image
+print('SART (1 iteration) rms reconstruction error: %.3g'
+      % np.sqrt(np.mean(error**2)))
+
+plt.figure(figsize=(8, 8.5))
+
+plt.subplot(221)
+plt.title("Reconstruction\nSART")
+plt.imshow(reconstruction_sart, cmap=plt.cm.Greys_r)
+plt.subplot(222)
+plt.title("Reconstruction error\nSART")
+plt.imshow(reconstruction_sart - image, cmap=plt.cm.Greys_r, **imkwargs)
+
+# Run a second iteration of SART by supplying the reconstruction
+# from the first iteration as an initial estimate
+reconstruction_sart2 = iradon_sart(sinogram, theta=theta,
+                                   image=reconstruction_sart)
+error = reconstruction_sart2 - image
+print('SART (2 iterations) rms reconstruction error: %.3g'
+      % np.sqrt(np.mean(error**2)))
+
+plt.subplot(223)
+plt.title("Reconstruction\nSART, 2 iterations")
+plt.imshow(reconstruction_sart2, cmap=plt.cm.Greys_r)
+plt.subplot(224)
+plt.title("Reconstruction error\nSART, 2 iterations")
+plt.imshow(reconstruction_sart2 - image, cmap=plt.cm.Greys_r, **imkwargs)
+plt.show()
+
+"""
+.. image:: PLOT2RST.current_figure
+
+
+.. [1]  AC Kak, M Slaney, "Principles of Computerized Tomographic Imaging",
+        IEEE Press 1988. http://www.slaney.org/pct/pct-toc.html
+.. [2]  Wikipedia, Radon transform,
+        http://en.wikipedia.org/wiki/Radon_transform#Relationship_with_the_Fourier_transform
+.. [3]  S Kaczmarz, "Angenäherte auflösung von systemen linearer
+        gleichungen", Bulletin International de l’Academie Polonaise des
+        Sciences et des Lettres 35 pp 355--357 (1937)
+.. [4]  AH Andersen, AC Kak, "Simultaneous algebraic reconstruction technique
+        (SART): a superior implementation of the ART algorithm", Ultrasonic
+        Imaging 6 pp 81--94 (1984)
+
+"""

skimage/transform/__init__.py

@@ -2,7 +2,7 @@ from ._hough_transform import (hough_circle, hough_ellipse, hough_line,
                                probabilistic_hough_line)
 from .hough_transform import (hough, probabilistic_hough, hough_peaks,
                               hough_line_peaks)
-from .radon_transform import radon, iradon
+from .radon_transform import radon, iradon, iradon_sart
 from .finite_radon_transform import frt2, ifrt2
 from .integral import integral_image, integrate
 from ._geometric import (warp, warp_coords, estimate_transform,
@@ -24,6 +24,7 @@ __all__ = ['hough_circle',
            'hough_line_peaks',
            'radon',
            'iradon',
+           'iradon_sart',
            'frt2',
            'ifrt2',
            'integral_image',

skimage/transform/_radon_transform.pyx

@@ -0,0 +1,202 @@
+#cython: cdivision=True
+#cython: boundscheck=False
+#cython: nonecheck=False
+#cython: wraparound=False
+import numpy as np
+
+cimport numpy as cnp
+cimport cython
+from libc.math cimport cos, sin, floor, ceil, sqrt, abs, M_PI
+
+
+cpdef bilinear_ray_sum(cnp.double_t[:, :] image, cnp.double_t theta,
+                       cnp.double_t ray_position):
+    """
+    Compute the projection of an image along a ray.
+
+    Parameters
+    ----------
+    image : 2D array, dtype=float
+        Image to project.
+    theta : float
+        Angle of the projection
+    ray_position : float
+        Position of the ray within the projection
+
+    Returns
+    -------
+    projected_value : float
+        Ray sum along the projection
+    norm_of_weights :
+        A measure of how long the ray's path through the reconstruction
+        circle was
+    """
+    theta = theta / 180. * M_PI
+    cdef cnp.double_t radius = image.shape[0] // 2 - 1
+    cdef cnp.double_t projection_center = image.shape[0] // 2
+    cdef cnp.double_t rotation_center = image.shape[0] // 2
+    # (s, t) is the (x, y) system rotated by theta
+    cdef cnp.double_t t = ray_position - projection_center
+    # s0 is the half-length of the ray's path in the reconstruction circle
+    cdef cnp.double_t s0
+    s0 = sqrt(radius**2 - t**2) if radius**2 >= t**2 else 0.
+    cdef Py_ssize_t Ns = 2 * (<Py_ssize_t> ceil(2 * s0))  # number of steps
+                                                          # along the ray
+    cdef cnp.double_t ray_sum = 0.
+    cdef cnp.double_t weight_norm = 0.
+    cdef cnp.double_t ds, dx, dy, x0, y0, x, y, di, dj,
+    cdef cnp.double_t index_i, index_j, weight
+    cdef Py_ssize_t k, i, j
+    if Ns > 0:
+        # step length between samples
+        ds = 2 * s0 / Ns
+        dx = -ds * cos(theta)
+        dy = -ds * sin(theta)
+        # point of entry of the ray into the reconstruction circle
+        x0 = s0 * cos(theta) - t * sin(theta)
+        y0 = s0 * sin(theta) + t * cos(theta)
+        for k in range(Ns+1):
+            x = x0 + k * dx
+            y = y0 + k * dy
+            index_i = x + rotation_center
+            index_j = y + rotation_center
+            i = <Py_ssize_t> floor(index_i)
+            j = <Py_ssize_t> floor(index_j)
+            di = index_i - floor(index_i)
+            dj = index_j - floor(index_j)
+            # Use linear interpolation between values
+            # Where values fall outside the array, assume zero
+            if i > 0 and j > 0:
+                weight = (1. - di) * (1. - dj) * ds
+                ray_sum += weight * image[i, j]
+                weight_norm += weight**2
+            if i > 0 and j < image.shape[1] - 1:
+                weight = (1. - di) * dj * ds
+                ray_sum += weight * image[i, j+1]
+                weight_norm += weight**2
+            if i < image.shape[0] - 1 and j > 0:
+                weight = di * (1 - dj) * ds
+                ray_sum += weight * image[i+1, j]
+                weight_norm += weight**2
+            if i < image.shape[0] - 1 and j < image.shape[1] - 1:
+                weight = di * dj * ds
+                ray_sum += weight * image[i+1, j+1]
+                weight_norm += weight**2
+    return ray_sum, weight_norm
+
+
+cpdef bilinear_ray_update(cnp.double_t[:, :] image,
+                          cnp.double_t[:, :] image_update,
+                          cnp.double_t theta, cnp.double_t ray_position,
+                          cnp.double_t projected_value):
+    """
+    Compute the update along a ray using bilinear interpolation.
+
+    Parameters
+    ----------
+    image : 2D array, dtype=float
+        Current reconstruction estimate
+    image_update : 2D array, dtype=float
+        Array of same shape as ``image``. Updates will be added to this array.
+    theta : float
+        Angle of the projection
+    ray_position : float
+        Position of the ray within the projection
+    projected_value : float
+        Projected value (from the sinogram)
+
+    Returns
+    -------
+    deviation :
+        Deviation before updating the image
+    """
+    cdef cnp.double_t ray_sum, weight_norm, deviation
+    ray_sum, weight_norm = bilinear_ray_sum(image, theta, ray_position)
+    if weight_norm > 0.:
+        deviation = -(ray_sum - projected_value) / weight_norm
+    else:
+        deviation = 0.
+    theta = theta / 180. * M_PI
+    cdef cnp.double_t radius = image.shape[0] // 2 - 1
+    cdef cnp.double_t projection_center = image.shape[0] // 2
+    cdef cnp.double_t rotation_center = image.shape[0] // 2
+    # (s, t) is the (x, y) system rotated by theta
+    cdef cnp.double_t t = ray_position - projection_center
+    # s0 is the half-length of the ray's path in the reconstruction circle
+    cdef cnp.double_t s0
+    s0 = sqrt(radius*radius - t*t) if radius**2 >= t**2 else 0.
+    cdef Py_ssize_t Ns = 2 * (<Py_ssize_t> ceil(2 * s0))
+    cdef cnp.double_t hamming_beta = 0.46164    # beta for equiripple Hamming window
+
+    cdef cnp.double_t ds, dx, dy, x0, y0, x, y, di, dj, index_i, index_j
+    cdef cnp.double_t hamming_window
+    cdef Py_ssize_t k, i, j
+    if Ns > 0:
+        # Step length between samples
+        ds = 2 * s0 / Ns
+        dx = -ds * cos(theta)
+        dy = -ds * sin(theta)
+        # Point of entry of the ray into the reconstruction circle
+        x0 = s0 * cos(theta) - t * sin(theta)
+        y0 = s0 * sin(theta) + t * cos(theta)
+        for k in range(Ns+1):
+            x = x0 + k * dx
+            y = y0 + k * dy
+            index_i = x + rotation_center
+            index_j = y + rotation_center
+            i = <Py_ssize_t> floor(index_i)
+            j = <Py_ssize_t> floor(index_j)
+            di = index_i - floor(index_i)
+            dj = index_j - floor(index_j)
+            hamming_window = ((1 - hamming_beta)
+                              - hamming_beta * cos(2 * M_PI * k / (Ns - 1)))
+            if i > 0 and j > 0:
+                image_update[i, j] += (deviation * (1. - di) * (1. - dj)
+                                       * ds * hamming_window)
+            if i > 0 and j < image.shape[1] - 1:
+                image_update[i, j+1] += (deviation * (1. - di) * dj
+                                         * ds * hamming_window)
+            if i < image.shape[0] - 1 and j > 0:
+                image_update[i+1, j] += (deviation * di * (1 - dj)
+                                         * ds * hamming_window)
+            if i < image.shape[0] - 1 and j < image.shape[1] - 1:
+                image_update[i+1, j+1] += (deviation * di * dj
+                                           * ds * hamming_window)
+    return deviation
+
+
+@cython.boundscheck(True)
+def sart_projection_update(cnp.double_t[:, :] image not None,
+                           cnp.double_t theta,
+                           cnp.double_t[:] projection not None,
+                           cnp.double_t projection_shift=0.):
+    """
+    Compute update to a reconstruction estimate from a single projection
+    using bilinear interpolation.
+
+    Parameters
+    ----------
+    image : 2D array, dtype=float
+        Current reconstruction estimate
+    theta : float
+        Angle of the projection
+    projection : 1D array, dtype=float
+        Projected values, taken from the sinogram
+    projection_shift : float
+        Shift the position of the projection by this many pixels before
+        using it to compute an update to the reconstruction estimate
+
+    Returns
+    -------
+    image_update : 2D array, dtype=float
+        Array of same shape as ``image`` containing updates that should be
+        added to ``image`` to improve the reconstruction estimate
+    """
+    cdef cnp.ndarray[cnp.double_t, ndim=2] image_update = np.zeros_like(image)
+    cdef cnp.double_t ray_position
+    cdef Py_ssize_t i
+    for i in range(projection.shape[0]):
+        ray_position = i + projection_shift
+        bilinear_ray_update(image, image_update, theta, ray_position,
+                            projection[i])
+    return image_update

skimage/transform/radon_transform.py

@@ -1,3 +1,4 @@
+# -*- coding: utf-8 -*-
 """
 radon.py - Radon and inverse radon transforms
 
@@ -16,8 +17,9 @@ from __future__ import division
 import numpy as np
 from scipy.fftpack import fftshift, fft, ifft
 from ._warps_cy import _warp_fast
+from ._radon_transform import sart_projection_update
 
-__all__ = ["radon", "iradon"]
+__all__ = ["radon", "iradon", "iradon_sart"]
 
 
 def radon(image, theta=None, circle=False):
@@ -254,3 +256,162 @@ def iradon(radon_image, theta=None, output_size=None,
         raise ValueError("Unknown interpolation: %s" % interpolation)
 
     return reconstructed * np.pi / (2 * len(th))
+
+
+def order_angles_golden_ratio(theta):
+    """
+    Order angles to reduce the amount of correlated information
+    in subsequent projections.
+
+    Parameters
+    ----------
+    theta : 1D array of floats
+        Projection angles in degrees. Duplicate angles are not allowed.
+
+    Returns
+    -------
+    indices : 1D array of unsigned integers
+        Indices into ``theta`` such that ``theta[indices]`` gives the
+        approximate golden ratio ordering of the projections.
+
+    Notes
+    -----
+    The method used here is that of the golden ratio introduced
+    by T. Kohler.
+
+    References:
+        -Kohler, T. "A projection access scheme for iterative
+        reconstruction based on the golden section." Nuclear Science
+        Symposium Conference Record, 2004 IEEE. Vol. 6. IEEE, 2004.
+        -Winkelmann, Stefanie, et al. "An optimal radial profile order
+        based on the Golden Ratio for time-resolved MRI."
+        Medical Imaging, IEEE Transactions on 26.1 (2007): 68-76.
+    """
+    interval = 180
+
+    def angle_distance(a, b):
+        difference = a - b
+        return min(abs(difference % interval), abs(difference % -interval))
+
+    remaining = list(np.argsort(theta))   # indices into theta
+    # yield an arbitrary angle to start things off
+    index = remaining.pop(0)
+    angle = theta[index]
+    yield index
+    # determine subsequent angles using the golden ratio method
+    angle_increment = interval * (1 - (np.sqrt(5) - 1) / 2)
+    while remaining:
+        angle = (angle + angle_increment) % interval
+        insert_point = np.searchsorted(theta[remaining], angle)
+        index_below = insert_point - 1
+        index_above = 0 if insert_point == len(remaining) else insert_point
+        distance_below = angle_distance(angle, theta[remaining[index_below]])
+        distance_above = angle_distance(angle, theta[remaining[index_above]])
+        if distance_below < distance_above:
+            yield remaining.pop(index_below)
+        else:
+            yield remaining.pop(index_above)
+
+
+def iradon_sart(radon_image, theta=None, image=None, projection_shifts=None,
+                clip=None, relaxation=0.15):
+    """
+    Inverse radon transform
+
+    Reconstruct an image from the radon transform, using a single iteration of
+    the Simultaneous Algebraic Reconstruction Technique (SART) algorithm.
+
+    Parameters
+    ----------
+    radon_image : 2D array, dtype=float
+        Image containing radon transform (sinogram). Each column of
+        the image corresponds to a projection along a different angle.
+    theta : 1D array, dtype=float, optional
+        Reconstruction angles (in degrees). Default: m angles evenly spaced
+        between 0 and 180 (if the shape of `radon_image` is (N, M)).
+    image : 2D array, dtype=float, optional
+        Image containing an initial reconstruction estimate. Shape of this
+        array should be ``(radon_image.shape[0], radon_image.shape[0])``. The
+        default is an array of zeros.
+    projection_shifts : 1D array, dtype=float
+        Shift the projections contained in ``radon_image`` (the sinogram) by
+        this many pixels before reconstructing the image. The i'th value
+        defines the shift of the i'th column of ``radon_image``.
+    clip : length-2 sequence of floats
+        Force all values in the reconstructed tomogram to lie in the range
+        ``[clip[0], clip[1]]``
+    relaxation : float
+        Relaxation parameter for the update step. A higher value can
+        improve the convergence rate, but one runs the risk of instabilities.
+        Values close to or higher than 1 are not recommended.
+
+    Returns
+    -------
+    output : ndarray
+      Reconstructed image.
+
+    Notes
+    -----
+    Algebraic Reconstruction Techniques are based on formulating the tomography
+    reconstruction problem as a set of linear equations. Along each ray,
+    the projected value is the sum of all the values of the cross section along
+    the ray. A typical feature of SART (and a few other variants of algebraic
+    techniques) is that it samples the cross section at equidistant points
+    along the ray, using linear interpolation between the pixel values of the
+    cross section. The resulting set of linear equations are then solved using
+    a slightly modified Kaczmarz method.
+
+    When using SART, a single iteration is usually sufficient to obtain a good
+    reconstruction. Further iterations will tend to enhance high-frequency
+    information, but will also often increase the noise.
+
+    References:
+        -AC Kak, M Slaney, "Principles of Computerized Tomographic
+        Imaging", IEEE Press 1988.
+        -AH Andersen, AC Kak, "Simultaneous algebraic reconstruction technique
+        (SART): a superior implementation of the ART algorithm", Ultrasonic
+        Imaging 6 pp 81--94 (1984)
+        -S Kaczmarz, "Angenäherte auflösung von systemen linearer
+        gleichungen", Bulletin International de l’Academie Polonaise des
+        Sciences et des Lettres 35 pp 355--357 (1937)
+        -Kohler, T. "A projection access scheme for iterative
+        reconstruction based on the golden section." Nuclear Science
+        Symposium Conference Record, 2004 IEEE. Vol. 6. IEEE, 2004.
+        -Kaczmarz' method, Wikipedia,
+        http://en.wikipedia.org/wiki/Kaczmarz_method
+    """
+    if radon_image.ndim != 2:
+        raise ValueError('radon_image must be two dimensional')
+    reconstructed_shape = (radon_image.shape[0], radon_image.shape[0])
+    if theta is None:
+        theta = np.linspace(0, 180, radon_image.shape[1], endpoint=False)
+    elif theta.shape != (radon_image.shape[1],):
+        raise ValueError('Shape of theta (%s) does not match the '
+                         'number of projections (%d)'
+                         % (projection_shifts.shape, radon_image.shape[1]))
+    if image is None:
+        image = np.zeros(reconstructed_shape, dtype=np.float)
+    elif image.shape != reconstructed_shape:
+        raise ValueError('Shape of image (%s) does not match first dimension '
+                         'of radon_image (%s)'
+                         % (image.shape, reconstructed_shape))
+    if projection_shifts is None:
+        projection_shifts = np.zeros((radon_image.shape[1],), dtype=np.float)
+    elif projection_shifts.shape != (radon_image.shape[1],):
+        raise ValueError('Shape of projection_shifts (%s) does not match the '
+                         'number of projections (%d)'
+                         % (projection_shifts.shape, radon_image.shape[1]))
+    if not clip is None:
+        if len(clip) != 2:
+            raise ValueError('clip must be a length-2 sequence')
+        clip = (float(clip[0]), float(clip[1]))
+    relaxation = float(relaxation)
+
+    for angle_index in order_angles_golden_ratio(theta):
+        image_update = sart_projection_update(image, theta[angle_index],
+                                              radon_image[:, angle_index],
+                                              projection_shifts[angle_index])
+        image += relaxation * image_update
+        if not clip is None:
+            image = np.clip(image, clip[0], clip[1])
+    return image

skimage/transform/setup.py

@@ -15,6 +15,7 @@ def configuration(parent_package='', top_path=None):
 
     cython(['_hough_transform.pyx'], working_path=base_path)
     cython(['_warps_cy.pyx'], working_path=base_path)
+    cython(['_radon_transform.pyx'], working_path=base_path)
 
     config.add_extension('_hough_transform', sources=['_hough_transform.c'],
                          include_dirs=[get_numpy_include_dirs()])
@@ -22,6 +23,10 @@ def configuration(parent_package='', top_path=None):
     config.add_extension('_warps_cy', sources=['_warps_cy.c'],
                          include_dirs=[get_numpy_include_dirs(), '../_shared'])
 
+    config.add_extension('_radon_transform',
+                         sources=['_radon_transform.c'],
+                         include_dirs=[get_numpy_include_dirs()])
+
     return config
 
 if __name__ == '__main__':

skimage/transform/tests/test_radon_transform.py

@@ -1,15 +1,17 @@
-from __future__ import print_function
-from __future__ import division
+from __future__ import print_function, division
 
 import numpy as np
 from numpy.testing import assert_raises
 import itertools
+import os.path
+
 from skimage.transform import radon, iradon
 from skimage.io import imread
 from skimage import data_dir
 
 
-__PHANTOM = imread(data_dir + "/phantom.png", as_grey=True)[::2, ::2]
+__PHANTOM = imread(os.path.join(data_dir, "phantom.png"),
+                   as_grey=True)[::2, ::2]
 
 
 def _get_phantom():
@@ -296,6 +298,85 @@ def test_radon_iradon_circle():
         yield check_radon_iradon_circle, interpolation, shape, output_size
 
 
+def test_order_angles_golden_ratio():
+    from skimage.transform.radon_transform import order_angles_golden_ratio
+    np.random.seed(1231)
+    lengths = [1, 4, 10, 180]
+    for l in lengths:
+        theta_ordered = np.linspace(0, 180, l, endpoint=False)
+        theta_random = np.random.uniform(0, 180, l)
+        for theta in (theta_random, theta_ordered):
+            indices = [x for x in order_angles_golden_ratio(theta)]
+            # no duplicate indices allowed
+            assert len(indices) == len(set(indices))
+
+
+def test_iradon_sart():
+    from skimage.io import imread
+    from skimage import data_dir
+    from skimage.transform import rescale, radon, iradon_sart
+
+    debug = False
+
+    shepp_logan = imread(os.path.join(data_dir, "phantom.png"), as_grey=True)
+    image = rescale(shepp_logan, scale=0.4)
+    theta_ordered = np.linspace(0., 180., image.shape[0], endpoint=False)
+    theta_missing_wedge = np.linspace(0., 150., image.shape[0], endpoint=True)
+    for theta, error_factor in ((theta_ordered, 1.),
+                                (theta_missing_wedge, 2.)):
+        sinogram = radon(image, theta, circle=True)
+        reconstructed = iradon_sart(sinogram, theta)
+
+        if debug:
+            from matplotlib import pyplot as plt
+            plt.figure()
+            plt.subplot(221)
+            plt.imshow(image, interpolation='nearest')
+            plt.subplot(222)
+            plt.imshow(sinogram, interpolation='nearest')
+            plt.subplot(223)
+            plt.imshow(reconstructed, interpolation='nearest')
+            plt.subplot(224)
+            plt.imshow(reconstructed - image, interpolation='nearest')
+            plt.show()
+
+        delta = np.mean(np.abs(reconstructed - image))
+        print('delta (1 iteration) =', delta)
+        assert delta < 0.016 * error_factor
+        reconstructed = iradon_sart(sinogram, theta, reconstructed)
+        delta = np.mean(np.abs(reconstructed - image))
+        print('delta (2 iterations) =', delta)
+        assert delta < 0.013 * error_factor
+        reconstructed = iradon_sart(sinogram, theta, clip=(0, 1))
+        delta = np.mean(np.abs(reconstructed - image))
+        print('delta (1 iteration, clip) =', delta)
+        assert delta < 0.015 * error_factor
+
+        np.random.seed(1239867)
+        shifts = np.random.uniform(-3, 3, sinogram.shape[1])
+        x = np.arange(sinogram.shape[0])
+        sinogram_shifted = np.vstack(np.interp(x + shifts[i], x,
+                                               sinogram[:, i])
+                                     for i in range(sinogram.shape[1])).T
+        reconstructed = iradon_sart(sinogram_shifted, theta,
+                                    projection_shifts=shifts)
+        if debug:
+            from matplotlib import pyplot as plt
+            plt.figure()
+            plt.subplot(221)
+            plt.imshow(image, interpolation='nearest')
+            plt.subplot(222)
+            plt.imshow(sinogram_shifted, interpolation='nearest')
+            plt.subplot(223)
+            plt.imshow(reconstructed, interpolation='nearest')
+            plt.subplot(224)
+            plt.imshow(reconstructed - image, interpolation='nearest')
+            plt.show()
+
+        delta = np.mean(np.abs(reconstructed - image))
+        print('delta (1 iteration, shifted sinogram) =', delta)
+        assert delta < 0.018 * error_factor
+
 if __name__ == "__main__":
     from numpy.testing import run_module_suite
     run_module_suite()