From 4e0cbf97feb81258663b5a8b59cb94f3afb1dff6 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Jostein=20B=C3=B8=20Fl=C3=B8ystad?=
 <jostein.floystad@gmail.com>
Date: Tue, 18 Jun 2013 11:46:50 +0200
Subject: [PATCH 01/21] transform.radon: Add test to verify the projection
 center.

This test is designed for issue gh-592.
---
 .../transform/tests/test_radon_transform.py   | 25 +++++++++++++++++++
 1 file changed, 25 insertions(+)

diff --git a/skimage/transform/tests/test_radon_transform.py b/skimage/transform/tests/test_radon_transform.py
index 79cf80c5..8a85d58c 100644
--- a/skimage/transform/tests/test_radon_transform.py
+++ b/skimage/transform/tests/test_radon_transform.py
@@ -3,6 +3,7 @@ from __future__ import division
 
 import numpy as np
 from numpy.testing import *
+from numpy.fft import ifftshift, ifftn
 import itertools
 from skimage.transform import *
 
@@ -14,6 +15,30 @@ def rescale(x):
     return x
 
 
+def check_radon_center(shape, circle):
+    # Determine the center of an array as defined by the fft
+    ft_image = np.abs(ifftshift(ifftn(np.ones(shape))))**2
+    fft_center = np.unravel_index(np.argmax(ft_image), shape)
+    print('fft_center =', fft_center)
+    # Create a test image with only a single non-zero pixel at the origin
+    image = np.zeros(shape, dtype=np.float)
+    image[fft_center] = 1.
+    # Calculate the sinogram
+    theta = np.linspace(0., 180., max(shape), endpoint=False)
+    sinogram = radon(image, theta=theta, circle=circle)
+    # The sinogram should be a straight, horizontal line
+    sinogram_max = np.argmax(sinogram, axis=0)
+    print(sinogram_max)
+    assert np.std(sinogram_max) < 1e-6
+
+
+def test_radon_center():
+    shapes = [(16, 16), (17, 17)]
+    circles = [False, True]
+    for shape, circle in itertools.product(shapes, circles):
+        yield check_radon_center, shape, circle
+
+
 def test_radon_iradon():
     size = 100
     debug = False

From b8a20bcb59874f55ddcd6c3a8c9901e5de06a7a5 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Jostein=20B=C3=B8=20Fl=C3=B8ystad?=
 <jostein.floystad@gmail.com>
Date: Tue, 18 Jun 2013 12:03:31 +0200
Subject: [PATCH 02/21] transform.radon: Add testcases for rectangular input
 arrays.

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

diff --git a/skimage/transform/tests/test_radon_transform.py b/skimage/transform/tests/test_radon_transform.py
index 8a85d58c..44d3ac77 100644
--- a/skimage/transform/tests/test_radon_transform.py
+++ b/skimage/transform/tests/test_radon_transform.py
@@ -37,6 +37,9 @@ def test_radon_center():
     circles = [False, True]
     for shape, circle in itertools.product(shapes, circles):
         yield check_radon_center, shape, circle
+    rectangular_shapes = [(32, 16), (33, 17)]
+    for shape in rectangular_shapes:
+        yield check_radon_center, shape, False
 
 
 def test_radon_iradon():

From 1d64eb59eb535d84dd19f52d2d33d110a4a35c1c Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Jostein=20B=C3=B8=20Fl=C3=B8ystad?=
 <jostein.floystad@gmail.com>
Date: Tue, 18 Jun 2013 13:36:21 +0200
Subject: [PATCH 03/21] transform.iradon: Add tests for center of projection.

This is a test designed for resolving gh-592.
---
 .../transform/tests/test_radon_transform.py   | 46 +++++++++++++++++++
 1 file changed, 46 insertions(+)

diff --git a/skimage/transform/tests/test_radon_transform.py b/skimage/transform/tests/test_radon_transform.py
index 44d3ac77..8e31c658 100644
--- a/skimage/transform/tests/test_radon_transform.py
+++ b/skimage/transform/tests/test_radon_transform.py
@@ -42,6 +42,52 @@ def test_radon_center():
         yield check_radon_center, shape, False
 
 
+def check_iradon_center(size, theta, circle):
+    debug = False
+    # Create a test sinogram corresponding to a single projection
+    # with a single non-zero pixel at the rotation center
+    if circle:
+        sinogram = np.zeros((size, 1), dtype=np.float)
+        sinogram[size // 2, 0] = 1.
+    else:
+        diagonal = int(np.ceil(np.sqrt(2) * size))
+        sinogram = np.zeros((diagonal, 1), dtype=np.float)
+        sinogram[sinogram.shape[0] // 2, 0] = 1.
+    maxpoint = np.unravel_index(np.argmax(sinogram), sinogram.shape)
+    print('shape of generated sinogram', sinogram.shape)
+    print('maximum in generated sinogram', maxpoint)
+    # Compare reconstructions for theta=angle and theta=angle + 180;
+    # these should be exactly equal
+    reconstruction = iradon(sinogram, theta=[theta], circle=circle)
+    reconstruction_opposite = iradon(sinogram, theta=[theta + 180],
+                                     circle=circle)
+    print('rms deviance:',
+          np.sqrt(np.mean((reconstruction_opposite - reconstruction)**2)))
+    if debug:
+        import matplotlib.pyplot as plt
+        imkwargs = dict(cmap='gray', interpolation='nearest')
+        plt.figure()
+        plt.subplot(221)
+        plt.imshow(sinogram, **imkwargs)
+        plt.subplot(222)
+        plt.imshow(reconstruction_opposite - reconstruction, **imkwargs)
+        plt.subplot(223)
+        plt.imshow(reconstruction, **imkwargs)
+        plt.subplot(224)
+        plt.imshow(reconstruction_opposite, **imkwargs)
+        plt.show()
+
+    assert np.allclose(reconstruction, reconstruction_opposite)
+
+
+def test_iradon_center():
+    sizes = [16, 17]
+    thetas = [0, 90]
+    circles = [False, True]
+    for size, theta, circle in itertools.product(sizes, thetas, circles):
+        yield check_iradon_center, size, theta, circle
+
+
 def test_radon_iradon():
     size = 100
     debug = False

From 28de2f978abf5fa7f9a1c4d71548d36325f1da44 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Jostein=20B=C3=B8=20Fl=C3=B8ystad?=
 <jostein.floystad@gmail.com>
Date: Tue, 18 Jun 2013 18:23:09 +0200
Subject: [PATCH 04/21] transform.radon: Remove unneccesary matrix inverse.

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

diff --git a/skimage/transform/radon_transform.py b/skimage/transform/radon_transform.py
index 59e6c338..e3f46f30 100644
--- a/skimage/transform/radon_transform.py
+++ b/skimage/transform/radon_transform.py
@@ -98,17 +98,16 @@ def radon(image, theta=None, circle=False):
                        [0, 0, 1]])
 
     def build_rotation(theta):
-        T = -np.deg2rad(theta)
+        T = np.deg2rad(theta)
 
-        R = np.array([[np.cos(T), -np.sin(T), 0],
-                      [np.sin(T), np.cos(T), 0],
+        R = np.array([[np.cos(T), np.sin(T), 0],
+                      [-np.sin(T), np.cos(T), 0],
                       [0, 0, 1]])
 
         return shift1.dot(R).dot(shift0)
 
     for i in range(len(theta)):
-        rotated = _warp_fast(padded_image,
-                             np.linalg.inv(build_rotation(-theta[i])))
+        rotated = _warp_fast(padded_image, build_rotation(theta[i]))
 
         out[:, i] = rotated.sum(0)[::-1]
 

From 380c916c9221b144686de878eb0822fd38d06525 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Jostein=20B=C3=B8=20Fl=C3=B8ystad?=
 <jostein.floystad@gmail.com>
Date: Tue, 18 Jun 2013 18:36:23 +0200
Subject: [PATCH 05/21] transform.radon: Use correct padding for rectangular
 images.

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

diff --git a/skimage/transform/radon_transform.py b/skimage/transform/radon_transform.py
index e3f46f30..fffe4e6f 100644
--- a/skimage/transform/radon_transform.py
+++ b/skimage/transform/radon_transform.py
@@ -74,9 +74,9 @@ def radon(image, theta=None, circle=False):
         out = np.zeros((min(padded_image.shape), len(theta)))
     else:
         height, width = image.shape
-        diagonal = np.sqrt(height**2 + width**2)
-        heightpad = np.ceil(diagonal - height)
-        widthpad = np.ceil(diagonal - width)
+        diagonal = np.sqrt(2) * max(image.shape)
+        heightpad = int(np.ceil(diagonal - height))
+        widthpad = int(np.ceil(diagonal - width))
         padded_image = np.zeros((int(height + heightpad),
                                 int(width + widthpad)))
         y0 = int(np.ceil(heightpad / 2))

From a44f1d4ef92342027f2901aa77f70e86b98ad1ee Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Jostein=20B=C3=B8=20Fl=C3=B8ystad?=
 <jostein.floystad@gmail.com>
Date: Tue, 18 Jun 2013 18:49:24 +0200
Subject: [PATCH 06/21] transform.radon: Consistent definition of center of
 array.

The center is now defined as shape[i] // 2.
---
 skimage/transform/radon_transform.py            | 11 +++++------
 skimage/transform/tests/test_radon_transform.py |  6 +-----
 2 files changed, 6 insertions(+), 11 deletions(-)

diff --git a/skimage/transform/radon_transform.py b/skimage/transform/radon_transform.py
index fffe4e6f..2733e243 100644
--- a/skimage/transform/radon_transform.py
+++ b/skimage/transform/radon_transform.py
@@ -79,11 +79,10 @@ def radon(image, theta=None, circle=False):
         widthpad = int(np.ceil(diagonal - width))
         padded_image = np.zeros((int(height + heightpad),
                                 int(width + widthpad)))
-        y0 = int(np.ceil(heightpad / 2))
-        y1 = int((np.ceil(heightpad / 2) + height))
-        x0 = int((np.ceil(widthpad / 2)))
-        x1 = int((np.ceil(widthpad / 2) + width))
-
+        y0 = heightpad // 2
+        y1 = y0 + height
+        x0 = widthpad // 2
+        x1 = x0 + width
         padded_image[y0:y1, x0:x1] = image
         out = np.zeros((max(padded_image.shape), len(theta)))
 
@@ -109,7 +108,7 @@ def radon(image, theta=None, circle=False):
     for i in range(len(theta)):
         rotated = _warp_fast(padded_image, build_rotation(theta[i]))
 
-        out[:, i] = rotated.sum(0)[::-1]
+        out[:, i] = rotated.sum(0)
 
     return out
 
diff --git a/skimage/transform/tests/test_radon_transform.py b/skimage/transform/tests/test_radon_transform.py
index 8e31c658..9a03b1f2 100644
--- a/skimage/transform/tests/test_radon_transform.py
+++ b/skimage/transform/tests/test_radon_transform.py
@@ -16,13 +16,9 @@ def rescale(x):
 
 
 def check_radon_center(shape, circle):
-    # Determine the center of an array as defined by the fft
-    ft_image = np.abs(ifftshift(ifftn(np.ones(shape))))**2
-    fft_center = np.unravel_index(np.argmax(ft_image), shape)
-    print('fft_center =', fft_center)
     # Create a test image with only a single non-zero pixel at the origin
     image = np.zeros(shape, dtype=np.float)
-    image[fft_center] = 1.
+    image[(shape[0] // 2, shape[1] // 2)] = 1.
     # Calculate the sinogram
     theta = np.linspace(0., 180., max(shape), endpoint=False)
     sinogram = radon(image, theta=theta, circle=circle)

From 364e82176f6aa40533b4f4116eb49dd0c658d3c7 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Jostein=20B=C3=B8=20Fl=C3=B8ystad?=
 <jostein.floystad@gmail.com>
Date: Tue, 18 Jun 2013 22:03:56 +0200
Subject: [PATCH 07/21] radon tests: Refactor test for circle reconstructions.

---
 .../transform/tests/test_radon_transform.py   | 57 ++++++++++---------
 1 file changed, 29 insertions(+), 28 deletions(-)

diff --git a/skimage/transform/tests/test_radon_transform.py b/skimage/transform/tests/test_radon_transform.py
index 9a03b1f2..2f4884ea 100644
--- a/skimage/transform/tests/test_radon_transform.py
+++ b/skimage/transform/tests/test_radon_transform.py
@@ -214,40 +214,41 @@ def test_radon_circle():
     assert np.all(relative_error < 3e-3)
 
 
+def check_radon_iradon_circle(interpolation, shape, output_size):
+    # Forward and inverse radon on synthetic data
+    image = _random_circle(shape)
+    radius = min(shape) // 2
+    sinogram_rectangle = radon(image, circle=False)
+    reconstruction_rectangle = iradon(sinogram_rectangle,
+                                        output_size=output_size,
+                                        interpolation=interpolation,
+                                        circle=False)
+    sinogram_circle = radon(image, circle=True)
+    reconstruction_circle = iradon(sinogram_circle,
+                                    output_size=output_size,
+                                    interpolation=interpolation,
+                                    circle=True)
+    # Crop rectangular reconstruction to match circle=True reconstruction
+    width = reconstruction_circle.shape[0]
+    excess = int(np.ceil((reconstruction_rectangle.shape[0] - width) / 2))
+    s = np.s_[excess:width + excess, excess:width + excess]
+    reconstruction_rectangle = reconstruction_rectangle[s]
+    # Find the reconstruction circle, set reconstruction to zero outside
+    c0, c1 = np.ogrid[0:width, 0:width]
+    r = np.sqrt((c0 - width // 2)**2 + (c1 - width // 2)**2)
+    reconstruction_rectangle[r >= radius] = 0.
+    print(reconstruction_circle.shape)
+    print(reconstruction_rectangle.shape)
+    np.allclose(reconstruction_rectangle, reconstruction_circle)
+
+
 def test_radon_iradon_circle():
     shape = (61, 79)
-    radius = min(shape) // 2
-    image = _random_circle(shape)
     interpolations = ('nearest', 'linear')
     output_sizes = (None, min(shape), max(shape), 97)
-
     for interpolation, output_size in itertools.product(interpolations,
                                                         output_sizes):
-        print('interpolation =', interpolation)
-        print('output_size =', output_size)
-        # Forward and inverse radon on synthetic data
-        sinogram_rectangle = radon(image, circle=False)
-        reconstruction_rectangle = iradon(sinogram_rectangle,
-                                          output_size=output_size,
-                                          interpolation=interpolation,
-                                          circle=False)
-        sinogram_circle = radon(image, circle=True)
-        reconstruction_circle = iradon(sinogram_circle,
-                                       output_size=output_size,
-                                       interpolation=interpolation,
-                                       circle=True)
-        # Crop rectangular reconstruction to match circle=True reconstruction
-        width = reconstruction_circle.shape[0]
-        excess = int(np.ceil((reconstruction_rectangle.shape[0] - width) / 2))
-        s = np.s_[excess:width + excess, excess:width + excess]
-        reconstruction_rectangle = reconstruction_rectangle[s]
-        # Find the reconstruction circle, set reconstruction to zero outside
-        c0, c1 = np.ogrid[0:width, 0:width]
-        r = np.sqrt((c0 - width // 2)**2 + (c1 - width // 2)**2)
-        reconstruction_rectangle[r >= radius] = 0.
-        print(reconstruction_circle.shape)
-        print(reconstruction_rectangle.shape)
-        np.allclose(reconstruction_rectangle, reconstruction_circle)
+        yield check_radon_iradon_circle, interpolation, shape, output_size
 
 
 if __name__ == "__main__":

From 131cfc73edf64ec5e25ef946cb677625356e5658 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Jostein=20B=C3=B8=20Fl=C3=B8ystad?=
 <jostein.floystad@gmail.com>
Date: Tue, 18 Jun 2013 22:22:22 +0200
Subject: [PATCH 08/21] transform.iradon: Redefine slice and projection center.

These changes should match those made to radon.
---
 skimage/transform/radon_transform.py | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/skimage/transform/radon_transform.py b/skimage/transform/radon_transform.py
index 2733e243..4bd03fdb 100644
--- a/skimage/transform/radon_transform.py
+++ b/skimage/transform/radon_transform.py
@@ -179,7 +179,7 @@ def iradon(radon_image, theta=None, output_size=None,
     if circle:
         radon_size = int(np.ceil(np.sqrt(2) * radon_image.shape[0]))
         radon_image_padded = np.zeros((radon_size, radon_image.shape[1]))
-        radon_pad = (radon_size - radon_image.shape[0]) // 2
+        radon_pad = int(np.floor((radon_size - radon_image.shape[0] - 1) / 2.))
         radon_image_padded[radon_pad:radon_pad + radon_image.shape[0], :] \
             = radon_image
         radon_image = radon_image_padded
@@ -221,7 +221,7 @@ def iradon(radon_image, theta=None, output_size=None,
     # resize filtered image back to original size
     radon_filtered = radon_filtered[:radon_image.shape[0], :]
     reconstructed = np.zeros((output_size, output_size))
-    mid_index = np.ceil(n / 2.0)
+    mid_index = n // 2 + 1
 
     x = output_size
     y = output_size
@@ -236,7 +236,7 @@ def iradon(radon_image, theta=None, output_size=None,
     # reconstruct image by interpolation
     if interpolation == "nearest":
         for i in range(len(theta)):
-            k = np.round(mid_index + xpr * np.sin(th[i]) - ypr * np.cos(th[i]))
+            k = np.round(mid_index + ypr * np.cos(th[i]) - xpr * np.sin(th[i]))
             backprojected = radon_filtered[
                 ((((k > 0) & (k < n)) * k) - 1).astype(np.int), i]
             if circle:
@@ -245,7 +245,7 @@ def iradon(radon_image, theta=None, output_size=None,
 
     elif interpolation == "linear":
         for i in range(len(theta)):
-            t = xpr * np.sin(th[i]) - ypr * np.cos(th[i])
+            t = ypr * np.cos(th[i]) - xpr * np.sin(th[i])
             a = np.floor(t)
             b = mid_index + a
             b0 = ((((b + 1 > 0) & (b + 1 < n)) * (b + 1)) - 1).astype(np.int)

From da423931b5f10bb11d6662834833cc5a7188607a Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Jostein=20B=C3=B8=20Fl=C3=B8ystad?=
 <jostein.floystad@gmail.com>
Date: Tue, 18 Jun 2013 22:30:09 +0200
Subject: [PATCH 09/21] test_radon_transform: Add helper functions.

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

diff --git a/skimage/transform/tests/test_radon_transform.py b/skimage/transform/tests/test_radon_transform.py
index 2f4884ea..7884ef83 100644
--- a/skimage/transform/tests/test_radon_transform.py
+++ b/skimage/transform/tests/test_radon_transform.py
@@ -6,6 +6,35 @@ from numpy.testing import *
 from numpy.fft import ifftshift, ifftn
 import itertools
 from skimage.transform import *
+from skimage.io import imread
+from skimage import data_dir
+
+
+__PHANTOM = imread(data_dir + "/phantom.png", as_grey=True)[::2, ::2]
+
+def _get_phantom():
+    return __PHANTOM
+
+
+def _debug_plot(original, result, sinogram=None):
+    from matplotlib import pyplot as plt
+    imkwargs = dict(cmap='gray', interpolation='nearest')
+    if sinogram is None:
+        plt.figure(figsize=(15, 6))
+        sp = 130
+    else:
+        plt.figure(figsize=(11, 11))
+        sp = 221
+        plt.subplot(sp + 0)
+        plt.imshow(sinogram, aspect='auto', **imkwargs)
+    plt.subplot(sp + 1)
+    plt.imshow(original, **imkwargs)
+    plt.subplot(sp + 2)
+    plt.imshow(result, vmin=original.min(), vmax=original.max(), **imkwargs)
+    plt.subplot(sp + 3)
+    plt.imshow(result - original, **imkwargs)
+    plt.colorbar()
+    plt.show()
 
 
 def rescale(x):

From b61ff7513efd53d1fdc330b8b84ee75777a8b0cb Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Jostein=20B=C3=B8=20Fl=C3=B8ystad?=
 <jostein.floystad@gmail.com>
Date: Wed, 19 Jun 2013 00:12:14 +0200
Subject: [PATCH 10/21] transform.iradon: Refactoring for shorter functions.

Will facilitate testing.
---
 skimage/transform/radon_transform.py | 15 +++++++++------
 1 file changed, 9 insertions(+), 6 deletions(-)

diff --git a/skimage/transform/radon_transform.py b/skimage/transform/radon_transform.py
index 4bd03fdb..8ed11cc7 100644
--- a/skimage/transform/radon_transform.py
+++ b/skimage/transform/radon_transform.py
@@ -113,6 +113,14 @@ def radon(image, theta=None, circle=False):
     return out
 
 
+def _sinogram_circle_to_square(sinogram):
+    size = int(np.ceil(np.sqrt(2) * sinogram.shape[0]))
+    sinogram_padded = np.zeros((size, sinogram.shape[1]))
+    pad = int(np.ceil((size - sinogram.shape[0] - 1) / 2))
+    sinogram_padded[pad:pad + sinogram.shape[0], :] = sinogram
+    return sinogram_padded
+
+
 def iradon(radon_image, theta=None, output_size=None,
            filter="ramp", interpolation="linear", circle=False):
     """
@@ -177,12 +185,7 @@ def iradon(radon_image, theta=None, output_size=None,
             output_size = int(np.floor(np.sqrt((radon_image.shape[0])**2
                                                / 2.0)))
     if circle:
-        radon_size = int(np.ceil(np.sqrt(2) * radon_image.shape[0]))
-        radon_image_padded = np.zeros((radon_size, radon_image.shape[1]))
-        radon_pad = int(np.floor((radon_size - radon_image.shape[0] - 1) / 2.))
-        radon_image_padded[radon_pad:radon_pad + radon_image.shape[0], :] \
-            = radon_image
-        radon_image = radon_image_padded
+        radon_image = _sinogram_circle_to_square(radon_image)
 
     n = radon_image.shape[0]
 

From 1caafd44516714cf5c4fd21b920a345036c5bab8 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Jostein=20B=C3=B8=20Fl=C3=B8ystad?=
 <jostein.floystad@gmail.com>
Date: Wed, 19 Jun 2013 00:13:20 +0200
Subject: [PATCH 11/21] test_radon_transform: Test sinogram conversions.

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

diff --git a/skimage/transform/tests/test_radon_transform.py b/skimage/transform/tests/test_radon_transform.py
index 7884ef83..5a7ca963 100644
--- a/skimage/transform/tests/test_radon_transform.py
+++ b/skimage/transform/tests/test_radon_transform.py
@@ -243,6 +243,23 @@ def test_radon_circle():
     assert np.all(relative_error < 3e-3)
 
 
+def check_sinogram_circle_to_square(size):
+    from skimage.transform.radon_transform import _sinogram_circle_to_square
+    image = _random_circle((size, size))
+    theta = np.linspace(0., 180., size, False)
+    sinogram_circle = radon(image, theta, circle=True)
+    sinogram_square = radon(image, theta, circle=False)
+    sinogram_circle_to_square = _sinogram_circle_to_square(sinogram_circle)
+    error = abs(sinogram_square - sinogram_circle_to_square)
+    print(np.mean(error), np.max(error))
+    assert np.allclose(sinogram_square, sinogram_circle_to_square)
+
+
+def test_sinogram_circle_to_square():
+    for size in (50, 51):
+        yield check_sinogram_circle_to_square, size
+
+
 def check_radon_iradon_circle(interpolation, shape, output_size):
     # Forward and inverse radon on synthetic data
     image = _random_circle(shape)

From cf51de6b37b8ed535fd97a45afe102e5128d77b6 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Jostein=20B=C3=B8=20Fl=C3=B8ystad?=
 <jostein.floystad@gmail.com>
Date: Wed, 19 Jun 2013 00:17:58 +0200
Subject: [PATCH 12/21] test_radon_transform: Refactor and improve
 test_radon_iradon.

Aside from refactoring, the Shepp-Logan phantom is now used as it is a
more challenging test object.
---
 .../transform/tests/test_radon_transform.py   | 48 +++++++++----------
 1 file changed, 22 insertions(+), 26 deletions(-)

diff --git a/skimage/transform/tests/test_radon_transform.py b/skimage/transform/tests/test_radon_transform.py
index 5a7ca963..bc938da8 100644
--- a/skimage/transform/tests/test_radon_transform.py
+++ b/skimage/transform/tests/test_radon_transform.py
@@ -113,35 +113,31 @@ def test_iradon_center():
         yield check_iradon_center, size, theta, circle
 
 
-def test_radon_iradon():
-    size = 100
+def check_radon_iradon(interpolation_type, filter_type):
     debug = False
-    image = np.tri(size) + np.tri(size)[::-1]
-    for filter_type in ["ramp", "shepp-logan", "cosine", "hamming", "hann"]:
-        reconstructed = iradon(radon(image), filter=filter_type)
+    image = _get_phantom()
+    reconstructed = iradon(radon(image), filter=filter_type,
+                           interpolation=interpolation_type)
+    delta = np.mean(np.abs(image - reconstructed))
+    print('\n\tmean error:', delta)
+    if debug:
+        _debug_plot(image, reconstructed)
+    if filter_type == 'ramp':
+        if interpolation_type == 'linear':
+            allowed_delta = 0.02
+        else:
+            allowed_delta = 0.03
+    else:
+        allowed_delta = 0.05
+    assert delta < allowed_delta
 
-        image = rescale(image)
-        reconstructed = rescale(reconstructed)
-        delta = np.mean(np.abs(image - reconstructed))
 
-        if debug:
-            print(delta)
-            import matplotlib.pyplot as plt
-            f, (ax1, ax2) = plt.subplots(1, 2)
-            ax1.imshow(image, cmap=plt.cm.gray)
-            ax2.imshow(reconstructed, cmap=plt.cm.gray)
-            plt.show()
-
-        assert delta < 0.05
-
-    reconstructed = iradon(radon(image), filter="ramp",
-                           interpolation="nearest")
-    delta = np.mean(abs(image - reconstructed))
-    assert delta < 0.05
-    size = 20
-    image = np.tri(size) + np.tri(size)[::-1]
-    reconstructed = iradon(radon(image), filter="ramp",
-                           interpolation="nearest")
+def test_radon_iradon():
+    filter_types = ["ramp", "shepp-logan", "cosine", "hamming", "hann"]
+    interpolation_types = ["linear", "nearest"]
+    for interpolation_type, filter_type in \
+            itertools.product(interpolation_types, filter_types):
+        yield check_radon_iradon, interpolation_type, filter_type
 
 
 def test_iradon_angles():

From e63a1fb341052019991b62aff0ea5b1cf414b4cb Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Jostein=20B=C3=B8=20Fl=C3=B8ystad?=
 <jostein.floystad@gmail.com>
Date: Wed, 19 Jun 2013 00:19:31 +0200
Subject: [PATCH 13/21] test_radon_transform: debug option for
 test_iradon_minimal.

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

diff --git a/skimage/transform/tests/test_radon_transform.py b/skimage/transform/tests/test_radon_transform.py
index bc938da8..1c580a4a 100644
--- a/skimage/transform/tests/test_radon_transform.py
+++ b/skimage/transform/tests/test_radon_transform.py
@@ -172,6 +172,7 @@ def test_radon_minimal():
     """
     Test for small images for various angles
     """
+    debug = False
     thetas = [np.arange(180)]
     for theta in thetas:
         a = np.zeros((3, 3))
@@ -179,6 +180,8 @@ def test_radon_minimal():
         p = radon(a, theta)
         reconstructed = iradon(p, theta)
         reconstructed /= np.max(reconstructed)
+        if debug:
+            _debug_plot(a, reconstructed)
         assert np.all(abs(a - reconstructed) < 0.4)
 
         b = np.zeros((4, 4))
@@ -186,6 +189,8 @@ def test_radon_minimal():
         p = radon(b, theta)
         reconstructed = iradon(p, theta)
         reconstructed /= np.max(reconstructed)
+        if debug:
+            _debug_plot(b, reconstructed)
         assert np.all(abs(b - reconstructed) < 0.4)
 
         c = np.zeros((5, 5))
@@ -193,6 +198,8 @@ def test_radon_minimal():
         p = radon(c, theta)
         reconstructed = iradon(p, theta)
         reconstructed /= np.max(reconstructed)
+        if debug:
+            _debug_plot(c, reconstructed)
         assert np.all(abs(c - reconstructed) < 0.4)
 
 

From b8a6b4fa00249646212569ffbdb93134938e92f2 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Jostein=20B=C3=B8=20Fl=C3=B8ystad?=
 <jostein.floystad@gmail.com>
Date: Wed, 19 Jun 2013 00:23:55 +0200
Subject: [PATCH 14/21] radon_transform: Stylistic changes.

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

diff --git a/skimage/transform/radon_transform.py b/skimage/transform/radon_transform.py
index 8ed11cc7..2cde2232 100644
--- a/skimage/transform/radon_transform.py
+++ b/skimage/transform/radon_transform.py
@@ -91,25 +91,20 @@ def radon(image, theta=None, circle=False):
     shift0 = np.array([[1, 0, -dw],
                        [0, 1, -dh],
                        [0, 0, 1]])
-
     shift1 = np.array([[1, 0, dw],
                        [0, 1, dh],
                        [0, 0, 1]])
 
     def build_rotation(theta):
         T = np.deg2rad(theta)
-
         R = np.array([[np.cos(T), np.sin(T), 0],
                       [-np.sin(T), np.cos(T), 0],
                       [0, 0, 1]])
-
         return shift1.dot(R).dot(shift0)
 
     for i in range(len(theta)):
         rotated = _warp_fast(padded_image, build_rotation(theta[i]))
-
         out[:, i] = rotated.sum(0)
-
     return out
 
 
@@ -158,27 +153,23 @@ def iradon(radon_image, theta=None, output_size=None,
 
     Notes
     -----
-    It applies the fourier slice theorem to reconstruct an image by
+    It applies the Fourier slice theorem to reconstruct an image by
     multiplying the frequency domain of the filter with the FFT of the
     projection data. This algorithm is called filtered back projection.
 
     """
     if radon_image.ndim != 2:
         raise ValueError('The input image must be 2-D')
-
     if theta is None:
         m, n = radon_image.shape
         theta = np.linspace(0, 180, n, endpoint=False)
     else:
         theta = np.asarray(theta)
-
     if len(theta) != radon_image.shape[1]:
         raise ValueError("The given ``theta`` does not match the number of "
                          "projections in ``radon_image``.")
-
-    th = (np.pi / 180.0) * theta
-    # if output size not specified, estimate from input radon image
     if not output_size:
+        # If output size not specified, estimate from input radon image
         if circle:
             output_size = radon_image.shape[0]
         else:
@@ -187,19 +178,19 @@ def iradon(radon_image, theta=None, output_size=None,
     if circle:
         radon_image = _sinogram_circle_to_square(radon_image)
 
+    th = (np.pi / 180.0) * theta
     n = radon_image.shape[0]
-
     img = radon_image.copy()
     # resize image to next power of two for fourier analysis
     # speeds up fourier and lessens artifacts
     order = max(64., 2**np.ceil(np.log(2 * n) / np.log(2)))
     # zero pad input image
     img.resize((order, img.shape[1]))
-    # construct the fourier filter
 
+    # Construct the Fourier filter
     f = fftshift(abs(np.mgrid[-1:1:2 / order])).reshape(-1, 1)
     w = 2 * np.pi * f
-    # start from first element to avoid divide by zero
+    # Start from first element to avoid divide by zero
     if filter == "ramp":
         pass
     elif filter == "shepp-logan":
@@ -214,14 +205,12 @@ def iradon(radon_image, theta=None, output_size=None,
         f[1:] = 1
     else:
         raise ValueError("Unknown filter: %s" % filter)
-
     filter_ft = np.tile(f, (1, len(theta)))
-
-    # apply filter in fourier domain
+    # Apply filter in Fourier domain
     projection = fft(img, axis=0) * filter_ft
     radon_filtered = np.real(ifft(projection, axis=0))
 
-    # resize filtered image back to original size
+    # Resize filtered image back to original size
     radon_filtered = radon_filtered[:radon_image.shape[0], :]
     reconstructed = np.zeros((output_size, output_size))
     mid_index = n // 2 + 1
@@ -231,12 +220,11 @@ def iradon(radon_image, theta=None, output_size=None,
     [X, Y] = np.mgrid[0.0:x, 0.0:y]
     xpr = X - int(output_size) // 2
     ypr = Y - int(output_size) // 2
-
     if circle:
         radius = (output_size - 1) // 2
         reconstruction_circle = (xpr**2 + ypr**2) < radius**2
 
-    # reconstruct image by interpolation
+    # Reconstruct image by interpolation
     if interpolation == "nearest":
         for i in range(len(theta)):
             k = np.round(mid_index + ypr * np.cos(th[i]) - xpr * np.sin(th[i]))
@@ -245,7 +233,6 @@ def iradon(radon_image, theta=None, output_size=None,
             if circle:
                 backprojected[~reconstruction_circle] = 0.
             reconstructed += backprojected
-
     elif interpolation == "linear":
         for i in range(len(theta)):
             t = ypr * np.cos(th[i]) - xpr * np.sin(th[i])

From afaab4fea7a6bc766a06a8d9df296f3f6ec63e49 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Jostein=20B=C3=B8=20Fl=C3=B8ystad?=
 <jostein.floystad@gmail.com>
Date: Wed, 19 Jun 2013 00:51:25 +0200
Subject: [PATCH 15/21] test_radon_transform: Refactor tests and make them
 stricter.

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

diff --git a/skimage/transform/tests/test_radon_transform.py b/skimage/transform/tests/test_radon_transform.py
index 1c580a4a..d2e0a490 100644
--- a/skimage/transform/tests/test_radon_transform.py
+++ b/skimage/transform/tests/test_radon_transform.py
@@ -168,39 +168,26 @@ def test_iradon_angles():
     assert delta_80 > delta_200
 
 
-def test_radon_minimal():
-    """
-    Test for small images for various angles
-    """
+def check_radon_iradon_minimal(shape, slices):
     debug = False
-    thetas = [np.arange(180)]
-    for theta in thetas:
-        a = np.zeros((3, 3))
-        a[1, 1] = 1
-        p = radon(a, theta)
-        reconstructed = iradon(p, theta)
-        reconstructed /= np.max(reconstructed)
-        if debug:
-            _debug_plot(a, reconstructed)
-        assert np.all(abs(a - reconstructed) < 0.4)
+    theta = np.arange(180)
+    image = np.zeros(shape, dtype=np.float)
+    image[slices] = 1.
+    sinogram = radon(image, theta)
+    reconstructed = iradon(sinogram, theta)
+    print('\n\tMaximum deviation:', np.max(np.abs(image - reconstructed)))
+    if debug:
+        _debug_plot(image, reconstructed)
+    assert np.allclose(image, reconstructed, rtol=1e-1, atol=1e-1)
 
-        b = np.zeros((4, 4))
-        b[1:3, 1:3] = 1
-        p = radon(b, theta)
-        reconstructed = iradon(p, theta)
-        reconstructed /= np.max(reconstructed)
-        if debug:
-            _debug_plot(b, reconstructed)
-        assert np.all(abs(b - reconstructed) < 0.4)
 
-        c = np.zeros((5, 5))
-        c[1:3, 1:3] = 1
-        p = radon(c, theta)
-        reconstructed = iradon(p, theta)
-        reconstructed /= np.max(reconstructed)
-        if debug:
-            _debug_plot(c, reconstructed)
-        assert np.all(abs(c - reconstructed) < 0.4)
+def test_radon_iradon_minimal():
+    # Each testcase is a (image shape, slice of image to set to one) tuple
+    testcases = [((3, 3), np.s_[1, 1]),
+                 ((4, 4), np.s_[1:3, 1:3]),
+                 ((5, 5), np.s_[1:3, 1:3])]
+    for shape, slices in testcases:
+        yield check_radon_iradon_minimal, shape, slices
 
 
 def test_reconstruct_with_wrong_angles():

From 934f1040ad67466c5ede195a01e656eb31bf5699 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Jostein=20B=C3=B8=20Fl=C3=B8ystad?=
 <jostein.floystad@gmail.com>
Date: Wed, 19 Jun 2013 11:31:41 +0200
Subject: [PATCH 16/21] test_radon_transform: Change test criterion for
 iradon_minimal.

---
 .../transform/tests/test_radon_transform.py   | 19 +++++++++++--------
 1 file changed, 11 insertions(+), 8 deletions(-)

diff --git a/skimage/transform/tests/test_radon_transform.py b/skimage/transform/tests/test_radon_transform.py
index d2e0a490..7f517801 100644
--- a/skimage/transform/tests/test_radon_transform.py
+++ b/skimage/transform/tests/test_radon_transform.py
@@ -177,17 +177,20 @@ def check_radon_iradon_minimal(shape, slices):
     reconstructed = iradon(sinogram, theta)
     print('\n\tMaximum deviation:', np.max(np.abs(image - reconstructed)))
     if debug:
-        _debug_plot(image, reconstructed)
-    assert np.allclose(image, reconstructed, rtol=1e-1, atol=1e-1)
+        _debug_plot(image, reconstructed, sinogram)
+    if image.sum() == 1:
+        assert (np.unravel_index(np.argmax(reconstructed), image.shape)
+                == np.unravel_index(np.argmax(image), image.shape))
 
 
 def test_radon_iradon_minimal():
-    # Each testcase is a (image shape, slice of image to set to one) tuple
-    testcases = [((3, 3), np.s_[1, 1]),
-                 ((4, 4), np.s_[1:3, 1:3]),
-                 ((5, 5), np.s_[1:3, 1:3])]
-    for shape, slices in testcases:
-        yield check_radon_iradon_minimal, shape, slices
+    shapes = [(3, 3), (4, 4), (5, 5)]
+    for shape in shapes:
+        c0, c1 = shape[0] // 2, shape[1] // 2
+        coordinates = itertools.product((c0 - 1, c0, c0 + 1),
+                                        (c1 - 1, c1, c1 + 1))
+        for coordinate in coordinates:
+            yield check_radon_iradon_minimal, shape, coordinate
 
 
 def test_reconstruct_with_wrong_angles():

From df0b060c696da4fe03322b40931833a3beabb5ca Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Jostein=20B=C3=B8=20Fl=C3=B8ystad?=
 <jostein.floystad@gmail.com>
Date: Wed, 19 Jun 2013 11:32:38 +0200
Subject: [PATCH 17/21] test_radon_transform: Change test criterion for
 sinogram_circle.

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

diff --git a/skimage/transform/tests/test_radon_transform.py b/skimage/transform/tests/test_radon_transform.py
index 7f517801..5eef1e1c 100644
--- a/skimage/transform/tests/test_radon_transform.py
+++ b/skimage/transform/tests/test_radon_transform.py
@@ -241,11 +241,17 @@ def check_sinogram_circle_to_square(size):
     image = _random_circle((size, size))
     theta = np.linspace(0., 180., size, False)
     sinogram_circle = radon(image, theta, circle=True)
+    argmax_shape = lambda a: np.unravel_index(np.argmax(a), a.shape)
+    print('\n\targmax of circle:', argmax_shape(sinogram_circle))
     sinogram_square = radon(image, theta, circle=False)
+    print('\targmax of square:', argmax_shape(sinogram_square))
     sinogram_circle_to_square = _sinogram_circle_to_square(sinogram_circle)
+    print('\targmax of circle to square:',
+          argmax_shape(sinogram_circle_to_square))
     error = abs(sinogram_square - sinogram_circle_to_square)
     print(np.mean(error), np.max(error))
-    assert np.allclose(sinogram_square, sinogram_circle_to_square)
+    assert (argmax_shape(sinogram_square)
+            == argmax_shape(sinogram_circle_to_square))
 
 
 def test_sinogram_circle_to_square():

From cca66a04ef50b877d219b1f719f6f16c3dfe3891 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Jostein=20B=C3=B8=20Fl=C3=B8ystad?=
 <jostein.floystad@gmail.com>
Date: Wed, 19 Jun 2013 12:25:39 +0200
Subject: [PATCH 18/21] transform.radon: Robust determination of center of
 projection.

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

diff --git a/skimage/transform/radon_transform.py b/skimage/transform/radon_transform.py
index 2cde2232..6bd0f628 100644
--- a/skimage/transform/radon_transform.py
+++ b/skimage/transform/radon_transform.py
@@ -72,6 +72,8 @@ def radon(image, theta=None, circle=False):
         slices = tuple(slices)
         padded_image = image[slices]
         out = np.zeros((min(padded_image.shape), len(theta)))
+        dh = padded_image.shape[0] // 2
+        dw = padded_image.shape[1] // 2
     else:
         height, width = image.shape
         diagonal = np.sqrt(2) * max(image.shape)
@@ -85,9 +87,9 @@ def radon(image, theta=None, circle=False):
         x1 = x0 + width
         padded_image[y0:y1, x0:x1] = image
         out = np.zeros((max(padded_image.shape), len(theta)))
+        dh = y0 + height // 2
+        dw = x0 + width // 2
 
-    h, w = padded_image.shape
-    dh, dw = h // 2, w // 2
     shift0 = np.array([[1, 0, -dw],
                        [0, 1, -dh],
                        [0, 0, 1]])
@@ -111,7 +113,7 @@ def radon(image, theta=None, circle=False):
 def _sinogram_circle_to_square(sinogram):
     size = int(np.ceil(np.sqrt(2) * sinogram.shape[0]))
     sinogram_padded = np.zeros((size, sinogram.shape[1]))
-    pad = int(np.ceil((size - sinogram.shape[0] - 1) / 2))
+    pad = (size - sinogram.shape[0]) // 2
     sinogram_padded[pad:pad + sinogram.shape[0], :] = sinogram
     return sinogram_padded
 

From 4b25c482452da3c04c5803e48f64bb8813d99174 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Jostein=20B=C3=B8=20Fl=C3=B8ystad?=
 <jostein.floystad@gmail.com>
Date: Wed, 19 Jun 2013 12:26:47 +0200
Subject: [PATCH 19/21] transform.iradon: Correct determination of center of
 projection.

---
 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 6bd0f628..b114adad 100644
--- a/skimage/transform/radon_transform.py
+++ b/skimage/transform/radon_transform.py
@@ -215,7 +215,10 @@ def iradon(radon_image, theta=None, output_size=None,
     # Resize filtered image back to original size
     radon_filtered = radon_filtered[:radon_image.shape[0], :]
     reconstructed = np.zeros((output_size, output_size))
-    mid_index = n // 2 + 1
+    # Determine the center of the projections (= center of sinogram)
+    circle_size = int(np.floor(radon_image.shape[0] / np.sqrt(2)))
+    square_size = radon_image.shape[0]
+    mid_index = (square_size - circle_size) // 2 + circle_size // 2 + 1
 
     x = output_size
     y = output_size

From 186a238e4814fcb86bc8351f59a52e50f8a92583 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Jostein=20B=C3=B8=20Fl=C3=B8ystad?=
 <jostein.floystad@gmail.com>
Date: Wed, 19 Jun 2013 12:37:46 +0200
Subject: [PATCH 20/21] test_radon_transform: Clean up imports.

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

diff --git a/skimage/transform/tests/test_radon_transform.py b/skimage/transform/tests/test_radon_transform.py
index 5eef1e1c..c385c84c 100644
--- a/skimage/transform/tests/test_radon_transform.py
+++ b/skimage/transform/tests/test_radon_transform.py
@@ -2,10 +2,9 @@ from __future__ import print_function
 from __future__ import division
 
 import numpy as np
-from numpy.testing import *
-from numpy.fft import ifftshift, ifftn
+from numpy.testing import assert_raises
 import itertools
-from skimage.transform import *
+from skimage.transform import radon, iradon
 from skimage.io import imread
 from skimage import data_dir
 

From 6d2f082c1196aae972e85319a67fe9cc60537fce Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Jostein=20B=C3=B8=20Fl=C3=B8ystad?=
 <jostein.floystad@gmail.com>
Date: Wed, 19 Jun 2013 12:38:50 +0200
Subject: [PATCH 21/21] test_radon_transform: Style fixes, PEP8.

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

diff --git a/skimage/transform/tests/test_radon_transform.py b/skimage/transform/tests/test_radon_transform.py
index c385c84c..f7127aa5 100644
--- a/skimage/transform/tests/test_radon_transform.py
+++ b/skimage/transform/tests/test_radon_transform.py
@@ -11,6 +11,7 @@ from skimage import data_dir
 
 __PHANTOM = imread(data_dir + "/phantom.png", as_grey=True)[::2, ::2]
 
+
 def _get_phantom():
     return __PHANTOM
 
@@ -264,14 +265,14 @@ def check_radon_iradon_circle(interpolation, shape, output_size):
     radius = min(shape) // 2
     sinogram_rectangle = radon(image, circle=False)
     reconstruction_rectangle = iradon(sinogram_rectangle,
-                                        output_size=output_size,
-                                        interpolation=interpolation,
-                                        circle=False)
+                                      output_size=output_size,
+                                      interpolation=interpolation,
+                                      circle=False)
     sinogram_circle = radon(image, circle=True)
     reconstruction_circle = iradon(sinogram_circle,
-                                    output_size=output_size,
-                                    interpolation=interpolation,
-                                    circle=True)
+                                   output_size=output_size,
+                                   interpolation=interpolation,
+                                   circle=True)
     # Crop rectangular reconstruction to match circle=True reconstruction
     width = reconstruction_circle.shape[0]
     excess = int(np.ceil((reconstruction_rectangle.shape[0] - width) / 2))