Merge radon transform by Pieter Holtzhausen.

doc/source/tutorials/radon_transform.txt

@@ -0,0 +1,116 @@
+***************
+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 of these projections (the 
+default).
+
+.. 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()
+
+

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 *
 

scikits/image/transform/radon_transform.py

@@ -0,0 +1,157 @@
+"""
+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, fft, ifft
+import math
+
+
+def radon(image, theta=None):
+    """
+    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)    
+    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
+
+
+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)))
+    n = radon_image.shape[0]
+  
+    img = radon_image.copy()
+    # resize image to next power of two for fourier analysis
+    # speeds up fourier and lessens artifacts
+    order = max(64, 2 ** np.ceil(np.log(2 * n) / np.log(2)))
+    # zero pad input image
+    img.resize((order, img.shape[1]))
+    #construct the fourier filter
+    freqs = np.zeros((order, 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)
+    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(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))
+    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
+    
+    # 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 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)
+          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))
+
+

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()
+