diff --git a/doc/source/tutorials/radon_transform.txt b/doc/source/tutorials/radon_transform.txt new file mode 100644 index 00000000..f7b341f7 --- /dev/null +++ b/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() + + diff --git a/scikits/image/data/phantom.png b/scikits/image/data/phantom.png new file mode 100644 index 00000000..f9fd2329 Binary files /dev/null and b/scikits/image/data/phantom.png differ 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..21101e39 --- /dev/null +++ b/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)) + + 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() +