From be28bb9fba1162322a3e49e0ce9c01d6f277e1d2 Mon Sep 17 00:00:00 2001
From: Pieter Holtzhausen <holtzhau@gmail.com>
Date: Wed, 17 Aug 2011 17:46:13 +0200
Subject: [PATCH 01/17] First working version of radon and iradon

---
 scikits/image/transform/__init__.py        |   1 +
 scikits/image/transform/radon_transform.py | 140 +++++++++++++++++++++
 2 files changed, 141 insertions(+)
 create mode 100644 scikits/image/transform/radon_transform.py

diff --git a/scikits/image/transform/__init__.py b/scikits/image/transform/__init__.py
index 91711c6e..d8adede9 100644
--- a/scikits/image/transform/__init__.py
+++ b/scikits/image/transform/__init__.py
@@ -1,4 +1,5 @@
 from hough_transform import *
 from finite_radon_transform import *
+from radon_transform import *
 from project import *
 
diff --git a/scikits/image/transform/radon_transform.py b/scikits/image/transform/radon_transform.py
new file mode 100644
index 00000000..b3aab162
--- /dev/null
+++ b/scikits/image/transform/radon_transform.py
@@ -0,0 +1,140 @@
+import numpy as np
+from scipy.misc import imrotate
+from scipy.interpolate import interp1d
+from scipy.fftpack import fftshift, ifftshift, fft, ifft
+import math
+
+def radon(image, theta=None):
+    """
+    Calculates the projections given the current object and projection angle
+    Justin K. Romberg
+    """
+    if theta == None:
+        theta = np.arange(180)
+    height, width = image.shape
+    diagonal = np.sqrt(height**2 + width**2)
+    heightpad = np.ceil(diagonal - height) + 2
+    widthpad = np.ceil(diagonal - width) + 2
+    padded_image = np.zeros((int(height+heightpad), int(width+widthpad)))
+    y0, y1 = int(np.ceil(heightpad/2)), int((np.ceil(heightpad/2)+height))
+    x0, x1 = int((np.ceil(widthpad/2))), int((np.ceil(widthpad/2)+width))
+    padded_image[y0:y1, x0:x1] = image
+    out = np.zeros((max(padded_image.shape), len(theta)))
+    for i in range(len(theta)):
+        rotated = imrotate(padded_image, -theta[i])
+        out[:,i] = rotated.sum(0)[::-1]
+    return out
+
+"""
+  if 0:
+        # filter the projections
+        freqs = np.zeros((n, 1))
+        freqs[:, 0] = np.linspace(-1, 1, n).T;
+        filter_ft = np.tile(np.abs(freqs), (1, len(theta)))
+        # fourier domain filtering
+        radon_ft = fft(radon_image, axis=0)
+        projection = radon_ft * fftshift(filter_ft)
+        radon_filtered = np.real(ifft(projection, axis=0))        
+    #    print np.max(projection)
+    #    print projection
+        #projection = ifftshift(projection, axes=1);
+    if 0:
+        height, width = radon_image.shape
+        w = np.mgrid[-math.pi:math.pi:(2*math.pi)/height]
+        f = fftshift(abs(w))
+        g = np.array([np.real(ifft(fft(i)*f)) for i in radon_image.T])
+        radon_filtered = np.transpose(g)
+    if 0:
+        img = radon_image.copy()
+        order = 1024
+        filt = np.zeros((order/2, 1))
+        filt[:, 0] = 2.0*np.arange(0, order/2)/order;
+        filt = np.vstack((filt, filt[ ::-1])).T
+        #filt = fftshift(abs(filt))
+    #    order = radon_image.shape[0]
+        w = np.mgrid[-math.pi:math.pi:(2*math.pi)/order]
+        filt = fftshift(abs(w))
+        img.resize((order, img.shape[1]))
+        radon_filtered = np.array([np.real(ifft(fft(column)*filt)) for column in img.T]).T
+        radon_filtered = radon_filtered[:radon_image.shape[0], :]
+    if 0:
+        ### bestest
+        img = radon_image.copy()
+        order = max(64, 2 ** np.ceil(np.log(2*n)/np.log(2)))
+#        filt = np.zeros((order/2, 1))
+#        filt[:, 0] = 2.0*np.arange(0, order/2)/order;
+#        filt = np.vstack((filt, filt[ ::-1])).T
+        #filt = fftshift(abs(filt))
+    #    order = radon_image.shape[0]
+        w = np.mgrid[-math.pi:math.pi:(2*math.pi)/order]
+        filt = fftshift(abs(w))
+        img.resize((order, img.shape[1]))
+        img = fft(img, axis=0)
+        #radon_filtered = np.array([np.real(ifft(column*filt)) for column in img.T]).T
+        radon_filtered = np.array([column*filt for column in img.T]).T
+        
+        radon_filtered = np.real(ifft(radon_filtered, axis=0))
+        radon_filtered = radon_filtered[:radon_image.shape[0], :]
+"""
+
+def iradon(radon_image, theta=None, output_size=None, filter="ramp", interpolate="nearest"):
+    if theta == None:
+        theta = np.mgrid[0:180]
+    th = (math.pi/180.0)*theta        
+    # if output size not specified, estimate from input radon image
+    if not output_size:
+        output_size = 2*np.floor(radon_image.shape[0] / (2*np.sqrt(2)))
+    n = radon_image.shape[0]
+  
+    img = radon_image.copy()
+    # resize image to next power of two for fourier analysis
+    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
+    freqs = np.zeros((order, 1))
+    
+    #w = np.sqrt(np.sum((np.mgrid[-pi:pi:(2*pi)/Length1, -pi:pi:(2*pi)/Length2])**2, 0))
+
+    w = fftshift(abs(np.mgrid[-1:1:2/order])).reshape(-1, 1)
+#    if filter == "ramp":
+#    elif filter == "shepp-logan":
+#        rn1 = abs(2/a*s.sin(a*w/2))
+#        rn2 = s.sin(a*w/2)
+#        rd = (a*w)/2
+#        r = rn1*(rn2/rd)**2
+#        r = where(w!=0, r, w!=0)
+#        f = fftshift(r)
+#    elif filter == 'cosine':
+#    elif filter == 'hamming':
+#    elif filter == 'hann':
+#    elif filter == None:
+        
+        
+    filter_ft = np.tile(w, (1, len(theta)))        
+    # 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
+    radon_filtered = radon_filtered[:radon_image.shape[0], :]
+    reconstructed = np.zeros((output_size, output_size))
+    midindex = (n + 1.0) / 2.0
+    x = output_size
+    y = output_size
+    [X, Y] = np.mgrid[0.0:x, 0.0:y]
+    xpr = X - (output_size+1.0)/2.0
+    ypr = Y - (output_size+1.0)/2.0
+    if interpolate == "nearest":        
+        for i in range(len(theta)):
+            filtIndex = np.round(midindex + xpr*np.sin(th[i]) - ypr*np.cos(th[i]))
+            reconstructed += radon_filtered[((((filtIndex > 0) & \
+                (filtIndex <= n))*filtIndex) - 1).astype('i'), i]
+    elif interpolate == "linear":
+        pass
+    elif interpolate == "spline":
+        pass
+        
+    return reconstructed * math.pi / (2*len(th))
+
+
+

From e26089bad09088f47a8f8e96327a4a34f1cce5fc Mon Sep 17 00:00:00 2001
From: Pieter Holtzhausen <holtzhau@gmail.com>
Date: Thu, 18 Aug 2011 15:52:59 +0200
Subject: [PATCH 02/17] Implemented various filters and interpolation methods.

---
 scikits/image/transform/radon_transform.py | 72 +++++++++++++---------
 1 file changed, 43 insertions(+), 29 deletions(-)

diff --git a/scikits/image/transform/radon_transform.py b/scikits/image/transform/radon_transform.py
index b3aab162..635e5823 100644
--- a/scikits/image/transform/radon_transform.py
+++ b/scikits/image/transform/radon_transform.py
@@ -77,9 +77,9 @@ def radon(image, theta=None):
         radon_filtered = radon_filtered[:radon_image.shape[0], :]
 """
 
-def iradon(radon_image, theta=None, output_size=None, filter="ramp", interpolate="nearest"):
+def iradon(radon_image, theta=None, output_size=None, filter="ramp", interpolate="linear"):
     if theta == None:
-        theta = np.mgrid[0:180]
+        theta = np.arange(180)
     th = (math.pi/180.0)*theta        
     # if output size not specified, estimate from input radon image
     if not output_size:
@@ -94,45 +94,59 @@ def iradon(radon_image, theta=None, output_size=None, filter="ramp", interpolate
     #construct the fourier filter
     freqs = np.zeros((order, 1))
     
-    #w = np.sqrt(np.sum((np.mgrid[-pi:pi:(2*pi)/Length1, -pi:pi:(2*pi)/Length2])**2, 0))
-
-    w = fftshift(abs(np.mgrid[-1:1:2/order])).reshape(-1, 1)
-#    if filter == "ramp":
-#    elif filter == "shepp-logan":
-#        rn1 = abs(2/a*s.sin(a*w/2))
-#        rn2 = s.sin(a*w/2)
-#        rd = (a*w)/2
-#        r = rn1*(rn2/rd)**2
-#        r = where(w!=0, r, w!=0)
-#        f = fftshift(r)
-#    elif filter == 'cosine':
-#    elif filter == 'hamming':
-#    elif filter == 'hann':
-#    elif filter == None:
+    f = fftshift(abs(np.mgrid[-1:1:2/order])).reshape(-1, 1)
+    w = 2 * math.pi * f
+    # start from first element to avoid divide by zero
+    if filter == "ramp":
+        pass
+    elif filter == "shepp-logan":
+        f[1:] = f[1:] * np.sin(w[1:] / 2) / (w[1:]/2)
+    elif filter == "cosine":
+       f[1:] = f[1:] * np.cos(w[1:] / 2)
+    elif filter == "hamming":
+       f[1:] = f[1:] * (0.54 + 0.46 * np.cos(w[1:]))
+    elif filter == "hann":
+       f[1:] = f[1:] * (1 + np.cos(w[1:])) / 2
+    elif filter == None:
+        f[1:] = 1
+    else:
+        raise ValueError("Unknown filter: %s" % filter)
         
-        
-    filter_ft = np.tile(w, (1, len(theta)))        
+    filter_ft = np.tile(f, (1, len(theta)))        
     # 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
     radon_filtered = radon_filtered[:radon_image.shape[0], :]
     reconstructed = np.zeros((output_size, output_size))
-    midindex = (n + 1.0) / 2.0
+    mid_index = np.ceil(n/2);
     x = output_size
     y = output_size
     [X, Y] = np.mgrid[0.0:x, 0.0:y]
-    xpr = X - (output_size+1.0)/2.0
-    ypr = Y - (output_size+1.0)/2.0
-    if interpolate == "nearest":        
+    xpr = X - (output_size + 1.0) / 2.0
+    ypr = Y - (output_size + 1.0) / 2.0
+    if interpolate == "nearest":   
         for i in range(len(theta)):
-            filtIndex = np.round(midindex + xpr*np.sin(th[i]) - ypr*np.cos(th[i]))
-            reconstructed += radon_filtered[((((filtIndex > 0) & \
-                (filtIndex <= n))*filtIndex) - 1).astype('i'), i]
+            k = np.round(mid_index + xpr*np.sin(th[i]) - ypr*np.cos(th[i]))
+            reconstructed += radon_filtered[((((k > 0) & (k < n))*k) - 1).astype(np.int), i]             
     elif interpolate == "linear":
-        pass
-    elif interpolate == "spline":
-        pass
+        for i in range(len(theta)):
+          t = xpr*np.sin(th[i]) - ypr*np.cos(th[i])
+          a = np.floor(t)
+          b = mid_index + a 
+          reconstructed += (t - a) * radon_filtered[((((b+1 > 0) & (b+1 < n))*(b+1)) - 1).astype(np.int), i] \
+           + (a - t + 1) * radon_filtered[((((b > 0) & (b < n))*b) - 1).astype(np.int), i]
+# XXX slow and inaccurate           
+#    elif interpolate == "spline":
+#        axis = np.arange(0, radon_filtered.shape[0]) - mid_index
+#        for i in range(len(theta)):
+#            print i
+#            t = xpr*np.sin(th[i]) - ypr*np.cos(th[i])
+#            #f = interp1d(axis, radon_filtered[:, i], kind="cubic", bounds_error=False, fill_value=0) 
+#            f = interp1d(axis, radon_filtered[:, i], kind="linear", bounds_error=False, fill_value=0) # cubic
+#            reconstructed += f(t).reshape(output_size, output_size)
+    else:
+        raise ValueError("Unknown interpolation: %s" % interpolate)
         
     return reconstructed * math.pi / (2*len(th))
 

From 494be9f96765643a1f7fb443dd7d1294b78b59e9 Mon Sep 17 00:00:00 2001
From: Pieter Holtzhausen <holtzhau@gmail.com>
Date: Fri, 19 Aug 2011 00:20:14 +0200
Subject: [PATCH 03/17] Test and image dimension checking

---
 scikits/image/transform/radon_transform.py    | 150 ++++++++++--------
 .../transform/tests/test_radon_transform.py   |  20 +++
 2 files changed, 101 insertions(+), 69 deletions(-)
 create mode 100644 scikits/image/transform/tests/test_radon_transform.py

diff --git a/scikits/image/transform/radon_transform.py b/scikits/image/transform/radon_transform.py
index 635e5823..77f84f0d 100644
--- a/scikits/image/transform/radon_transform.py
+++ b/scikits/image/transform/radon_transform.py
@@ -1,16 +1,45 @@
+"""
+radon.py - Radon and inverse radon transforms
+
+Based on code of Justin K. Romberg 
+(http://www.clear.rice.edu/elec431/projects96/DSP/bpanalysis.html)
+J. Gillam and Chris Griffin.
+
+References: 
+    -B.R. Ramesh, N. Srinivasa, K. Rajgopal, "An Algorithm for Computing
+    the Discrete Radon Transform With Some Applications", Proceedings of
+    the Fourth IEEE Region 10 International Conference, TENCON '89, 1989.
+    -A. C. Kak, Malcolm Slaney, "Principles of Computerized Tomographic
+    Imaging", IEEE Press 1988.
+"""
+
 import numpy as np
 from scipy.misc import imrotate
 from scipy.interpolate import interp1d
-from scipy.fftpack import fftshift, ifftshift, fft, ifft
+from scipy.fftpack import fftshift, fft, ifft
 import math
 
+
 def radon(image, theta=None):
     """
-    Calculates the projections given the current object and projection angle
-    Justin K. Romberg
+    Calculates the radon transform of an image given specified projection angles.
+
+    Parameters
+    ----------
+    image : array_like, dtype=float
+        Input image.
+    theta : array_like, dtype=float, optional (default np.arange(180))
+        Projection angles (in degrees).
+    
+    Returns
+    -------
+    output : ndarray
+        Radon transform.
     """
+    if image.ndim != 2:
+        raise ValueError('The input image must be 2-D')
     if theta == None:
-        theta = np.arange(180)
+        theta = np.arange(180)    
     height, width = image.shape
     diagonal = np.sqrt(height**2 + width**2)
     heightpad = np.ceil(diagonal - height) + 2
@@ -25,82 +54,63 @@ def radon(image, theta=None):
         out[:,i] = rotated.sum(0)[::-1]
     return out
 
-"""
-  if 0:
-        # filter the projections
-        freqs = np.zeros((n, 1))
-        freqs[:, 0] = np.linspace(-1, 1, n).T;
-        filter_ft = np.tile(np.abs(freqs), (1, len(theta)))
-        # fourier domain filtering
-        radon_ft = fft(radon_image, axis=0)
-        projection = radon_ft * fftshift(filter_ft)
-        radon_filtered = np.real(ifft(projection, axis=0))        
-    #    print np.max(projection)
-    #    print projection
-        #projection = ifftshift(projection, axes=1);
-    if 0:
-        height, width = radon_image.shape
-        w = np.mgrid[-math.pi:math.pi:(2*math.pi)/height]
-        f = fftshift(abs(w))
-        g = np.array([np.real(ifft(fft(i)*f)) for i in radon_image.T])
-        radon_filtered = np.transpose(g)
-    if 0:
-        img = radon_image.copy()
-        order = 1024
-        filt = np.zeros((order/2, 1))
-        filt[:, 0] = 2.0*np.arange(0, order/2)/order;
-        filt = np.vstack((filt, filt[ ::-1])).T
-        #filt = fftshift(abs(filt))
-    #    order = radon_image.shape[0]
-        w = np.mgrid[-math.pi:math.pi:(2*math.pi)/order]
-        filt = fftshift(abs(w))
-        img.resize((order, img.shape[1]))
-        radon_filtered = np.array([np.real(ifft(fft(column)*filt)) for column in img.T]).T
-        radon_filtered = radon_filtered[:radon_image.shape[0], :]
-    if 0:
-        ### bestest
-        img = radon_image.copy()
-        order = max(64, 2 ** np.ceil(np.log(2*n)/np.log(2)))
-#        filt = np.zeros((order/2, 1))
-#        filt[:, 0] = 2.0*np.arange(0, order/2)/order;
-#        filt = np.vstack((filt, filt[ ::-1])).T
-        #filt = fftshift(abs(filt))
-    #    order = radon_image.shape[0]
-        w = np.mgrid[-math.pi:math.pi:(2*math.pi)/order]
-        filt = fftshift(abs(w))
-        img.resize((order, img.shape[1]))
-        img = fft(img, axis=0)
-        #radon_filtered = np.array([np.real(ifft(column*filt)) for column in img.T]).T
-        radon_filtered = np.array([column*filt for column in img.T]).T
-        
-        radon_filtered = np.real(ifft(radon_filtered, axis=0))
-        radon_filtered = radon_filtered[:radon_image.shape[0], :]
-"""
 
-def iradon(radon_image, theta=None, output_size=None, filter="ramp", interpolate="linear"):
+def iradon(radon_image, theta=None, output_size=None, filter="ramp", interpolation="linear"):
+    """
+    Reconstructs an image from radon transformed data.
+        
+    Parameters
+    ----------
+    radon_image : array_like, dtype=float
+        Image containing radon transform.
+    theta : array_like, dtype=float, optional (default np.arange(180))
+        Reconstruction angles (in degrees).
+    output_size : int
+        Number of rows and columns in the reconstruction.
+    filter : str, optional (default ramp)
+        Filter used in frequency domain filtering. Ramp filter used by default.
+        Filters available: ramp, shepp-logan, cosine, hamming, hann
+        Assign None to use no filter.
+    interpolation : str, optional (default linear)
+        Interpolation method used in reconstruction.
+        Methods available: nearest, linear.
+    
+    Returns
+    -------
+    output : ndarray
+      Reconstructed image.
+
+    Notes
+    -----
+    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. 
+    """
+    if radon_image.ndim != 2:
+        raise ValueError('The input image must be 2-D')
     if theta == None:
         theta = np.arange(180)
     th = (math.pi/180.0)*theta        
     # if output size not specified, estimate from input radon image
     if not output_size:
-        output_size = 2*np.floor(radon_image.shape[0] / (2*np.sqrt(2)))
+        output_size = 2*np.floor(radon_image.shape[0] / (2 * np.sqrt(2)))
     n = radon_image.shape[0]
   
     img = radon_image.copy()
     # resize image to next power of two for fourier analysis
-    order = max(64, 2 ** np.ceil(np.log(2*n)/np.log(2)))
+    # 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
     freqs = np.zeros((order, 1))
     
-    f = fftshift(abs(np.mgrid[-1:1:2/order])).reshape(-1, 1)
+    f = fftshift(abs(np.mgrid[-1:1:2 / order])).reshape(-1, 1)
     w = 2 * math.pi * f
     # start from first element to avoid divide by zero
     if filter == "ramp":
         pass
     elif filter == "shepp-logan":
-        f[1:] = f[1:] * np.sin(w[1:] / 2) / (w[1:]/2)
+        f[1:] = f[1:] * np.sin(w[1:] / 2) / (w[1:] / 2)
     elif filter == "cosine":
        f[1:] = f[1:] * np.cos(w[1:] / 2)
     elif filter == "hamming":
@@ -125,30 +135,32 @@ def iradon(radon_image, theta=None, output_size=None, filter="ramp", interpolate
     [X, Y] = np.mgrid[0.0:x, 0.0:y]
     xpr = X - (output_size + 1.0) / 2.0
     ypr = Y - (output_size + 1.0) / 2.0
-    if interpolate == "nearest":   
+    
+    # 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]))
             reconstructed += radon_filtered[((((k > 0) & (k < n))*k) - 1).astype(np.int), i]             
-    elif interpolate == "linear":
+    elif interpolation == "linear":
         for i in range(len(theta)):
           t = xpr*np.sin(th[i]) - ypr*np.cos(th[i])
           a = np.floor(t)
           b = mid_index + a 
-          reconstructed += (t - a) * radon_filtered[((((b+1 > 0) & (b+1 < n))*(b+1)) - 1).astype(np.int), i] \
-           + (a - t + 1) * radon_filtered[((((b > 0) & (b < n))*b) - 1).astype(np.int), i]
-# XXX slow and inaccurate           
-#    elif interpolate == "spline":
+          b0 = ((((b + 1 > 0) & (b + 1 < n))*(b + 1)) - 1).astype(np.int)
+          b1 = ((((b > 0) & (b < n))*b) - 1).astype(np.int)
+          reconstructed += (t - a) * radon_filtered[b0, i] + (a - t + 1) * radon_filtered[b1, i]
+# XXX slow with some artifacts
+#    elif interpolation == "spline":
 #        axis = np.arange(0, radon_filtered.shape[0]) - mid_index
 #        for i in range(len(theta)):
 #            print i
 #            t = xpr*np.sin(th[i]) - ypr*np.cos(th[i])
 #            #f = interp1d(axis, radon_filtered[:, i], kind="cubic", bounds_error=False, fill_value=0) 
-#            f = interp1d(axis, radon_filtered[:, i], kind="linear", bounds_error=False, fill_value=0) # cubic
+#            f = interp1d(axis, radon_filtered[:, i], kind="linear", bounds_error=False, fill_value=0)
 #            reconstructed += f(t).reshape(output_size, output_size)
     else:
-        raise ValueError("Unknown interpolation: %s" % interpolate)
+        raise ValueError("Unknown interpolation: %s" % interpolation)
         
     return reconstructed * math.pi / (2*len(th))
 
 
-
diff --git a/scikits/image/transform/tests/test_radon_transform.py b/scikits/image/transform/tests/test_radon_transform.py
new file mode 100644
index 00000000..d4bce280
--- /dev/null
+++ b/scikits/image/transform/tests/test_radon_transform.py
@@ -0,0 +1,20 @@
+import numpy as np
+from numpy.testing import *
+from scikits.image.transform import *
+
+
+def test_radon_iradon():
+    size = 100
+    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)
+        delta = np.sum(abs(image/np.max(image) - reconstructed/np.max(reconstructed)))/(size*size)
+        assert delta < 0.1
+    reconstructed = iradon(radon(image), filter="ramp", interpolation="nearest")
+    delta = np.sum(abs(image/np.max(image) - reconstructed/np.max(reconstructed)))/(size*size)
+    assert delta < 0.1
+    
+        
+if __name__ == "__main__":
+    run_module_suite()
+

From b492295575a4bce8ff2d09ccd0ef7c8700a448ce Mon Sep 17 00:00:00 2001
From: Pieter Holtzhausen <holtzhau@gmail.com>
Date: Fri, 19 Aug 2011 00:20:34 +0200
Subject: [PATCH 04/17] Tutorial

---
 doc/source/tutorials/radon_transform.txt | 109 +++++++++++++++++++++++
 1 file changed, 109 insertions(+)
 create mode 100644 doc/source/tutorials/radon_transform.txt

diff --git a/doc/source/tutorials/radon_transform.txt b/doc/source/tutorials/radon_transform.txt
new file mode 100644
index 00000000..03eafa40
--- /dev/null
+++ b/doc/source/tutorials/radon_transform.txt
@@ -0,0 +1,109 @@
+***************
+Radon transform
+***************
+
+The radon transform is a technique widely used in tomography, where you reconstruct an object from its different projections. A projection for example the scattering data obtained as the output of a tomographic scan.
+
+For more information:
+    http://en.wikipedia.org/wiki/Radon_transform
+    http://www.clear.rice.edu/elec431/projects96/DSP/bpanalysis.html
+    
+Forward transform
+=================
+
+First we load the Schepp-Logan phantom, a classic test image representing a tomographic scan.
+
+.. ipython::
+
+    In [1]: from scikits.image.io import imread
+
+    In [1]: from scikits.image import data_dir
+
+    In [2]: from scikits.image.transform import radon, iradon
+
+    In [3]: from scikits.image.color import rgb2gray
+
+    In [4]: import matplotlib.pyplot as plt
+
+    In [5]: import matplotlib.cm as cm
+
+    In [6]: image = rgb2gray(imread(data_dir + "/phantom.png"))
+
+    In [7]: plt.title("original image");
+
+    In [8]: plt.imshow(image, cmap=cm.Greys_r)
+
+    @savefig radon_original_image.png width=4in
+    In [9]: plt.show()
+
+
+Let us illustrate the transform by looking at projections taken at specific angles.
+
+.. ipython::
+
+    In [10]: projections = radon(image, theta=[0, 45, 90])
+
+    In [11]: plt.plot(projections);
+
+    In [12]: plt.title("radon projections");
+
+    In [13]: plt.xlabel("projection axis");
+
+    In [14]: plt.ylabel("intensity");
+
+    @savefig radon_projection_plot1.png width=4in
+    In [15]: plt.show()
+
+We are going to reconstruct an image from 180 (the default) of these projections.
+
+.. ipython::
+
+    In [16]: projections = radon(image)
+
+    In [17]: plt.figure()
+
+    In [18]: plt.title("radon projections");
+
+    In [19]: plt.xlabel("projection axis");
+
+    In [20]: plt.ylabel("intensity");
+
+    In [21]: plt.plot(projections)
+
+    @savefig radon_projection_plot2.png width=4in
+    In [22]: plt.show()
+    
+
+We have now constructed various projections, line integrals of an image, at specific angles. This image is called a sinogram.
+
+.. ipython::
+
+    In [23]: plt.figure()
+    
+    In [24]: plt.title("sinogram");
+
+    In [25]: plt.xlabel("projection axis");
+
+    In [26]: plt.ylabel("intensity");
+
+    In [27]: plt.imshow(projections)
+
+    @savefig radon_sinogram.png width=4in
+    In [28]: plt.show()
+
+Inverse transform
+=================
+To reconstruct the image from this sinogram, we apply the inverse transform.
+
+.. ipython::
+
+    In [29]: reconstruction = iradon(projections)
+
+    In [30]: plt.title("reconstructed image");
+
+    In [31]: plt.imshow(reconstruction, cmap=cm.Greys_r)
+
+    @savefig radon_reconstructed_image.png width=4in
+    In [32]: plt.show()
+
+

From ed6f452c7c1a4aa46282dc73a3903b8cc86b22ef Mon Sep 17 00:00:00 2001
From: Pieter Holtzhausen <holtzhau@gmail.com>
Date: Fri, 19 Aug 2011 01:44:39 +0200
Subject: [PATCH 05/17] Tutorial edit

---
 doc/source/tutorials/radon_transform.txt | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/doc/source/tutorials/radon_transform.txt b/doc/source/tutorials/radon_transform.txt
index 03eafa40..7983b414 100644
--- a/doc/source/tutorials/radon_transform.txt
+++ b/doc/source/tutorials/radon_transform.txt
@@ -54,7 +54,7 @@ Let us illustrate the transform by looking at projections taken at specific angl
     @savefig radon_projection_plot1.png width=4in
     In [15]: plt.show()
 
-We are going to reconstruct an image from 180 (the default) of these projections.
+We are going to reconstruct an image from 180 of these projections (the default).
 
 .. ipython::
 

From 39da9f95c5bfbf6cdc42009adeb5873004fdcc64 Mon Sep 17 00:00:00 2001
From: Pieter Holtzhausen <holtzhau@gmail.com>
Date: Fri, 19 Aug 2011 01:48:55 +0200
Subject: [PATCH 06/17] Phantom image example for radon transform

---
 scikits/image/data/phantom.png | Bin 0 -> 3386 bytes
 1 file changed, 0 insertions(+), 0 deletions(-)
 create mode 100644 scikits/image/data/phantom.png

diff --git a/scikits/image/data/phantom.png b/scikits/image/data/phantom.png
new file mode 100644
index 0000000000000000000000000000000000000000..f9fd232986e73148b9c25bbf87bd27545401620e
GIT binary patch
literal 3386
zcmZu!dt4Iv7T1Cm#0N^SYMGA~Hf7jM&9t>t%zS34fS_iBDW%!1EK^g{NZi)P%`HVm
z&0KYT06F3MDtioOmSm}mrS{e}vxlyaEVS$n?Y)28`?<qs4rk8qd%kmi=XcI;W`5r3
z@2RV8qOGB!q3i9nHBduCGf#bDmLV4XClw{gpXQN3PrAm%d#fgpi*F))J-2FnQU5Mn
ztvQESu!p=jM-cai>Qj@xwj>cTY9007?xFPvV@O<QlKs?{sG*^A*?X%SBXRiMlhD9X
zFT87jGfesO!1L-xoxyF%#&<nEd-Yr*y;<G`#p5YGbFTUh)nb5C)UpvF`=4X*J{)u-
z>%@r@7lQ|36L;v9`8KhX9)Kv#I~N3?LDN_17&xNVgSrLa1S5+m3F8a=t-j?87<Wi~
zk9Yo7g`fpfCiq);6J=z!*8rjBkBNz{l+;uNLS~vWql@2#QbPG2^^7J(kPG%44YoWa
zozGSNFx?NpmVxyPIm%npXKC=SZy24gl@4QSwIRe~it#2cS_ocBtDlLlk8X?fg$8R_
z#J~Z|P!X9O8zPRW3)!qz?RPN8{F~Mm#$rj^W5~>0Sv7~>DQY@f|6mB=mxEDCI;6ug
zYC>h7K!$(h`8J|_rB&8}&Gx>K(no5KgU?zc<U6T72t1Hlr!A!Q7y?a!)E+uysr0G0
z{M(}!K-_u(m@c(P7^q1I1?MY;HiA~Fh1%sSg}u1;i!=JlPU^EQer(NK)cseV*GAn&
z0Q5!_PAz<;&cPxi%xcy?7BQS0>n}~IT+i(^Y5*dR2rU#DM4lanbOuYk#kUJLmfigd
z6w)l-1W<7|&Vll*_2C?eUwNa>e+FU%WUjm$v5J86h?OFdZaf8^EQTTg$|!N29_J-M
zc(6^WoqaQdxTjr)=N!Ig8;NCgfOP}aT%*Y)l`3WaPrq0pn4A1zMZ0adm6EUjN5eF^
zhC=~LVAD<S8AEYHM%uz2un1Z3+zMJn4|togYw2@ymP4i*axsXRG%Q;mb|aPDSjW8f
zrA`nj?{*3&aw0b~DsR`0MF_FkI(GhL0XTwi+S;39eV3DogtVaSE{u%EvWljCq$yls
zHNtyg?cmB`BN23agS_>dGgw6;5kR$#1MWiX8)nQHu6Yney_;Cgx)zILLZ#EI=0)nf
z?xg|2ux+?;gEI<uFrj*Pf|1kiIwt1l-S(@JI-wbbfCAjeh(ac2EJKp>U%Cby11FAE
zaPs&>o(o2J33Mtcc5`-B1VW{SHj!B2v^L-m$8D4Cnz5lp!*Snq7!rZ>6Duxr=KKD7
zyFJgYPW9VV{N(YCTunYPvVEB%6p}+`b9uC2co|;8)Z&d(`*k=T$B+nIp)frS--H8G
zQ6En2QY|g5RyCNJnf){HH{ku;#>lIS4j)SE#b=6{7iHJ`Q;jkox11R^G^t;d4h(eN
z|GWYKVVT+E3RXj!CGnVmKlWA`Ak7J3td>1~fNx-Jqjh($l5FW3ji**}@`s7S5@_DE
zDs--MzV~8phDYB43PuPgX8*R2G?5@^qDbb}{LUL9ilEc^hDed;+#EUW?h5JL7YPJ3
z_KYUq!>GKxe1lnQOAKD}*Od<<m&7QUxSi^kWUR=j=9sgOp3~{E)SYt7vg<CpkQ^FA
zct8rfVH??VJY@A|yBFs#Sr%T3#oAp-xhnB$q+HwTI7a)ST{jct-g});Plq;{7Hpu)
zW>dAOVGCYF$%qSm?S<Uc$*;br4Y_fBiypt*P<FG&q%A#4Nd#(bgI||{Lj#YDjzgrB
z0kAE7!~V+rDVrZ=Y_iQvnlp-|eldK(#GWaw)8Nqkpu_F5w`iPaQ<ujbux%tf>5rrD
zdzIGHnh|o_NT>GA!F-QW3+mkAeeH%4ov_0V(0YbACi_{2z}*y<=S8*Cp?@U(+?!X;
z^3Wp?kUm6`a%?PXQMP<bcWHzjdmeAYjb^=F7b^8*E5|;1LN<}c4S2xBUi?VAwb%b6
z*j7{Fh4K09psVFufmQfssG37~R4qQuV*Rhrdxll-C#{}a>vq@xODq4`H9AoKNj2LA
z3b5JDL!uw$tIl*i3d#MXqBhxvBO6&tN(y6e^_K1Zk9oyXJn)h{Iy-x$Af#sjnpHW6
zQ%4IxSe|m%r0kZ<$B8ev2(35s*9<?{OWsU5HoAFIme%D2H$Gmi*LSdU#@b?*AEz=I
z;fp7_MyF@uJJ@TOz0|YC7!oyK+~<tucfQ^A;>n#6+f6~+q|e}(!_v1}1`^u5!Ni&0
zgD3R_6YV};MQ-=crGpkE3@IJ#yPwzed`2asPHxUgoQ%c;@)DwXf}J`0%hBYLE9Yik
zcbN9+Yr_%chx}B=t}fr71`YNcT+5Vn=2ph`NnS0wl^iFKl8y<?**k6?uhgm8x5|+M
z=o~GAT(|ml@zLgNf8{nO6wVd0SuDyn0PGE5`Gg4yN5Er9MW@lYuJmbLq_=>(Zal8K
z$L2~eQ}zFWA(^vx#+%rpaSyunemXU(>0^qp=EP>aXn6^V)O%b~>+ad8r48I)KGh~N
z#M56;{LbG>TDeB&Dbd|sj=z|*)4>;!Thgv<T!JsQu3^swiaX0kR`5Yhz98uG6q@hR
z5|&JdN(WCW;e~~T8Pm?_HPdcfx%`?T{OgU0R9GK7uR0<)HpXuf^d8LD)H6m3X#Y<Z
z`$sbGB_a>dJ$sBMxAX8omV@|Fh;aidJ4@@nD&<;ewP%_hF0m`%*x>3-?jTa8Cqcmo
zFQIg$p_$%`Qw38jz@4N=3Z)-n*>GPRS_Us0Uyf3^A>kYdcxra5*cDnh`?5OhCz4)9
z20TNtbJSYl<QI!I1<(7bQ+OAhXoAKmwc%lR*+b7&2UobDWM~CWa6477{HCxspire*
zUx39$RJln}xSz^vaIoAB6VH#2=kz6$gYO-lmhGl3`K|LVhIBaREm1hFKUeE-N1#=!
z;E>QV1v?XHI{w`>8XyRH<Rok7l@tKe8&PKn=x}~aBL^+fPchxFX}pi)<Kr{Cchm$;
z9^~p@|6H-ladOM1vg6-<KucO$S`wa(qK>fgpRozYvu+CUuFz{;PNq4ES|xH})S6UX
zPmhxK19-q5hsu5`BL9MUAa37C8ep^{cH6K0mvumgTiFvd(#3kW_vmCy0->EUtS^G{
zh<Ko<-z^sVP8&C=JJ5ZiH$!mM0z>+?XJ&siME;#Bk-0Odi*Rbv^6ULS-kcaZYmLg*
zK&09OK5@w5+b5dvGBq$mwCF}a3asy*Kd6ag?FfC=0b6&?kbrb)f7t83%X9wMK7(~n
zPGNC)N%~%EXdhZ!qYL>?vHSUG$m0TT3$O6$ol9mw>mD1fCIwkgqXs{=D~BIH{unK+
z1>gTpIF`1Y1c6ScPIg+c?*)JpBB-<r0?>+LZ6E}WbLil&0LH0M7+T?^3wXeB59D+r
zP&x!*(26!35CF%W=$Ih@du<@c<Ya4SbYc4r($(DH7$iZ(d+uO`Ln{HVP`PP>z#)^S
z)A85h{qfl$!VIzrBI|Y}#Ix!0=U{|2MAfb(ET&cDa@(f6Z}svH2@`g}mC>J?IU}~~
ze@Jub>$ECyQx|oM|8}c~@5aMZbCAdLd+mo1uXM~hEUG9To6<dY<e_rF2l9Kl`0(}n
z8xtb(_w3MeeS<X$ZzxvhlJQn3F%~YEdbYz@_1G#+TDZahaA!T@^+<E9jT^!@8*N4t
zm{%0sPx3!1=USz1Qc>LKs~jh>Qd*_Zj>rs&%IFdv%AWm>?a-Xzka=fgN&HLO5;wSC
zGT0MHjm7G!c9q20Wq8PG-#!`FuYP))JHTUxdkIe>ld^C)#TW1@So~^EEs@an>S@Z!
zi1ncSb;51iz0<X<TU5e>r*~RVml%HRKK1Jc@<xL4ZQMzP>5V(ZlWmzz?`Jw&RuWj(
z-Gxr#IB8C>ny3J*P>X0Gua5Ztulrg#vVt|t7&2QV2(eHLEJ;+&!mui0qgfY1yw%RJ
z#uB7b=}224vp{YdDZav}D*(AlTebyz8%>9W1di|vV!i);lS>B?%&vN+50po%qeJ$-
z%`Be+gjhyqVuZ1GgbXanNkk3}3amz4)n1*rA?v`4R%K_cVS00wHvH@M=KlQs5j9ml
z(Cs|wg)z}nE2L=+dEzgLJz4pPyP$i(yPdSvrv{aMd1$zGw?doD<O|}Oq=VT9`mJDk
z8=V$bW(ef_*5HAgzEya@cHXxF4=lVZq+76$8f;jv{_Odm<LpPXe5pwjdw%a0^@z8J
K|JI9ikpC|b9t?H>

literal 0
HcmV?d00001


From 3e8cc62e675bffe5eb46bd6d3497b4d693b0ff75 Mon Sep 17 00:00:00 2001
From: Pieter Holtzhausen <holtzhau@gmail.com>
Date: Mon, 22 Aug 2011 14:08:45 +0200
Subject: [PATCH 07/17] Tutorial 80 line length

---
 doc/source/tutorials/radon_transform.txt | 17 ++++++++++++-----
 1 file changed, 12 insertions(+), 5 deletions(-)

diff --git a/doc/source/tutorials/radon_transform.txt b/doc/source/tutorials/radon_transform.txt
index 7983b414..f7b341f7 100644
--- a/doc/source/tutorials/radon_transform.txt
+++ b/doc/source/tutorials/radon_transform.txt
@@ -2,7 +2,9 @@
 Radon transform
 ***************
 
-The radon transform is a technique widely used in tomography, where you reconstruct an object from its different projections. A projection for example the scattering data obtained as the output of a tomographic scan.
+The radon transform is a technique widely used in tomography, where you 
+reconstruct an object from its different projections. A projection for example 
+the scattering data obtained as the output of a tomographic scan.
 
 For more information:
     http://en.wikipedia.org/wiki/Radon_transform
@@ -11,7 +13,8 @@ For more information:
 Forward transform
 =================
 
-First we load the Schepp-Logan phantom, a classic test image representing a tomographic scan.
+First we load the Schepp-Logan phantom, a classic test image representing a 
+tomographic scan.
 
 .. ipython::
 
@@ -37,7 +40,8 @@ First we load the Schepp-Logan phantom, a classic test image representing a tomo
     In [9]: plt.show()
 
 
-Let us illustrate the transform by looking at projections taken at specific angles.
+Let us illustrate the transform by looking at projections taken at specific 
+angles.
 
 .. ipython::
 
@@ -54,7 +58,8 @@ Let us illustrate the transform by looking at projections taken at specific angl
     @savefig radon_projection_plot1.png width=4in
     In [15]: plt.show()
 
-We are going to reconstruct an image from 180 of these projections (the default).
+We are going to reconstruct an image from 180 of these projections (the 
+default).
 
 .. ipython::
 
@@ -74,7 +79,8 @@ We are going to reconstruct an image from 180 of these projections (the default)
     In [22]: plt.show()
     
 
-We have now constructed various projections, line integrals of an image, at specific angles. This image is called a sinogram.
+We have now constructed various projections, line integrals of an image, at 
+specific angles. This image is called a sinogram.
 
 .. ipython::
 
@@ -91,6 +97,7 @@ We have now constructed various projections, line integrals of an image, at spec
     @savefig radon_sinogram.png width=4in
     In [28]: plt.show()
 
+
 Inverse transform
 =================
 To reconstruct the image from this sinogram, we apply the inverse transform.

From 72fc24fc90142d1e4d5116cc8176675b5593bec0 Mon Sep 17 00:00:00 2001
From: Pieter Holtzhausen <holtzhau@gmail.com>
Date: Mon, 22 Aug 2011 14:40:13 +0200
Subject: [PATCH 08/17] Fixes

---
 scikits/image/transform/radon_transform.py | 23 +++++++---------------
 1 file changed, 7 insertions(+), 16 deletions(-)

diff --git a/scikits/image/transform/radon_transform.py b/scikits/image/transform/radon_transform.py
index 77f84f0d..21101e39 100644
--- a/scikits/image/transform/radon_transform.py
+++ b/scikits/image/transform/radon_transform.py
@@ -41,12 +41,12 @@ def radon(image, theta=None):
     if theta == None:
         theta = np.arange(180)    
     height, width = image.shape
-    diagonal = np.sqrt(height**2 + width**2)
+    diagonal = np.sqrt(height ** 2 + width ** 2)
     heightpad = np.ceil(diagonal - height) + 2
     widthpad = np.ceil(diagonal - width) + 2
-    padded_image = np.zeros((int(height+heightpad), int(width+widthpad)))
-    y0, y1 = int(np.ceil(heightpad/2)), int((np.ceil(heightpad/2)+height))
-    x0, x1 = int((np.ceil(widthpad/2))), int((np.ceil(widthpad/2)+width))
+    padded_image = np.zeros((int(height + heightpad), int(width + widthpad)))
+    y0, y1 = int(np.ceil(heightpad / 2)), int((np.ceil(heightpad / 2) + height))
+    x0, x1 = int((np.ceil(widthpad / 2))), int((np.ceil(widthpad / 2) + width))
     padded_image[y0:y1, x0:x1] = image
     out = np.zeros((max(padded_image.shape), len(theta)))
     for i in range(len(theta)):
@@ -140,24 +140,15 @@ def iradon(radon_image, theta=None, output_size=None, filter="ramp", interpolati
     if interpolation == "nearest":   
         for i in range(len(theta)):
             k = np.round(mid_index + xpr*np.sin(th[i]) - ypr*np.cos(th[i]))
-            reconstructed += radon_filtered[((((k > 0) & (k < n))*k) - 1).astype(np.int), i]             
+            reconstructed += radon_filtered[((((k > 0) & (k < n)) * k) - 1).astype(np.int), i]             
     elif interpolation == "linear":
         for i in range(len(theta)):
           t = xpr*np.sin(th[i]) - ypr*np.cos(th[i])
           a = np.floor(t)
           b = mid_index + a 
-          b0 = ((((b + 1 > 0) & (b + 1 < n))*(b + 1)) - 1).astype(np.int)
-          b1 = ((((b > 0) & (b < n))*b) - 1).astype(np.int)
+          b0 = ((((b + 1 > 0) & (b + 1 < n)) * (b + 1)) - 1).astype(np.int)
+          b1 = ((((b > 0) & (b < n)) * b) - 1).astype(np.int)
           reconstructed += (t - a) * radon_filtered[b0, i] + (a - t + 1) * radon_filtered[b1, i]
-# XXX slow with some artifacts
-#    elif interpolation == "spline":
-#        axis = np.arange(0, radon_filtered.shape[0]) - mid_index
-#        for i in range(len(theta)):
-#            print i
-#            t = xpr*np.sin(th[i]) - ypr*np.cos(th[i])
-#            #f = interp1d(axis, radon_filtered[:, i], kind="cubic", bounds_error=False, fill_value=0) 
-#            f = interp1d(axis, radon_filtered[:, i], kind="linear", bounds_error=False, fill_value=0)
-#            reconstructed += f(t).reshape(output_size, output_size)
     else:
         raise ValueError("Unknown interpolation: %s" % interpolation)
         

From 9fb46e1e441e375497cd5c8c9148ef95e52705e2 Mon Sep 17 00:00:00 2001
From: Stefan van der Walt <stefan@sun.ac.za>
Date: Sun, 25 Sep 2011 15:03:01 -0700
Subject: [PATCH 09/17] Add Cython homography implementation.  Not well
 optimised yet.

---
 scikits/image/transform/_project.pyx | 150 +++++++++++++++++++++++++++
 scikits/image/transform/setup.py     |   4 +
 2 files changed, 154 insertions(+)
 create mode 100644 scikits/image/transform/_project.pyx

diff --git a/scikits/image/transform/_project.pyx b/scikits/image/transform/_project.pyx
new file mode 100644
index 00000000..5cbbb608
--- /dev/null
+++ b/scikits/image/transform/_project.pyx
@@ -0,0 +1,150 @@
+#cython: cdivison=True boundscheck=False
+
+cimport cython
+cimport numpy as np
+
+import numpy as np
+import cython
+
+np.import_array()
+
+cdef extern from "math.h":
+    double floor(double)
+    double fmod(double, double)
+
+cdef double get_pixel(np.ndarray image, int r, int c, char mode,
+                      double cval=0):
+    cdef np.ndarray[dtype=np.double_t, ndim=2] img = image
+    cdef int rows = img.shape[0]
+    cdef int cols = img.shape[1]
+
+    if mode == 'C':
+        if (r < 0) or (r >= cols) or (c < 0) or (c >= cols):
+            return cval
+        else:
+            return img[r, c]
+    else:
+        return img[coord_map(rows, r, mode),
+                   coord_map(cols, c, mode)]
+
+cdef int coord_map(int dim, int coord, char mode):
+    dim = dim - 1
+    if mode == 'M': # mirror
+        if (coord < 0):
+            return <int>(-coord % dim)
+        elif (coord > dim):
+            return <int>(dim  - (coord % dim))
+    elif mode == 'W': # wrap
+        if (coord < 0):
+            return <int>(dim - (-coord % dim))
+        elif (coord > dim):
+            return <int>(coord % dim)
+
+    return coord
+
+cdef tf(double x, double y, H):
+    cdef np.ndarray[np.double_t, ndim=2] M = H
+    cdef double xx, yy, zz
+
+    xx = M[0, 0] * x + M[0, 1] * y + M[0, 2]
+    yy = M[1, 0] * x + M[1, 1] * y + M[1, 2]
+    zz = M[2, 0] * x + M[2, 1] * y + M[2, 2]
+
+    xx /= zz
+    yy /= zz
+
+    return xx, yy
+
+@cython.boundscheck(False)
+def homography(np.ndarray image, np.ndarray H, output_shape=None,
+               mode='C', double cval=0):
+    """
+    Projective transformation (homography).
+    
+    Perform a projective transformation (homography) of a
+    floating point image, using bi-linear interpolation.
+
+    For each pixel, given its homogeneous coordinate :math:`\mathbf{x}
+    = [x, y, 1]^T`, its target position is calculated by multiplying
+    with the given matrix, :math:`H`, to give :math:`H \mathbf{x}`.
+    E.g., to rotate by theta degrees clockwise, the matrix should be
+
+    ::
+
+      [[cos(theta) -sin(theta) 0]
+       [sin(theta)  cos(theta) 0]
+       [0            0         1]]
+
+    or, to translate x by 10 and y by 20,
+
+    ::
+
+      [[1 0 10]
+       [0 1 20]
+       [0 0 1 ]].
+
+    Parameters
+    ----------
+    image : 2-D array
+        Input image.
+    H : array of shape ``(3, 3)``
+        Transformation matrix H that defines the homography.
+    output_shape : tuple (rows, cols)
+        Shape of the output image generated.
+    order : int
+        Order of splines used in interpolation.
+    mode : {'C', 'M', 'W'}
+        How to handle values outside the image borders.
+        Constant, Mirror or Wrap.
+    cval : string
+        Used in conjunction with mode 'C' (constant), the value
+        outside the image boundaries.
+
+    """
+
+    cdef np.ndarray[dtype=np.double_t, ndim=2] img = image
+    cdef np.ndarray[dtype=np.double_t, ndim=2] M = np.linalg.inv(H)
+
+    if mode not in ('C', 'W', 'M'):
+        raise ValueError("Invalid mode specified.  Please use "
+                         "C [constant], W [wrap] or M [mirror].")
+    cdef char mode_c = ord(mode[0])
+
+    cdef int out_r, out_c, columns, rows
+    if output_shape is None:
+        out_r = img.shape[0]
+        out_c = img.shape[1]
+    else:
+        out_r = output_shape[0]
+        out_c = output_shape[1]
+
+    rows = img.shape[0]
+    columns = img.shape[1]
+
+    cdef np.ndarray[dtype=np.double_t, ndim=2] out = \
+         np.zeros((out_r, out_c), dtype=np.float64)
+    
+    cdef int tfr, tfc, r_int, c_int
+    cdef double y0, y1, y2, y3
+    cdef double r, c, z, t, u
+
+    for tfr in range(out_r):
+        for tfc in range(out_c):
+            c, r = tf(tfc, tfr, M)
+            r_int = <int>floor(r)
+            c_int = <int>floor(c)
+
+            t = r - r_int
+            u = c - c_int
+
+            y0 = get_pixel(img, r_int, c_int, mode_c)
+            y1 = get_pixel(img, r_int + 1, c_int, mode_c)
+            y2 = get_pixel(img, r_int + 1, c_int + 1, mode_c)
+            y3 = get_pixel(img, r_int, c_int + 1, mode_c)
+
+            out[tfr, tfc] = \
+                (1 - t) * (1 - u) * y0 + \
+                t * (1 - u) * y1 + \
+                t * u * y2 + (1 - t) * u * y3;
+
+    return out
diff --git a/scikits/image/transform/setup.py b/scikits/image/transform/setup.py
index f15dd08c..c5629eb0 100644
--- a/scikits/image/transform/setup.py
+++ b/scikits/image/transform/setup.py
@@ -15,10 +15,14 @@ def configuration(parent_package='', top_path=None):
     config.add_data_dir('tests')
 
     cython(['_hough_transform.pyx'], working_path=base_path)
+    cython(['_project.pyx'], working_path=base_path)
 
     config.add_extension('_hough_transform', sources=['_hough_transform.c'],
                          include_dirs=[get_numpy_include_dirs()])
 
+    config.add_extension('_project', sources=['_project.c'],
+                         include_dirs=[get_numpy_include_dirs()])
+
     return config
 
 if __name__ == '__main__':

From 167b43fc33a6fb1f5cf25cd3ac1f2707c6073247 Mon Sep 17 00:00:00 2001
From: Stefan van der Walt <stefan@sun.ac.za>
Date: Sun, 25 Sep 2011 16:45:52 -0700
Subject: [PATCH 10/17] ENH: fast_homography: 20x speed improvement.

---
 scikits/image/transform/_project.pyx | 71 ++++++++++++++++------------
 1 file changed, 42 insertions(+), 29 deletions(-)

diff --git a/scikits/image/transform/_project.pyx b/scikits/image/transform/_project.pyx
index 5cbbb608..06bc53cb 100644
--- a/scikits/image/transform/_project.pyx
+++ b/scikits/image/transform/_project.pyx
@@ -6,26 +6,24 @@ cimport numpy as np
 import numpy as np
 import cython
 
+from cython.operator import dereference
+
 np.import_array()
 
 cdef extern from "math.h":
     double floor(double)
     double fmod(double, double)
 
-cdef double get_pixel(np.ndarray image, int r, int c, char mode,
-                      double cval=0):
-    cdef np.ndarray[dtype=np.double_t, ndim=2] img = image
-    cdef int rows = img.shape[0]
-    cdef int cols = img.shape[1]
-
+cdef double get_pixel(double *image, int rows, int cols,
+                      int r, int c, char mode, double cval=0):
     if mode == 'C':
         if (r < 0) or (r >= cols) or (c < 0) or (c >= cols):
             return cval
         else:
-            return img[r, c]
+            return image[r * cols + c]
     else:
-        return img[coord_map(rows, r, mode),
-                   coord_map(cols, c, mode)]
+        return image[coord_map(rows, r, mode) * cols +
+                     coord_map(cols, c, mode)]
 
 cdef int coord_map(int dim, int coord, char mode):
     dim = dim - 1
@@ -42,22 +40,28 @@ cdef int coord_map(int dim, int coord, char mode):
 
     return coord
 
-cdef tf(double x, double y, H):
-    cdef np.ndarray[np.double_t, ndim=2] M = H
+cdef tf(double x, double y, double* H, double *x_, double *y_):
     cdef double xx, yy, zz
 
-    xx = M[0, 0] * x + M[0, 1] * y + M[0, 2]
-    yy = M[1, 0] * x + M[1, 1] * y + M[1, 2]
-    zz = M[2, 0] * x + M[2, 1] * y + M[2, 2]
+    ## print
+    ## print H[0], H[1], H[2]
+    ## print H[3], H[4], H[5]
+    ## print H[6], H[7], H[8]
+    ## print
 
-    xx /= zz
-    yy /= zz
+    xx = H[0] * x + H[1] * y + H[2]
+    yy = H[3] * x + H[4] * y + H[5]
+    zz =  H[6] * x + H[7] * y + H[8]
 
-    return xx, yy
+    xx = xx / zz
+    yy = yy / zz
+
+    x_[0] = xx
+    y_[0] = yy
 
 @cython.boundscheck(False)
 def homography(np.ndarray image, np.ndarray H, output_shape=None,
-               mode='C', double cval=0):
+               mode='constant', double cval=0):
     """
     Projective transformation (homography).
     
@@ -93,9 +97,8 @@ def homography(np.ndarray image, np.ndarray H, output_shape=None,
         Shape of the output image generated.
     order : int
         Order of splines used in interpolation.
-    mode : {'C', 'M', 'W'}
+    mode : {'constant', 'mirror', 'wrap'}
         How to handle values outside the image borders.
-        Constant, Mirror or Wrap.
     cval : string
         Used in conjunction with mode 'C' (constant), the value
         outside the image boundaries.
@@ -103,12 +106,18 @@ def homography(np.ndarray image, np.ndarray H, output_shape=None,
     """
 
     cdef np.ndarray[dtype=np.double_t, ndim=2] img = image
-    cdef np.ndarray[dtype=np.double_t, ndim=2] M = np.linalg.inv(H)
+    cdef np.ndarray[dtype=np.double_t, ndim=2] M = \
+         np.ascontiguousarray(np.linalg.inv(H))
 
-    if mode not in ('C', 'W', 'M'):
+    if mode not in ('constant', 'wrap', 'mirror'):
         raise ValueError("Invalid mode specified.  Please use "
-                         "C [constant], W [wrap] or M [mirror].")
-    cdef char mode_c = ord(mode[0])
+                         "`constant`, `wrap` or `mirror`.")
+    if mode == 'constant':
+        mode_c = ord('C')
+    elif mode == 'wrap':
+        mode_c = ord('W')
+    elif mode == 'mirror':
+        mode_c = ord('M')
 
     cdef int out_r, out_c, columns, rows
     if output_shape is None:
@@ -130,17 +139,21 @@ def homography(np.ndarray image, np.ndarray H, output_shape=None,
 
     for tfr in range(out_r):
         for tfc in range(out_c):
-            c, r = tf(tfc, tfr, M)
+            tf(tfc, tfr, <double*>M.data, &c, &r)
             r_int = <int>floor(r)
             c_int = <int>floor(c)
 
             t = r - r_int
             u = c - c_int
 
-            y0 = get_pixel(img, r_int, c_int, mode_c)
-            y1 = get_pixel(img, r_int + 1, c_int, mode_c)
-            y2 = get_pixel(img, r_int + 1, c_int + 1, mode_c)
-            y3 = get_pixel(img, r_int, c_int + 1, mode_c)
+            y0 = get_pixel(<double*>img.data, rows, columns,
+                           r_int, c_int, mode_c)
+            y1 = get_pixel(<double*>img.data, rows, columns,
+                           r_int + 1, c_int, mode_c)
+            y2 = get_pixel(<double*>img.data, rows, columns,
+                           r_int + 1, c_int + 1, mode_c)
+            y3 = get_pixel(<double*>img.data, rows, columns,
+                           r_int, c_int + 1, mode_c)
 
             out[tfr, tfc] = \
                 (1 - t) * (1 - u) * y0 + \

From 69ae2f022b4a1f8ae48301e68689094b10b06954 Mon Sep 17 00:00:00 2001
From: Stefan van der Walt <stefan@sun.ac.za>
Date: Sun, 25 Sep 2011 16:53:55 -0700
Subject: [PATCH 11/17] BUG: fast_homography: Fix mirror mode.

---
 scikits/image/transform/_project.pyx | 7 +++++--
 1 file changed, 5 insertions(+), 2 deletions(-)

diff --git a/scikits/image/transform/_project.pyx b/scikits/image/transform/_project.pyx
index 06bc53cb..e2d08a02 100644
--- a/scikits/image/transform/_project.pyx
+++ b/scikits/image/transform/_project.pyx
@@ -17,7 +17,7 @@ cdef extern from "math.h":
 cdef double get_pixel(double *image, int rows, int cols,
                       int r, int c, char mode, double cval=0):
     if mode == 'C':
-        if (r < 0) or (r >= cols) or (c < 0) or (c >= cols):
+        if (r < 0) or (r > cols - 1) or (c < 0) or (c > cols - 1):
             return cval
         else:
             return image[r * cols + c]
@@ -31,7 +31,10 @@ cdef int coord_map(int dim, int coord, char mode):
         if (coord < 0):
             return <int>(-coord % dim)
         elif (coord > dim):
-            return <int>(dim  - (coord % dim))
+            if (<int>(coord / dim) % 2 != 0):
+                return <int>(dim - (coord % dim))
+            else:
+                return <int>(coord % dim)
     elif mode == 'W': # wrap
         if (coord < 0):
             return <int>(dim - (-coord % dim))

From 2d94a273f6c77aa0c2d0f70092fd8966f141b068 Mon Sep 17 00:00:00 2001
From: Stefan van der Walt <stefan@sun.ac.za>
Date: Sun, 25 Sep 2011 17:22:34 -0700
Subject: [PATCH 12/17] Fix mirror mode. Update docstrings. Update tests

---
 scikits/image/transform/_project.pyx          | 61 ++++++++++++++++---
 scikits/image/transform/tests/test_project.py | 36 +++++++++++
 2 files changed, 88 insertions(+), 9 deletions(-)

diff --git a/scikits/image/transform/_project.pyx b/scikits/image/transform/_project.pyx
index e2d08a02..d1788ab8 100644
--- a/scikits/image/transform/_project.pyx
+++ b/scikits/image/transform/_project.pyx
@@ -1,5 +1,7 @@
 #cython: cdivison=True boundscheck=False
 
+__all__ = ['homography']
+
 cimport cython
 cimport numpy as np
 
@@ -16,6 +18,22 @@ cdef extern from "math.h":
 
 cdef double get_pixel(double *image, int rows, int cols,
                       int r, int c, char mode, double cval=0):
+    """Get a pixel from the image, taking wrapping mode into consideration.
+
+    Parameters
+    ----------
+    image : *double
+        Input image.
+    rows, cols : int
+        Dimensions of image.
+    r, c : int
+        Position at which to get the pixel.
+    mode : {'C', 'W', 'M'}
+        Wrapping mode.  Constant, Wrap or Mirror.
+    cval : double
+        Constant value to use for mode constant.
+    
+    """
     if mode == 'C':
         if (r < 0) or (r > cols - 1) or (c < 0) or (c > cols - 1):
             return cval
@@ -26,10 +44,28 @@ cdef double get_pixel(double *image, int rows, int cols,
                      coord_map(cols, c, mode)]
 
 cdef int coord_map(int dim, int coord, char mode):
+    """
+    Wrap a coordinate, according to a given dimension and mode.
+    
+    Parameters
+    ----------
+    dim : int
+        Maximum coordinate.
+    coord : int
+        Coord provided by user.  May be < 0 or > dim.
+    mode : {'W', 'M'}
+        Whether to wrap or mirror the coordinate if it
+        falls outside [0, dim).
+    
+    """
     dim = dim - 1
     if mode == 'M': # mirror
         if (coord < 0):
-            return <int>(-coord % dim)
+            # How many times times does the coordinate wrap?
+            if (<int>(-coord / dim) % 2 != 0):
+                return dim - <int>(-coord % dim)
+            else:
+                return <int>(-coord % dim)
         elif (coord > dim):
             if (<int>(coord / dim) % 2 != 0):
                 return <int>(dim - (coord % dim))
@@ -44,13 +80,19 @@ cdef int coord_map(int dim, int coord, char mode):
     return coord
 
 cdef tf(double x, double y, double* H, double *x_, double *y_):
-    cdef double xx, yy, zz
+    """Apply a homography to a coordinate.
 
-    ## print
-    ## print H[0], H[1], H[2]
-    ## print H[3], H[4], H[5]
-    ## print H[6], H[7], H[8]
-    ## print
+    Parameters
+    ----------
+    x, y : double
+        Input coordinate.
+    H : (3,3) *double
+        Transformation matrix.
+    x_, y_ : *double
+        Output coordinate.
+
+    """
+    cdef double xx, yy, zz
 
     xx = H[0] * x + H[1] * y + H[2]
     yy = H[3] * x + H[4] * y + H[5]
@@ -108,7 +150,8 @@ def homography(np.ndarray image, np.ndarray H, output_shape=None,
 
     """
 
-    cdef np.ndarray[dtype=np.double_t, ndim=2] img = image
+    cdef np.ndarray[dtype=np.double_t, ndim=2] img = \
+         np.asarray(image, dtype=np.double)
     cdef np.ndarray[dtype=np.double_t, ndim=2] M = \
          np.ascontiguousarray(np.linalg.inv(H))
 
@@ -134,7 +177,7 @@ def homography(np.ndarray image, np.ndarray H, output_shape=None,
     columns = img.shape[1]
 
     cdef np.ndarray[dtype=np.double_t, ndim=2] out = \
-         np.zeros((out_r, out_c), dtype=np.float64)
+         np.zeros((out_r, out_c), dtype=np.double)
     
     cdef int tfr, tfc, r_int, c_int
     cdef double y0, y1, y2, y3
diff --git a/scikits/image/transform/tests/test_project.py b/scikits/image/transform/tests/test_project.py
index 44c0e557..7ae24049 100644
--- a/scikits/image/transform/tests/test_project.py
+++ b/scikits/image/transform/tests/test_project.py
@@ -2,6 +2,8 @@ import numpy as np
 from numpy.testing import assert_array_almost_equal
 
 from scikits.image.transform.project import _stackcopy, homography
+from scikits.image.transform._project import homography as fast_homography
+from scikits.image import data
 
 def test_stackcopy():
     layers = 4
@@ -19,3 +21,37 @@ def test_homography():
                   [0,             0,            1]])
     x90 = homography(x, M, order=1)
     assert_array_almost_equal(x90, np.rot90(x))
+
+def test_fast_homography():
+    img = data.lena()
+    img = img[:, :100]
+    
+    theta = np.deg2rad(30)
+    scale = 0.5
+    tx, ty = 50, 50
+
+    H = np.eye(3)
+    S = scale * np.sin(theta)
+    C = scale * np.cos(theta)
+
+    H[:2, :2] = [[C, -S], [S, C]]
+    H[:2, 2] = [tx, ty]
+
+    for mode in ('constant', 'mirror', 'wrap'):
+        print 'Transform mode:', mode
+
+        p0 = homography(img, H, mode=mode)
+        p1 = fast_homography(img, H, mode=mode)
+        p1 = np.round(p1)
+
+        ## import matplotlib.pyplot as plt
+        ## plt.imshow(np.abs(p0 - p1), cmap=plt.cm.gray)
+        ## plt.show()
+
+        d = np.mean(np.abs(p0 - p1))
+        assert d < 0.1
+    
+
+if __name__ == "__main__":
+    from numpy.testing import run_module_suite
+    run_module_suite()

From 2b71fd9a0eebb038c41fe61432bc4019dcbe52e6 Mon Sep 17 00:00:00 2001
From: Stefan van der Walt <stefan@sun.ac.za>
Date: Sun, 25 Sep 2011 17:29:38 -0700
Subject: [PATCH 13/17] BUG: fast_homography: Fix constant mode.

---
 scikits/image/transform/_project.pyx          |  2 +-
 scikits/image/transform/tests/test_project.py | 10 +++++++---
 2 files changed, 8 insertions(+), 4 deletions(-)

diff --git a/scikits/image/transform/_project.pyx b/scikits/image/transform/_project.pyx
index d1788ab8..3e2b9b7d 100644
--- a/scikits/image/transform/_project.pyx
+++ b/scikits/image/transform/_project.pyx
@@ -35,7 +35,7 @@ cdef double get_pixel(double *image, int rows, int cols,
     
     """
     if mode == 'C':
-        if (r < 0) or (r > cols - 1) or (c < 0) or (c > cols - 1):
+        if (r < 0) or (r > rows - 1) or (c < 0) or (c > cols - 1):
             return cval
         else:
             return image[r * cols + c]
diff --git a/scikits/image/transform/tests/test_project.py b/scikits/image/transform/tests/test_project.py
index 7ae24049..b7c25fdf 100644
--- a/scikits/image/transform/tests/test_project.py
+++ b/scikits/image/transform/tests/test_project.py
@@ -40,16 +40,20 @@ def test_fast_homography():
     for mode in ('constant', 'mirror', 'wrap'):
         print 'Transform mode:', mode
 
-        p0 = homography(img, H, mode=mode)
+        p0 = homography(img, H, mode=mode, order=1)
         p1 = fast_homography(img, H, mode=mode)
         p1 = np.round(p1)
 
         ## import matplotlib.pyplot as plt
-        ## plt.imshow(np.abs(p0 - p1), cmap=plt.cm.gray)
+        ## f, (ax0, ax1, ax2, ax3) = plt.subplots(1, 4)
+        ## ax0.imshow(img)
+        ## ax1.imshow(p0, cmap=plt.cm.gray)
+        ## ax2.imshow(p1, cmap=plt.cm.gray)
+        ## ax3.imshow(np.abs(p0 - p1), cmap=plt.cm.gray)
         ## plt.show()
 
         d = np.mean(np.abs(p0 - p1))
-        assert d < 0.1
+        assert d < 0.2
     
 
 if __name__ == "__main__":

From 3df7df263ea02d3589a26a395c955c33798bb2b8 Mon Sep 17 00:00:00 2001
From: Stefan van der Walt <stefan@sun.ac.za>
Date: Sun, 25 Sep 2011 17:37:04 -0700
Subject: [PATCH 14/17] Import fast homography as transform.fast_homography.

---
 scikits/image/transform/__init__.py           | 1 +
 scikits/image/transform/tests/test_project.py | 4 ++--
 2 files changed, 3 insertions(+), 2 deletions(-)

diff --git a/scikits/image/transform/__init__.py b/scikits/image/transform/__init__.py
index 341671d1..e13d3ac8 100644
--- a/scikits/image/transform/__init__.py
+++ b/scikits/image/transform/__init__.py
@@ -1,4 +1,5 @@
 from .hough_transform import *
 from .finite_radon_transform import *
 from .project import *
+from ._project import homography as fast_homography
 from .integral import *
diff --git a/scikits/image/transform/tests/test_project.py b/scikits/image/transform/tests/test_project.py
index b7c25fdf..ba8e277e 100644
--- a/scikits/image/transform/tests/test_project.py
+++ b/scikits/image/transform/tests/test_project.py
@@ -1,8 +1,8 @@
 import numpy as np
 from numpy.testing import assert_array_almost_equal
 
-from scikits.image.transform.project import _stackcopy, homography
-from scikits.image.transform._project import homography as fast_homography
+from scikits.image.transform.project import _stackcopy
+from scikits.image.transform import homography, fast_homography
 from scikits.image import data
 
 def test_stackcopy():

From b65af2727afc11075d9504e3f1b60759e2a3e40a Mon Sep 17 00:00:00 2001
From: Stefan van der Walt <stefan@sun.ac.za>
Date: Sun, 25 Sep 2011 18:04:27 -0700
Subject: [PATCH 15/17] ENH: Use fast homography instead of PIL's rotate. Minor
 cleanup.

---
 scikits/image/transform/__init__.py           |  2 +-
 scikits/image/transform/radon_transform.py    | 78 +++++++++++++------
 .../transform/tests/test_radon_transform.py   | 25 +++++-
 3 files changed, 76 insertions(+), 29 deletions(-)

diff --git a/scikits/image/transform/__init__.py b/scikits/image/transform/__init__.py
index 01b8fafb..42945fbe 100644
--- a/scikits/image/transform/__init__.py
+++ b/scikits/image/transform/__init__.py
@@ -1,5 +1,5 @@
 from .hough_transform import *
-from .radon import *
+from .radon_transform import *
 from .finite_radon_transform import *
 from .project import *
 from ._project import homography as fast_homography
diff --git a/scikits/image/transform/radon_transform.py b/scikits/image/transform/radon_transform.py
index 21101e39..dfe83eb6 100644
--- a/scikits/image/transform/radon_transform.py
+++ b/scikits/image/transform/radon_transform.py
@@ -14,15 +14,13 @@ References:
 """
 
 import numpy as np
-from scipy.misc import imrotate
-from scipy.interpolate import interp1d
 from scipy.fftpack import fftshift, fft, ifft
-import math
-
+from ._project import homography
 
 def radon(image, theta=None):
     """
-    Calculates the radon transform of an image given specified projection angles.
+    Calculates the radon transform of an image given specified
+    projection angles.
 
     Parameters
     ----------
@@ -35,6 +33,7 @@ def radon(image, theta=None):
     -------
     output : ndarray
         Radon transform.
+
     """
     if image.ndim != 2:
         raise ValueError('The input image must be 2-D')
@@ -44,25 +43,54 @@ def radon(image, theta=None):
     diagonal = np.sqrt(height ** 2 + width ** 2)
     heightpad = np.ceil(diagonal - height) + 2
     widthpad = np.ceil(diagonal - width) + 2
-    padded_image = np.zeros((int(height + heightpad), int(width + widthpad)))
-    y0, y1 = int(np.ceil(heightpad / 2)), int((np.ceil(heightpad / 2) + height))
-    x0, x1 = int((np.ceil(widthpad / 2))), int((np.ceil(widthpad / 2) + width))
+    padded_image = np.zeros((int(height + heightpad),
+                             int(width + widthpad)))
+    y0, y1 = int(np.ceil(heightpad / 2)), \
+             int((np.ceil(heightpad / 2) + height))
+    x0, x1 = int((np.ceil(widthpad / 2))), \
+             int((np.ceil(widthpad / 2) + width))
+
     padded_image[y0:y1, x0:x1] = image
     out = np.zeros((max(padded_image.shape), len(theta)))
+
+    h, w = padded_image.shape
+    shift0 = np.array([[1, 0, -w/2.],
+                       [0, 1, -h/2.],
+                       [0, 0, 1]])
+
+    shift1 = np.array([[1, 0, w/2.],
+                       [0, 1, h/2.],
+                       [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 = imrotate(padded_image, -theta[i])
+        rotated = homography(padded_image,
+                             build_rotation(-theta[i]))
+
         out[:,i] = rotated.sum(0)[::-1]
+
     return out
 
 
-def iradon(radon_image, theta=None, output_size=None, filter="ramp", interpolation="linear"):
+def iradon(radon_image, theta=None, output_size=None,
+           filter="ramp", interpolation="linear"):
     """
-    Reconstructs an image from radon transformed data.
+    Inverse radon transform.
+
+    Reconstruct an image from the radon transform.
         
     Parameters
     ----------
     radon_image : array_like, dtype=float
-        Image containing radon transform.
+        Image containing radon transform (sinogram).
     theta : array_like, dtype=float, optional (default np.arange(180))
         Reconstruction angles (in degrees).
     output_size : int
@@ -82,17 +110,19 @@ def iradon(radon_image, theta=None, output_size=None, filter="ramp", interpolati
 
     Notes
     -----
-    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. 
+    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.
+    
     """
     if radon_image.ndim != 2:
         raise ValueError('The input image must be 2-D')
     if theta == None:
         theta = np.arange(180)
-    th = (math.pi/180.0)*theta        
+    th = (np.pi / 180.0) * theta        
     # if output size not specified, estimate from input radon image
     if not output_size:
-        output_size = 2*np.floor(radon_image.shape[0] / (2 * np.sqrt(2)))
+        output_size = 2 * np.floor(radon_image.shape[0] / (2 * np.sqrt(2)))
     n = radon_image.shape[0]
   
     img = radon_image.copy()
@@ -101,11 +131,11 @@ def iradon(radon_image, theta=None, output_size=None, filter="ramp", interpolati
     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
     freqs = np.zeros((order, 1))
     
     f = fftshift(abs(np.mgrid[-1:1:2 / order])).reshape(-1, 1)
-    w = 2 * math.pi * f
+    w = 2 * np.pi * f
     # start from first element to avoid divide by zero
     if filter == "ramp":
         pass
@@ -139,8 +169,9 @@ def iradon(radon_image, theta=None, output_size=None, filter="ramp", interpolati
     # 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]))
-            reconstructed += radon_filtered[((((k > 0) & (k < n)) * k) - 1).astype(np.int), i]             
+            k = np.round(mid_index + xpr * np.sin(th[i]) - ypr * np.cos(th[i]))
+            reconstructed += radon_filtered[
+                ((((k > 0) & (k < n)) * k) - 1).astype(np.int), i]
     elif interpolation == "linear":
         for i in range(len(theta)):
           t = xpr*np.sin(th[i]) - ypr*np.cos(th[i])
@@ -148,10 +179,9 @@ def iradon(radon_image, theta=None, output_size=None, filter="ramp", interpolati
           b = mid_index + a 
           b0 = ((((b + 1 > 0) & (b + 1 < n)) * (b + 1)) - 1).astype(np.int)
           b1 = ((((b > 0) & (b < n)) * b) - 1).astype(np.int)
-          reconstructed += (t - a) * radon_filtered[b0, i] + (a - t + 1) * radon_filtered[b1, i]
+          reconstructed += (t - a) * radon_filtered[b0, i] + \
+                           (a - t + 1) * radon_filtered[b1, i]
     else:
         raise ValueError("Unknown interpolation: %s" % interpolation)
         
-    return reconstructed * math.pi / (2*len(th))
-
-
+    return reconstructed * np.pi / (2 * len(th))
diff --git a/scikits/image/transform/tests/test_radon_transform.py b/scikits/image/transform/tests/test_radon_transform.py
index d4bce280..a1bb64f8 100644
--- a/scikits/image/transform/tests/test_radon_transform.py
+++ b/scikits/image/transform/tests/test_radon_transform.py
@@ -2,17 +2,34 @@ import numpy as np
 from numpy.testing import *
 from scikits.image.transform import *
 
+def rescale(x):
+    x = x.astype(float, copy=True)
+    x -= x.min()
+    x /= x.max()
+    return x
 
 def test_radon_iradon():
     size = 100
     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)
-        delta = np.sum(abs(image/np.max(image) - reconstructed/np.max(reconstructed)))/(size*size)
-        assert delta < 0.1
+
+        image = rescale(image)
+        reconstructed = rescale(reconstructed)
+        delta = np.mean(np.abs(image - reconstructed))
+
+        ## 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.sum(abs(image/np.max(image) - reconstructed/np.max(reconstructed)))/(size*size)
-    assert delta < 0.1
+    delta = np.mean(abs(image - reconstructed))
+    assert delta < 0.05
     
         
 if __name__ == "__main__":

From 55349987aa7c1769396eb106f9726dded0dcabde Mon Sep 17 00:00:00 2001
From: Emmanuelle Gouillart <emma@aleph.(none)>
Date: Mon, 26 Sep 2011 23:47:35 +0200
Subject: [PATCH 16/17] Change default angles in iradon (as many angles as
 projections) + a few doc fixes

---
 scikits/image/transform/radon_transform.py | 18 +++++++++++-------
 1 file changed, 11 insertions(+), 7 deletions(-)

diff --git a/scikits/image/transform/radon_transform.py b/scikits/image/transform/radon_transform.py
index dfe83eb6..7e3dbe02 100644
--- a/scikits/image/transform/radon_transform.py
+++ b/scikits/image/transform/radon_transform.py
@@ -32,7 +32,7 @@ def radon(image, theta=None):
     Returns
     -------
     output : ndarray
-        Radon transform.
+        Radon transform (sinogram).
 
     """
     if image.ndim != 2:
@@ -85,14 +85,17 @@ def iradon(radon_image, theta=None, output_size=None,
     """
     Inverse radon transform.
 
-    Reconstruct an image from the radon transform.
+    Reconstruct an image from the radon transform, using the filtered
+    back projection algorithm.
         
     Parameters
     ----------
     radon_image : array_like, dtype=float
-        Image containing radon transform (sinogram).
-    theta : array_like, dtype=float, optional (default np.arange(180))
-        Reconstruction angles (in degrees).
+        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 nxm)
     output_size : int
         Number of rows and columns in the reconstruction.
     filter : str, optional (default ramp)
@@ -112,13 +115,14 @@ def iradon(radon_image, theta=None, output_size=None,
     -----
     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.
+    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 == None:
-        theta = np.arange(180)
+        m, n = radon_image.shape
+        theta = np.linspace(0, 180, n, endpoint=False)
     th = (np.pi / 180.0) * theta        
     # if output size not specified, estimate from input radon image
     if not output_size:

From 1e3ad55c84ed2cd997d6540f48bd6fc269b7f4e9 Mon Sep 17 00:00:00 2001
From: Emmanuelle Gouillart <emma@aleph.(none)>
Date: Mon, 26 Sep 2011 23:48:58 +0200
Subject: [PATCH 17/17] More tests on radon and iradon

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

diff --git a/scikits/image/transform/tests/test_radon_transform.py b/scikits/image/transform/tests/test_radon_transform.py
index a1bb64f8..6679528f 100644
--- a/scikits/image/transform/tests/test_radon_transform.py
+++ b/scikits/image/transform/tests/test_radon_transform.py
@@ -3,7 +3,7 @@ from numpy.testing import *
 from scikits.image.transform import *
 
 def rescale(x):
-    x = x.astype(float, copy=True)
+    x = x.astype(float)
     x -= x.min()
     x /= x.max()
     return x
@@ -30,7 +30,33 @@ def test_radon_iradon():
     reconstructed = iradon(radon(image), filter="ramp", interpolation="nearest")
     delta = np.mean(abs(image - reconstructed))
     assert delta < 0.05
-    
+   
+def test_iradon_angles():
+    """
+    Test with different number of projections
+    """
+    size = 100
+    # Synthetic data
+    image = np.tri(size) + np.tri(size)[::-1]
+    # Large number of projections: a good quality is expected
+    nb_angles = 200
+    radon_image_200 = radon(image, theta=np.linspace(0, 180, nb_angles, 
+                    endpoint=False))
+    reconstructed = iradon(radon_image_200)
+    delta_200 = np.mean(abs(rescale(image) - rescale(reconstructed)))
+    assert delta_200 < 0.03
+    # Lower number of projections
+    nb_angles = 80
+    radon_image_80 = radon(image, theta=np.linspace(0, 180, nb_angles,
+                    endpoint=False))
+    # Test whether the sum of all projections is approximately the same
+    s = radon_image_80.sum(axis=0)
+    assert np.allclose(s, s[0], rtol=0.01)
+    reconstructed = iradon(radon_image_80)
+    delta_80 = np.mean(abs(image/np.max(image) - reconstructed/np.max(reconstructed)))
+    # Loss of quality when the number of projections is reduced
+    assert delta_80 > delta_200
+
         
 if __name__ == "__main__":
     run_module_suite()