Merge pull request #2134 from jmetz/felzenszwalb_multichannel

Switched felzenszwalb gray to multichannel version

skimage/segmentation/_felzenszwalb.py

@@ -1,7 +1,7 @@
 import numpy as np
 
 from .._shared.utils import warn
-from ._felzenszwalb_cy import _felzenszwalb_grey
+from ._felzenszwalb_cy import _felzenszwalb_cython
 
 
 def felzenszwalb(image, scale=1, sigma=0.8, min_size=20):
@@ -17,9 +17,8 @@ def felzenszwalb(image, scale=1, sigma=0.8, min_size=20):
     controlled indirectly through ``scale``. Segment size within an image can
     vary greatly depending on local contrast.
 
-    For RGB images, the algorithm computes a separate segmentation for each
-    channel and then combines these. The combined segmentation is the
-    intersection of the separate segmentations on the color channels.
+    For RGB images, the algorithm uses the euclidean distance between pixels in
+    color space.
 
     Parameters
     ----------
@@ -50,35 +49,6 @@ def felzenszwalb(image, scale=1, sigma=0.8, min_size=20):
     >>> segments = felzenszwalb(img, scale=3.0, sigma=0.95, min_size=5)
     """
 
-    if image.ndim == 2:
-        # assume single channel image
-        return _felzenszwalb_grey(image, scale=scale, sigma=sigma,
-                                  min_size=min_size)
-
-    elif image.ndim != 3:
-        raise ValueError("Felzenswalb segmentation can only operate on RGB and"
-                         " grey images, but input array of ndim %d given."
-                         % image.ndim)
-
-    # assume we got 2d image with multiple channels
-    n_channels = image.shape[2]
-    if n_channels != 3:
-        warn("Got image with %d channels. Is that really what you"
-             " wanted?" % image.shape[2])
-    segmentations = []
-    # compute quickshift for each channel
-    for c in range(n_channels):
-        channel = np.ascontiguousarray(image[:, :, c])
-        s = _felzenszwalb_grey(channel, scale=scale, sigma=sigma,
-                               min_size=min_size)
-        segmentations.append(s)
-
-    # put pixels in same segment only if in the same segment in all images
-    # we do this by combining the channels to one number
-    n0 = segmentations[0].max() + 1
-    n1 = segmentations[1].max() + 1
-    segmentation = (segmentations[0] + segmentations[1] * n0
-                    + segmentations[2] * n0 * n1)
-    # make segment labels consecutive numbers starting at 0
-    labels = np.unique(segmentation, return_inverse=True)[1]
-    return labels.reshape(image.shape[:2])
+    image = np.atleast_3d(image)
+    return _felzenszwalb_cython(image, scale=scale, sigma=sigma,
+                                min_size=min_size)

skimage/segmentation/_felzenszwalb_cy.pyx

@@ -12,50 +12,56 @@ from ..measure._ccomp cimport find_root, join_trees
 from ..util import img_as_float
 
 
-def _felzenszwalb_grey(image, double scale=1, sigma=0.8,
-                       Py_ssize_t min_size=20):
-    """Felzenszwalb's efficient graph based segmentation for a single channel.
+def _felzenszwalb_cython(image, double scale=1, sigma=0.8,
+                         Py_ssize_t min_size=20):
+    """Felzenszwalb's efficient graph based segmentation for 
+    single or multiple channels.
 
-    Produces an oversegmentation of a 2d image using a fast, minimum spanning
-    tree based clustering on the image grid.
+    Produces an oversegmentation of a single or multi-channel image 
+    using a fast, minimum spanning tree based clustering on the image grid.
     The number of produced segments as well as their size can only be
     controlled indirectly through ``scale``. Segment size within an image can
     vary greatly depending on local contrast.
 
     Parameters
     ----------
-    image: ndarray
+    image : (N, M, C) ndarray
         Input image.
-    scale: float, optional (default 1)
+    scale : float, optional (default 1)
         Sets the obervation level. Higher means larger clusters.
-    sigma: float, optional (default 0.8)
+    sigma : float, optional (default 0.8)
         Width of Gaussian smoothing kernel used in preprocessing.
         Larger sigma gives smother segment boundaries.
-    min_size: int, optional (default 20)
+    min_size : int, optional (default 20)
         Minimum component size. Enforced using postprocessing.
 
     Returns
     -------
-    segment_mask: (height, width) ndarray
+    segment_mask : (N, M) ndarray
         Integer mask indicating segment labels.
     """
-    if image.ndim != 2:
-        raise ValueError("This algorithm works only on single-channel 2d"
-                "images. Got image of shape %s" % str(image.shape))
+    if image.ndim != 3:
+        raise ValueError("This algorithm works only on single or " 
+                         "multi-channel 2d images. "
+                         "Got image of shape %s" % str(image.shape))
 
     image = img_as_float(image)
 
     # rescale scale to behave like in reference implementation
     scale = float(scale) / 255.
-    image = ndi.gaussian_filter(image, sigma=sigma)
+    image = ndi.gaussian_filter(image, sigma=[sigma, sigma, 0])
 
     # compute edge weights in 8 connectivity:
-    right_cost = np.abs((image[1:, :] - image[:-1, :]))
-    down_cost = np.abs((image[:, 1:] - image[:, :-1]))
-    dright_cost = np.abs((image[1:, 1:] - image[:-1, :-1]))
-    uright_cost = np.abs((image[1:, :-1] - image[:-1, 1:]))
-    cdef cnp.ndarray[cnp.float_t, ndim=1] costs = np.hstack([right_cost.ravel(),
-        down_cost.ravel(), dright_cost.ravel(),
+    right_cost = np.sqrt(np.sum((image[1:, :, :] - image[:-1, :, :])
+    	*(image[1:, :, :] - image[:-1, :, :]), axis=-1))
+    down_cost = np.sqrt(np.sum((image[:, 1:, :] - image[:, :-1, :])
+    	*(image[:, 1:, :] - image[:, :-1, :]), axis=-1))
+    dright_cost = np.sqrt(np.sum((image[1:, 1:, :] - image[:-1, :-1, :])
+	*(image[1:, 1:, :] - image[:-1, :-1, :]), axis=-1))
+    uright_cost = np.sqrt(np.sum((image[1:, :-1, :] - image[:-1, 1:, :])
+    	*(image[1:, :-1, :] - image[:-1, 1:, :]), axis=-1))
+    cdef cnp.ndarray[cnp.float_t, ndim=1] costs = np.hstack([
+    	right_cost.ravel(), down_cost.ravel(), dright_cost.ravel(),
         uright_cost.ravel()]).astype(np.float)
 
     # compute edges between pixels:

skimage/segmentation/tests/test_felzenszwalb.py

@@ -34,9 +34,7 @@ def test_minsize():
         segments = felzenszwalb(coffee, min_size=min_size, sigma=3)
         counts = np.bincount(segments.ravel())
         # actually want to test greater or equal.
-        # the construction doesn't guarantee min_size is respected
-        # after intersecting the sementations for the colors
-        assert_greater(np.mean(counts) + 1, min_size)
+        assert_greater(counts.min() + 1, min_size)
 
 
 def test_color():