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 341671d1..42945fbe 100644
--- a/scikits/image/transform/__init__.py
+++ b/scikits/image/transform/__init__.py
@@ -1,4 +1,6 @@
 from .hough_transform import *
+from .radon_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/_project.pyx b/scikits/image/transform/_project.pyx
new file mode 100644
index 00000000..3e2b9b7d
--- /dev/null
+++ b/scikits/image/transform/_project.pyx
@@ -0,0 +1,209 @@
+#cython: cdivison=True boundscheck=False
+
+__all__ = ['homography']
+
+cimport cython
+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(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 > rows - 1) or (c < 0) or (c > cols - 1):
+            return cval
+        else:
+            return image[r * cols + c]
+    else:
+        return image[coord_map(rows, r, mode) * 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):
+            # 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))
+            else:
+                return <int>(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, double* H, double *x_, double *y_):
+    """Apply a homography to a coordinate.
+
+    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]
+    zz =  H[6] * x + H[7] * y + H[8]
+
+    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='constant', 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 : {'constant', 'mirror', 'wrap'}
+        How to handle values outside the image borders.
+    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 = \
+         np.asarray(image, dtype=np.double)
+    cdef np.ndarray[dtype=np.double_t, ndim=2] M = \
+         np.ascontiguousarray(np.linalg.inv(H))
+
+    if mode not in ('constant', 'wrap', 'mirror'):
+        raise ValueError("Invalid mode specified.  Please use "
+                         "`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:
+        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.double)
+    
+    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):
+            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(<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 + \
+                t * (1 - u) * y1 + \
+                t * u * y2 + (1 - t) * u * y3;
+
+    return out
diff --git a/scikits/image/transform/radon_transform.py b/scikits/image/transform/radon_transform.py
new file mode 100644
index 00000000..7e3dbe02
--- /dev/null
+++ b/scikits/image/transform/radon_transform.py
@@ -0,0 +1,191 @@
+"""
+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.fftpack import fftshift, fft, ifft
+from ._project import homography
+
+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 (sinogram).
+
+    """
+    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)))
+
+    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 = 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"):
+    """
+    Inverse 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). 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)
+        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. This algorithm is called filtered back projection.
+    
+    """
+    if radon_image.ndim != 2:
+        raise ValueError('The input image must be 2-D')
+    if theta == None:
+        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:
+        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 * np.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 * np.pi / (2 * len(th))
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__':
diff --git a/scikits/image/transform/tests/test_project.py b/scikits/image/transform/tests/test_project.py
index 44c0e557..ba8e277e 100644
--- a/scikits/image/transform/tests/test_project.py
+++ b/scikits/image/transform/tests/test_project.py
@@ -1,7 +1,9 @@
 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 _stackcopy
+from scikits.image.transform import homography, fast_homography
+from scikits.image import data
 
 def test_stackcopy():
     layers = 4
@@ -19,3 +21,41 @@ 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, order=1)
+        p1 = fast_homography(img, H, mode=mode)
+        p1 = np.round(p1)
+
+        ## import matplotlib.pyplot as plt
+        ## 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.2
+    
+
+if __name__ == "__main__":
+    from numpy.testing import run_module_suite
+    run_module_suite()
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..6679528f
--- /dev/null
+++ b/scikits/image/transform/tests/test_radon_transform.py
@@ -0,0 +1,63 @@
+import numpy as np
+from numpy.testing import *
+from scikits.image.transform import *
+
+def rescale(x):
+    x = x.astype(float)
+    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)
+
+        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.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()
+