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 1/8] 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 2/8] 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 3/8] 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 4/8] 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 5/8] 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 6/8] 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 7/8] 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 8/8] 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)