Merge commit 'v0.5-100-gfeb3e92' into debian

* commit 'v0.5-100-gfeb3e92':
  TST: Increase radon transform threshold to make test more robust.
  BUG: Fix plugin import on systems without PIL or FreeImage.
  BUG: Ensure that the appropriate I/O plugin is used in the test suite.
  BUG: Fix Python 3 syntax error.
  Rename test module to match module.
  Add `__all__` to grey module.
  Remove unused variable.
  Replace `import *` with `import grey`.
  Fix dtype compatibility for functions in morphology.grey
  Make test module runnable.
  Rename greyscale morphology functions.
  Remove unused imports
  Add "page.png" and use for threshold example.
  Fixed padding in radon, iradon. Tests for small images.
  DOC: Add CSS for LaTeX math.
  Change "text" image to grayscale.
  DOC: Add example of adaptive thresholding.

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

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;
+}

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")

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'))

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

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

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

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)):

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