Merge branches 'holtzhau_edge_filters' and 'ccomp'.

CONTRIBUTORS.txt

@@ -38,3 +38,5 @@
 - Dan Farmer
   Incorporating CellProfiler's Canny edge detector
 
+- Pieter Holtzhausen
+  Incorporating CellProfiler's Sobel edge detector

scikits/image/filter/__init__.py

@@ -1,3 +1,4 @@
 from lpi_filter import *
 from ctmf import median_filter
 from canny import canny
+from edges import sobel, hsobel, vsobel, hprewitt, vprewitt, prewitt

scikits/image/filter/edges.py

@@ -0,0 +1,204 @@
+"""edges.py - Sobel edge filter
+
+Originally part of CellProfiler, code licensed under both GPL and BSD licenses.
+Website: http://www.cellprofiler.org
+Copyright (c) 2003-2009 Massachusetts Institute of Technology
+Copyright (c) 2009-2011 Broad Institute
+All rights reserved.
+Original author: Lee Kamentsky
+
+"""
+import numpy as np
+from scipy.ndimage import convolve, binary_erosion, generate_binary_structure
+
+def sobel(image, mask=None):
+    """Calculate the absolute magnitude Sobel to find edges.
+
+    Parameters
+    ----------
+    image : array_like, dtype=float
+        Image to process.
+    mask : array_like, dtype=bool, optional
+        An optional mask to limit the application to a certain area.
+
+    Returns
+    -------
+    output : ndarray
+      The Sobel edge map.
+
+    Notes
+    -----
+    Take the square root of the sum of the squares of the horizontal and
+    vertical Sobels to get a magnitude that's somewhat insensitive to
+    direction.
+
+    Note that ``scipy.ndimage.sobel`` returns a directional Sobel which
+    has to be further processed to perform edge detection.
+    """
+    return np.sqrt(hsobel(image, mask)**2 + vsobel(image, mask)**2)
+
+def hsobel(image, mask=None):
+    """Find the horizontal edges of an image using the Sobel transform.
+
+    Parameters
+    ----------
+    image : array_like, dtype=float
+        Image to process.
+    mask : array_like, dtype=bool, optional
+        An optional mask to limit the application to a certain area.
+
+    Returns
+    -------
+    output : ndarray
+      The Sobel edge map.
+
+    Notes
+    -----
+    We use the following kernel and return the absolute value of the
+    result at each point::
+
+      1   2   1
+      0   0   0
+     -1  -2  -1
+
+    """
+    if mask is None:
+        mask = np.ones(image.shape, bool)
+    big_mask = binary_erosion(mask,
+                              generate_binary_structure(2, 2),
+                              border_value = 0)
+    result = np.abs(convolve(image, np.array([[ 1, 2, 1],
+                                              [ 0, 0, 0],
+                                              [-1,-2,-1]]).astype(float) / 4.0))
+    result[big_mask == False] = 0
+    return result
+
+def vsobel(image, mask=None):
+    """Find the vertical edges of an image using the Sobel transform.
+
+    Parameters
+    ----------
+    image : array_like, dtype=float
+        Image to process
+    mask : array_like, dtype=bool, optional
+        An optional mask to limit the application to a certain area
+
+    Returns
+    -------
+    output : ndarray
+      The Sobel edge map.
+
+    Notes
+    -----
+    We use the following kernel and return the absolute value of the
+    result at each point::
+
+      1   0  -1
+      2   0  -2
+      1   0  -1
+
+    """
+    if mask is None:
+        mask = np.ones(image.shape, bool)
+    big_mask = binary_erosion(mask,
+                              generate_binary_structure(2, 2),
+                              border_value=0)
+    result = np.abs(convolve(image, np.array([[1, 0, -1],
+                                              [2, 0, -2],
+                                              [1, 0, -1]]).astype(float) / 4.0))
+    result[big_mask == False] = 0
+    return result
+
+def prewitt(image, mask=None):
+    """Find the edge magnitude using the Prewitt transform.
+
+    Parameters
+    ----------
+    image : array_like, dtype=float
+        Image to process.
+    mask : array_like, dtype=bool, optional
+        An optional mask to limit the application to a certain area.
+
+    Returns
+    -------
+    output : ndarray
+      The Prewitt edge map.
+
+    Notes
+    -----
+    Return the square root of the sum of squares of the horizontal
+    and vertical Prewitt transforms.
+    """
+    return np.sqrt(hprewitt(image, mask) ** 2 + vprewitt(image, mask) ** 2)
+
+def hprewitt(image, mask=None):
+    """Find the horizontal edges of an image using the Prewitt transform.
+
+    Parameters
+    ----------
+    image : array_like, dtype=float
+        Image to process.
+    mask : array_like, dtype=bool, optional
+        An optional mask to limit the application to a certain area.
+
+    Returns
+    -------
+    output : ndarray
+      The Prewitt edge map.
+
+    Notes
+    -----
+    We use the following kernel and return the absolute value of the
+    result at each point::
+
+      1   1   1
+      0   0   0
+     -1  -1  -1
+
+    """
+    if mask is None:
+        mask = np.ones(image.shape, bool)
+    big_mask = binary_erosion(mask,
+                              generate_binary_structure(2, 2),
+                              border_value=0)
+    result = np.abs(convolve(image, np.array([[ 1, 1, 1],
+                                              [ 0, 0, 0],
+                                              [-1,-1,-1]]).astype(float) / 3.0))
+    result[big_mask == False] = 0
+    return result
+
+def vprewitt(image, mask=None):
+    """Find the vertical edges of an image using the Prewitt transform.
+
+    Parameters
+    ----------
+    image : array_like, dtype=float
+        Image to process.
+    mask : array_like, dtype=bool, optional
+        An optional mask to limit the application to a certain area.
+
+    Returns
+    -------
+    output : ndarray
+      The Prewitt edge map.
+
+    Notes
+    -----
+    We use the following kernel and return the absolute value of the
+    result at each point::
+
+      1   0  -1
+      1   0  -1
+      1   0  -1
+
+    """
+    if mask is None:
+        mask = np.ones(image.shape, bool)
+    big_mask = binary_erosion(mask,
+                              generate_binary_structure(2, 2),
+                              border_value=0)
+    result = np.abs(convolve(image, np.array([[1, 0, -1],
+                                              [1, 0, -1],
+                                              [1, 0, -1]]).astype(float) / 3.0))
+    result[big_mask == False] = 0
+    return result

scikits/image/filter/tests/test_edges.py

@@ -0,0 +1,201 @@
+from numpy.testing import *
+import numpy as np
+from scipy.ndimage import binary_dilation, binary_erosion
+import scikits.image.filter as F
+
+class TestSobel():
+    def test_00_00_zeros(self):
+        """Sobel on an array of all zeros"""
+        result = F.sobel(np.zeros((10, 10)), np.ones((10, 10), bool))
+        assert (np.all(result == 0))
+
+    def test_00_01_mask(self):
+        """Sobel on a masked array should be zero"""
+        np.random.seed(0)
+        result = F.sobel(np.random.uniform(size=(10, 10)),
+                         np.zeros((10, 10), bool))
+        assert (np.all(result == 0))
+
+    def test_01_01_horizontal(self):
+        """Sobel on an edge should be a horizontal line"""
+        i, j = np.mgrid[-5:6, -5:6]
+        image = (i >= 0).astype(float)
+        result = F.sobel(image)
+        # Fudge the eroded points
+        i[np.abs(j) == 5] = 10000
+        assert (np.all(result[i == 0] == 1))
+        assert (np.all(result[np.abs(i) > 1] == 0))
+
+    def test_01_02_vertical(self):
+        """Sobel on a vertical edge should be a vertical line"""
+        i, j = np.mgrid[-5:6, -5:6]
+        image = (j >= 0).astype(float)
+        result = F.sobel(image)
+        j[np.abs(i) == 5] = 10000
+        assert (np.all(result[j == 0] == 1))
+        assert (np.all(result[np.abs(j) > 1] == 0))
+
+class TestHSobel():
+    def test_00_00_zeros(self):
+        """Horizontal sobel on an array of all zeros"""
+        result = F.hsobel(np.zeros((10, 10)), np.ones((10, 10), bool))
+        assert (np.all(result == 0))
+
+    def test_00_01_mask(self):
+        """Horizontal Sobel on a masked array should be zero"""
+        np.random.seed(0)
+        result = F.hsobel(np.random.uniform(size=(10, 10)),
+                          np.zeros((10, 10), bool))
+        assert (np.all(result == 0))
+
+    def test_01_01_horizontal(self):
+        """Horizontal Sobel on an edge should be a horizontal line"""
+        i, j = np.mgrid[-5:6, -5:6]
+        image = (i >= 0).astype(float)
+        result = F.hsobel(image)
+        # Fudge the eroded points
+        i[np.abs(j) == 5] = 10000
+        assert (np.all(result[i == 0] == 1))
+        assert (np.all(result[np.abs(i) > 1] == 0))
+
+    def test_01_02_vertical(self):
+        """Horizontal Sobel on a vertical edge should be zero"""
+        i, j = np.mgrid[-5:6, -5:6]
+        image = (j >= 0).astype(float)
+        result = F.hsobel(image)
+        assert (np.all(result == 0))
+
+class TestVSobel():
+    def test_00_00_zeros(self):
+        """Vertical sobel on an array of all zeros"""
+        result = F.vsobel(np.zeros((10, 10)), np.ones((10, 10), bool))
+        assert (np.all(result == 0))
+
+    def test_00_01_mask(self):
+        """Vertical Sobel on a masked array should be zero"""
+        np.random.seed(0)
+        result = F.vsobel(np.random.uniform(size=(10, 10)),
+                          np.zeros((10, 10), bool))
+        assert (np.all(result == 0))
+
+    def test_01_01_vertical(self):
+        """Vertical Sobel on an edge should be a vertical line"""
+        i, j = np.mgrid[-5:6, -5:6]
+        image = (j >= 0).astype(float)
+        result = F.vsobel(image)
+        # Fudge the eroded points
+        j[np.abs(i) == 5] = 10000
+        assert (np.all(result[j == 0] == 1))
+        assert (np.all(result[np.abs(j) > 1] == 0))
+
+    def test_01_02_horizontal(self):
+        """vertical Sobel on a horizontal edge should be zero"""
+        i, j = np.mgrid[-5:6, -5:6]
+        image = (i >= 0).astype(float)
+        result = F.vsobel(image)
+        eps = .000001
+        assert (np.all(np.abs(result) < eps))
+
+class TestPrewitt():
+    def test_00_00_zeros(self):
+        """Prewitt on an array of all zeros"""
+        result = F.prewitt(np.zeros((10, 10)), np.ones((10, 10), bool))
+        assert (np.all(result == 0))
+
+    def test_00_01_mask(self):
+        """Prewitt on a masked array should be zero"""
+        np.random.seed(0)
+        result = F.prewitt(np.random.uniform(size=(10, 10)),
+                           np.zeros((10, 10), bool))
+        eps = .000001
+        assert (np.all(np.abs(result) < eps))
+
+    def test_01_01_horizontal(self):
+        """Prewitt on an edge should be a horizontal line"""
+        i, j = np.mgrid[-5:6, -5:6]
+        image = (i >= 0).astype(float)
+        result = F.prewitt(image)
+        # Fudge the eroded points
+        i[np.abs(j) == 5] = 10000
+        eps = .000001
+        assert (np.all(result[i == 0] == 1))
+        assert (np.all(np.abs(result[np.abs(i) > 1]) < eps))
+
+    def test_01_02_vertical(self):
+        """Prewitt on a vertical edge should be a vertical line"""
+        i, j = np.mgrid[-5:6, -5:6]
+        image = (j >= 0).astype(float)
+        result = F.prewitt(image)
+        eps = .000001
+        j[np.abs(i)==5] = 10000
+        assert (np.all(result[j == 0] == 1))
+        assert (np.all(np.abs(result[np.abs(j) > 1]) < eps))
+
+class TestHPrewitt():
+    def test_00_00_zeros(self):
+        """Horizontal sobel on an array of all zeros"""
+        result = F.hprewitt(np.zeros((10, 10)), np.ones((10, 10), bool))
+        assert (np.all(result == 0))
+
+    def test_00_01_mask(self):
+        """Horizontal prewitt on a masked array should be zero"""
+        np.random.seed(0)
+        result = F.hprewitt(np.random.uniform(size=(10, 10)),
+                            np.zeros((10, 10), bool))
+        eps = .000001
+        assert (np.all(np.abs(result) < eps))
+
+    def test_01_01_horizontal(self):
+        """Horizontal prewitt on an edge should be a horizontal line"""
+        i, j = np.mgrid[-5:6, -5:6]
+        image = (i >= 0).astype(float)
+        result = F.hprewitt(image)
+        # Fudge the eroded points
+        i[np.abs(j) == 5] = 10000
+        eps = .000001
+        assert (np.all(result[i == 0] == 1))
+        assert (np.all(np.abs(result[np.abs(i) > 1]) < eps))
+
+    def test_01_02_vertical(self):
+        """Horizontal prewitt on a vertical edge should be zero"""
+        i, j = np.mgrid[-5:6, -5:6]
+        image = (j >= 0).astype(float)
+        result = F.hprewitt(image)
+        eps = .000001
+        assert (np.all(np.abs(result) < eps))
+
+class TestVPrewitt():
+    def test_00_00_zeros(self):
+        """Vertical prewitt on an array of all zeros"""
+        result = F.vprewitt(np.zeros((10, 10)), np.ones((10, 10), bool))
+        assert (np.all(result == 0))
+
+    def test_00_01_mask(self):
+        """Vertical prewitt on a masked array should be zero"""
+        np.random.seed(0)
+        result = F.vprewitt(np.random.uniform(size=(10, 10)),
+                            np.zeros((10, 10), bool))
+        assert (np.all(result == 0))
+
+    def test_01_01_vertical(self):
+        """Vertical prewitt on an edge should be a vertical line"""
+        i, j = np.mgrid[-5:6, -5:6]
+        image = (j >= 0).astype(float)
+        result = F.vprewitt(image)
+        # Fudge the eroded points
+        j[np.abs(i) == 5] = 10000
+        assert (np.all(result[j == 0] == 1))
+        eps = .000001
+        assert (np.all(np.abs(result[np.abs(j) > 1]) < eps))
+
+    def test_01_02_horizontal(self):
+        """Vertical prewitt on a horizontal edge should be zero"""
+        i, j = np.mgrid[-5:6, -5:6]
+        image = (i >= 0).astype(float)
+        result = F.vprewitt(image)
+        eps = .000001
+        assert (np.all(np.abs(result) < eps))
+
+
+if __name__ == "__main__":
+    run_module_suite()

scikits/image/io/collection.py

@@ -40,6 +40,8 @@ class MultiImage(object):
     The last accessed frame is cached, all other frames will have to be read
     from file.
 
+    The current implementation makes use of PIL.
+
     Examples
     --------
     >>> from scikits.image import data_dir
@@ -92,6 +94,7 @@ class MultiImage(object):
 
     def _getframe(self, framenum):
         """Open the image and extract the frame."""
+        from PIL import Image
         img = Image.open(self.filename)
         img.seek(framenum)
         return np.asarray(img, dtype=self._dtype)

scikits/image/io/tests/test_collection.py

@@ -3,11 +3,19 @@ import os.path
 
 import numpy as np
 from numpy.testing import *
+from numpy.testing.decorators import skipif
 
 from scikits.image import data_dir
 from scikits.image.io import ImageCollection, MultiImage
 
 
+try:
+    from PIL import Image
+except ImportError:
+    PIL_available = False
+else:
+    PIL_available = True
+
 if sys.version_info[0] > 2:
     basestring = str
 
@@ -58,9 +66,11 @@ class TestMultiImage():
         # convert im1.tif im2.tif -adjoin multipage.tif
         self.img = MultiImage(os.path.join(data_dir, 'multipage.tif'))
 
+    @skipif(not PIL_available)
     def test_len(self):
         assert len(self.img) == 2
 
+    @skipif(not PIL_available)
     def test_getitem(self):
         num = len(self.img)
         for i in range(-num, num):
@@ -74,6 +84,7 @@ class TestMultiImage():
         assert_raises(IndexError, return_img, num)
         assert_raises(IndexError, return_img, -num-1)
 
+    @skipif(not PIL_available)
     def test_files_property(self):
         assert isinstance(self.img.filename, basestring)
 
@@ -81,6 +92,7 @@ class TestMultiImage():
             self.img.filename = f
         assert_raises(AttributeError, set_filename, 'newfile')
 
+    @skipif(not PIL_available)
     def test_conserve_memory_property(self):
         assert isinstance(self.img.conserve_memory, bool)
 

scikits/image/morphology/__init__.py

@@ -1,2 +1,3 @@
 from grey import *
 from selem import *
+from ccomp import label

scikits/image/morphology/ccomp.pyx

@@ -0,0 +1,138 @@
+# -*- python -*-
+#cython: cdivision=True
+
+import numpy as np
+cimport numpy as np
+
+"""
+See also:
+
+  Christophe Fiorio and Jens Gustedt,
+  "Two linear time Union-Find strategies for image processing",
+  Theoretical Computer Science 154 (1996), pp. 165-181.
+
+  Kensheng Wu, Ekow Otoo and Arie Shoshani,
+  "Optimizing connected component labeling algorithms",
+  Paper LBNL-56864, 2005,
+  Lawrence Berkeley National Laboratory
+  (University of California),
+  http://repositories.cdlib.org/lbnl/LBNL-56864.
+
+"""
+
+# Tree operations implemented by an array as described in Wu et al.
+
+DTYPE = np.int
+ctypedef np.int_t DTYPE_t
+
+cdef DTYPE_t find_root(np.int_t *work, np.int_t n):
+    """Find the root of node n.
+
+    """
+    cdef np.int_t root = n
+    while (work[root] < root):
+        root = work[root]
+    return root
+
+cdef set_root(np.int_t *work, np.int_t n, np.int_t root):
+    """
+    Set all nodes on a path to point to new_root.
+
+    """
+    cdef np.int_t j
+    while (work[n] < n):
+        j = work[n]
+        work[n] = root
+        n = j
+
+    work[n] = root
+
+
+cdef join_trees(np.int_t *work, np.int_t n, np.int_t m):
+    """Join two trees containing nodes n and m.
+
+    """
+    cdef np.int_t root = find_root(work, n)
+    cdef np.int_t root_m
+
+    if (n != m):
+        root_m = find_root(work, m)
+
+        if (root > root_m):
+            root = root_m
+
+        set_root(work, n, root)
+        set_root(work, m, root)
+
+# Connected components search as described in Fiorio et al.
+
+def label(np.ndarray[DTYPE_t, ndim=2] input):
+    """Label connected regions of an integer array.
+
+    Connectivity is defined as two (8-connected) neighboring entries
+    having equal value.
+
+    Parameters
+    ----------
+    input : ndarray of dtype int
+        Image to label.
+
+    Returns
+    -------
+    labels : ndarray of dtype int
+        Labeled array, where all connected regions are assigned the
+        same integer value.
+
+    """
+    cdef np.int_t rows = input.shape[0]
+    cdef np.int_t cols = input.shape[1]
+
+    cdef np.ndarray[DTYPE_t, ndim=2] data = input.copy()
+    cdef np.ndarray[DTYPE_t, ndim=2] work
+
+    work = np.arange(data.size, dtype=DTYPE).reshape((rows, cols))
+
+    cdef np.int_t *work_p = <np.int_t*>work.data
+    cdef np.int_t *data_p = <np.int_t*>data.data
+
+    cdef np.int_t i, j
+
+    # Initialize the first row
+    for j in range(1, cols):
+        if data[0, j] == data[0, j-1]:
+            join_trees(work_p, j, j-1)
+
+    for i in range(1, rows):
+        # Handle the first column
+        if data[i, 0] == data[i-1, 0]:
+            join_trees(work_p, i*cols, (i-1)*cols)
+
+        if data[i, 0] == data[i-1, 1]:
+            join_trees(work_p, i*cols, (i-1)*cols + 1)
+
+        for j in range(1, cols):
+            if data[i, j] == data[i-1, j-1]:
+                join_trees(work_p, i*cols + j, (i-1)*cols + j - 1)
+
+            if data[i, j] == data[i-1, j]:
+                join_trees(work_p, i*cols + j, (i-1)*cols + j)
+
+            if j < cols - 1:
+                if data[i, j] == data[i - 1, j + 1]:
+                    join_trees(work_p, i*cols + j, (i-1)*cols + j + 1)
+
+            if data[i, j] == data[i, j-1]:
+                join_trees(work_p, i*cols + j, i*cols + j - 1)
+
+    # Label output
+
+    cdef np.int_t ctr = 0
+    for i in range(rows):
+        for j in range(cols):
+            if (i*cols + j) == work[i, j]:
+                data[i, j] = ctr
+                ctr = ctr + 1
+            else:
+                data[i, j] = data_p[work[i, j]]
+
+    return data

scikits/image/morphology/setup.py

@@ -14,8 +14,11 @@ def configuration(parent_package='', top_path=None):
     config = Configuration('morphology', parent_package, top_path)
     config.add_data_dir('tests')
 
+    cython(['ccomp.pyx'], working_path=base_path)
     cython(['cmorph.pyx'], working_path=base_path)
 
+    config.add_extension('ccomp', sources=['ccomp.c'],
+                         include_dirs=[get_numpy_include_dirs()])
     config.add_extension('cmorph', sources=['cmorph.c'],
                          include_dirs=[get_numpy_include_dirs()])
 

scikits/image/morphology/tests/test_ccomp.py

@@ -0,0 +1,41 @@
+import numpy as np
+from numpy.testing import assert_array_equal, run_module_suite
+
+from scikits.image.morphology import label
+
+class TestConnectedComponents:
+    def setup(self):
+        self.x = np.array([[0, 0, 3, 2, 1, 9],
+                           [0, 1, 1, 9, 2, 9],
+                           [0, 0, 1, 9, 9, 9],
+                           [3, 1, 1, 5, 3, 0]])
+
+        self.labels = np.array([[0, 0, 1, 2, 3, 4],
+                                [0, 5, 5, 4, 2, 4],
+                                [0, 0, 5, 4, 4, 4],
+                                [6, 5, 5, 7, 8, 9]])
+
+    def test_basic(self):
+        assert_array_equal(label(self.x), self.labels)
+
+        # Make sure data wasn't modified
+        assert self.x[0, 2] == 3
+
+    def test_random(self):
+        x = (np.random.random((20, 30)) * 5).astype(np.int)
+
+        labels = label(x)
+        n = labels.max()
+        for i in range(n):
+            values = x[labels == i]
+            assert np.all(values == values[0])
+
+    def test_diag(self):
+        x = np.array([[0, 0, 1],
+                      [0, 1, 0],
+                      [1, 0, 0]])
+        assert_array_equal(label(x),
+                           x)
+
+if __name__ == "__main__":
+    run_module_suite()