Normalized Cut on RAGs (WIP)

bento.info

@@ -154,6 +154,9 @@ Library:
     Extension: skimage.feature._hessian_det_appx
         Sources:
             skimage/exposure/_hessian_det_appx.pyx
+    Extension: skimage.graph._ncut_cy
+        Sources:
+            skimage/graph/_ncut_cy.pyx
 
 Executable: skivi
     Module: skimage.scripts.skivi

doc/examples/plot_ncut.py

@@ -0,0 +1,32 @@
+"""
+==============
+Normalized Cut
+==============
+
+This example constructs a Region Adjacency Graph (RAG) and recursively performs
+a Normalized Cut on it.
+
+References
+----------
+.. [1] Shi, J.; Malik, J., "Normalized cuts and image segmentation",
+       Pattern Analysis and Machine Intelligence,
+       IEEE Transactions on , vol.22, no.8, pp.888,905, Aug 2000
+"""
+from skimage import graph, data, io, segmentation, color
+from matplotlib import pyplot as plt
+
+
+img = data.coffee()
+
+labels1 = segmentation.slic(img, compactness=30, n_segments=400)
+out1 = color.label2rgb(labels1, img, kind='avg')
+
+g = graph.rag_mean_color(img, labels1, mode='similarity')
+labels2 = graph.cut_normalized(labels1, g)
+out2 = color.label2rgb(labels2, img, kind='avg')
+
+plt.figure()
+io.imshow(out1)
+plt.figure()
+io.imshow(out2)
+io.show()

skimage/graph/__init__.py

@@ -1,7 +1,8 @@
 from .spath import shortest_path
 from .mcp import MCP, MCP_Geometric, MCP_Connect, MCP_Flexible, route_through_array
 from .rag import rag_mean_color, RAG
-from .graph_cut import cut_threshold
+from .graph_cut import cut_threshold, cut_normalized
+ncut = cut_normalized
 
 __all__ = ['shortest_path',
            'MCP',
@@ -11,4 +12,6 @@
            'route_through_array',
            'rag_mean_color',
            'cut_threshold',
+           'cut_normalized',
+           'ncut',
            'RAG']

skimage/graph/_ncut.py

@@ -0,0 +1,85 @@
+try:
+    import networkx as nx
+except ImportError:
+    import warnings
+    warnings.warn('RAGs require networkx')
+import numpy as np
+from scipy import sparse
+from . import _ncut_cy
+
+
+def DW_matrices(graph):
+    """Returns the diagonal and weight matrices of a graph.
+
+    Parameters
+    ----------
+    graph : RAG
+        A Region Adjacency Graph.
+
+    Returns
+    -------
+    D : csc_matrix
+        The diagonal matrix of the graph. ``D[i, i]`` is the sum of weights of
+        all edges incident on `i`. All other enteries are `0`.
+    W : csc_matrix
+        The weight matrix of the graph. ``W[i, j]`` is the weight of the edge
+        joining `i` to `j`.
+    """
+    # sparse.eighsh is most efficient with CSC-formatted input
+    W = nx.to_scipy_sparse_matrix(graph, format='csc')
+    entries = W.sum(axis=0)
+    D = sparse.dia_matrix((entries, 0), shape=W.shape).tocsc()
+    return D, W
+
+
+def ncut_cost(cut, D, W):
+    """Returns the N-cut cost of a bi-partition of a graph.
+
+    Parameters
+    ----------
+    cut : ndarray
+        The mask for the nodes in the graph. Nodes corressponding to a `True`
+        value are in one set.
+    D : csc_matrix
+        The diagonal matrix of the graph.
+    W : csc_matrix
+        The weight matrix of the graph.
+
+    Returns
+    -------
+    cost : float
+        The cost of performing the N-cut.
+
+    References
+    ----------
+    .. [1] Normalized Cuts and Image Segmentation, Jianbo Shi and
+           Jitendra Malik, IEEE Transactions on Pattern Analysis and Machine
+           Intelligence, Page 889, Equation 2.
+    """
+    cut = np.array(cut)
+    cut_cost = _ncut_cy.cut_cost(cut, W)
+
+    # D has elements only along the diagonal, one per node, so we can directly
+    # index the data attribute with cut.
+    assoc_a = D.data[cut].sum()
+    assoc_b = D.data[~cut].sum()
+
+    return (cut_cost / assoc_a) + (cut_cost / assoc_b)
+
+
+def normalize(a):
+    """Normalize values in an array between `0` and `1`.
+
+    Parameters
+    ----------
+    a : ndarray
+        The array to be normalized.
+
+    Returns
+    -------
+    out : ndarray
+        The normalized array.
+    """
+    mi = a.min()
+    mx = a.max()
+    return (a - mi) / (mx - mi)

skimage/graph/_ncut_cy.pyx

@@ -0,0 +1,78 @@
+# cython: cdivision=True
+# cython: boundscheck=False
+# cython: nonecheck=False
+# cython: wraparound=False
+cimport numpy as cnp
+import numpy as np
+
+
+def argmin2(cnp.double_t[:] array):
+    """Return the index of the 2nd smallest value in an array.
+
+    Parameters
+    ----------
+    a : array
+        The array to process.
+
+    Returns
+    -------
+    min_idx2 : int
+        The index of the second smallest value.
+    """
+    cdef cnp.float64_t min1 = np.inf
+    cdef cnp.float64_t min2 = np.inf
+    cdef Py_ssize_t min_idx1 = 0
+    cdef Py_ssize_t min_idx2 = 0
+    cdef Py_ssize_t i = 0
+    cdef Py_ssize_t n = array.shape[0]
+
+    for i in range(n):
+        x = array[i]
+        if x < min1:
+            min2 = min1
+            min_idx2 = min_idx1
+            min1 = x
+            min_idx1 = i
+        elif x > min1 and x < min2:
+            min2 = x
+            min_idx2 = i
+        i += 1
+
+    return min_idx2
+
+
+def cut_cost(cut, W):
+    """Return the total weight of crossing edges in a bi-partition.
+
+    Parameters
+    ----------
+    cut : array
+        A array of booleans. Elements set to `True` belong to one
+        set.
+    W : array
+        The weight matrix of the graph.
+
+    Returns
+    -------
+    cost : float
+        The total weight of crossing edges.
+    """
+    cdef cnp.ndarray[cnp.uint8_t, cast = True] cut_mask = np.array(cut)
+    cdef Py_ssize_t num_rows, num_cols
+    cdef cnp.int32_t row, col
+    cdef cnp.int32_t[:] indices = W.indices
+    cdef cnp.int32_t[:] indptr = W.indptr
+    cdef cnp.double_t[:] data = W.data.astype(np.double)
+    cdef cnp.int32_t row_index
+    cdef cnp.double_t cost = 0
+
+    num_rows = W.shape[0]
+    num_cols = W.shape[1]
+
+    for col in range(num_cols):
+        for row_index in range(indptr[col], indptr[col + 1]):
+            row = indices[row_index]
+            if cut_mask[row] != cut_mask[col]:
+                cost += data[row_index]
+
+    return cost * 0.5