diff --git a/doc/examples/plot_threshold_adaptive.py b/doc/examples/plot_threshold_adaptive.py
new file mode 100644
index 00000000..e1e3cbee
--- /dev/null
+++ b/doc/examples/plot_threshold_adaptive.py
@@ -0,0 +1,48 @@
+"""
+=====================
+Adaptive Thresholding
+=====================
+
+Thresholding is the simplest way to segment objects from a background. If that
+background is relatively uniform, then you can use a global threshold value to
+binarize the image by pixel-intensity. If there's large variation in the
+background intensity, however, adaptive thresholding (a.k.a. local or dynamic
+thresholding) may produce better results.
+
+Here, we binarize an image using the `threshold_adaptive` function, which
+calculates thresholds in regions of size `block_size` surrounding each pixel
+(i.e. local neighborhoods). Each threshold value is the weighted mean of the
+local neighborhood minus an offset value.
+
+"""
+import matplotlib.pyplot as plt
+
+from skimage import data
+from skimage.filter import threshold_otsu, threshold_adaptive
+
+
+image = data.page()
+
+global_thresh = threshold_otsu(image)
+binary_global = image > global_thresh
+
+block_size = 40
+binary_adaptive = threshold_adaptive(image, block_size, offset=10)
+
+fig, axes = plt.subplots(nrows=3, figsize=(7, 8))
+ax0, ax1, ax2 = axes
+plt.gray()
+
+ax0.imshow(image)
+ax0.set_title('Image')
+
+ax1.imshow(binary_global)
+ax1.set_title('Global thresholding')
+
+ax2.imshow(binary_adaptive)
+ax2.set_title('Adaptive thresholding')
+
+for ax in axes:
+    ax.axis('off')
+
+plt.show()
diff --git a/doc/source/themes/agogo/static/agogo.css_t b/doc/source/themes/agogo/static/agogo.css_t
index 26280cd3..e35202c2 100644
--- a/doc/source/themes/agogo/static/agogo.css_t
+++ b/doc/source/themes/agogo/static/agogo.css_t
@@ -737,3 +737,16 @@ p.rubric {
     font-weight: bold;
     font-size: 120%;
 }
+
+/* Math */
+img.math {
+    vertical-align: middle;
+}
+
+div.body div.math p {
+    text-align: center;
+}
+
+span.eqno {
+    float: right;
+}
diff --git a/skimage/data/__init__.py b/skimage/data/__init__.py
index cb3a5b71..c4467678 100644
--- a/skimage/data/__init__.py
+++ b/skimage/data/__init__.py
@@ -96,3 +96,12 @@ def moon():
 
     """
     return load("moon.png")
+
+def page():
+    """Scanned page.
+
+    This image of printed text is useful for demonstrations requiring uneven
+    background illumination.
+
+    """
+    return load("page.png")
diff --git a/skimage/data/page.png b/skimage/data/page.png
new file mode 100644
index 00000000..6c9554ea
Binary files /dev/null and b/skimage/data/page.png differ
diff --git a/skimage/data/text.png b/skimage/data/text.png
index 7135075b..9f6726af 100644
Binary files a/skimage/data/text.png and b/skimage/data/text.png differ
diff --git a/skimage/io/tests/test_freeimage.py b/skimage/io/tests/test_freeimage.py
index 50fef9ee..0f598fe6 100644
--- a/skimage/io/tests/test_freeimage.py
+++ b/skimage/io/tests/test_freeimage.py
@@ -14,6 +14,19 @@ try:
 except OSError:
     FI_available = False
 
+
+def setup_module(self):
+    """The effect of the `plugin.use` call may be overridden by later imports.
+    Call `use_plugin` directly before the tests to ensure that freeimage is
+    used.
+
+    """
+    try:
+        sio.use_plugin('freeimage')
+    except OSError:
+        pass
+
+
 @skipif(not FI_available)
 def test_imread():
     img = sio.imread(os.path.join(si.data_dir, 'color.png'))
diff --git a/skimage/io/tests/test_pil.py b/skimage/io/tests/test_pil.py
index 021c74ea..f54fa9e3 100644
--- a/skimage/io/tests/test_pil.py
+++ b/skimage/io/tests/test_pil.py
@@ -18,6 +18,16 @@ else:
     PIL_available = True
 
 
+def setup_module(self):
+    """The effect of the `plugin.use` call may be overridden by later imports.
+    Call `use_plugin` directly before the tests to ensure that PIL is used.
+
+    """
+    try:
+        use_plugin('pil')
+    except ImportError:
+        pass
+
 @skipif(not PIL_available)
 def test_imread_flatten():
     # a color image is flattened
diff --git a/skimage/morphology/grey.py b/skimage/morphology/grey.py
index c82e2585..ee8542dd 100644
--- a/skimage/morphology/grey.py
+++ b/skimage/morphology/grey.py
@@ -5,11 +5,20 @@
 
 __docformat__ = 'restructuredtext en'
 
+import warnings
+
 import numpy as np
 
-eps = np.finfo(float).eps
+import skimage
 
-def greyscale_erode(image, selem, out=None, shift_x=False, shift_y=False):
+
+__all__ = ['erosion', 'dilation', 'opening', 'closing', 'white_tophat',
+           'black_tophat', 'greyscale_erode', 'greyscale_dilate',
+           'greyscale_open', 'greyscale_close', 'greyscale_white_top_hat',
+           'greyscale_black_top_hat']
+
+
+def erosion(image, selem, out=None, shift_x=False, shift_y=False):
     """Return greyscale morphological erosion of an image.
 
     Morphological erosion sets a pixel at (i,j) to the minimum over all pixels
@@ -19,7 +28,7 @@ def greyscale_erode(image, selem, out=None, shift_x=False, shift_y=False):
     Parameters
     ----------
     image : ndarray
-       The image as a uint8 ndarray.
+       Image array.
 
     selem : ndarray
        The neighborhood expressed as a 2-D array of 1's and 0's.
@@ -34,7 +43,7 @@ def greyscale_erode(image, selem, out=None, shift_x=False, shift_y=False):
 
     Returns
     -------
-    eroded : ndarray
+    eroded : uint8 array
        The result of the morphological erosion.
 
     Examples
@@ -46,7 +55,7 @@ def greyscale_erode(image, selem, out=None, shift_x=False, shift_y=False):
     ...                           [0, 1, 1, 1, 0],
     ...                           [0, 1, 1, 1, 0],
     ...                           [0, 0, 0, 0, 0]], dtype=np.uint8)
-    >>> greyscale_erode(bright_square, square(3))
+    >>> erosion(bright_square, square(3))
     array([[0, 0, 0, 0, 0],
            [0, 0, 0, 0, 0],
            [0, 0, 1, 0, 0],
@@ -56,6 +65,8 @@ def greyscale_erode(image, selem, out=None, shift_x=False, shift_y=False):
     """
     if image is out:
         raise NotImplementedError("In-place erosion not supported!")
+    image = skimage.img_as_ubyte(image)
+
     try:
         import skimage.morphology.cmorph as cmorph
         out = cmorph.erode(image, selem, out=out,
@@ -64,7 +75,8 @@ def greyscale_erode(image, selem, out=None, shift_x=False, shift_y=False):
     except ImportError:
         raise ImportError("cmorph extension not available.")
 
-def greyscale_dilate(image, selem, out=None, shift_x=False, shift_y=False):
+
+def dilation(image, selem, out=None, shift_x=False, shift_y=False):
     """Return greyscale morphological dilation of an image.
 
     Morphological dilation sets a pixel at (i,j) to the maximum over all pixels
@@ -75,7 +87,7 @@ def greyscale_dilate(image, selem, out=None, shift_x=False, shift_y=False):
     ----------
 
     image : ndarray
-       The image as a uint8 ndarray.
+       Image array.
 
     selem : ndarray
        The neighborhood expressed as a 2-D array of 1's and 0's.
@@ -90,7 +102,7 @@ def greyscale_dilate(image, selem, out=None, shift_x=False, shift_y=False):
 
     Returns
     -------
-    dilated : ndarray
+    dilated : uint8 array
        The result of the morphological dilation.
 
     Examples
@@ -102,7 +114,7 @@ def greyscale_dilate(image, selem, out=None, shift_x=False, shift_y=False):
     ...                          [0, 0, 1, 0, 0],
     ...                          [0, 0, 0, 0, 0],
     ...                          [0, 0, 0, 0, 0]], dtype=np.uint8)
-    >>> greyscale_dilate(bright_pixel, square(3))
+    >>> dilation(bright_pixel, square(3))
     array([[0, 0, 0, 0, 0],
            [0, 1, 1, 1, 0],
            [0, 1, 1, 1, 0],
@@ -112,6 +124,8 @@ def greyscale_dilate(image, selem, out=None, shift_x=False, shift_y=False):
     """
     if image is out:
         raise NotImplementedError("In-place dilation not supported!")
+    image = skimage.img_as_ubyte(image)
+
     try:
         from . import cmorph
         out = cmorph.dilate(image, selem, out=out,
@@ -120,7 +134,8 @@ def greyscale_dilate(image, selem, out=None, shift_x=False, shift_y=False):
     except ImportError:
         raise ImportError("cmorph extension not available.")
 
-def greyscale_open(image, selem, out=None):
+
+def opening(image, selem, out=None):
     """Return greyscale morphological opening of an image.
 
     The morphological opening on an image is defined as an erosion followed by
@@ -131,7 +146,7 @@ def greyscale_open(image, selem, out=None):
     Parameters
     ----------
     image : ndarray
-       The image as a uint8 ndarray.
+       Image array.
 
     selem : ndarray
        The neighborhood expressed as a 2-D array of 1's and 0's.
@@ -142,7 +157,7 @@ def greyscale_open(image, selem, out=None):
 
     Returns
     -------
-    opening : ndarray
+    opening : uint8 array
        The result of the morphological opening.
 
     Examples
@@ -154,7 +169,7 @@ def greyscale_open(image, selem, out=None):
     ...                            [1, 1, 1, 1, 1],
     ...                            [1, 1, 0, 1, 1],
     ...                            [1, 0, 0, 0, 1]], dtype=np.uint8)
-    >>> greyscale_open(bad_connection, square(3))
+    >>> opening(bad_connection, square(3))
     array([[0, 0, 0, 0, 0],
            [1, 1, 0, 1, 1],
            [1, 1, 0, 1, 1],
@@ -166,12 +181,12 @@ def greyscale_open(image, selem, out=None):
     shift_x = True if (w % 2) == 0 else False
     shift_y = True if (h % 2) == 0 else False
 
-    eroded = greyscale_erode(image, selem)
-    out = greyscale_dilate(eroded, selem, out=out,
-                           shift_x=shift_x, shift_y=shift_y)
+    eroded = erosion(image, selem)
+    out = dilation(eroded, selem, out=out, shift_x=shift_x, shift_y=shift_y)
     return out
 
-def greyscale_close(image, selem, out=None):
+
+def closing(image, selem, out=None):
     """Return greyscale morphological closing of an image.
 
     The morphological closing on an image is defined as a dilation followed by
@@ -182,7 +197,7 @@ def greyscale_close(image, selem, out=None):
     Parameters
     ----------
     image : ndarray
-       The image as a uint8 ndarray.
+       Image array.
 
     selem : ndarray
        The neighborhood expressed as a 2-D array of 1's and 0's.
@@ -193,8 +208,8 @@ def greyscale_close(image, selem, out=None):
 
     Returns
     -------
-    opening : ndarray
-       The result of the morphological opening.
+    closing : uint8 array
+       The result of the morphological closing.
 
     Examples
     --------
@@ -205,7 +220,7 @@ def greyscale_close(image, selem, out=None):
     ...                         [1, 1, 0, 1, 1],
     ...                         [0, 0, 0, 0, 0],
     ...                         [0, 0, 0, 0, 0]], dtype=np.uint8)
-    >>> greyscale_close(broken_line, square(3))
+    >>> closing(broken_line, square(3))
     array([[0, 0, 0, 0, 0],
            [0, 0, 0, 0, 0],
            [1, 1, 1, 1, 1],
@@ -217,12 +232,12 @@ def greyscale_close(image, selem, out=None):
     shift_x = True if (w % 2) == 0 else False
     shift_y = True if (h % 2) == 0 else False
 
-    dilated = greyscale_dilate(image, selem)
-    out = greyscale_erode(dilated, selem, out=out,
-                          shift_x=shift_x, shift_y=shift_y)
+    dilated = dilation(image, selem)
+    out = erosion(dilated, selem, out=out, shift_x=shift_x, shift_y=shift_y)
     return out
 
-def greyscale_white_top_hat(image, selem, out=None):
+
+def white_tophat(image, selem, out=None):
     """Return white top hat of an image.
 
     The white top hat of an image is defined as the image minus its
@@ -232,7 +247,7 @@ def greyscale_white_top_hat(image, selem, out=None):
     Parameters
     ----------
     image : ndarray
-       The image as a uint8 ndarray.
+       Image array.
 
     selem : ndarray
        The neighborhood expressed as a 2-D array of 1's and 0's.
@@ -243,7 +258,7 @@ def greyscale_white_top_hat(image, selem, out=None):
 
     Returns
     -------
-    opening : ndarray
+    opening : uint8 array
        The result of the morphological white top hat.
 
     Examples
@@ -255,7 +270,7 @@ def greyscale_white_top_hat(image, selem, out=None):
     ...                            [3, 5, 9, 5, 3],
     ...                            [3, 4, 5, 4, 3],
     ...                            [2, 3, 3, 3, 2]], dtype=np.uint8)
-    >>> greyscale_white_top_hat(bright_on_grey, square(3))
+    >>> white_tophat(bright_on_grey, square(3))
     array([[0, 0, 0, 0, 0],
            [0, 0, 1, 0, 0],
            [0, 1, 5, 1, 0],
@@ -265,12 +280,14 @@ def greyscale_white_top_hat(image, selem, out=None):
    """
     if image is out:
         raise NotImplementedError("Cannot perform white top hat in place.")
+    image = skimage.img_as_ubyte(image)
 
-    out = greyscale_open(image, selem, out=out)
+    out = opening(image, selem, out=out)
     out = image - out
     return out
 
-def greyscale_black_top_hat(image, selem, out=None):
+
+def black_tophat(image, selem, out=None):
     """Return black top hat of an image.
 
     The black top hat of an image is defined as its morphological closing minus
@@ -281,7 +298,7 @@ def greyscale_black_top_hat(image, selem, out=None):
     Parameters
     ----------
     image : ndarray
-       The image as a uint8 ndarray.
+       Image array.
 
     selem : ndarray
        The neighborhood expressed as a 2-D array of 1's and 0's.
@@ -292,7 +309,7 @@ def greyscale_black_top_hat(image, selem, out=None):
 
     Returns
     -------
-    opening : ndarray
+    opening : uint8 array
        The result of the black top filter.
 
     Examples
@@ -304,7 +321,7 @@ def greyscale_black_top_hat(image, selem, out=None):
     ...                          [6, 4, 0, 4, 6],
     ...                          [6, 5, 4, 5, 6],
     ...                          [7, 6, 6, 6, 7]], dtype=np.uint8)
-    >>> greyscale_black_top_hat(dark_on_grey, square(3))
+    >>> black_tophat(dark_on_grey, square(3))
     array([[0, 0, 0, 0, 0],
            [0, 0, 1, 0, 0],
            [0, 1, 5, 1, 0],
@@ -314,7 +331,34 @@ def greyscale_black_top_hat(image, selem, out=None):
     """
     if image is out:
         raise NotImplementedError("Cannot perform white top hat in place.")
-    out = greyscale_close(image, selem, out=out)
+    image = skimage.img_as_ubyte(image)
+
+    out = closing(image, selem, out=out)
     out = out - image
     return out
 
+
+def greyscale_erode(*args, **kwargs):
+    warnings.warn("`greyscale_erode` renamed `erosion`.")
+    return erosion(*args, **kwargs)
+
+def greyscale_dilate(*args, **kwargs):
+    warnings.warn("`greyscale_dilate` renamed `dilation`.")
+    return dilation(*args, **kwargs)
+
+def greyscale_open(*args, **kwargs):
+    warnings.warn("`greyscale_open` renamed `opening`.")
+    return opening(*args, **kwargs)
+
+def greyscale_close(*args, **kwargs):
+    warnings.warn("`greyscale_close` renamed `closing`.")
+    return closing(*args, **kwargs)
+
+def greyscale_white_top_hat(*args, **kwargs):
+    warnings.warn("`greyscale_white_top_hat` renamed `white_tophat`.")
+    return white_tophat(*args, **kwargs)
+
+def greyscale_black_top_hat(*args, **kwargs):
+    warnings.warn("`greyscale_black_top_hat` renamed `black_tophat`.")
+    return black_tophat(*args, **kwargs)
+
diff --git a/skimage/morphology/tests/test_morphology.py b/skimage/morphology/tests/test_grey.py
similarity index 51%
rename from skimage/morphology/tests/test_morphology.py
rename to skimage/morphology/tests/test_grey.py
index 2d089c58..1103b20d 100644
--- a/skimage/morphology/tests/test_morphology.py
+++ b/skimage/morphology/tests/test_grey.py
@@ -1,12 +1,13 @@
 import os.path
 
 import numpy as np
-from numpy.testing import *
+from numpy import testing
 
+import skimage
 from skimage import data_dir
-from skimage.io import imread
-from skimage import data_dir
-from skimage.morphology import *
+from skimage.morphology import grey
+from skimage.morphology import selem
+
 
 lena = np.load(os.path.join(data_dir, 'lena_GRAY_U8.npy'))
 
@@ -19,48 +20,48 @@ class TestMorphology():
             expected_result = matlab_results[arrname]
             mask = strel_func(k)
             actual_result = morph_func(lena, mask)
-            assert_equal(expected_result, actual_result)
+            testing.assert_equal(expected_result, actual_result)
             k = k + 1
 
     def test_erode_diamond(self):
         self.morph_worker(lena, "diamond-erode-matlab-output.npz",
-                          greyscale_erode, diamond)
+                          grey.erosion, selem.diamond)
 
     def test_dilate_diamond(self):
         self.morph_worker(lena, "diamond-dilate-matlab-output.npz",
-                          greyscale_dilate, diamond)
+                          grey.dilation, selem.diamond)
 
     def test_open_diamond(self):
         self.morph_worker(lena, "diamond-open-matlab-output.npz",
-                          greyscale_open, diamond)
+                          grey.opening, selem.diamond)
 
     def test_close_diamond(self):
         self.morph_worker(lena, "diamond-close-matlab-output.npz",
-                          greyscale_close, diamond)
+                          grey.closing, selem.diamond)
 
     def test_tophat_diamond(self):
         self.morph_worker(lena, "diamond-tophat-matlab-output.npz",
-                          greyscale_white_top_hat, diamond)
+                          grey.white_tophat, selem.diamond)
 
     def test_bothat_diamond(self):
         self.morph_worker(lena, "diamond-bothat-matlab-output.npz",
-                          greyscale_black_top_hat, diamond)
+                          grey.black_tophat, selem.diamond)
 
     def test_erode_disk(self):
         self.morph_worker(lena, "disk-erode-matlab-output.npz",
-                          greyscale_erode, disk)
+                          grey.erosion, selem.disk)
 
     def test_dilate_disk(self):
         self.morph_worker(lena, "disk-dilate-matlab-output.npz",
-                          greyscale_dilate, disk)
+                          grey.dilation, selem.disk)
 
     def test_open_disk(self):
         self.morph_worker(lena, "disk-open-matlab-output.npz",
-                          greyscale_open, disk)
+                          grey.opening, selem.disk)
 
     def test_close_disk(self):
         self.morph_worker(lena, "disk-close-matlab-output.npz",
-                          greyscale_close, disk)
+                          grey.closing, selem.disk)
 
 
 class TestEccentricStructuringElements():
@@ -69,50 +70,89 @@ class TestEccentricStructuringElements():
         self.black_pixel = 255 * np.ones((4, 4), dtype=np.uint8)
         self.black_pixel[1, 1] = 0
         self.white_pixel = 255 - self.black_pixel
-        self.selems = [square(2), rectangle(2, 2),
-                       rectangle(2, 1), rectangle(1, 2)]
+        self.selems = [selem.square(2), selem.rectangle(2, 2),
+                       selem.rectangle(2, 1), selem.rectangle(1, 2)]
 
     def test_dilate_erode_symmetry(self):
         for s in self.selems:
-            c = greyscale_erode(self.black_pixel, s)
-            d = greyscale_dilate(self.white_pixel, s)
+            c = grey.erosion(self.black_pixel, s)
+            d = grey.dilation(self.white_pixel, s)
             assert np.all(c == (255 - d))
 
     def test_open_black_pixel(self):
         for s in self.selems:
-            grey_open = greyscale_open(self.black_pixel, s)
+            grey_open = grey.opening(self.black_pixel, s)
             assert np.all(grey_open == self.black_pixel)
 
     def test_close_white_pixel(self):
         for s in self.selems:
-            grey_close = greyscale_close(self.white_pixel, s)
+            grey_close = grey.closing(self.white_pixel, s)
             assert np.all(grey_close == self.white_pixel)
 
     def test_open_white_pixel(self):
         for s in self.selems:
-            assert np.all(greyscale_open(self.white_pixel, s) == 0)
+            assert np.all(grey.opening(self.white_pixel, s) == 0)
 
     def test_close_black_pixel(self):
         for s in self.selems:
-            assert np.all(greyscale_close(self.black_pixel, s) == 255)
+            assert np.all(grey.closing(self.black_pixel, s) == 255)
 
     def test_white_tophat_white_pixel(self):
         for s in self.selems:
-            tophat = greyscale_white_top_hat(self.white_pixel, s)
+            tophat = grey.white_tophat(self.white_pixel, s)
             assert np.all(tophat == self.white_pixel)
 
     def test_black_tophat_black_pixel(self):
         for s in self.selems:
-            tophat = greyscale_black_top_hat(self.black_pixel, s)
+            tophat = grey.black_tophat(self.black_pixel, s)
             assert np.all(tophat == (255 - self.black_pixel))
 
     def test_white_tophat_black_pixel(self):
         for s in self.selems:
-            tophat = greyscale_white_top_hat(self.black_pixel, s)
+            tophat = grey.white_tophat(self.black_pixel, s)
             assert np.all(tophat == 0)
 
     def test_black_tophat_white_pixel(self):
         for s in self.selems:
-            tophat = greyscale_black_top_hat(self.white_pixel, s)
+            tophat = grey.black_tophat(self.white_pixel, s)
             assert np.all(tophat == 0)
 
+
+class TestDTypes():
+
+    def setUp(self):
+        k = 5
+        arrname = '%03i' % k
+
+        self.disk = selem.disk(k)
+
+        fname_opening = os.path.join(data_dir, "disk-open-matlab-output.npz")
+        self.expected_opening = np.load(fname_opening)[arrname]
+
+        fname_closing = os.path.join(data_dir, "disk-close-matlab-output.npz")
+        self.expected_closing = np.load(fname_closing)[arrname]
+
+    def _test_image(self, image):
+        result_opening = grey.opening(image, self.disk)
+        testing.assert_equal(result_opening, self.expected_opening)
+
+        result_closing = grey.closing(image, self.disk)
+        testing.assert_equal(result_closing, self.expected_closing)
+
+    def test_float(self):
+        image = skimage.img_as_float(lena)
+        self._test_image(image)
+
+    @testing.decorators.skipif(True)
+    def test_int(self):
+        image = skimage.img_as_int(lena)
+        self._test_image(image)
+
+    def test_uint(self):
+        image = skimage.img_as_uint(lena)
+        self._test_image(image)
+
+
+if __name__ == '__main__':
+    testing.run_module_suite()
+
diff --git a/skimage/transform/radon_transform.py b/skimage/transform/radon_transform.py
index 759a2038..ef19319c 100644
--- a/skimage/transform/radon_transform.py
+++ b/skimage/transform/radon_transform.py
@@ -44,8 +44,8 @@ def radon(image, 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
+    heightpad = np.ceil(diagonal - height)
+    widthpad = np.ceil(diagonal - width)
     padded_image = np.zeros((int(height + heightpad),
                              int(width + widthpad)))
     y0, y1 = int(np.ceil(heightpad / 2)), \
@@ -57,14 +57,16 @@ def radon(image, theta=None):
     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.],
+    dh, dw = h / 2, w / 2
+    shift0 = np.array([[1, 0, -dw],
+                       [0, 1, -dh],
                        [0, 0, 1]])
 
-    shift1 = np.array([[1, 0, w/2.],
-                       [0, 1, h/2.],
+    shift1 = np.array([[1, 0, dw],
+                       [0, 1, dh],
                        [0, 0, 1]])
 
+
     def build_rotation(theta):
         T = -np.deg2rad(theta)
 
@@ -129,7 +131,7 @@ def iradon(radon_image, theta=None, output_size=None,
     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)))
+        output_size = int(np.floor(np.sqrt((radon_image.shape[0]) ** 2 / 2.0)))
     n = radon_image.shape[0]
   
     img = radon_image.copy()
@@ -166,13 +168,14 @@ def iradon(radon_image, theta=None, output_size=None,
     # 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);
+    mid_index = np.ceil(n / 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
-    
+    xpr = X - int(output_size) / 2
+    ypr = Y - int(output_size) / 2 
+
     # reconstruct image by interpolation
     if interpolation == "nearest":   
         for i in range(len(theta)):
diff --git a/skimage/transform/tests/test_radon_transform.py b/skimage/transform/tests/test_radon_transform.py
index b33a8887..b52ff69b 100644
--- a/skimage/transform/tests/test_radon_transform.py
+++ b/skimage/transform/tests/test_radon_transform.py
@@ -1,3 +1,5 @@
+from __future__ import print_function
+
 import numpy as np
 from numpy.testing import *
 from skimage.transform import *
@@ -10,6 +12,7 @@ def rescale(x):
 
 def test_radon_iradon():
     size = 100
+    debug = False
     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)
@@ -18,12 +21,13 @@ def test_radon_iradon():
         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()
+        if debug:
+            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
 
@@ -33,7 +37,7 @@ def test_radon_iradon():
     size = 20
     image = np.tri(size) + np.tri(size)[::-1]
     reconstructed = iradon(radon(image), filter="ramp", interpolation="nearest")
-   
+
 def test_iradon_angles():
     """
     Test with different number of projections
@@ -43,7 +47,7 @@ def test_iradon_angles():
     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, 
+    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)))
@@ -60,7 +64,33 @@ def test_iradon_angles():
     # Loss of quality when the number of projections is reduced
     assert delta_80 > delta_200
 
-        
+def test_radon_minimal():
+    """
+    Test for small images for various angles
+    """
+    thetas = [np.arange(180)]
+    for theta in thetas:
+        a = np.zeros((3, 3))
+        a[1, 1] = 1
+        p = radon(a, theta)
+        reconstructed = iradon(p, theta)
+        reconstructed /= np.max(reconstructed)
+        assert np.all(abs(a - reconstructed) < 0.4)
+
+        b = np.zeros((4, 4))
+        b[1:3, 1:3] = 1
+        p = radon(b, theta)
+        reconstructed = iradon(p, theta)
+        reconstructed /= np.max(reconstructed)
+        assert np.all(abs(b - reconstructed) < 0.4)
+
+        c = np.zeros((5, 5))
+        c[1:3, 1:3] = 1
+        p = radon(c, theta)
+        reconstructed = iradon(p, theta)
+        reconstructed /= np.max(reconstructed)
+        assert np.all(abs(c - reconstructed) < 0.4)
+
+
 if __name__ == "__main__":
     run_module_suite()
-