DOC: Move radon example to gallery.

This commit is contained in:
Stefan van der Walt
2011-09-26 21:04:32 -07:00
parent 2c0f48cd35
commit 15e65fa7ca
2 changed files with 57 additions and 116 deletions
+57
View File
@@ -0,0 +1,57 @@
"""
===============
Radon transform
===============
The radon transform is a technique widely used in tomography to
reconstruct an object from different projections. A projection is, for
example, the scattering data obtained as the output of a tomographic
scan.
For more information see:
- http://en.wikipedia.org/wiki/Radon_transform
- http://www.clear.rice.edu/elec431/projects96/DSP/bpanalysis.html
This script performs the radon transform, and reconstructs the
input image based on the resulting sinogram.
"""
import matplotlib.pyplot as plt
from scikits.image.io import imread
from scikits.image import data_dir
from scikits.image.transform import radon, iradon
from scipy.ndimage import zoom
image = imread(data_dir + "/phantom.png", as_grey=True)
image = zoom(image, 0.4)
plt.figure(figsize=(9, 8.5), dpi=75)
plt.subplot(221)
plt.title("Original");
plt.imshow(image, cmap=plt.cm.Greys_r)
plt.subplot(222)
projections = radon(image, theta=[0, 45, 90])
plt.plot(projections);
plt.title("Projections at\n0, 45 and 90 degrees")
plt.xlabel("Projection axis");
plt.ylabel("Intensity");
projections = radon(image)
plt.subplot(223)
plt.title("Radon transform\n(Sinogram)");
plt.xlabel("Projection axis");
plt.ylabel("Intensity");
plt.imshow(projections)
reconstruction = iradon(projections)
plt.subplot(224)
plt.title("Reconstruction\nfrom sinogram")
plt.imshow(reconstruction, cmap=plt.cm.Greys_r)
plt.subplots_adjust(hspace=0.4, wspace=0.5)
plt.show()
-116
View File
@@ -1,116 +0,0 @@
***************
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()