From 9484afeed17531dbfb747776d05df6c662e5d366 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Jostein=20B=C3=B8=20Fl=C3=B8ystad?=
 <jostein.floystad@gmail.com>
Date: Sun, 2 Jun 2013 21:05:53 +0200
Subject: [PATCH 01/34] Add SART tomography reconstruction to radon_transform.

---
 skimage/transform/_radon_transform.pyx | 170 +++++++++++++++++++++++++
 skimage/transform/radon_transform.py   |  91 ++++++++++++-
 skimage/transform/setup.py             |   5 +
 3 files changed, 265 insertions(+), 1 deletion(-)
 create mode 100644 skimage/transform/_radon_transform.pyx

diff --git a/skimage/transform/_radon_transform.pyx b/skimage/transform/_radon_transform.pyx
new file mode 100644
index 00000000..b1bcc300
--- /dev/null
+++ b/skimage/transform/_radon_transform.pyx
@@ -0,0 +1,170 @@
+#cython: cdivision=True
+#cython: boundscheck=True
+#cython: nonecheck=True
+#cython: wraparound=False
+import numpy as np
+from numpy import pi
+
+cimport numpy as cnp
+cimport cython
+from libc.math cimport cos, sin, floor, ceil, sqrt, abs
+
+
+cpdef bilinear_ray_sum(cnp.ndarray[cnp.double_t, ndim=2] image, double theta,
+                       double ray_position):
+    '''Compute the projection of an image along a ray.
+
+    Parameters
+    ----------
+    image : 2D array, dtype=float
+        Image to project.
+    :param theta: Angle of the projection.
+    :param ray_position: 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. * pi
+    cdef double radius = image.shape[0] // 2 - 1
+    cdef double projection_center = image.shape[0] // 2 - 1
+    cdef double rotation_center = image.shape[0] // 2
+    # (s, t) is the (x, y) system rotated by theta
+    cdef double t = ray_position - projection_center
+    # s0 is the half-length of the ray's path in the reconstruction circle
+    cdef double s0
+    s0 = sqrt(radius**2 - t**2) if radius**2 >= t**2 else 0.
+    cdef Py_ssize_t Ns = 2 * int(ceil(2 * s0))  # number of steps along the ray
+    cdef double ray_sum = 0.
+    cdef double weight_norm = 0.
+    cdef double ds, dx, dy, x0, y0, x, y, di, dj, index_i, index_j
+    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:
+                ray_sum += (1. - di) * (1. - dj) * image[i, j] * ds
+                weight_norm += ((1 - di) * (1 - dj) * ds)**2
+            if i > 0 and j < image.shape[1] - 1:
+                ray_sum += (1. - di) * dj * image[i, j+1] * ds
+                weight_norm += ((1 - di) * dj * ds)**2
+            if i < image.shape[0] - 1 and j > 0:
+                ray_sum += di * (1 - dj) * image[i+1, j] * ds
+                weight_norm += (di * (1 - dj) * ds)**2
+            if i < image.shape[0] - 1 and j < image.shape[1] - 1:
+                ray_sum += di * dj * image[i+1, j+1] * ds
+                weight_norm += (di * dj * ds)**2
+    return ray_sum, weight_norm
+
+
+cpdef bilinear_ray_update(cnp.ndarray[cnp.double_t, ndim=2] image,
+        cnp.ndarray[cnp.double_t, ndim=2] image_update,
+        double theta, double ray_position, double projected_value):
+    """Compute the update along a ray using bilinear interpolation.
+
+    Parameters
+    ----------
+    image :
+        Current reconstruction estimate
+    image_update :
+        Array of same shape as ``image``. Updates will be added to this array.
+    theta :
+        Angle of the projection
+    ray_position :
+        Position of the ray within the projection
+    projected_value :
+        Projected value (from the sinogram)
+
+    Returns
+    -------
+    deviation :
+        Deviation before updating the image
+    """
+    cdef double 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. * pi
+    cdef double radius = image.shape[0] // 2 - 1
+    cdef double projection_center = image.shape[0] // 2 - 1
+    cdef double rotation_center = image.shape[0] // 2
+    # (s, t) is the (x, y) system rotated by theta
+    cdef double t = ray_position - projection_center
+    # s0 is the half-length of the ray's path in the reconstruction circle
+    cdef double s0
+    s0 = sqrt(radius*radius - t*t) if radius**2 >= t**2 else 0.
+    cdef unsigned int Ns = 2 * int(ceil(2 * s0))
+    cdef double hamming_beta = 0.46164
+
+    cdef double ds, dx, dy, x0, y0, x, y, di, dj, index_i, index_j
+    cdef double hamming_window
+    cdef unsigned int 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*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
+
+
+def sart_projection_update(cnp.ndarray[cnp.double_t, ndim=2] image, \
+                           double theta, \
+                           cnp.ndarray[cnp.double_t, ndim=1] projection):
+    cdef cnp.ndarray[cnp.double_t, ndim=2] image_update = np.zeros_like(image)
+    cdef unsigned int ray_position
+    cdef Py_ssize_t i
+    for i in range(projection.shape[0]):
+        # TODO:
+        # ip may differ from i in the future (for alignment of projections)
+        ray_position = i
+        bilinear_ray_update(image, image_update, theta, ray_position,
+                            projection[i])
+    return image_update
diff --git a/skimage/transform/radon_transform.py b/skimage/transform/radon_transform.py
index b114adad..55bd1d59 100644
--- a/skimage/transform/radon_transform.py
+++ b/skimage/transform/radon_transform.py
@@ -16,8 +16,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 +255,91 @@ def iradon(radon_image, theta=None, output_size=None,
         raise ValueError("Unknown interpolation: %s" % interpolation)
 
     return reconstructed * np.pi / (2 * len(th))
+
+
+def _sart_next_angle(remaining, used):
+    used = np.array(used)
+    used.shape = (-1, 1)
+    remaining = np.array(remaining)
+    remaining.shape = (1, -1)
+    time = np.arange(used.shape[0]) + 1
+    time.shape = (-1, 1)
+    tau = 3.
+    difference = used - remaining
+    distance = np.minimum(abs(difference % 180), abs(difference % -180))
+    #print distance
+    cost = np.exp(-distance * time / tau).sum(axis=0).squeeze()
+    next_angle_index = np.argmin(cost)
+    return remaining[0, next_angle_index]
+
+
+def iradon_sart(radon_image, theta=None, image=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 : array_like, dtype=float
+        Image containing radon transform (sinogram). Each column of
+        the image corresponds to a projection along a different angle.
+    theta : array_like, 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 : array_like, 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.
+
+    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:
+        -A. C. Kak, Malcolm 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)
+    """
+    if theta is None:
+        theta = np.linspace(0, 180, radon_image.shape[1], endpoint=False)
+    angle_indices = {theta[i]: i for i in range(theta.shape[0])}
+    reconstructed_shape = (radon_image.shape[0], radon_image.shape[0])
+    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))
+    used_angles = []
+    while angle_indices:
+        angle_index = angle_indices.pop(_sart_next_angle(angle_indices.keys(),
+                                                         used_angles))
+        image_update = sart_projection_update(image, theta[angle_index],
+                                                  radon_image[:, angle_index])
+        image += relaxation * image_update
+        if not clip is None:
+            image = clip(image, clip[0], clip[1])
+        used_angles.append(theta[angle_index])
+
+    return image
diff --git a/skimage/transform/setup.py b/skimage/transform/setup.py
index b0093d87..22f31696 100644
--- a/skimage/transform/setup.py
+++ b/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__':

From d4b33059cb15af905252da2edb01db5bf7177e1c Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Jostein=20B=C3=B8=20Fl=C3=B8ystad?=
 <jostein.floystad@gmail.com>
Date: Sun, 2 Jun 2013 21:07:58 +0200
Subject: [PATCH 02/34] Make iradon_sart visible in the transform package.

---
 skimage/transform/__init__.py | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/skimage/transform/__init__.py b/skimage/transform/__init__.py
index f079ada4..2cbc7007 100644
--- a/skimage/transform/__init__.py
+++ b/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',

From 8e6468eef5df2bf3ef79f030ceb2e78a5946cc51 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Jostein=20B=C3=B8=20Fl=C3=B8ystad?=
 <jostein.floystad@gmail.com>
Date: Sun, 9 Jun 2013 15:01:48 +0200
Subject: [PATCH 03/34] Add tests for transform.iradon_sart.

---
 .../transform/tests/test_radon_transform.py   | 35 +++++++++++++++++++
 1 file changed, 35 insertions(+)

diff --git a/skimage/transform/tests/test_radon_transform.py b/skimage/transform/tests/test_radon_transform.py
index f7127aa5..fd0d47a3 100644
--- a/skimage/transform/tests/test_radon_transform.py
+++ b/skimage/transform/tests/test_radon_transform.py
@@ -296,6 +296,41 @@ def test_radon_iradon_circle():
         yield check_radon_iradon_circle, interpolation, shape, output_size
 
 
+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(data_dir + "/phantom.png", as_grey=True)
+    image = rescale(shepp_logan, scale=0.4)
+    theta = np.linspace(0., 180., image.shape[0], endpoint=False)
+    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.025
+    reconstructed = iradon_sart(sinogram, theta, reconstructed)
+    delta = np.mean(np.abs(reconstructed - image))
+    print('delta (2 iterations) =', delta)
+    assert delta < 0.015
+
+
 if __name__ == "__main__":
     from numpy.testing import run_module_suite
     run_module_suite()

From 1e1dd180d2252a448d1771d0f06dca8a1beea228 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Jostein=20B=C3=B8=20Fl=C3=B8ystad?=
 <jostein.floystad@gmail.com>
Date: Sun, 9 Jun 2013 15:32:25 +0200
Subject: [PATCH 04/34] Implement projection shifts.

---
 skimage/transform/_radon_transform.pyx |  9 ++++-----
 skimage/transform/radon_transform.py   | 11 +++++++++--
 2 files changed, 13 insertions(+), 7 deletions(-)

diff --git a/skimage/transform/_radon_transform.pyx b/skimage/transform/_radon_transform.pyx
index b1bcc300..e198c4a5 100644
--- a/skimage/transform/_radon_transform.pyx
+++ b/skimage/transform/_radon_transform.pyx
@@ -157,14 +157,13 @@ cpdef bilinear_ray_update(cnp.ndarray[cnp.double_t, ndim=2] image,
 
 def sart_projection_update(cnp.ndarray[cnp.double_t, ndim=2] image, \
                            double theta, \
-                           cnp.ndarray[cnp.double_t, ndim=1] projection):
+                           cnp.ndarray[cnp.double_t, ndim=1] projection,
+                           double projection_shift=0.):
     cdef cnp.ndarray[cnp.double_t, ndim=2] image_update = np.zeros_like(image)
-    cdef unsigned int ray_position
+    cdef double ray_position
     cdef Py_ssize_t i
     for i in range(projection.shape[0]):
-        # TODO:
-        # ip may differ from i in the future (for alignment of projections)
-        ray_position = i
+        ray_position = i + projection_shift
         bilinear_ray_update(image, image_update, theta, ray_position,
                             projection[i])
     return image_update
diff --git a/skimage/transform/radon_transform.py b/skimage/transform/radon_transform.py
index 55bd1d59..cc1b76c9 100644
--- a/skimage/transform/radon_transform.py
+++ b/skimage/transform/radon_transform.py
@@ -273,7 +273,7 @@ def _sart_next_angle(remaining, used):
     return remaining[0, next_angle_index]
 
 
-def iradon_sart(radon_image, theta=None, image=None,
+def iradon_sart(radon_image, theta=None, image=None, projection_shifts=None,
                 clip=None, relaxation=0.15):
     """
     Inverse radon transform
@@ -293,6 +293,10 @@ def iradon_sart(radon_image, theta=None, image=None,
         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``.
 
     Returns
     -------
@@ -327,6 +331,8 @@ def iradon_sart(radon_image, theta=None, image=None,
     reconstructed_shape = (radon_image.shape[0], radon_image.shape[0])
     if image is None:
         image = np.zeros(reconstructed_shape, dtype=np.float)
+    if projection_shifts is None:
+        projection_shifts = np.zeros((radon_image.shape[1],), dtype=np.float)
     elif image.shape != reconstructed_shape:
         raise ValueError('Shape of image (%s) does not match first dimension '
                          'of radon_image (%s)'
@@ -336,7 +342,8 @@ def iradon_sart(radon_image, theta=None, image=None,
         angle_index = angle_indices.pop(_sart_next_angle(angle_indices.keys(),
                                                          used_angles))
         image_update = sart_projection_update(image, theta[angle_index],
-                                                  radon_image[:, angle_index])
+                                              radon_image[:, angle_index],
+                                              projection_shifts[angle_index])
         image += relaxation * image_update
         if not clip is None:
             image = clip(image, clip[0], clip[1])

From d5f323eec037194111a0f2ef9788347a1093e9ae Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Jostein=20B=C3=B8=20Fl=C3=B8ystad?=
 <jostein.floystad@gmail.com>
Date: Sun, 9 Jun 2013 17:06:39 +0200
Subject: [PATCH 05/34] Tests for the projection_shifts functionality of
 iradon_sart.

---
 .../transform/tests/test_radon_transform.py   | 23 +++++++++++++++++++
 1 file changed, 23 insertions(+)

diff --git a/skimage/transform/tests/test_radon_transform.py b/skimage/transform/tests/test_radon_transform.py
index fd0d47a3..650d6961 100644
--- a/skimage/transform/tests/test_radon_transform.py
+++ b/skimage/transform/tests/test_radon_transform.py
@@ -330,6 +330,29 @@ def test_iradon_sart():
     print('delta (2 iterations) =', delta)
     assert delta < 0.015
 
+    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.025
 
 if __name__ == "__main__":
     from numpy.testing import run_module_suite

From b474804654cae9652931eac0b5f94fab4828936c Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Jostein=20B=C3=B8=20Fl=C3=B8ystad?=
 <jostein.floystad@gmail.com>
Date: Sun, 9 Jun 2013 15:33:08 +0200
Subject: [PATCH 06/34] Fix docstrings for transform.iradon_sart with
 subroutines.

---
 skimage/transform/_radon_transform.pyx | 46 ++++++++++++++++++++------
 skimage/transform/radon_transform.py   | 13 ++++++--
 2 files changed, 46 insertions(+), 13 deletions(-)

diff --git a/skimage/transform/_radon_transform.pyx b/skimage/transform/_radon_transform.pyx
index e198c4a5..65cf1023 100644
--- a/skimage/transform/_radon_transform.pyx
+++ b/skimage/transform/_radon_transform.pyx
@@ -12,14 +12,17 @@ from libc.math cimport cos, sin, floor, ceil, sqrt, abs
 
 cpdef bilinear_ray_sum(cnp.ndarray[cnp.double_t, ndim=2] image, double theta,
                        double ray_position):
-    '''Compute the projection of an image along a ray.
+    """
+    Compute the projection of an image along a ray.
 
     Parameters
     ----------
     image : 2D array, dtype=float
         Image to project.
-    :param theta: Angle of the projection.
-    :param ray_position: Position of the ray within the projection
+    theta : float
+        Angle of the projection
+    ray_position : float
+        Position of the ray within the projection
 
     Returns
     -------
@@ -28,7 +31,7 @@ cpdef bilinear_ray_sum(cnp.ndarray[cnp.double_t, ndim=2] image, double theta,
     norm_of_weights :
         A measure of how long the ray's path through the reconstruction
         circle was
-    '''
+    """
     theta = theta / 180. * pi
     cdef double radius = image.shape[0] // 2 - 1
     cdef double projection_center = image.shape[0] // 2 - 1
@@ -80,19 +83,20 @@ cpdef bilinear_ray_sum(cnp.ndarray[cnp.double_t, ndim=2] image, double theta,
 cpdef bilinear_ray_update(cnp.ndarray[cnp.double_t, ndim=2] image,
         cnp.ndarray[cnp.double_t, ndim=2] image_update,
         double theta, double ray_position, double projected_value):
-    """Compute the update along a ray using bilinear interpolation.
+    """
+    Compute the update along a ray using bilinear interpolation.
 
     Parameters
     ----------
-    image :
+    image : 2D array, dtype=float
         Current reconstruction estimate
-    image_update :
+    image_update : 2D array, dtype=float
         Array of same shape as ``image``. Updates will be added to this array.
-    theta :
+    theta : float
         Angle of the projection
-    ray_position :
+    ray_position : float
         Position of the ray within the projection
-    projected_value :
+    projected_value : float
         Projected value (from the sinogram)
 
     Returns
@@ -159,6 +163,28 @@ def sart_projection_update(cnp.ndarray[cnp.double_t, ndim=2] image, \
                            double theta, \
                            cnp.ndarray[cnp.double_t, ndim=1] projection,
                            double 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 double ray_position
     cdef Py_ssize_t i
diff --git a/skimage/transform/radon_transform.py b/skimage/transform/radon_transform.py
index cc1b76c9..042d37e7 100644
--- a/skimage/transform/radon_transform.py
+++ b/skimage/transform/radon_transform.py
@@ -283,13 +283,13 @@ def iradon_sart(radon_image, theta=None, image=None, projection_shifts=None,
 
     Parameters
     ----------
-    radon_image : array_like, dtype=float
+    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 : array_like, dtype=float, optional
+    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 : array_like, dtype=float, optional
+    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.
@@ -297,6 +297,13 @@ def iradon_sart(radon_image, theta=None, image=None, projection_shifts=None,
         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
     -------

From 372f0127f97bc75abb9786e012b9f4495522cb80 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Jostein=20B=C3=B8=20Fl=C3=B8ystad?=
 <jostein.floystad@gmail.com>
Date: Sun, 9 Jun 2013 18:21:16 +0200
Subject: [PATCH 07/34] transform.iradon_sart: Clean up code for ordering
 projections.

---
 skimage/transform/radon_transform.py | 46 +++++++++++++++-------------
 1 file changed, 25 insertions(+), 21 deletions(-)

diff --git a/skimage/transform/radon_transform.py b/skimage/transform/radon_transform.py
index 042d37e7..445cd24c 100644
--- a/skimage/transform/radon_transform.py
+++ b/skimage/transform/radon_transform.py
@@ -257,20 +257,30 @@ def iradon(radon_image, theta=None, output_size=None,
     return reconstructed * np.pi / (2 * len(th))
 
 
-def _sart_next_angle(remaining, used):
-    used = np.array(used)
-    used.shape = (-1, 1)
-    remaining = np.array(remaining)
-    remaining.shape = (1, -1)
-    time = np.arange(used.shape[0]) + 1
-    time.shape = (-1, 1)
-    tau = 3.
-    difference = used - remaining
-    distance = np.minimum(abs(difference % 180), abs(difference % -180))
-    #print distance
-    cost = np.exp(-distance * time / tau).sum(axis=0).squeeze()
-    next_angle_index = np.argmin(cost)
-    return remaining[0, next_angle_index]
+def _sart_order_angles(theta, tau=3.):
+    """
+    Order angles to reduce the amount of correlated information
+    in subsequent projections, i.e. make sure subsequent angles
+    are as far away from each other mod 180 degrees as possible.
+    Indices into the ``theta`` array are yielded.
+    """
+    used_indices = [0]
+    remaining_indices = range(1, len(theta))
+    yield 0
+    while remaining_indices:
+        used = np.array(theta[used_indices])
+        used.shape = (-1, 1)
+        remaining = np.array(theta[remaining_indices])
+        remaining.shape = (1, -1)
+        time = (np.arange(used.shape[0]) + 1)[::-1]
+        time.shape = (-1, 1)
+        difference = used - remaining
+        distance = np.minimum(abs(difference % 180), abs(difference % -180))
+        cost = np.exp(-distance * time / tau).sum(axis=0).squeeze()
+        next_angle_remaining_index = np.argmin(cost)
+        next_angle_index = remaining_indices.pop(next_angle_remaining_index)
+        used_indices.append(next_angle_index)
+        yield next_angle_index
 
 
 def iradon_sart(radon_image, theta=None, image=None, projection_shifts=None,
@@ -334,7 +344,6 @@ def iradon_sart(radon_image, theta=None, image=None, projection_shifts=None,
     """
     if theta is None:
         theta = np.linspace(0, 180, radon_image.shape[1], endpoint=False)
-    angle_indices = {theta[i]: i for i in range(theta.shape[0])}
     reconstructed_shape = (radon_image.shape[0], radon_image.shape[0])
     if image is None:
         image = np.zeros(reconstructed_shape, dtype=np.float)
@@ -344,16 +353,11 @@ def iradon_sart(radon_image, theta=None, image=None, projection_shifts=None,
         raise ValueError('Shape of image (%s) does not match first dimension '
                          'of radon_image (%s)'
                          % (image.shape, reconstructed_shape))
-    used_angles = []
-    while angle_indices:
-        angle_index = angle_indices.pop(_sart_next_angle(angle_indices.keys(),
-                                                         used_angles))
+    for angle_index in _sart_order_angles(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 = clip(image, clip[0], clip[1])
-        used_angles.append(theta[angle_index])
-
     return image

From 496902145f88169f53ae40bec3f5ef5ffd6fe260 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Jostein=20B=C3=B8=20Fl=C3=B8ystad?=
 <jostein.floystad@gmail.com>
Date: Mon, 10 Jun 2013 09:38:28 +0200
Subject: [PATCH 08/34] iradon_sart: Add wikipedia reference to Kaczmarz
 method.

---
 skimage/transform/radon_transform.py | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/skimage/transform/radon_transform.py b/skimage/transform/radon_transform.py
index 445cd24c..ebf06d35 100644
--- a/skimage/transform/radon_transform.py
+++ b/skimage/transform/radon_transform.py
@@ -341,6 +341,8 @@ def iradon_sart(radon_image, theta=None, image=None, projection_shifts=None,
         -AH Andersen, AC Kak, "Simultaneous algebraic reconstruction technique
         (SART): a superior implementation of the ART algorithm", Ultrasonic
         Imaging 6 pp 81--94 (1984)
+        -Kaczmarz' method, Wikipedia,
+        http://en.wikipedia.org/wiki/Kaczmarz_method
     """
     if theta is None:
         theta = np.linspace(0, 180, radon_image.shape[1], endpoint=False)

From 94eba1c9116eeb7cdf87447c98a8005141b4cd4d Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Jostein=20B=C3=B8=20Fl=C3=B8ystad?=
 <jostein.floystad@gmail.com>
Date: Mon, 10 Jun 2013 10:56:12 +0200
Subject: [PATCH 09/34] iradon_sart: Clarify how constants are chosen.

---
 skimage/transform/_radon_transform.pyx | 2 +-
 skimage/transform/radon_transform.py   | 3 ++-
 2 files changed, 3 insertions(+), 2 deletions(-)

diff --git a/skimage/transform/_radon_transform.pyx b/skimage/transform/_radon_transform.pyx
index 65cf1023..321d9556 100644
--- a/skimage/transform/_radon_transform.pyx
+++ b/skimage/transform/_radon_transform.pyx
@@ -120,7 +120,7 @@ cpdef bilinear_ray_update(cnp.ndarray[cnp.double_t, ndim=2] image,
     cdef double s0
     s0 = sqrt(radius*radius - t*t) if radius**2 >= t**2 else 0.
     cdef unsigned int Ns = 2 * int(ceil(2 * s0))
-    cdef double hamming_beta = 0.46164
+    cdef double hamming_beta = 0.46164    # beta for equiripple Hamming window
 
     cdef double ds, dx, dy, x0, y0, x, y, di, dj, index_i, index_j
     cdef double hamming_window
diff --git a/skimage/transform/radon_transform.py b/skimage/transform/radon_transform.py
index ebf06d35..bba73c3f 100644
--- a/skimage/transform/radon_transform.py
+++ b/skimage/transform/radon_transform.py
@@ -257,13 +257,14 @@ def iradon(radon_image, theta=None, output_size=None,
     return reconstructed * np.pi / (2 * len(th))
 
 
-def _sart_order_angles(theta, tau=3.):
+def _sart_order_angles(theta):
     """
     Order angles to reduce the amount of correlated information
     in subsequent projections, i.e. make sure subsequent angles
     are as far away from each other mod 180 degrees as possible.
     Indices into the ``theta`` array are yielded.
     """
+    tau = 3.    # time constant for correlations; 0.1 < tau < 100 works well
     used_indices = [0]
     remaining_indices = range(1, len(theta))
     yield 0

From 1b620a0a1294cbffb5529f28e5974b35804ac5e3 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Jostein=20B=C3=B8=20Fl=C3=B8ystad?=
 <jostein.floystad@gmail.com>
Date: Mon, 10 Jun 2013 11:00:19 +0200
Subject: [PATCH 10/34] Tests for iradon_sart: Robust path handling.

---
 skimage/transform/tests/test_radon_transform.py | 10 ++++++----
 1 file changed, 6 insertions(+), 4 deletions(-)

diff --git a/skimage/transform/tests/test_radon_transform.py b/skimage/transform/tests/test_radon_transform.py
index 650d6961..5ad87877 100644
--- a/skimage/transform/tests/test_radon_transform.py
+++ b/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():
@@ -303,7 +305,7 @@ def test_iradon_sart():
 
     debug = False
 
-    shepp_logan = imread(data_dir + "/phantom.png", as_grey=True)
+    shepp_logan = imread(os.path.join(data_dir, "phantom.png"), as_grey=True)
     image = rescale(shepp_logan, scale=0.4)
     theta = np.linspace(0., 180., image.shape[0], endpoint=False)
     sinogram = radon(image, theta, circle=True)

From 59bdb24c92cc7b12d6c04087f51300c1c003b49e Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Jostein=20B=C3=B8=20Fl=C3=B8ystad?=
 <jostein.floystad@gmail.com>
Date: Mon, 10 Jun 2013 11:25:23 +0200
Subject: [PATCH 11/34] iradon_sart: Clean up cython code to minimize python
 calls.

---
 skimage/transform/_radon_transform.pyx | 26 +++++++++++++-------------
 1 file changed, 13 insertions(+), 13 deletions(-)

diff --git a/skimage/transform/_radon_transform.pyx b/skimage/transform/_radon_transform.pyx
index 321d9556..d17640ad 100644
--- a/skimage/transform/_radon_transform.pyx
+++ b/skimage/transform/_radon_transform.pyx
@@ -1,13 +1,12 @@
 #cython: cdivision=True
-#cython: boundscheck=True
-#cython: nonecheck=True
+#cython: boundscheck=False
+#cython: nonecheck=False
 #cython: wraparound=False
 import numpy as np
-from numpy import pi
 
 cimport numpy as cnp
 cimport cython
-from libc.math cimport cos, sin, floor, ceil, sqrt, abs
+from libc.math cimport cos, sin, floor, ceil, sqrt, abs, M_PI
 
 
 cpdef bilinear_ray_sum(cnp.ndarray[cnp.double_t, ndim=2] image, double theta,
@@ -32,7 +31,7 @@ cpdef bilinear_ray_sum(cnp.ndarray[cnp.double_t, ndim=2] image, double theta,
         A measure of how long the ray's path through the reconstruction
         circle was
     """
-    theta = theta / 180. * pi
+    theta = theta / 180. * M_PI
     cdef double radius = image.shape[0] // 2 - 1
     cdef double projection_center = image.shape[0] // 2 - 1
     cdef double rotation_center = image.shape[0] // 2
@@ -41,7 +40,7 @@ cpdef bilinear_ray_sum(cnp.ndarray[cnp.double_t, ndim=2] image, double theta,
     # s0 is the half-length of the ray's path in the reconstruction circle
     cdef double s0
     s0 = sqrt(radius**2 - t**2) if radius**2 >= t**2 else 0.
-    cdef Py_ssize_t Ns = 2 * int(ceil(2 * s0))  # number of steps along the ray
+    cdef Py_ssize_t Ns = 2 * (<Py_ssize_t> ceil(2 * s0))  # number of steps along the ray
     cdef double ray_sum = 0.
     cdef double weight_norm = 0.
     cdef double ds, dx, dy, x0, y0, x, y, di, dj, index_i, index_j
@@ -110,7 +109,7 @@ cpdef bilinear_ray_update(cnp.ndarray[cnp.double_t, ndim=2] image,
         deviation = -(ray_sum - projected_value) / weight_norm
     else:
         deviation = 0.
-    theta = theta / 180. * pi
+    theta = theta / 180. * M_PI
     cdef double radius = image.shape[0] // 2 - 1
     cdef double projection_center = image.shape[0] // 2 - 1
     cdef double rotation_center = image.shape[0] // 2
@@ -119,12 +118,12 @@ cpdef bilinear_ray_update(cnp.ndarray[cnp.double_t, ndim=2] image,
     # s0 is the half-length of the ray's path in the reconstruction circle
     cdef double s0
     s0 = sqrt(radius*radius - t*t) if radius**2 >= t**2 else 0.
-    cdef unsigned int Ns = 2 * int(ceil(2 * s0))
+    cdef Py_ssize_t Ns = 2 * (<Py_ssize_t> ceil(2 * s0))
     cdef double hamming_beta = 0.46164    # beta for equiripple Hamming window
 
     cdef double ds, dx, dy, x0, y0, x, y, di, dj, index_i, index_j
     cdef double hamming_window
-    cdef unsigned int k, i, j
+    cdef Py_ssize_t k, i, j
     if Ns > 0:
         # Step length between samples
         ds = 2 * s0 / Ns
@@ -143,7 +142,7 @@ cpdef bilinear_ray_update(cnp.ndarray[cnp.double_t, ndim=2] image,
             di = index_i - floor(index_i)
             dj = index_j - floor(index_j)
             hamming_window = ((1 - hamming_beta)
-                              - hamming_beta * cos(2*pi*k / (Ns - 1)))
+                              - 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)
@@ -159,9 +158,10 @@ cpdef bilinear_ray_update(cnp.ndarray[cnp.double_t, ndim=2] image,
     return deviation
 
 
-def sart_projection_update(cnp.ndarray[cnp.double_t, ndim=2] image, \
-                           double theta, \
-                           cnp.ndarray[cnp.double_t, ndim=1] projection,
+@cython.boundscheck(True)
+def sart_projection_update(cnp.ndarray[cnp.double_t, ndim=2] image not None,
+                           double theta,
+                           cnp.ndarray[cnp.double_t, ndim=1] projection not None,
                            double projection_shift=0.):
     """
     Compute update to a reconstruction estimate from a single projection

From c977c5974dd6a1f31cbd402c5bb6c93627f028c6 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Jostein=20B=C3=B8=20Fl=C3=B8ystad?=
 <jostein.floystad@gmail.com>
Date: Mon, 10 Jun 2013 11:45:47 +0200
Subject: [PATCH 12/34] iradon_sart: Improve argument checking.

---
 skimage/transform/radon_transform.py | 22 +++++++++++++++++++---
 1 file changed, 19 insertions(+), 3 deletions(-)

diff --git a/skimage/transform/radon_transform.py b/skimage/transform/radon_transform.py
index bba73c3f..1ca56feb 100644
--- a/skimage/transform/radon_transform.py
+++ b/skimage/transform/radon_transform.py
@@ -345,17 +345,33 @@ def iradon_sart(radon_image, theta=None, image=None, projection_shifts=None,
         -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)
-    reconstructed_shape = (radon_image.shape[0], radon_image.shape[0])
+    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)
-    if projection_shifts is None:
-        projection_shifts = np.zeros((radon_image.shape[1],), 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 _sart_order_angles(theta):
         image_update = sart_projection_update(image, theta[angle_index],
                                               radon_image[:, angle_index],

From 91aeb658ddb7d6b9867fcd0e69cda4fb5380c71a Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Jostein=20B=C3=B8=20Fl=C3=B8ystad?=
 <jostein.floystad@gmail.com>
Date: Mon, 10 Jun 2013 13:22:55 +0200
Subject: [PATCH 13/34] iradon_sart's cython exts: Use typed memoryviews, not
 ndarrays.

---
 skimage/transform/_radon_transform.pyx | 59 +++++++++++++-------------
 1 file changed, 30 insertions(+), 29 deletions(-)

diff --git a/skimage/transform/_radon_transform.pyx b/skimage/transform/_radon_transform.pyx
index d17640ad..8622666a 100644
--- a/skimage/transform/_radon_transform.pyx
+++ b/skimage/transform/_radon_transform.pyx
@@ -9,8 +9,8 @@ cimport cython
 from libc.math cimport cos, sin, floor, ceil, sqrt, abs, M_PI
 
 
-cpdef bilinear_ray_sum(cnp.ndarray[cnp.double_t, ndim=2] image, double theta,
-                       double ray_position):
+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.
 
@@ -32,18 +32,18 @@ cpdef bilinear_ray_sum(cnp.ndarray[cnp.double_t, ndim=2] image, double theta,
         circle was
     """
     theta = theta / 180. * M_PI
-    cdef double radius = image.shape[0] // 2 - 1
-    cdef double projection_center = image.shape[0] // 2 - 1
-    cdef double rotation_center = image.shape[0] // 2
+    cdef cnp.double_t radius = image.shape[0] // 2 - 1
+    cdef cnp.double_t projection_center = image.shape[0] // 2 - 1
+    cdef cnp.double_t rotation_center = image.shape[0] // 2
     # (s, t) is the (x, y) system rotated by theta
-    cdef double t = ray_position - projection_center
+    cdef cnp.double_t t = ray_position - projection_center
     # s0 is the half-length of the ray's path in the reconstruction circle
-    cdef double s0
+    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 double ray_sum = 0.
-    cdef double weight_norm = 0.
-    cdef double ds, dx, dy, x0, y0, x, y, di, dj, index_i, index_j
+    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, index_i, index_j
     cdef Py_ssize_t k, i, j
     if Ns > 0:
         # step length between samples
@@ -79,9 +79,10 @@ cpdef bilinear_ray_sum(cnp.ndarray[cnp.double_t, ndim=2] image, double theta,
     return ray_sum, weight_norm
 
 
-cpdef bilinear_ray_update(cnp.ndarray[cnp.double_t, ndim=2] image,
-        cnp.ndarray[cnp.double_t, ndim=2] image_update,
-        double theta, double ray_position, double projected_value):
+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.
 
@@ -103,26 +104,26 @@ cpdef bilinear_ray_update(cnp.ndarray[cnp.double_t, ndim=2] image,
     deviation :
         Deviation before updating the image
     """
-    cdef double ray_sum, weight_norm, deviation
+    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 double radius = image.shape[0] // 2 - 1
-    cdef double projection_center = image.shape[0] // 2 - 1
-    cdef double rotation_center = image.shape[0] // 2
+    cdef cnp.double_t radius = image.shape[0] // 2 - 1
+    cdef cnp.double_t projection_center = image.shape[0] // 2 - 1
+    cdef cnp.double_t rotation_center = image.shape[0] // 2
     # (s, t) is the (x, y) system rotated by theta
-    cdef double t = ray_position - projection_center
+    cdef cnp.double_t t = ray_position - projection_center
     # s0 is the half-length of the ray's path in the reconstruction circle
-    cdef double s0
+    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 double hamming_beta = 0.46164    # beta for equiripple Hamming window
+    cdef cnp.double_t hamming_beta = 0.46164    # beta for equiripple Hamming window
 
-    cdef double ds, dx, dy, x0, y0, x, y, di, dj, index_i, index_j
-    cdef double 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
@@ -159,10 +160,10 @@ cpdef bilinear_ray_update(cnp.ndarray[cnp.double_t, ndim=2] image,
 
 
 @cython.boundscheck(True)
-def sart_projection_update(cnp.ndarray[cnp.double_t, ndim=2] image not None,
-                           double theta,
-                           cnp.ndarray[cnp.double_t, ndim=1] projection not None,
-                           double projection_shift=0.):
+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.
@@ -185,11 +186,11 @@ def sart_projection_update(cnp.ndarray[cnp.double_t, ndim=2] image not None,
         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 double ray_position
+    cdef cnp.double_t[:, :] 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
+    return np.asarray(image_update)

From 5baaf785646cefbabd86bcb1d61e85f56fd8caa3 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Jostein=20B=C3=B8=20Fl=C3=B8ystad?=
 <jostein.floystad@gmail.com>
Date: Mon, 10 Jun 2013 22:13:43 +0200
Subject: [PATCH 14/34] iradon_sart: Reduce code duplication in interpolation.

---
 skimage/transform/_radon_transform.pyx | 26 ++++++++++++++++----------
 1 file changed, 16 insertions(+), 10 deletions(-)

diff --git a/skimage/transform/_radon_transform.pyx b/skimage/transform/_radon_transform.pyx
index 8622666a..bd2f356c 100644
--- a/skimage/transform/_radon_transform.pyx
+++ b/skimage/transform/_radon_transform.pyx
@@ -40,10 +40,12 @@ cpdef bilinear_ray_sum(cnp.double_t[:, :] image, cnp.double_t theta,
     # 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 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, index_i, index_j
+    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
@@ -65,17 +67,21 @@ cpdef bilinear_ray_sum(cnp.double_t[:, :] image, cnp.double_t theta,
             # Use linear interpolation between values
             # Where values fall outside the array, assume zero
             if i > 0 and j > 0:
-                ray_sum += (1. - di) * (1. - dj) * image[i, j] * ds
-                weight_norm += ((1 - di) * (1 - dj) * ds)**2
+                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:
-                ray_sum += (1. - di) * dj * image[i, j+1] * ds
-                weight_norm += ((1 - di) * dj * ds)**2
+                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:
-                ray_sum += di * (1 - dj) * image[i+1, j] * ds
-                weight_norm += (di * (1 - dj) * ds)**2
+                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:
-                ray_sum += di * dj * image[i+1, j+1] * ds
-                weight_norm += (di * dj * ds)**2
+                weight = di * dj * ds
+                ray_sum += weight * image[i+1, j+1]
+                weight_norm += weight**2
     return ray_sum, weight_norm
 
 

From 89f227c5ffbb7a5c405c1a81d8dd3e6009a4e70e Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Jostein=20B=C3=B8=20Fl=C3=B8ystad?=
 <jostein.floystad@gmail.com>
Date: Sat, 15 Jun 2013 17:10:03 +0200
Subject: [PATCH 15/34] iradon_sart: Remove needless memoryview/ndarray
 conversion.

image_update is not manipulated in sart_projection_update and the
conversion to memoryview will happen when it is passed to the
subroutines anyway.
---
 skimage/transform/_radon_transform.pyx | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/skimage/transform/_radon_transform.pyx b/skimage/transform/_radon_transform.pyx
index bd2f356c..d6112a57 100644
--- a/skimage/transform/_radon_transform.pyx
+++ b/skimage/transform/_radon_transform.pyx
@@ -192,11 +192,11 @@ def sart_projection_update(cnp.double_t[:, :] image not None,
         Array of same shape as ``image`` containing updates that should be
         added to ``image`` to improve the reconstruction estimate
     """
-    cdef cnp.double_t[:, :] image_update = np.zeros_like(image)
+    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 np.asarray(image_update)
+    return image_update

From 8d4b1e671014693742f82bf7967fc85f248b5e10 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Jostein=20B=C3=B8=20Fl=C3=B8ystad?=
 <jostein.floystad@gmail.com>
Date: Sun, 16 Jun 2013 09:37:45 +0200
Subject: [PATCH 16/34] iradon_sart: Add Kaczmarz reference and reformat
 citations.

---
 skimage/transform/radon_transform.py | 5 ++++-
 1 file changed, 4 insertions(+), 1 deletion(-)

diff --git a/skimage/transform/radon_transform.py b/skimage/transform/radon_transform.py
index 1ca56feb..b69439a2 100644
--- a/skimage/transform/radon_transform.py
+++ b/skimage/transform/radon_transform.py
@@ -337,11 +337,14 @@ def iradon_sart(radon_image, theta=None, image=None, projection_shifts=None,
     information, but will also often increase the noise.
 
     References:
-        -A. C. Kak, Malcolm Slaney, "Principles of Computerized Tomographic
+        -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)
         -Kaczmarz' method, Wikipedia,
         http://en.wikipedia.org/wiki/Kaczmarz_method
     """

From dcd31c18826ba491ee4442601e3ea66696cbcba3 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Jostein=20B=C3=B8=20Fl=C3=B8ystad?=
 <jostein.floystad@gmail.com>
Date: Sun, 16 Jun 2013 10:18:48 +0200
Subject: [PATCH 17/34] radon_transform: Declare encoding.

---
 skimage/transform/radon_transform.py | 1 +
 1 file changed, 1 insertion(+)

diff --git a/skimage/transform/radon_transform.py b/skimage/transform/radon_transform.py
index b69439a2..32572f68 100644
--- a/skimage/transform/radon_transform.py
+++ b/skimage/transform/radon_transform.py
@@ -1,3 +1,4 @@
+# -*- coding: utf-8 -*-
 """
 radon.py - Radon and inverse radon transforms
 

From 2082b57e4171225caa40ab94990e8ae168f7fb49 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Jostein=20B=C3=B8=20Fl=C3=B8ystad?=
 <jostein.floystad@gmail.com>
Date: Sun, 16 Jun 2013 14:24:15 +0200
Subject: [PATCH 18/34] iradon_sart: Order angles using a golden ratio
 approach.

---
 skimage/transform/radon_transform.py | 76 +++++++++++++++++++---------
 1 file changed, 53 insertions(+), 23 deletions(-)

diff --git a/skimage/transform/radon_transform.py b/skimage/transform/radon_transform.py
index 32572f68..82f6bd21 100644
--- a/skimage/transform/radon_transform.py
+++ b/skimage/transform/radon_transform.py
@@ -258,31 +258,58 @@ def iradon(radon_image, theta=None, output_size=None,
     return reconstructed * np.pi / (2 * len(th))
 
 
-def _sart_order_angles(theta):
+def order_angles_golden_ratio(theta):
     """
     Order angles to reduce the amount of correlated information
-    in subsequent projections, i.e. make sure subsequent angles
-    are as far away from each other mod 180 degrees as possible.
-    Indices into the ``theta`` array are yielded.
+    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.
     """
-    tau = 3.    # time constant for correlations; 0.1 < tau < 100 works well
-    used_indices = [0]
-    remaining_indices = range(1, len(theta))
-    yield 0
-    while remaining_indices:
-        used = np.array(theta[used_indices])
-        used.shape = (-1, 1)
-        remaining = np.array(theta[remaining_indices])
-        remaining.shape = (1, -1)
-        time = (np.arange(used.shape[0]) + 1)[::-1]
-        time.shape = (-1, 1)
-        difference = used - remaining
-        distance = np.minimum(abs(difference % 180), abs(difference % -180))
-        cost = np.exp(-distance * time / tau).sum(axis=0).squeeze()
-        next_angle_remaining_index = np.argmin(cost)
-        next_angle_index = remaining_indices.pop(next_angle_remaining_index)
-        used_indices.append(next_angle_index)
-        yield next_angle_index
+    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 ration 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,
@@ -346,6 +373,9 @@ def iradon_sart(radon_image, theta=None, image=None, projection_shifts=None,
         -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
     """
@@ -376,7 +406,7 @@ def iradon_sart(radon_image, theta=None, image=None, projection_shifts=None,
         clip = (float(clip[0]), float(clip[1]))
     relaxation = float(relaxation)
 
-    for angle_index in _sart_order_angles(theta):
+    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])

From ace07a0b18987e734511e8f865b6140525a0f584 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Jostein=20B=C3=B8=20Fl=C3=B8ystad?=
 <jostein.floystad@gmail.com>
Date: Sun, 16 Jun 2013 14:24:47 +0200
Subject: [PATCH 19/34] iradon_sart: Add test for order_angles_golden_ratio.

---
 skimage/transform/tests/test_radon_transform.py | 13 +++++++++++++
 1 file changed, 13 insertions(+)

diff --git a/skimage/transform/tests/test_radon_transform.py b/skimage/transform/tests/test_radon_transform.py
index 5ad87877..e5df9b8d 100644
--- a/skimage/transform/tests/test_radon_transform.py
+++ b/skimage/transform/tests/test_radon_transform.py
@@ -298,6 +298,19 @@ 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 116e1dd57140fe99c9ff7456a71c50a0595c3fc3 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Jostein=20B=C3=B8=20Fl=C3=B8ystad?=
 <jostein.floystad@gmail.com>
Date: Sun, 16 Jun 2013 14:25:23 +0200
Subject: [PATCH 20/34] iradon_sart: Also test accuracy with a missing wedge.

---
 .../transform/tests/test_radon_transform.py   | 91 ++++++++++---------
 1 file changed, 47 insertions(+), 44 deletions(-)

diff --git a/skimage/transform/tests/test_radon_transform.py b/skimage/transform/tests/test_radon_transform.py
index e5df9b8d..ae25c970 100644
--- a/skimage/transform/tests/test_radon_transform.py
+++ b/skimage/transform/tests/test_radon_transform.py
@@ -320,54 +320,57 @@ def test_iradon_sart():
 
     shepp_logan = imread(os.path.join(data_dir, "phantom.png"), as_grey=True)
     image = rescale(shepp_logan, scale=0.4)
-    theta = np.linspace(0., 180., image.shape[0], endpoint=False)
-    sinogram = radon(image, theta, circle=True)
-    reconstructed = iradon_sart(sinogram, theta)
+    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()
+        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.025
-    reconstructed = iradon_sart(sinogram, theta, reconstructed)
-    delta = np.mean(np.abs(reconstructed - image))
-    print('delta (2 iterations) =', delta)
-    assert delta < 0.015
+        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
 
-    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()
+        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.025
+        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

From 266b2022262a8cbcad9af6962030cc8fd4826c00 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Jostein=20B=C3=B8=20Fl=C3=B8ystad?=
 <jostein.floystad@gmail.com>
Date: Sun, 16 Jun 2013 14:34:53 +0200
Subject: [PATCH 21/34] iradon_sart: Style fixes.

---
 skimage/transform/radon_transform.py | 1 +
 1 file changed, 1 insertion(+)

diff --git a/skimage/transform/radon_transform.py b/skimage/transform/radon_transform.py
index 82f6bd21..0f26e343 100644
--- a/skimage/transform/radon_transform.py
+++ b/skimage/transform/radon_transform.py
@@ -288,6 +288,7 @@ def order_angles_golden_ratio(theta):
         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))

From dde4288865f48beca0f81b7ee91a4ddec6c5a664 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Jostein=20B=C3=B8=20Fl=C3=B8ystad?=
 <jostein.floystad@gmail.com>
Date: Sun, 16 Jun 2013 14:37:31 +0200
Subject: [PATCH 22/34] test_radon_transform: Style fixes.

---
 skimage/transform/tests/test_radon_transform.py | 5 +++--
 1 file changed, 3 insertions(+), 2 deletions(-)

diff --git a/skimage/transform/tests/test_radon_transform.py b/skimage/transform/tests/test_radon_transform.py
index ae25c970..74d20a24 100644
--- a/skimage/transform/tests/test_radon_transform.py
+++ b/skimage/transform/tests/test_radon_transform.py
@@ -351,8 +351,9 @@ def test_iradon_sart():
         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
+        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 a5df8d463085e31dbaf8fecdf11945df4112c599 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Jostein=20B=C3=B8=20Fl=C3=B8ystad?=
 <jostein.floystad@gmail.com>
Date: Sun, 16 Jun 2013 18:43:49 +0200
Subject: [PATCH 23/34] Rewrite Radon example; include SART.

The old example had some flaws. The new example corrects these, expands
on the topic and adds content relating to the newly implemented SART
algorithm.
---
 doc/examples/plot_radon_transform.py | 177 ++++++++++++++++++++++-----
 1 file changed, 147 insertions(+), 30 deletions(-)

diff --git a/doc/examples/plot_radon_transform.py b/doc/examples/plot_radon_transform.py
index b4b9721d..df3b1577 100644
--- a/doc/examples/plot_radon_transform.py
+++ b/doc/examples/plot_radon_transform.py
@@ -3,55 +3,172 @@
 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. 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 with the rays making a different angle with
+the axes of the object. A collection of projections at several angles is
+called a sinogram.
+
+The inverse Radon transform is used in computed tomography to reconstruct
+a 2D image from can hence be used to reconstruct an object 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.
+
+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.
 
 For more information see:
 
+  - AC Kak, M Slaney, "Principles of Computerized Tomographic Imaging",
+    http://www.slaney.org/pct/pct-toc.html
   - http://en.wikipedia.org/wiki/Radon_transform
-  - http://www.clear.rice.edu/elec431/projects96/DSP/bpanalysis.html
 
-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 will be used: The Filtered Back
+Projection (FBP) and the Simultaneous Algebraic Reconstruction
+Technique (SART).
 
+
+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_:
 """
 
+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)
-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)
+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 (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 (http://en.wikipedia.org/wiki/Projection-slice_theorem). 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 is needed 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)
+
+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)
+
+"""
+.. 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, 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). It uses Kaczmarz' method 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
+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)
+
+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)
+
+d = reconstruction_sart - image
+print(d.max(), d.min())
+
+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)
+
+"""
+.. image:: PLOT2RST.current_figure
+"""

From a4870242b63458ed9ee79db6a951007c94deebf4 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Jostein=20B=C3=B8=20Fl=C3=B8ystad?=
 <jostein.floystad@gmail.com>
Date: Mon, 24 Jun 2013 10:18:54 +0200
Subject: [PATCH 24/34] iradon_sart: Test clip functionality and add a test for
 it.

---
 skimage/transform/radon_transform.py            | 2 +-
 skimage/transform/tests/test_radon_transform.py | 3 +++
 2 files changed, 4 insertions(+), 1 deletion(-)

diff --git a/skimage/transform/radon_transform.py b/skimage/transform/radon_transform.py
index 0f26e343..ff3af596 100644
--- a/skimage/transform/radon_transform.py
+++ b/skimage/transform/radon_transform.py
@@ -413,5 +413,5 @@ def iradon_sart(radon_image, theta=None, image=None, projection_shifts=None,
                                               projection_shifts[angle_index])
         image += relaxation * image_update
         if not clip is None:
-            image = clip(image, clip[0], clip[1])
+            image = np.clip(image, clip[0], clip[1])
     return image
diff --git a/skimage/transform/tests/test_radon_transform.py b/skimage/transform/tests/test_radon_transform.py
index 74d20a24..11fced56 100644
--- a/skimage/transform/tests/test_radon_transform.py
+++ b/skimage/transform/tests/test_radon_transform.py
@@ -347,6 +347,9 @@ def test_iradon_sart():
         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))
+        print('delta (1 iteration, clip) =', delta)
+        assert delta < 0.013 * error_factor
 
         np.random.seed(1239867)
         shifts = np.random.uniform(-3, 3, sinogram.shape[1])

From 9e11e453b64424182e9221f32425df12612177f0 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Jostein=20B=C3=B8=20Fl=C3=B8ystad?=
 <jostein.floystad@gmail.com>
Date: Mon, 24 Jun 2013 10:20:16 +0200
Subject: [PATCH 25/34] radon example: Display plots when run from command
 line.

Call to matplotlib's show() was missing.
---
 doc/examples/plot_radon_transform.py | 3 +++
 1 file changed, 3 insertions(+)

diff --git a/doc/examples/plot_radon_transform.py b/doc/examples/plot_radon_transform.py
index df3b1577..fb2d9dbd 100644
--- a/doc/examples/plot_radon_transform.py
+++ b/doc/examples/plot_radon_transform.py
@@ -76,6 +76,7 @@ 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
@@ -107,6 +108,7 @@ 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
@@ -168,6 +170,7 @@ 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

From e5dfcb0ab45dabdfc670b8badc7352857ffdc094 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Jostein=20B=C3=B8=20Fl=C3=B8ystad?=
 <jostein.floystad@gmail.com>
Date: Thu, 4 Jul 2013 15:23:51 +0200
Subject: [PATCH 26/34] Improve the text in and add references to the Radon
 example.

Based on comments from Emmanuelle.
---
 doc/examples/plot_radon_transform.py | 75 +++++++++++++++++-----------
 1 file changed, 45 insertions(+), 30 deletions(-)

diff --git a/doc/examples/plot_radon_transform.py b/doc/examples/plot_radon_transform.py
index fb2d9dbd..4a2c01a5 100644
--- a/doc/examples/plot_radon_transform.py
+++ b/doc/examples/plot_radon_transform.py
@@ -4,38 +4,36 @@ Radon transform
 ===============
 
 In computed tomography, the tomography reconstruction problem is to obtain
-a tomographic slice image from a set of projections. A projection is formed
+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 with the rays making a different angle with
-the axes of the object. A collection of projections at several angles is
-called a sinogram.
+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.
 
 The inverse Radon transform is used in computed tomography to reconstruct
-a 2D image from can hence be used to reconstruct an object 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.
+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.
 
 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.
 
-For more information see:
-
-  - AC Kak, M Slaney, "Principles of Computerized Tomographic Imaging",
-    http://www.slaney.org/pct/pct-toc.html
-  - http://en.wikipedia.org/wiki/Radon_transform
-
 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 will be used: The Filtered Back
+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
 =====================
@@ -85,14 +83,14 @@ Reconstruction with the Filtered Back Projection (FBP)
 ======================================================
 
 The mathematical foundation of the filtered back projection is the Fourier
-slice theorem (http://en.wikipedia.org/wiki/Projection-slice_theorem). 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 is needed to suppress high frequency noise in the
-reconstruction. ``skimage`` provides a few different options for the filter.
+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.
 
 """
 
@@ -117,15 +115,16 @@ 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
+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, which has the property
-that the solution will approach a least-squares solution of the equation set.
+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
@@ -134,12 +133,15 @@ with relative ease.
 
 ``skimage`` provides one of the more popular variations of the algebraic
 reconstruction techniques: the Simultaneous Algebraic Reconstruction Technique
-(SART). It uses Kaczmarz' method as the iterative solver. A good
+(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
-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.
+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.
 
 """
 
@@ -174,4 +176,17 @@ 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)
+
 """

From 1d611878a6c3cbd56cbc55a7abf6220643d3459a Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Jostein=20B=C3=B8=20Fl=C3=B8ystad?=
 <jostein.floystad@gmail.com>
Date: Thu, 4 Jul 2013 15:24:44 +0200
Subject: [PATCH 27/34] Radon example: Style fixes.

---
 doc/examples/plot_radon_transform.py | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/doc/examples/plot_radon_transform.py b/doc/examples/plot_radon_transform.py
index 4a2c01a5..4b5b8d8d 100644
--- a/doc/examples/plot_radon_transform.py
+++ b/doc/examples/plot_radon_transform.py
@@ -67,9 +67,9 @@ plt.imshow(image, cmap=plt.cm.Greys_r)
 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 (deg)");
-plt.ylabel("Projection position (pixels)");
+plt.title("Radon transform\n(Sinogram)")
+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')
 

From 789c19923d7b60b471e0f7e240d8d2957c0adb10 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Jostein=20B=C3=B8=20Fl=C3=B8ystad?=
 <jostein.floystad@gmail.com>
Date: Thu, 4 Jul 2013 15:29:53 +0200
Subject: [PATCH 28/34] Radon example: print RMS errors of the reconstructions.

---
 doc/examples/plot_radon_transform.py | 13 ++++++++++---
 1 file changed, 10 insertions(+), 3 deletions(-)

diff --git a/doc/examples/plot_radon_transform.py b/doc/examples/plot_radon_transform.py
index 4b5b8d8d..d457d5a6 100644
--- a/doc/examples/plot_radon_transform.py
+++ b/doc/examples/plot_radon_transform.py
@@ -48,6 +48,8 @@ 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
 
@@ -97,6 +99,8 @@ 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))
@@ -148,6 +152,9 @@ 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))
 
@@ -162,9 +169,9 @@ plt.imshow(reconstruction_sart - image, cmap=plt.cm.Greys_r, **imkwargs)
 # from the first iteration as an initial estimate
 reconstruction_sart2 = iradon_sart(sinogram, theta=theta,
                                    image=reconstruction_sart)
-
-d = reconstruction_sart - image
-print(d.max(), d.min())
+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")

From e2a0f7fff810c08e71db632782392fd533a949f9 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Jostein=20B=C3=B8=20Fl=C3=B8ystad?=
 <jostein.floystad@gmail.com>
Date: Thu, 4 Jul 2013 16:39:35 +0200
Subject: [PATCH 29/34] Radon example: Declare UTF-8 encoding.

This is needed for non-ASCII characters in the references to papers
cited in the text.
---
 doc/examples/plot_radon_transform.py | 1 +
 1 file changed, 1 insertion(+)

diff --git a/doc/examples/plot_radon_transform.py b/doc/examples/plot_radon_transform.py
index d457d5a6..61175f8a 100644
--- a/doc/examples/plot_radon_transform.py
+++ b/doc/examples/plot_radon_transform.py
@@ -1,3 +1,4 @@
+# -*- coding: utf-8 -*-
 """
 ===============
 Radon transform

From 06857d3b92c65127b11853ab2bd41824d47693cf Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Jostein=20B=C3=B8=20Fl=C3=B8ystad?=
 <jostein.floystad@gmail.com>
Date: Thu, 4 Jul 2013 17:35:52 +0200
Subject: [PATCH 30/34] Radon tests: Check the correct error metric.

---
 skimage/transform/tests/test_radon_transform.py | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/skimage/transform/tests/test_radon_transform.py b/skimage/transform/tests/test_radon_transform.py
index 11fced56..333c298d 100644
--- a/skimage/transform/tests/test_radon_transform.py
+++ b/skimage/transform/tests/test_radon_transform.py
@@ -348,8 +348,9 @@ def test_iradon_sart():
         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.013 * error_factor
+        assert delta < 0.015 * error_factor
 
         np.random.seed(1239867)
         shifts = np.random.uniform(-3, 3, sinogram.shape[1])

From f2790d658c0e6f94bc1a479c4eed162907bf5a05 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Jostein=20B=C3=B8=20Fl=C3=B8ystad?=
 <jostein.floystad@gmail.com>
Date: Fri, 5 Jul 2013 13:02:03 +0200
Subject: [PATCH 31/34] iradon_sart: redefine projection center.

This reflects the bugfixes related to issue gh-592 implemented in
gh-596.
---
 skimage/transform/_radon_transform.pyx | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/skimage/transform/_radon_transform.pyx b/skimage/transform/_radon_transform.pyx
index d6112a57..3f865a93 100644
--- a/skimage/transform/_radon_transform.pyx
+++ b/skimage/transform/_radon_transform.pyx
@@ -33,7 +33,7 @@ cpdef bilinear_ray_sum(cnp.double_t[:, :] image, cnp.double_t theta,
     """
     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 - 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
@@ -118,7 +118,7 @@ cpdef bilinear_ray_update(cnp.double_t[:, :] image,
         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 - 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

From af479414353d6655497e0e9a55dcfe23a766f2d6 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Jostein=20B=C3=B8=20Fl=C3=B8ystad?=
 <jostein.floystad@gmail.com>
Date: Fri, 5 Jul 2013 13:03:13 +0200
Subject: [PATCH 32/34] iradon_sart: fix comment spelling.

---
 skimage/transform/radon_transform.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/skimage/transform/radon_transform.py b/skimage/transform/radon_transform.py
index ff3af596..944b4a6b 100644
--- a/skimage/transform/radon_transform.py
+++ b/skimage/transform/radon_transform.py
@@ -298,7 +298,7 @@ def order_angles_golden_ratio(theta):
     index = remaining.pop(0)
     angle = theta[index]
     yield index
-    # determine subsequent angles using the golden ration method
+    # determine subsequent angles using the golden ratio method
     angle_increment = interval * (1 - (np.sqrt(5) - 1) / 2)
     while remaining:
         angle = (angle + angle_increment) % interval

From 3cdfeccc7a80884676a3f364e406c4d7a47d1e98 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Jostein=20B=C3=B8=20Fl=C3=B8ystad?=
 <jostein.floystad@gmail.com>
Date: Fri, 5 Jul 2013 13:12:34 +0200
Subject: [PATCH 33/34] iradon_sart: flip signs to reflect changes to radon.

The changes to radon were introduced in gh-596.
---
 skimage/transform/_radon_transform.pyx | 16 ++++++++--------
 1 file changed, 8 insertions(+), 8 deletions(-)

diff --git a/skimage/transform/_radon_transform.pyx b/skimage/transform/_radon_transform.pyx
index 3f865a93..91b943b9 100644
--- a/skimage/transform/_radon_transform.pyx
+++ b/skimage/transform/_radon_transform.pyx
@@ -50,11 +50,11 @@ cpdef bilinear_ray_sum(cnp.double_t[:, :] image, cnp.double_t theta,
     if Ns > 0:
         # step length between samples
         ds = 2 * s0 / Ns
-        dx = ds * cos(theta)
-        dy = ds * sin(theta)
+        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)
+        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
@@ -134,11 +134,11 @@ cpdef bilinear_ray_update(cnp.double_t[:, :] image,
     if Ns > 0:
         # Step length between samples
         ds = 2 * s0 / Ns
-        dx = ds * cos(theta)
-        dy = ds * sin(theta)
+        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)
+        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

From 553767122bd18dec9dd1b75b25a272279d32b09d Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Jostein=20B=C3=B8=20Fl=C3=B8ystad?=
 <jostein.floystad@gmail.com>
Date: Fri, 5 Jul 2013 13:36:55 +0200
Subject: [PATCH 34/34] Add SART to contributors.

---
 CONTRIBUTORS.txt | 1 +
 1 file changed, 1 insertion(+)

diff --git a/CONTRIBUTORS.txt b/CONTRIBUTORS.txt
index e03271b8..ac11457c 100644
--- a/CONTRIBUTORS.txt
+++ b/CONTRIBUTORS.txt
@@ -144,3 +144,4 @@
 
 - Jostein Bø Fløystad
   Reconstruction circle mode for Radon transform
+  Simultaneous Algebraic Reconstruction Technique for inverse Radon transform