From fb893c277d572ec19d9b8f44d2214c5b1bbe9fbd Mon Sep 17 00:00:00 2001
From: Juan Nunez-Iglesias <juan.n@unimelb.edu.au>
Date: Sun, 11 Jan 2015 18:53:42 +1100
Subject: [PATCH] Add correct subpixel boundary estimation.

---
 skimage/segmentation/boundaries.py | 52 +++++++++++++++++++++++-------
 1 file changed, 41 insertions(+), 11 deletions(-)

diff --git a/skimage/segmentation/boundaries.py b/skimage/segmentation/boundaries.py
index 255f9816..306e213d 100644
--- a/skimage/segmentation/boundaries.py
+++ b/skimage/segmentation/boundaries.py
@@ -1,10 +1,49 @@
 import numpy as np
 from scipy import ndimage as nd
 from ..morphology import dilation, erosion, square
-from ..util import img_as_float
+from ..util import img_as_float, view_as_windows, pad
 from ..color import gray2rgb
 
 
+def _find_boundaries_subpixel(label_img):
+    """See ``find_boundaries(..., mode='subpixel')``.
+
+    Notes
+    -----
+    This function puts in an empty row and column between each *actual*
+    row and column of the image, for a corresponding shape of $2s - 1$
+    for every image dimension of size $s$. These "interstitial" rows
+    and columns are filled as ``True`` if they separate two labels in
+    `label_img`, ``False`` otherwise.
+
+    I used ``view_as_windows`` to get the neighborhood of each pixel.
+    Then I check whether there are two labels or more in that
+    neighborhood.
+    """
+    ndim = label_img.ndim
+    max_label = np.iinfo(label_img.dtype).max
+
+    label_img_expanded = np.zeros([(2 * s - 1) for s in label_img.shape],
+                                  label_img.dtype)
+    pixels = [slice(None, None, 2)] * ndim
+    label_img_expanded[pixels] = label_img
+
+    edges = np.ones(label_img_expanded.shape, dtype=bool)
+    edges[pixels] = False
+    label_img_expanded[edges] = max_label
+    windows = view_as_windows(pad(label_img_expanded, 1,
+                                  mode='constant', constant_values=0),
+                              (3,) * ndim)
+
+    boundaries = np.zeros_like(edges)
+    for index in np.ndindex(label_img_expanded.shape):
+        if edges[index]:
+            values = np.unique(windows[index].ravel())
+            if len(values) > 2: # single value and max_label
+                boundaries[index] = True
+    return boundaries
+
+
 def find_boundaries(label_img, connectivity=1, mode='thick', background=0):
     """Return bool array where boundaries between labeled regions are True.
 
@@ -124,16 +163,7 @@ def find_boundaries(label_img, connectivity=1, mode='thick', background=0):
             boundaries &= (background_image | adjacent_objects)
         return boundaries
     else:
-        label_img_expanded = np.zeros([(2 * s - 1) for s in label_img.shape],
-                                      label_img.dtype)
-        pixels = [slice(None, None, 2)] * ndim
-        selem = nd.generate_binary_structure(ndim, ndim)
-        label_img_expanded[pixels] = label_img
-        max_label = np.iinfo(label_img.dtype).max
-        label_img_edge_inverted = np.array(label_img_expanded, copy=True)
-        label_img_edge_inverted[label_img_expanded == 0] = max_label
-        boundaries = (dilation(label_img_expanded, selem) !=
-                      erosion(label_img_edge_inverted, selem))
+        boundaries = _find_boundaries_subpixel(label_img)
         return boundaries