From 7215936f5c39a85ba43cb0b93574141c7e6e69f3 Mon Sep 17 00:00:00 2001 From: Vighnesh Birodkar Date: Wed, 23 Jul 2014 02:14:46 +0530 Subject: [PATCH 01/36] First working copy --- skimage/graph/__init__.py | 3 +- skimage/graph/_ncut.py | 38 +++++++++++++++++++++ skimage/graph/_ncut_cy.pyx | 27 +++++++++++++++ skimage/graph/graph_cut.py | 70 +++++++++++++++++++++++++++++++++++++- skimage/graph/rag.py | 49 ++++++++++++++++++++++++++ skimage/graph/setup.py | 4 ++- 6 files changed, 188 insertions(+), 3 deletions(-) create mode 100644 skimage/graph/_ncut.py create mode 100644 skimage/graph/_ncut_cy.pyx diff --git a/skimage/graph/__init__.py b/skimage/graph/__init__.py index 68c1206c..9fe00d3c 100644 --- a/skimage/graph/__init__.py +++ b/skimage/graph/__init__.py @@ -1,7 +1,7 @@ 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_n __all__ = ['shortest_path', 'MCP', @@ -11,4 +11,5 @@ __all__ = ['shortest_path', 'route_through_array', 'rag_mean_color', 'cut_threshold', + 'cut_n', 'RAG'] diff --git a/skimage/graph/_ncut.py b/skimage/graph/_ncut.py new file mode 100644 index 00000000..9b321bd8 --- /dev/null +++ b/skimage/graph/_ncut.py @@ -0,0 +1,38 @@ +import networkx as nx +import numpy as np +from scipy import sparse + +def DW_matrix(graph): + #print graph[4][4]['weight'] + W = nx.to_scipy_sparse_matrix(graph, format='csc') + entries = W.sum(0) + D = sparse.dia_matrix( (entries,0),shape = W.shape).tocsc() + #print W[4,4] + + m,n = W.shape + #for i in range(n): + # W[i,i] = 1.0 + return D,W + +def ncut_cost(mask,D,W): + + mask = np.array(mask) + mask_list = [ np.logical_xor(mask[i], mask) for i in range(mask.shape[0])] + mask_array = np.array(mask_list) + + cut = float(W[mask_array].sum()/2.0) + #print W.todense() + #print mask_array.astype(int) + #print "cut=",cut + + assoc_a = D.data[mask].sum() + assoc_b = D.data[np.logical_not(mask)].sum() + + #print cut + #print assoc_a,assoc_b + return (cut/assoc_a) + (cut/assoc_b) + +def norml(a): + mi = a.min() + mx = a.max() + return (a-mi)/(mx-mi) diff --git a/skimage/graph/_ncut_cy.pyx b/skimage/graph/_ncut_cy.pyx new file mode 100644 index 00000000..fc3c0f3b --- /dev/null +++ b/skimage/graph/_ncut_cy.pyx @@ -0,0 +1,27 @@ +# cython: cdivision=True +# cython: boundscheck=False +# cython: nonecheck=False +# cython: wraparound=False +cimport numpy as cnp +import numpy as np + +def argmin2(cnp.float64_t[:] array): + cdef cnp.float64_t min1 = np.inf + cdef cnp.float64_t min2 = np.inf + cdef Py_ssize_t i1 = 0 + cdef Py_ssize_t i2 = 0 + cdef Py_ssize_t i = 0 + + while i < array.shape[0]: + x = array[i] + if x < min1 : + min2 = min1 + i2 = i1 + min1 = x + i1 = i + elif x > min1 and x < min2 : + min2 = x + i2 = i + i += 1 + + return i2 diff --git a/skimage/graph/graph_cut.py b/skimage/graph/graph_cut.py index d9fcfdef..15f6ff34 100644 --- a/skimage/graph/graph_cut.py +++ b/skimage/graph/graph_cut.py @@ -4,7 +4,11 @@ except ImportError: import warnings warnings.warn('"cut_threshold" requires networkx') import numpy as np - +import _ncut +import _ncut_cy +from scipy.sparse import linalg +from scipy.sparse.linalg.eigen.arpack.arpack import ArpackNoConvergence as ANC +from scipy.sparse.linalg.eigen.arpack.arpack import ArpackError as APE def cut_threshold(labels, rag, thresh): """Combine regions seperated by weight less than threshold. @@ -60,3 +64,67 @@ def cut_threshold(labels, rag, thresh): map_array[label] = i return map_array[labels] + + +def cut_n(labels, rag, thresh): + + _ncut_relabel(rag,thresh) + + from_ = range(labels.max()+1) + to = [ rag.node[x]['ncut label'] for x in from_ ] + map_array = np.array(to) + + return map_array[labels] + +def _ncut_relabel(rag, cut_thresh = 0.0001): + d, w = _ncut.DW_matrix(rag) + error = False + + try: + m = w.shape[0] + vals,vectors = linalg.eigsh(d-w,M=d,which='SM',k = min(100,m-2)) + except ANC as e: + vals = e.eigenvalues + vectors = e.eigenvectors + if len(vals) == 0: + error = True + except ValueError: + error = True + except APE: + error = True + + if not error : + vals,vectors = np.real(vals), np.real(vectors) + index2 = _ncut_cy.argmin2(vals) + + ev = np.real(vectors[:,index2]) + ev = _ncut.norml(ev) + + mcut = np.inf + thresh = None + for t in np.arange(0,1,0.1): + mask = ev > t + cost = _ncut.ncut_cost(mask,d,w) + if cost < mcut : + mcut = cost + thresh = t + + if ( mcut < cut_thresh ): + mask = ev > thresh + + nodes1 = [ n for i,n in enumerate(rag.nodes()) if mask[i]] + nodes2 = [ n for i,n in enumerate(rag.nodes()) if not mask[i]] + + sub1 = rag.subgraph(nodes1) + sub2 = rag.subgraph(nodes2) + + _ncut_relabel(sub1,cut_thresh) + _ncut_relabel(sub2, cut_thresh) + return + + node = rag.nodes()[0] + new_label = rag.node[node]['labels'][0] + for n in rag.nodes(): + rag.node[n]['ncut label'] = new_label + + diff --git a/skimage/graph/rag.py b/skimage/graph/rag.py index 7c1b7f07..4592ea36 100644 --- a/skimage/graph/rag.py +++ b/skimage/graph/rag.py @@ -11,6 +11,7 @@ except ImportError: import numpy as np from scipy.ndimage import filters from scipy import ndimage as nd +import math def min_weight(graph, src, dst, n): @@ -210,3 +211,51 @@ def rag_mean_color(image, labels, connectivity=2): graph[x][y]['weight'] = np.linalg.norm(diff) return graph + +def rag_similarity(image, labels, connectivity=2, sigma = 30.0): + graph = RAG() + + # The footprint is constructed in such a way that the first + # element in the array being passed to _add_edge_filter is + # the central value. + fp = nd.generate_binary_structure(labels.ndim, connectivity) + for d in range(fp.ndim): + fp = fp.swapaxes(0, d) + fp[0, ...] = 0 + fp = fp.swapaxes(0, d) + + + filters.generic_filter( + labels, + function=_add_edge_filter, + footprint=fp, + mode='nearest', + output=np.zeros(labels.shape, dtype=np.uint8), + extra_arguments=(graph,)) + + for n in graph: + graph.node[n].update({'labels': [n], + 'pixel count': 0, + 'total color': np.array([0, 0, 0], + dtype=np.double)}) + + for index in np.ndindex(labels.shape): + current = labels[index] + graph.node[current]['pixel count'] += 1 + graph.node[current]['total color'] += image[index] + + for n in graph: + graph.node[n]['mean color'] = (graph.node[n]['total color'] / + graph.node[n]['pixel count']) + + for x, y in graph.edges_iter(): + diff = graph.node[x]['mean color'] - graph.node[y]['mean color'] + diff = np.linalg.norm(diff) + + #if diff == 0: + # graph[x][y]['weight'] = 99999 + #else: + graph[x][y]['weight'] = math.e**(-(diff**2)/sigma) + #print graph[x][y]['weight'] + #print "diff",diff,"weight",math.e**(-(diff**2)/sigma) + return graph diff --git a/skimage/graph/setup.py b/skimage/graph/setup.py index 463d2739..edc7b653 100644 --- a/skimage/graph/setup.py +++ b/skimage/graph/setup.py @@ -17,6 +17,7 @@ def configuration(parent_package='', top_path=None): cython(['_spath.pyx'], working_path=base_path) cython(['_mcp.pyx'], working_path=base_path) cython(['heap.pyx'], working_path=base_path) + cython(['_ncut_cy.pyx'], working_path=base_path) config.add_extension('_spath', sources=['_spath.c'], include_dirs=[get_numpy_include_dirs()]) @@ -24,7 +25,8 @@ def configuration(parent_package='', top_path=None): include_dirs=[get_numpy_include_dirs()]) config.add_extension('heap', sources=['heap.c'], include_dirs=[get_numpy_include_dirs()]) - + config.add_extension('_ncut_cy', sources=['_ncut_cy.c'], + include_dirs=[get_numpy_include_dirs()]) return config if __name__ == '__main__': From ddf11bdfe85590a1319790ed97a1832bb3c789c2 Mon Sep 17 00:00:00 2001 From: Vighnesh Birodkar Date: Wed, 23 Jul 2014 02:19:47 +0530 Subject: [PATCH 02/36] pep8 changes --- skimage/graph/graph_cut.py | 53 +++++++++++++++++++------------------- 1 file changed, 27 insertions(+), 26 deletions(-) diff --git a/skimage/graph/graph_cut.py b/skimage/graph/graph_cut.py index 15f6ff34..f1971198 100644 --- a/skimage/graph/graph_cut.py +++ b/skimage/graph/graph_cut.py @@ -4,11 +4,12 @@ except ImportError: import warnings warnings.warn('"cut_threshold" requires networkx') import numpy as np -import _ncut -import _ncut_cy +from . import _ncut +from . import _ncut_cy from scipy.sparse import linalg -from scipy.sparse.linalg.eigen.arpack.arpack import ArpackNoConvergence as ANC -from scipy.sparse.linalg.eigen.arpack.arpack import ArpackError as APE +from scipy.sparse.linalg.eigen.arpack.arpack import ArpackNoConvergence +from scipy.sparse.linalg.eigen.arpack.arpack import ArpackError + def cut_threshold(labels, rag, thresh): """Combine regions seperated by weight less than threshold. @@ -68,57 +69,59 @@ def cut_threshold(labels, rag, thresh): def cut_n(labels, rag, thresh): - _ncut_relabel(rag,thresh) - - from_ = range(labels.max()+1) - to = [ rag.node[x]['ncut label'] for x in from_ ] + _ncut_relabel(rag, thresh) + + from_ = range(labels.max() + 1) + to = [rag.node[x]['ncut label'] for x in from_] map_array = np.array(to) - + return map_array[labels] -def _ncut_relabel(rag, cut_thresh = 0.0001): + +def _ncut_relabel(rag, cut_thresh=0.0001): d, w = _ncut.DW_matrix(rag) error = False try: m = w.shape[0] - vals,vectors = linalg.eigsh(d-w,M=d,which='SM',k = min(100,m-2)) - except ANC as e: + vals, vectors = linalg.eigsh(d - w, M=d, which='SM', + k=min(100, m - 2)) + except ArpackNoConvergence as e: vals = e.eigenvalues vectors = e.eigenvectors if len(vals) == 0: error = True except ValueError: error = True - except APE: + except ArpackError: error = True - - if not error : - vals,vectors = np.real(vals), np.real(vectors) + + if not error: + vals, vectors = np.real(vals), np.real(vectors) index2 = _ncut_cy.argmin2(vals) - ev = np.real(vectors[:,index2]) + ev = np.real(vectors[:, index2]) ev = _ncut.norml(ev) mcut = np.inf thresh = None - for t in np.arange(0,1,0.1): + for t in np.arange(0, 1, 0.1): mask = ev > t - cost = _ncut.ncut_cost(mask,d,w) - if cost < mcut : + cost = _ncut.ncut_cost(mask, d, w) + if cost < mcut: mcut = cost thresh = t - if ( mcut < cut_thresh ): + if (mcut < cut_thresh): mask = ev > thresh - nodes1 = [ n for i,n in enumerate(rag.nodes()) if mask[i]] - nodes2 = [ n for i,n in enumerate(rag.nodes()) if not mask[i]] + nodes1 = [n for i, n in enumerate(rag.nodes()) if mask[i]] + nodes2 = [n for i, n in enumerate(rag.nodes()) if not mask[i]] sub1 = rag.subgraph(nodes1) sub2 = rag.subgraph(nodes2) - _ncut_relabel(sub1,cut_thresh) + _ncut_relabel(sub1, cut_thresh) _ncut_relabel(sub2, cut_thresh) return @@ -126,5 +129,3 @@ def _ncut_relabel(rag, cut_thresh = 0.0001): new_label = rag.node[node]['labels'][0] for n in rag.nodes(): rag.node[n]['ncut label'] = new_label - - From 0971c11db6a662fd6be77e6843a51a700087337e Mon Sep 17 00:00:00 2001 From: Vighnesh Birodkar Date: Wed, 23 Jul 2014 02:27:33 +0530 Subject: [PATCH 03/36] Cleaned up rag.py --- skimage/graph/rag.py | 47 ++++---------------------------------------- 1 file changed, 4 insertions(+), 43 deletions(-) diff --git a/skimage/graph/rag.py b/skimage/graph/rag.py index 4592ea36..f4db630d 100644 --- a/skimage/graph/rag.py +++ b/skimage/graph/rag.py @@ -211,51 +211,12 @@ def rag_mean_color(image, labels, connectivity=2): graph[x][y]['weight'] = np.linalg.norm(diff) return graph - -def rag_similarity(image, labels, connectivity=2, sigma = 30.0): - graph = RAG() - - # The footprint is constructed in such a way that the first - # element in the array being passed to _add_edge_filter is - # the central value. - fp = nd.generate_binary_structure(labels.ndim, connectivity) - for d in range(fp.ndim): - fp = fp.swapaxes(0, d) - fp[0, ...] = 0 - fp = fp.swapaxes(0, d) - filters.generic_filter( - labels, - function=_add_edge_filter, - footprint=fp, - mode='nearest', - output=np.zeros(labels.shape, dtype=np.uint8), - extra_arguments=(graph,)) +def rag_similarity(image, labels, connectivity=2, sigma=30.0): + graph = rag_mean_color(image, labels, connectivity) - for n in graph: - graph.node[n].update({'labels': [n], - 'pixel count': 0, - 'total color': np.array([0, 0, 0], - dtype=np.double)}) + for x, w, d in graph.edges_iter(data=True): + d['weight'] = math.e ** (-(d['weight'] ** 2) / sigma) - for index in np.ndindex(labels.shape): - current = labels[index] - graph.node[current]['pixel count'] += 1 - graph.node[current]['total color'] += image[index] - - for n in graph: - graph.node[n]['mean color'] = (graph.node[n]['total color'] / - graph.node[n]['pixel count']) - - for x, y in graph.edges_iter(): - diff = graph.node[x]['mean color'] - graph.node[y]['mean color'] - diff = np.linalg.norm(diff) - - #if diff == 0: - # graph[x][y]['weight'] = 99999 - #else: - graph[x][y]['weight'] = math.e**(-(diff**2)/sigma) - #print graph[x][y]['weight'] - #print "diff",diff,"weight",math.e**(-(diff**2)/sigma) return graph From 997d4a3a682e4998539ecbfbfac197293b82914b Mon Sep 17 00:00:00 2001 From: Vighnesh Birodkar Date: Wed, 23 Jul 2014 02:37:15 +0530 Subject: [PATCH 04/36] cleaned up _ncut.py --- skimage/graph/_ncut.py | 36 ++++++++++++++---------------------- skimage/graph/_ncut_cy.pyx | 7 ++++--- 2 files changed, 18 insertions(+), 25 deletions(-) diff --git a/skimage/graph/_ncut.py b/skimage/graph/_ncut.py index 9b321bd8..e31b9431 100644 --- a/skimage/graph/_ncut.py +++ b/skimage/graph/_ncut.py @@ -2,37 +2,29 @@ import networkx as nx import numpy as np from scipy import sparse + def DW_matrix(graph): - #print graph[4][4]['weight'] + W = nx.to_scipy_sparse_matrix(graph, format='csc') entries = W.sum(0) - D = sparse.dia_matrix( (entries,0),shape = W.shape).tocsc() - #print W[4,4] - - m,n = W.shape - #for i in range(n): - # W[i,i] = 1.0 - return D,W + D = sparse.dia_matrix((entries, 0), shape=W.shape).tocsc() + return D, W -def ncut_cost(mask,D,W): + +def ncut_cost(mask, D, W): mask = np.array(mask) - mask_list = [ np.logical_xor(mask[i], mask) for i in range(mask.shape[0])] + mask_list = [np.logical_xor(mask[i], mask) for i in range(mask.shape[0])] mask_array = np.array(mask_list) - - cut = float(W[mask_array].sum()/2.0) - #print W.todense() - #print mask_array.astype(int) - #print "cut=",cut - + + cut = float(W[mask_array].sum() / 2.0) assoc_a = D.data[mask].sum() assoc_b = D.data[np.logical_not(mask)].sum() - - #print cut - #print assoc_a,assoc_b - return (cut/assoc_a) + (cut/assoc_b) -def norml(a): + return (cut / assoc_a) + (cut / assoc_b) + + +def normalize(a): mi = a.min() mx = a.max() - return (a-mi)/(mx-mi) + return (a - mi) / (mx - mi) diff --git a/skimage/graph/_ncut_cy.pyx b/skimage/graph/_ncut_cy.pyx index fc3c0f3b..77cb18ea 100644 --- a/skimage/graph/_ncut_cy.pyx +++ b/skimage/graph/_ncut_cy.pyx @@ -5,6 +5,7 @@ cimport numpy as cnp import numpy as np + def argmin2(cnp.float64_t[:] array): cdef cnp.float64_t min1 = np.inf cdef cnp.float64_t min2 = np.inf @@ -14,14 +15,14 @@ def argmin2(cnp.float64_t[:] array): while i < array.shape[0]: x = array[i] - if x < min1 : + if x < min1: min2 = min1 i2 = i1 min1 = x i1 = i - elif x > min1 and x < min2 : + elif x > min1 and x < min2: min2 = x i2 = i i += 1 - + return i2 From 68b087ca48ee2fbd9a1f0ef2e303ff44f03a7528 Mon Sep 17 00:00:00 2001 From: Vighnesh Birodkar Date: Wed, 23 Jul 2014 03:19:06 +0530 Subject: [PATCH 05/36] Docstrings --- skimage/graph/graph_cut.py | 2 +- skimage/graph/rag.py | 43 +++++++++++++++++++++++++------------- 2 files changed, 30 insertions(+), 15 deletions(-) diff --git a/skimage/graph/graph_cut.py b/skimage/graph/graph_cut.py index f1971198..e49d1182 100644 --- a/skimage/graph/graph_cut.py +++ b/skimage/graph/graph_cut.py @@ -101,7 +101,7 @@ def _ncut_relabel(rag, cut_thresh=0.0001): index2 = _ncut_cy.argmin2(vals) ev = np.real(vectors[:, index2]) - ev = _ncut.norml(ev) + ev = _ncut.normalize(ev) mcut = np.inf thresh = None diff --git a/skimage/graph/rag.py b/skimage/graph/rag.py index f4db630d..0acede62 100644 --- a/skimage/graph/rag.py +++ b/skimage/graph/rag.py @@ -119,14 +119,15 @@ def _add_edge_filter(values, graph): return 0 -def rag_mean_color(image, labels, connectivity=2): +def rag_mean_color(image, labels, connectivity=2, mode='dissimilarity', + sigma=30.0): """Compute the Region Adjacency Graph using mean colors. Given an image and its initial segmentation, this method constructs the corresponsing Region Adjacency Graph (RAG). Each node in the RAG represents a set of pixels within `image` with the same label in `labels`. - The weight between two adjacent regions is the difference in their mean - color. + The weight between two adjacent regions represents how similar or + dissimilar two regions are depending on the `mode` parameter. Parameters ---------- @@ -141,6 +142,23 @@ def rag_mean_color(image, labels, connectivity=2): are considered adjacent. It can range from 1 to `labels.ndim`. Its behavior is the same as `connectivity` parameter in `scipy.ndimage.filters.generate_binary_structure`. + mode : str['similarity' | 'dissimilarity'] + The strategy to assign edge weights. + + 'similarity' : The weight between two adjacent regions is the + :math:`|c_1 - c_2|`, where :math:`c1` and :math:`c2` are the mean + colors of the two regions. It represents how different two regions + are. + + 'dissimilarity' : The weight between two adjacent is + :math:`e^{-d^2/sigma}` where :math:`d=|c_1 - c_2|`, where + :math:`c1` and :math:`c2` are the mean colors of the two regions. + It represents how similar two regions are. + sigma : float, optional + Used for computation when `mode='dissimilarity'`. It governs how close + to each other, two colors should be, for their corresponding edge + weight to be significant. A very large value of `sigma` could make + any two colors behave as though they were similar. Returns ------- @@ -206,17 +224,14 @@ def rag_mean_color(image, labels, connectivity=2): graph.node[n]['mean color'] = (graph.node[n]['total color'] / graph.node[n]['pixel count']) - for x, y in graph.edges_iter(): + for x, y, d in graph.edges_iter(data=True): diff = graph.node[x]['mean color'] - graph.node[y]['mean color'] - graph[x][y]['weight'] = np.linalg.norm(diff) - - return graph - - -def rag_similarity(image, labels, connectivity=2, sigma=30.0): - graph = rag_mean_color(image, labels, connectivity) - - for x, w, d in graph.edges_iter(data=True): - d['weight'] = math.e ** (-(d['weight'] ** 2) / sigma) + diff = np.linalg.norm(diff) + if mode == 'similarity': + d['weight'] = math.e ** (-(diff ** 2) / sigma) + elif mode == 'dissimilarity': + d['weight'] = diff + else: + raise ValueError("The mode '%s' is not recodnized" % mode) return graph From b05646e2015f7b0d67eefc073534e2fcf80bb195 Mon Sep 17 00:00:00 2001 From: Vighnesh Birodkar Date: Wed, 23 Jul 2014 03:46:45 +0530 Subject: [PATCH 06/36] Docstring of graph_cut.py --- skimage/graph/graph_cut.py | 63 +++++++++++++++++++++++++++++++++++--- skimage/graph/rag.py | 2 +- 2 files changed, 59 insertions(+), 6 deletions(-) diff --git a/skimage/graph/graph_cut.py b/skimage/graph/graph_cut.py index e49d1182..3ae7b499 100644 --- a/skimage/graph/graph_cut.py +++ b/skimage/graph/graph_cut.py @@ -67,8 +67,44 @@ def cut_threshold(labels, rag, thresh): return map_array[labels] -def cut_n(labels, rag, thresh): +def cut_n(labels, rag, thresh=0.0001): + """Perform Normalized Graph cut on the Region Adjacency Graph. + Given an image's labels and its similarity RAG, recursively perform + a 2-wat Normalized cut on it. All nodes belonging to a subgraph + which cannot be cut further, are assigned a unique label in the + output. + + Parameters + ---------- + labels : ndarray + The array of labels. + rag : RAG + The region adjacency graph. + thresh : float + The threshold. A subgraph won't be further subdivided if the + value of the N-cut exceeds `thresh`. + + Returns + ------- + out : ndarray + The new labelled array. + + Examples + -------- + >>> from skimage import data, graph, segmentation, color, io + >>> img = data.lena() + >>> labels = segmentation.slic(img, compactness=30, n_segments=400) + >>> rag = graph.rag_mean_color(img, labels, mode='similarity') + >>> new_labels = graph.cut_n(labels, rag) + + 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 + + """ _ncut_relabel(rag, thresh) from_ = range(labels.max() + 1) @@ -78,7 +114,24 @@ def cut_n(labels, rag, thresh): return map_array[labels] -def _ncut_relabel(rag, cut_thresh=0.0001): +def _ncut_relabel(rag, thresh): + """Perform Normalized Graph cut on the Region Adjacency Graph. + + Recursively partition the graph into 2, untill further subdividing + yields a cut greather than `thresh` or such a cut cannot be computed. + For such a subgraph, assign a 'ncut label` attribute to all its nodes, + which is a their new unique label. + + Parameters + ---------- + labels : ndarray + The array of labels. + rag : RAG + The region adjacency graph. + thresh : float + The threshold. A subgraph won't be further subdivided if the + value of the N-cut exceeds `thresh`. + """ d, w = _ncut.DW_matrix(rag) error = False @@ -112,7 +165,7 @@ def _ncut_relabel(rag, cut_thresh=0.0001): mcut = cost thresh = t - if (mcut < cut_thresh): + if (mcut < thresh): mask = ev > thresh nodes1 = [n for i, n in enumerate(rag.nodes()) if mask[i]] @@ -121,8 +174,8 @@ def _ncut_relabel(rag, cut_thresh=0.0001): sub1 = rag.subgraph(nodes1) sub2 = rag.subgraph(nodes2) - _ncut_relabel(sub1, cut_thresh) - _ncut_relabel(sub2, cut_thresh) + _ncut_relabel(sub1, thresh) + _ncut_relabel(sub2, thresh) return node = rag.nodes()[0] diff --git a/skimage/graph/rag.py b/skimage/graph/rag.py index 0acede62..65dff3ba 100644 --- a/skimage/graph/rag.py +++ b/skimage/graph/rag.py @@ -120,7 +120,7 @@ def _add_edge_filter(values, graph): def rag_mean_color(image, labels, connectivity=2, mode='dissimilarity', - sigma=30.0): + sigma=255.0): """Compute the Region Adjacency Graph using mean colors. Given an image and its initial segmentation, this method constructs the From 1e39dfc7a1ae8e5942b9916749b44ce7af4fed33 Mon Sep 17 00:00:00 2001 From: Vighnesh Birodkar Date: Wed, 23 Jul 2014 03:54:57 +0530 Subject: [PATCH 07/36] rectified threshold mistake --- skimage/graph/graph_cut.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/skimage/graph/graph_cut.py b/skimage/graph/graph_cut.py index 3ae7b499..9f4f4789 100644 --- a/skimage/graph/graph_cut.py +++ b/skimage/graph/graph_cut.py @@ -157,16 +157,16 @@ def _ncut_relabel(rag, thresh): ev = _ncut.normalize(ev) mcut = np.inf - thresh = None + threshold = None for t in np.arange(0, 1, 0.1): mask = ev > t cost = _ncut.ncut_cost(mask, d, w) if cost < mcut: mcut = cost - thresh = t + threshold = t if (mcut < thresh): - mask = ev > thresh + mask = ev > threshold nodes1 = [n for i, n in enumerate(rag.nodes()) if mask[i]] nodes2 = [n for i, n in enumerate(rag.nodes()) if not mask[i]] From 273b9d081e3195f0acb5f5ae7e427c4f984864e7 Mon Sep 17 00:00:00 2001 From: Vighnesh Birodkar Date: Wed, 23 Jul 2014 04:08:51 +0530 Subject: [PATCH 08/36] docstrings to _ncut.py --- skimage/graph/_ncut.py | 44 ++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 44 insertions(+) diff --git a/skimage/graph/_ncut.py b/skimage/graph/_ncut.py index e31b9431..3e961ad4 100644 --- a/skimage/graph/_ncut.py +++ b/skimage/graph/_ncut.py @@ -4,7 +4,23 @@ from scipy import sparse def DW_matrix(graph): + """Returns the diagonal and weight matrix 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`. + """ + #Cause sparse.eigsh prefers CSC format W = nx.to_scipy_sparse_matrix(graph, format='csc') entries = W.sum(0) D = sparse.dia_matrix((entries, 0), shape=W.shape).tocsc() @@ -12,7 +28,23 @@ def DW_matrix(graph): def ncut_cost(mask, D, W): + """Returns the N-cut cost of a bi-partition of a graph. + Parameters + ---------- + mask : ndarray + The mask for the nodes in the graph. Nodes corrsesponding 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. + """ mask = np.array(mask) mask_list = [np.logical_xor(mask[i], mask) for i in range(mask.shape[0])] mask_array = np.array(mask_list) @@ -25,6 +57,18 @@ def ncut_cost(mask, D, W): 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) From 9c8b1efd2630a5e63c95bca614e00c0cc929d646 Mon Sep 17 00:00:00 2001 From: Vighnesh Birodkar Date: Wed, 23 Jul 2014 04:12:23 +0530 Subject: [PATCH 09/36] docstring to _ncut_cy.pys --- skimage/graph/_ncut_cy.pyx | 12 ++++++++++++ 1 file changed, 12 insertions(+) diff --git a/skimage/graph/_ncut_cy.pyx b/skimage/graph/_ncut_cy.pyx index 77cb18ea..d8156ba4 100644 --- a/skimage/graph/_ncut_cy.pyx +++ b/skimage/graph/_ncut_cy.pyx @@ -7,6 +7,18 @@ import numpy as np def argmin2(cnp.float64_t[:] array): + """Return the index of the 2nd smallest value in an array. + + Parameters + ---------- + a : array + The array to process. + + Returns + ------- + i : 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 i1 = 0 From 2017e350ff630106069adbfdc7f66a5f31f32086 Mon Sep 17 00:00:00 2001 From: Vighnesh Birodkar Date: Wed, 23 Jul 2014 04:35:54 +0530 Subject: [PATCH 10/36] Comments to graph_cut.py --- skimage/graph/graph_cut.py | 26 ++++++++++++++++++++------ 1 file changed, 20 insertions(+), 6 deletions(-) diff --git a/skimage/graph/graph_cut.py b/skimage/graph/graph_cut.py index 9f4f4789..5b7ff974 100644 --- a/skimage/graph/graph_cut.py +++ b/skimage/graph/graph_cut.py @@ -67,7 +67,7 @@ def cut_threshold(labels, rag, thresh): return map_array[labels] -def cut_n(labels, rag, thresh=0.0001): +def cut_n(labels, rag, thresh=0.0001, num_cuts=10): """Perform Normalized Graph cut on the Region Adjacency Graph. Given an image's labels and its similarity RAG, recursively perform @@ -84,6 +84,8 @@ def cut_n(labels, rag, thresh=0.0001): thresh : float The threshold. A subgraph won't be further subdivided if the value of the N-cut exceeds `thresh`. + num_cuts : int + The number or N-cuts to perform before determining the optimal one. Returns ------- @@ -105,7 +107,7 @@ def cut_n(labels, rag, thresh=0.0001): IEEE Transactions on , vol.22, no.8, pp.888,905, Aug 2000 """ - _ncut_relabel(rag, thresh) + _ncut_relabel(rag, thresh, num_cuts) from_ = range(labels.max() + 1) to = [rag.node[x]['ncut label'] for x in from_] @@ -114,7 +116,7 @@ def cut_n(labels, rag, thresh=0.0001): return map_array[labels] -def _ncut_relabel(rag, thresh): +def _ncut_relabel(rag, thresh, num_cuts): """Perform Normalized Graph cut on the Region Adjacency Graph. Recursively partition the graph into 2, untill further subdividing @@ -131,6 +133,8 @@ def _ncut_relabel(rag, thresh): thresh : float The threshold. A subgraph won't be further subdivided if the value of the N-cut exceeds `thresh`. + num_cuts : int + The number or N-cuts to perform before determining the optimal one. """ d, w = _ncut.DW_matrix(rag) error = False @@ -140,13 +144,17 @@ def _ncut_relabel(rag, thresh): vals, vectors = linalg.eigsh(d - w, M=d, which='SM', k=min(100, m - 2)) except ArpackNoConvergence as e: + # Not all eigenvectors converged, salvage the remaining. vals = e.eigenvalues vectors = e.eigenvectors if len(vals) == 0: + # No eigenvector converged. error = True except ValueError: + # k is too less, happens when the graph is of size 1 error = True except ArpackError: + # Arpack failing when two eigenvectors are same error = True if not error: @@ -158,7 +166,8 @@ def _ncut_relabel(rag, thresh): mcut = np.inf threshold = None - for t in np.arange(0, 1, 0.1): + # Perform evenly spaced n-cuts and determine the optimal one. + for t in np.linspace(0, 1, num_cuts, endpoint=False): mask = ev > t cost = _ncut.ncut_cost(mask, d, w) if cost < mcut: @@ -171,13 +180,18 @@ def _ncut_relabel(rag, thresh): nodes1 = [n for i, n in enumerate(rag.nodes()) if mask[i]] nodes2 = [n for i, n in enumerate(rag.nodes()) if not mask[i]] + # Sub divide and perform N-cut again sub1 = rag.subgraph(nodes1) sub2 = rag.subgraph(nodes2) - _ncut_relabel(sub1, thresh) - _ncut_relabel(sub2, thresh) + _ncut_relabel(sub1, thresh, num_cuts) + _ncut_relabel(sub2, thresh, num_cuts) return + # Either an errornous condition occurred, or N-cut wasn't small enough. + # The remaining graph is a region. + # Assign `ncut label` by picking any label from the existing nodes, since + # `labels` are unique, 'ncut label' is also unique. node = rag.nodes()[0] new_label = rag.node[node]['labels'][0] for n in rag.nodes(): From 22b5f4addf800fe959906978f1d5086faa0ddfcb Mon Sep 17 00:00:00 2001 From: Vighnesh Birodkar Date: Wed, 23 Jul 2014 04:47:17 +0530 Subject: [PATCH 11/36] Added example for ncut --- doc/examples/plot_ncut.py | 32 ++++++++++++++++++++++++++++++++ 1 file changed, 32 insertions(+) create mode 100644 doc/examples/plot_ncut.py diff --git a/doc/examples/plot_ncut.py b/doc/examples/plot_ncut.py new file mode 100644 index 00000000..941ef014 --- /dev/null +++ b/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_n(labels1, g) +out2 = color.label2rgb(labels2, img, kind='avg') + +plt.figure() +io.imshow(out1) +plt.figure() +io.imshow(out2) +io.show() From c915aef7591c809dc3c3540de72038e0e11a3f01 Mon Sep 17 00:00:00 2001 From: Vighnesh Birodkar Date: Wed, 23 Jul 2014 05:08:52 +0530 Subject: [PATCH 12/36] Updated bento file --- bento.info | 3 +++ skimage/graph/graph_cut.py | 2 +- 2 files changed, 4 insertions(+), 1 deletion(-) diff --git a/bento.info b/bento.info index 227e51e8..d0891272 100644 --- a/bento.info +++ b/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 diff --git a/skimage/graph/graph_cut.py b/skimage/graph/graph_cut.py index 5b7ff974..acb7d481 100644 --- a/skimage/graph/graph_cut.py +++ b/skimage/graph/graph_cut.py @@ -67,7 +67,7 @@ def cut_threshold(labels, rag, thresh): return map_array[labels] -def cut_n(labels, rag, thresh=0.0001, num_cuts=10): +def cut_n(labels, rag, thresh=0.001, num_cuts=10): """Perform Normalized Graph cut on the Region Adjacency Graph. Given an image's labels and its similarity RAG, recursively perform From d398b7fd9648b9d21cae3adc860d6d68650a35e1 Mon Sep 17 00:00:00 2001 From: Vighnesh Birodkar Date: Wed, 23 Jul 2014 23:35:31 +0530 Subject: [PATCH 13/36] used standard eignen solver --- skimage/graph/graph_cut.py | 9 ++++++++- 1 file changed, 8 insertions(+), 1 deletion(-) diff --git a/skimage/graph/graph_cut.py b/skimage/graph/graph_cut.py index acb7d481..c15eb7c0 100644 --- a/skimage/graph/graph_cut.py +++ b/skimage/graph/graph_cut.py @@ -138,10 +138,17 @@ def _ncut_relabel(rag, thresh, num_cuts): """ d, w = _ncut.DW_matrix(rag) error = False + try: m = w.shape[0] - vals, vectors = linalg.eigsh(d - w, M=d, which='SM', + d2 = d.copy() + # Since d is diagonal, we can directly operate on it's data + # the inverse + d2.data = 1.0/d2.data + # the square root + d2.data = np.sqrt(d2.data) + vals, vectors = linalg.eigsh(d2*(d - w)*d2, which='SM', k=min(100, m - 2)) except ArpackNoConvergence as e: # Not all eigenvectors converged, salvage the remaining. From fc9e8c4670404746e1b7cddf587d51fb82b3f264 Mon Sep 17 00:00:00 2001 From: Vighnesh Birodkar Date: Wed, 23 Jul 2014 23:41:02 +0530 Subject: [PATCH 14/36] doc string correction --- skimage/graph/graph_cut.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/skimage/graph/graph_cut.py b/skimage/graph/graph_cut.py index c15eb7c0..e703a111 100644 --- a/skimage/graph/graph_cut.py +++ b/skimage/graph/graph_cut.py @@ -71,7 +71,7 @@ def cut_n(labels, rag, thresh=0.001, num_cuts=10): """Perform Normalized Graph cut on the Region Adjacency Graph. Given an image's labels and its similarity RAG, recursively perform - a 2-wat Normalized cut on it. All nodes belonging to a subgraph + a 2-way normalized cut on it. All nodes belonging to a subgraph which cannot be cut further, are assigned a unique label in the output. From 62f090977662585a6a9cabdcbcfb07c23140c399 Mon Sep 17 00:00:00 2001 From: Vighnesh Birodkar Date: Wed, 23 Jul 2014 23:46:33 +0530 Subject: [PATCH 15/36] Remove unused except caluses --- skimage/graph/graph_cut.py | 11 +---------- 1 file changed, 1 insertion(+), 10 deletions(-) diff --git a/skimage/graph/graph_cut.py b/skimage/graph/graph_cut.py index e703a111..77e92bbd 100644 --- a/skimage/graph/graph_cut.py +++ b/skimage/graph/graph_cut.py @@ -150,19 +150,10 @@ def _ncut_relabel(rag, thresh, num_cuts): d2.data = np.sqrt(d2.data) vals, vectors = linalg.eigsh(d2*(d - w)*d2, which='SM', k=min(100, m - 2)) - except ArpackNoConvergence as e: - # Not all eigenvectors converged, salvage the remaining. - vals = e.eigenvalues - vectors = e.eigenvectors - if len(vals) == 0: - # No eigenvector converged. - error = True except ValueError: # k is too less, happens when the graph is of size 1 error = True - except ArpackError: - # Arpack failing when two eigenvectors are same - error = True + if not error: vals, vectors = np.real(vals), np.real(vectors) From 6cc0cf4a88f62d28b83aacbbe7e2c76dd31aa76b Mon Sep 17 00:00:00 2001 From: Vighnesh Birodkar Date: Thu, 24 Jul 2014 00:00:09 +0530 Subject: [PATCH 16/36] use map_array instead of modifying graph --- skimage/graph/graph_cut.py | 29 +++++++++++++---------------- skimage/graph/rag.py | 2 +- 2 files changed, 14 insertions(+), 17 deletions(-) diff --git a/skimage/graph/graph_cut.py b/skimage/graph/graph_cut.py index 77e92bbd..a5a54141 100644 --- a/skimage/graph/graph_cut.py +++ b/skimage/graph/graph_cut.py @@ -7,8 +7,6 @@ import numpy as np from . import _ncut from . import _ncut_cy from scipy.sparse import linalg -from scipy.sparse.linalg.eigen.arpack.arpack import ArpackNoConvergence -from scipy.sparse.linalg.eigen.arpack.arpack import ArpackError def cut_threshold(labels, rag, thresh): @@ -107,22 +105,19 @@ def cut_n(labels, rag, thresh=0.001, num_cuts=10): IEEE Transactions on , vol.22, no.8, pp.888,905, Aug 2000 """ - _ncut_relabel(rag, thresh, num_cuts) - - from_ = range(labels.max() + 1) - to = [rag.node[x]['ncut label'] for x in from_] - map_array = np.array(to) + map_array = np.arange(labels.max() + 1) + _ncut_relabel(rag, thresh, num_cuts, map_array) return map_array[labels] -def _ncut_relabel(rag, thresh, num_cuts): +def _ncut_relabel(rag, thresh, num_cuts, map_array): """Perform Normalized Graph cut on the Region Adjacency Graph. Recursively partition the graph into 2, untill further subdividing yields a cut greather than `thresh` or such a cut cannot be computed. - For such a subgraph, assign a 'ncut label` attribute to all its nodes, - which is a their new unique label. + For such a subgraph, indices to labels of all its nodes map to a single + unique value. Parameters ---------- @@ -135,10 +130,12 @@ def _ncut_relabel(rag, thresh, num_cuts): value of the N-cut exceeds `thresh`. num_cuts : int The number or N-cuts to perform before determining the optimal one. + map_array : array + The array which maps old labels to new ones. This is modified inside + the function. """ d, w = _ncut.DW_matrix(rag) error = False - try: m = w.shape[0] @@ -154,7 +151,6 @@ def _ncut_relabel(rag, thresh, num_cuts): # k is too less, happens when the graph is of size 1 error = True - if not error: vals, vectors = np.real(vals), np.real(vectors) index2 = _ncut_cy.argmin2(vals) @@ -182,8 +178,8 @@ def _ncut_relabel(rag, thresh, num_cuts): sub1 = rag.subgraph(nodes1) sub2 = rag.subgraph(nodes2) - _ncut_relabel(sub1, thresh, num_cuts) - _ncut_relabel(sub2, thresh, num_cuts) + _ncut_relabel(sub1, thresh, num_cuts, map_array) + _ncut_relabel(sub2, thresh, num_cuts, map_array) return # Either an errornous condition occurred, or N-cut wasn't small enough. @@ -192,5 +188,6 @@ def _ncut_relabel(rag, thresh, num_cuts): # `labels` are unique, 'ncut label' is also unique. node = rag.nodes()[0] new_label = rag.node[node]['labels'][0] - for n in rag.nodes(): - rag.node[n]['ncut label'] = new_label + for n, d in rag.nodes_iter(data=True): + for l in d['labels']: + map_array[l] = new_label diff --git a/skimage/graph/rag.py b/skimage/graph/rag.py index 65dff3ba..209ee602 100644 --- a/skimage/graph/rag.py +++ b/skimage/graph/rag.py @@ -156,7 +156,7 @@ def rag_mean_color(image, labels, connectivity=2, mode='dissimilarity', It represents how similar two regions are. sigma : float, optional Used for computation when `mode='dissimilarity'`. It governs how close - to each other, two colors should be, for their corresponding edge + to each other two colors should be, for their corresponding edge weight to be significant. A very large value of `sigma` could make any two colors behave as though they were similar. From f261e51faa0ff75d70832eb635481363e2a71aac Mon Sep 17 00:00:00 2001 From: Vighnesh Birodkar Date: Thu, 24 Jul 2014 00:01:56 +0530 Subject: [PATCH 17/36] rename method to cut_normalized --- doc/examples/plot_ncut.py | 2 +- skimage/graph/graph_cut.py | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/doc/examples/plot_ncut.py b/doc/examples/plot_ncut.py index 941ef014..3160f710 100644 --- a/doc/examples/plot_ncut.py +++ b/doc/examples/plot_ncut.py @@ -22,7 +22,7 @@ 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_n(labels1, g) +labels2 = graph.cut_normalized(labels1, g) out2 = color.label2rgb(labels2, img, kind='avg') plt.figure() diff --git a/skimage/graph/graph_cut.py b/skimage/graph/graph_cut.py index a5a54141..d02dde3a 100644 --- a/skimage/graph/graph_cut.py +++ b/skimage/graph/graph_cut.py @@ -65,7 +65,7 @@ def cut_threshold(labels, rag, thresh): return map_array[labels] -def cut_n(labels, rag, thresh=0.001, num_cuts=10): +def cut_normalized(labels, rag, thresh=0.001, num_cuts=10): """Perform Normalized Graph cut on the Region Adjacency Graph. Given an image's labels and its similarity RAG, recursively perform @@ -96,7 +96,7 @@ def cut_n(labels, rag, thresh=0.001, num_cuts=10): >>> img = data.lena() >>> labels = segmentation.slic(img, compactness=30, n_segments=400) >>> rag = graph.rag_mean_color(img, labels, mode='similarity') - >>> new_labels = graph.cut_n(labels, rag) + >>> new_labels = graph.cut_normalized(labels, rag) References ---------- From 7e2f33cb11b4de9e87e7e5ca55f46330b20865a1 Mon Sep 17 00:00:00 2001 From: Vighnesh Birodkar Date: Thu, 24 Jul 2014 00:11:11 +0530 Subject: [PATCH 18/36] Refer sections of paper in comments --- skimage/graph/__init__.py | 4 ++-- skimage/graph/graph_cut.py | 4 ++++ 2 files changed, 6 insertions(+), 2 deletions(-) diff --git a/skimage/graph/__init__.py b/skimage/graph/__init__.py index 9fe00d3c..2d774b44 100644 --- a/skimage/graph/__init__.py +++ b/skimage/graph/__init__.py @@ -1,7 +1,7 @@ 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, cut_n +from .graph_cut import cut_threshold, cut_normalized __all__ = ['shortest_path', 'MCP', @@ -11,5 +11,5 @@ __all__ = ['shortest_path', 'route_through_array', 'rag_mean_color', 'cut_threshold', - 'cut_n', + 'cut_normalized', 'RAG'] diff --git a/skimage/graph/graph_cut.py b/skimage/graph/graph_cut.py index d02dde3a..1010877f 100644 --- a/skimage/graph/graph_cut.py +++ b/skimage/graph/graph_cut.py @@ -145,6 +145,7 @@ def _ncut_relabel(rag, thresh, num_cuts, map_array): d2.data = 1.0/d2.data # the square root d2.data = np.sqrt(d2.data) + # Refer to Equation 7 vals, vectors = linalg.eigsh(d2*(d - w)*d2, which='SM', k=min(100, m - 2)) except ValueError: @@ -152,6 +153,7 @@ def _ncut_relabel(rag, thresh, num_cuts, map_array): error = True if not error: + # Refer Section 3.2.3 vals, vectors = np.real(vals), np.real(vectors) index2 = _ncut_cy.argmin2(vals) @@ -160,6 +162,7 @@ def _ncut_relabel(rag, thresh, num_cuts, map_array): mcut = np.inf threshold = None + # Refer Section 3.1.3 # Perform evenly spaced n-cuts and determine the optimal one. for t in np.linspace(0, 1, num_cuts, endpoint=False): mask = ev > t @@ -178,6 +181,7 @@ def _ncut_relabel(rag, thresh, num_cuts, map_array): sub1 = rag.subgraph(nodes1) sub2 = rag.subgraph(nodes2) + # Refer Section 3.2.5 _ncut_relabel(sub1, thresh, num_cuts, map_array) _ncut_relabel(sub2, thresh, num_cuts, map_array) return From d8c0b2e7dd87672e112af49c6e85d45edb8c52a1 Mon Sep 17 00:00:00 2001 From: Vighnesh Birodkar Date: Thu, 24 Jul 2014 00:58:33 +0530 Subject: [PATCH 19/36] Add unit test --- skimage/graph/tests/test_rag.py | 30 ++++++++++++++++++++++++++++++ 1 file changed, 30 insertions(+) diff --git a/skimage/graph/tests/test_rag.py b/skimage/graph/tests/test_rag.py index 753bcb74..8bdefee8 100644 --- a/skimage/graph/tests/test_rag.py +++ b/skimage/graph/tests/test_rag.py @@ -2,6 +2,8 @@ import numpy as np from skimage import graph from skimage._shared.version_requirements import is_installed from numpy.testing.decorators import skipif +from skimage import segmentation +from numpy import testing def max_edge(g, src, dst, n): @@ -66,3 +68,31 @@ def test_threshold_cut(): # Two labels assert new_labels.max() == 1 + + +def test_cut_normalized(): + + img = np.zeros((100, 100, 3), dtype='uint8') + img[:50, :50] = 255, 255, 255 + img[:50, 50:] = 254, 254, 254 + img[50:, :50] = 2, 2, 2 + img[50:, 50:] = 1, 1, 1 + + labels = np.zeros((100, 100), dtype='uint8') + labels[:50, :50] = 0 + labels[:50, 50:] = 1 + labels[50:, :50] = 2 + labels[50:, 50:] = 3 + + rag = graph.rag_mean_color(img, labels, mode='similarity') + new_labels = graph.cut_normalized(labels, rag) + new_labels, _, _ = segmentation.relabel_sequential(new_labels) + + # Two labels + assert new_labels.max() == 1 + + +def test_rag_error(): + img = np.zeros((10, 10, 3), dtype='uint8') + labels = np.zeros((10, 10), dtype='uint8') + testing.assert_raises(ValueError, graph.rag_mean_color,img, labels, 2, 'non existant mode') From 07cb79cd27838c2a9fbbe0e951209ea067eb9181 Mon Sep 17 00:00:00 2001 From: Vighnesh Birodkar Date: Thu, 31 Jul 2014 23:05:47 +0530 Subject: [PATCH 20/36] Cut cost is computed by Cython code --- skimage/graph/_ncut.py | 18 +++++++++--------- skimage/graph/_ncut_cy.pyx | 27 +++++++++++++++++++++++++++ skimage/graph/graph_cut.py | 2 +- 3 files changed, 37 insertions(+), 10 deletions(-) diff --git a/skimage/graph/_ncut.py b/skimage/graph/_ncut.py index 3e961ad4..4874b252 100644 --- a/skimage/graph/_ncut.py +++ b/skimage/graph/_ncut.py @@ -1,10 +1,11 @@ import networkx as nx import numpy as np from scipy import sparse +from . import _ncut_cy -def DW_matrix(graph): - """Returns the diagonal and weight matrix of a graph. +def DW_matrices(graph): + """Returns the diagonal and weight matrices of a graph. Parameters ---------- @@ -14,15 +15,15 @@ def DW_matrix(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`. + 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 + The weight matrix of the graph. `W[i, j]` is the weight of the edge joining `i` to `j`. """ #Cause sparse.eigsh prefers CSC format W = nx.to_scipy_sparse_matrix(graph, format='csc') - entries = W.sum(0) + entries = W.sum(axis=0) D = sparse.dia_matrix((entries, 0), shape=W.shape).tocsc() return D, W @@ -46,10 +47,9 @@ def ncut_cost(mask, D, W): The cost of performing the N-cut. """ mask = np.array(mask) - mask_list = [np.logical_xor(mask[i], mask) for i in range(mask.shape[0])] - mask_array = np.array(mask_list) + cut = _ncut_cy.cut_cost(mask, W) - cut = float(W[mask_array].sum() / 2.0) + # Cause D has elements only along diagonal assoc_a = D.data[mask].sum() assoc_b = D.data[np.logical_not(mask)].sum() diff --git a/skimage/graph/_ncut_cy.pyx b/skimage/graph/_ncut_cy.pyx index d8156ba4..8c416cbc 100644 --- a/skimage/graph/_ncut_cy.pyx +++ b/skimage/graph/_ncut_cy.pyx @@ -38,3 +38,30 @@ def argmin2(cnp.float64_t[:] array): i += 1 return i2 + + +def cut_cost(mask, W): + mask = np.array(mask) + + 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.float64_t[:] data = W.data + cdef cnp.int32_t row_index + cdef cnp.double_t cost = 0 + + num_rows = W.shape[0] + num_cols = W.shape[1] + + col = 0 + while col < num_cols: + row_index = indptr[col] + while row_index < indptr[col+1]: + row = indices[row_index] + if mask[row] != mask[col]: + cost += data[row_index] + row_index += 1 + col += 1 + + return cost*0.5 diff --git a/skimage/graph/graph_cut.py b/skimage/graph/graph_cut.py index 1010877f..af5048fd 100644 --- a/skimage/graph/graph_cut.py +++ b/skimage/graph/graph_cut.py @@ -134,7 +134,7 @@ def _ncut_relabel(rag, thresh, num_cuts, map_array): The array which maps old labels to new ones. This is modified inside the function. """ - d, w = _ncut.DW_matrix(rag) + d, w = _ncut.DW_matrices(rag) error = False try: From 452921d9f239423aee9e3aae35be9ecc6e40ab1e Mon Sep 17 00:00:00 2001 From: Vighnesh Birodkar Date: Wed, 6 Aug 2014 02:44:33 +0530 Subject: [PATCH 21/36] API changes, comments and docstrings --- skimage/graph/__init__.py | 2 + skimage/graph/_ncut.py | 27 +++++---- skimage/graph/_ncut_cy.pyx | 49 +++++++++++------ skimage/graph/graph_cut.py | 110 +++++++++++++++++++++++++------------ 4 files changed, 126 insertions(+), 62 deletions(-) diff --git a/skimage/graph/__init__.py b/skimage/graph/__init__.py index 2d774b44..aca73bc4 100644 --- a/skimage/graph/__init__.py +++ b/skimage/graph/__init__.py @@ -2,6 +2,7 @@ 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, cut_normalized +ncut = cut_normalized __all__ = ['shortest_path', 'MCP', @@ -12,4 +13,5 @@ __all__ = ['shortest_path', 'rag_mean_color', 'cut_threshold', 'cut_normalized', + 'ncut', 'RAG'] diff --git a/skimage/graph/_ncut.py b/skimage/graph/_ncut.py index 4874b252..fbefca27 100644 --- a/skimage/graph/_ncut.py +++ b/skimage/graph/_ncut.py @@ -21,20 +21,20 @@ def DW_matrices(graph): The weight matrix of the graph. `W[i, j]` is the weight of the edge joining `i` to `j`. """ - #Cause sparse.eigsh prefers CSC format + # 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(mask, D, W): +def ncut_cost(cut, D, W): """Returns the N-cut cost of a bi-partition of a graph. Parameters ---------- - mask : ndarray - The mask for the nodes in the graph. Nodes corrsesponding to a `True` + 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. @@ -45,15 +45,22 @@ def ncut_cost(mask, D, W): ------- 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. """ - mask = np.array(mask) - cut = _ncut_cy.cut_cost(mask, W) + cut = np.array(cut) + cut_cost = _ncut_cy.cut_cost(cut, W) - # Cause D has elements only along diagonal - assoc_a = D.data[mask].sum() - assoc_b = D.data[np.logical_not(mask)].sum() + # 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[np.logical_not(cut)].sum() - return (cut / assoc_a) + (cut / assoc_b) + return (cut_cost / assoc_a) + (cut_cost / assoc_b) def normalize(a): diff --git a/skimage/graph/_ncut_cy.pyx b/skimage/graph/_ncut_cy.pyx index 8c416cbc..0f859edf 100644 --- a/skimage/graph/_ncut_cy.pyx +++ b/skimage/graph/_ncut_cy.pyx @@ -16,52 +16,65 @@ def argmin2(cnp.float64_t[:] array): Returns ------- - i : int + 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 i1 = 0 - cdef Py_ssize_t i2 = 0 + 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] - while i < array.shape[0]: + for i in range(n): x = array[i] if x < min1: min2 = min1 - i2 = i1 + min_idx2 = min_idx1 min1 = x - i1 = i + min_idx1 = i elif x > min1 and x < min2: min2 = x - i2 = i + min_idx2 = i i += 1 - return i2 + return min_idx2 -def cut_cost(mask, W): - mask = np.array(mask) +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.float64_t[:] data = W.data + 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] - col = 0 - while col < num_cols: - row_index = indptr[col] - while row_index < indptr[col+1]: + for col in range(num_cols): + for row_index in range(indptr[col], indptr[col + 1]): row = indices[row_index] - if mask[row] != mask[col]: - cost += data[row_index] + if cut_mask[row] != cut_mask[col]: + cost += data[row_index] row_index += 1 col += 1 - return cost*0.5 + return cost * 0.5 diff --git a/skimage/graph/graph_cut.py b/skimage/graph/graph_cut.py index af5048fd..9c133bd3 100644 --- a/skimage/graph/graph_cut.py +++ b/skimage/graph/graph_cut.py @@ -111,6 +111,65 @@ def cut_normalized(labels, rag, thresh=0.001, num_cuts=10): return map_array[labels] +def partition_by_cut(cut, rag): + """Compute resulting subgraphs from given bi-parition. + + Parameters + ---------- + cut : array + A array of booleans. Elements set to `True` belong to one + set. + rag : RAG + The Region Adjacency Graph. + + Returns + ------- + sub1, sub2 : RAG + The two resulting subgraphs from the bi-partition. + """ + nodes1 = [n for i, n in enumerate(rag.nodes()) if cut[i]] + nodes2 = [n for i, n in enumerate(rag.nodes()) if not cut[i]] + + sub1 = rag.subgraph(nodes1) + sub2 = rag.subgraph(nodes2) + + return sub1, sub2 + + +def get_min_ncut(ev, d, w, num_cuts): + """Threshold an eigen vector evenly, to determine minimum ncut. + + Parameters + ---------- + ev : array + The eigenvector to threshold. + d : ndarray + The diagonal matrix of the graph. + w : ndarray + The weight matrix of the graph. + num_cuts : int + The number of evenly spaced thresholds to check for. + + Returns + ------- + threshold, mcut : float + The threshold which produced the minimum ncut, and the value of the + ncut itself. + """ + mcut = np.inf + + # Refer Shi & Malik 2001, Section 3.1.3, Page 892 + # Perform evenly spaced n-cuts and determine the optimal one. + for t in np.linspace(0, 1, num_cuts, endpoint=False): + mask = ev > t + cost = _ncut.ncut_cost(mask, d, w) + if cost < mcut: + mcut = cost + threshold = t + + return threshold, mcut + + def _ncut_relabel(rag, thresh, num_cuts, map_array): """Perform Normalized Graph cut on the Region Adjacency Graph. @@ -135,58 +194,41 @@ def _ncut_relabel(rag, thresh, num_cuts, map_array): the function. """ d, w = _ncut.DW_matrices(rag) - error = False + stop = False + m = w.shape[0] - try: - m = w.shape[0] + if m > 2: d2 = d.copy() # Since d is diagonal, we can directly operate on it's data # the inverse - d2.data = 1.0/d2.data + d2.data = 1.0 / d2.data # the square root d2.data = np.sqrt(d2.data) - # Refer to Equation 7 - vals, vectors = linalg.eigsh(d2*(d - w)*d2, which='SM', + # Refer Shi & Malik 2001, Equation 7, Page 891 + vals, vectors = linalg.eigsh(d2 * (d - w) * d2, which='SM', k=min(100, m - 2)) - except ValueError: - # k is too less, happens when the graph is of size 1 - error = True + else: + stop = True - if not error: - # Refer Section 3.2.3 + if not stop: + # Pick second smalles eigen vector. + # Refer Shi & Malik 2001, Section 3.2.3, Page 893 vals, vectors = np.real(vals), np.real(vectors) index2 = _ncut_cy.argmin2(vals) + ev = _ncut.normalize(vectors[:, index2]) - ev = np.real(vectors[:, index2]) - ev = _ncut.normalize(ev) - - mcut = np.inf - threshold = None - # Refer Section 3.1.3 - # Perform evenly spaced n-cuts and determine the optimal one. - for t in np.linspace(0, 1, num_cuts, endpoint=False): - mask = ev > t - cost = _ncut.ncut_cost(mask, d, w) - if cost < mcut: - mcut = cost - threshold = t - + threshold, mcut = get_min_ncut(ev, d, w, num_cuts) if (mcut < thresh): - mask = ev > threshold - - nodes1 = [n for i, n in enumerate(rag.nodes()) if mask[i]] - nodes2 = [n for i, n in enumerate(rag.nodes()) if not mask[i]] - + cut_mask = ev > threshold # Sub divide and perform N-cut again - sub1 = rag.subgraph(nodes1) - sub2 = rag.subgraph(nodes2) + # Refer Shi & Malik 2001, Section 3.2.5, Page 893 + sub1, sub2 = partition_by_cut(cut_mask, rag) - # Refer Section 3.2.5 _ncut_relabel(sub1, thresh, num_cuts, map_array) _ncut_relabel(sub2, thresh, num_cuts, map_array) return - # Either an errornous condition occurred, or N-cut wasn't small enough. + # The N-cut wasn't small enough, or could not be computed. # The remaining graph is a region. # Assign `ncut label` by picking any label from the existing nodes, since # `labels` are unique, 'ncut label' is also unique. From 69fec94326e56ec009f62c7626422be71fbe3f08 Mon Sep 17 00:00:00 2001 From: Vighnesh Birodkar Date: Wed, 6 Aug 2014 02:46:29 +0530 Subject: [PATCH 22/36] rag copy in threshold cut --- skimage/graph/graph_cut.py | 1 + 1 file changed, 1 insertion(+) diff --git a/skimage/graph/graph_cut.py b/skimage/graph/graph_cut.py index 9c133bd3..bc694409 100644 --- a/skimage/graph/graph_cut.py +++ b/skimage/graph/graph_cut.py @@ -47,6 +47,7 @@ def cut_threshold(labels, rag, thresh): """ # Because deleting edges while iterating through them produces an error. + rag = rag.copy() to_remove = [(x, y) for x, y, d in rag.edges_iter(data=True) if d['weight'] >= thresh] rag.remove_edges_from(to_remove) From 3ad2b682d6c96dff7ab232cee092710bd0adb2aa Mon Sep 17 00:00:00 2001 From: Vighnesh Birodkar Date: Wed, 6 Aug 2014 02:47:24 +0530 Subject: [PATCH 23/36] test skip --- skimage/graph/tests/test_rag.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/skimage/graph/tests/test_rag.py b/skimage/graph/tests/test_rag.py index 8bdefee8..4cd5fd36 100644 --- a/skimage/graph/tests/test_rag.py +++ b/skimage/graph/tests/test_rag.py @@ -69,7 +69,7 @@ def test_threshold_cut(): # Two labels assert new_labels.max() == 1 - +@skipif(not is_installed('networkx')) def test_cut_normalized(): img = np.zeros((100, 100, 3), dtype='uint8') From 35c9746b7407f495875dc8df116442eb287c58dc Mon Sep 17 00:00:00 2001 From: Vighnesh Birodkar Date: Wed, 6 Aug 2014 02:48:03 +0530 Subject: [PATCH 24/36] comment --- skimage/graph/graph_cut.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/skimage/graph/graph_cut.py b/skimage/graph/graph_cut.py index bc694409..65d14585 100644 --- a/skimage/graph/graph_cut.py +++ b/skimage/graph/graph_cut.py @@ -232,7 +232,7 @@ def _ncut_relabel(rag, thresh, num_cuts, map_array): # The N-cut wasn't small enough, or could not be computed. # The remaining graph is a region. # Assign `ncut label` by picking any label from the existing nodes, since - # `labels` are unique, 'ncut label' is also unique. + # `labels` are unique, `new_label` is also unique. node = rag.nodes()[0] new_label = rag.node[node]['labels'][0] for n, d in rag.nodes_iter(data=True): From 9371c4f0e738763238ede8e99788f19770b8cb92 Mon Sep 17 00:00:00 2001 From: Vighnesh Birodkar Date: Wed, 6 Aug 2014 22:38:18 +0530 Subject: [PATCH 25/36] Fixed potential import error --- skimage/graph/_ncut.py | 6 +++++- skimage/graph/graph_cut.py | 2 +- 2 files changed, 6 insertions(+), 2 deletions(-) diff --git a/skimage/graph/_ncut.py b/skimage/graph/_ncut.py index fbefca27..46ffb528 100644 --- a/skimage/graph/_ncut.py +++ b/skimage/graph/_ncut.py @@ -1,4 +1,8 @@ -import networkx as nx +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 diff --git a/skimage/graph/graph_cut.py b/skimage/graph/graph_cut.py index 65d14585..e44f8a5c 100644 --- a/skimage/graph/graph_cut.py +++ b/skimage/graph/graph_cut.py @@ -2,7 +2,7 @@ try: import networkx as nx except ImportError: import warnings - warnings.warn('"cut_threshold" requires networkx') + warnings.warn('RAGs require networkx') import numpy as np from . import _ncut from . import _ncut_cy From 1fa7bced8d10d08a82ed039b9959179e1ca19a88 Mon Sep 17 00:00:00 2001 From: Vighnesh Birodkar Date: Wed, 6 Aug 2014 23:45:35 +0530 Subject: [PATCH 26/36] skip error test --- skimage/graph/tests/test_rag.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/skimage/graph/tests/test_rag.py b/skimage/graph/tests/test_rag.py index 4cd5fd36..81889d6f 100644 --- a/skimage/graph/tests/test_rag.py +++ b/skimage/graph/tests/test_rag.py @@ -91,7 +91,7 @@ def test_cut_normalized(): # Two labels assert new_labels.max() == 1 - +@skipif(not is_installed('networkx')) def test_rag_error(): img = np.zeros((10, 10, 3), dtype='uint8') labels = np.zeros((10, 10), dtype='uint8') From 6d0bf723d700e943e977db27012821906a33bf86 Mon Sep 17 00:00:00 2001 From: Vighnesh Birodkar Date: Fri, 8 Aug 2014 23:16:38 +0530 Subject: [PATCH 27/36] Docstring changes and typos --- skimage/graph/_ncut.py | 6 +++--- skimage/graph/_ncut_cy.pyx | 2 +- skimage/graph/graph_cut.py | 12 ++++++------ skimage/graph/rag.py | 8 ++++---- skimage/graph/tests/test_rag.py | 2 +- 5 files changed, 15 insertions(+), 15 deletions(-) diff --git a/skimage/graph/_ncut.py b/skimage/graph/_ncut.py index 46ffb528..2735bbd2 100644 --- a/skimage/graph/_ncut.py +++ b/skimage/graph/_ncut.py @@ -19,10 +19,10 @@ def DW_matrices(graph): Returns ------- D : csc_matrix - The diagonal matrix of the graph. `D[i, i]` is the sum of weights of + 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 + 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 @@ -62,7 +62,7 @@ def ncut_cost(cut, D, 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[np.logical_not(cut)].sum() + assoc_b = D.data[~cut].sum() return (cut_cost / assoc_a) + (cut_cost / assoc_b) diff --git a/skimage/graph/_ncut_cy.pyx b/skimage/graph/_ncut_cy.pyx index 0f859edf..0888be72 100644 --- a/skimage/graph/_ncut_cy.pyx +++ b/skimage/graph/_ncut_cy.pyx @@ -6,7 +6,7 @@ cimport numpy as cnp import numpy as np -def argmin2(cnp.float64_t[:] array): +def argmin2(cnp.double_t[:] array): """Return the index of the 2nd smallest value in an array. Parameters diff --git a/skimage/graph/graph_cut.py b/skimage/graph/graph_cut.py index e44f8a5c..b86f6c48 100644 --- a/skimage/graph/graph_cut.py +++ b/skimage/graph/graph_cut.py @@ -89,11 +89,11 @@ def cut_normalized(labels, rag, thresh=0.001, num_cuts=10): Returns ------- out : ndarray - The new labelled array. + The new labeled array. Examples -------- - >>> from skimage import data, graph, segmentation, color, io + >>> from skimage import data, graph, segmentation >>> img = data.lena() >>> labels = segmentation.slic(img, compactness=30, n_segments=400) >>> rag = graph.rag_mean_color(img, labels, mode='similarity') @@ -103,7 +103,7 @@ def cut_normalized(labels, rag, thresh=0.001, num_cuts=10): ---------- .. [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 + IEEE Transactions on , vol.22, no.8, pp.888,905, August 2000 """ map_array = np.arange(labels.max() + 1) @@ -138,7 +138,7 @@ def partition_by_cut(cut, rag): def get_min_ncut(ev, d, w, num_cuts): - """Threshold an eigen vector evenly, to determine minimum ncut. + """Threshold an eigenvector evenly, to determine minimum ncut. Parameters ---------- @@ -174,7 +174,7 @@ def get_min_ncut(ev, d, w, num_cuts): def _ncut_relabel(rag, thresh, num_cuts, map_array): """Perform Normalized Graph cut on the Region Adjacency Graph. - Recursively partition the graph into 2, untill further subdividing + Recursively partition the graph into 2, until further subdivision yields a cut greather than `thresh` or such a cut cannot be computed. For such a subgraph, indices to labels of all its nodes map to a single unique value. @@ -212,7 +212,7 @@ def _ncut_relabel(rag, thresh, num_cuts, map_array): stop = True if not stop: - # Pick second smalles eigen vector. + # Pick second smalles eigenvector. # Refer Shi & Malik 2001, Section 3.2.3, Page 893 vals, vectors = np.real(vals), np.real(vectors) index2 = _ncut_cy.argmin2(vals) diff --git a/skimage/graph/rag.py b/skimage/graph/rag.py index 209ee602..f33936c2 100644 --- a/skimage/graph/rag.py +++ b/skimage/graph/rag.py @@ -146,17 +146,17 @@ def rag_mean_color(image, labels, connectivity=2, mode='dissimilarity', The strategy to assign edge weights. 'similarity' : The weight between two adjacent regions is the - :math:`|c_1 - c_2|`, where :math:`c1` and :math:`c2` are the mean + :math:`|c_1 - c_2|`, where :math:`c_1` and :math:`c_2` are the mean colors of the two regions. It represents how different two regions are. 'dissimilarity' : The weight between two adjacent is :math:`e^{-d^2/sigma}` where :math:`d=|c_1 - c_2|`, where - :math:`c1` and :math:`c2` are the mean colors of the two regions. + :math:`c_1` and :math:`c_2` are the mean colors of the two regions. It represents how similar two regions are. sigma : float, optional - Used for computation when `mode='dissimilarity'`. It governs how close - to each other two colors should be, for their corresponding edge + Used for computation when `mode` is "dissimilarity". It governs how + close to each other two colors should be, for their corresponding edge weight to be significant. A very large value of `sigma` could make any two colors behave as though they were similar. diff --git a/skimage/graph/tests/test_rag.py b/skimage/graph/tests/test_rag.py index 81889d6f..1e643fc4 100644 --- a/skimage/graph/tests/test_rag.py +++ b/skimage/graph/tests/test_rag.py @@ -95,4 +95,4 @@ def test_cut_normalized(): def test_rag_error(): img = np.zeros((10, 10, 3), dtype='uint8') labels = np.zeros((10, 10), dtype='uint8') - testing.assert_raises(ValueError, graph.rag_mean_color,img, labels, 2, 'non existant mode') + testing.assert_raises(ValueError, graph.rag_mean_color, img, labels, 2, 'non existant mode') From 3973b3030039c05bedefba5e70464dfb66af2250 Mon Sep 17 00:00:00 2001 From: Vighnesh Birodkar Date: Sat, 9 Aug 2014 21:29:41 +0530 Subject: [PATCH 28/36] docstring changes and code movement --- skimage/graph/_ncut_cy.pyx | 2 - skimage/graph/graph_cut.py | 84 +++++++++++++++++++++++++------------- skimage/graph/rag.py | 17 ++++---- 3 files changed, 64 insertions(+), 39 deletions(-) diff --git a/skimage/graph/_ncut_cy.pyx b/skimage/graph/_ncut_cy.pyx index 0888be72..de3b5c82 100644 --- a/skimage/graph/_ncut_cy.pyx +++ b/skimage/graph/_ncut_cy.pyx @@ -74,7 +74,5 @@ def cut_cost(cut, W): row = indices[row_index] if cut_mask[row] != cut_mask[col]: cost += data[row_index] - row_index += 1 - col += 1 return cost * 0.5 diff --git a/skimage/graph/graph_cut.py b/skimage/graph/graph_cut.py index b86f6c48..1e3d31dc 100644 --- a/skimage/graph/graph_cut.py +++ b/skimage/graph/graph_cut.py @@ -9,7 +9,7 @@ from . import _ncut_cy from scipy.sparse import linalg -def cut_threshold(labels, rag, thresh): +def cut_threshold(labels, rag, thresh, in_place=True): """Combine regions seperated by weight less than threshold. Given an image's labels and its RAG, output new labels by @@ -25,6 +25,10 @@ def cut_threshold(labels, rag, thresh): thresh : float The threshold. Regions connected by edges with smaller weights are combined. + in_place : bool + If set, modifies `rag` in place. The function will remove the edges + with weights less that `thresh`. If set to `False` the function + makes a copy of `rag` before proceeding. Returns ------- @@ -47,7 +51,9 @@ def cut_threshold(labels, rag, thresh): """ # Because deleting edges while iterating through them produces an error. - rag = rag.copy() + if not in_place: + rag = rag.copy() + to_remove = [(x, y) for x, y, d in rag.edges_iter(data=True) if d['weight'] >= thresh] rag.remove_edges_from(to_remove) @@ -66,7 +72,7 @@ def cut_threshold(labels, rag, thresh): return map_array[labels] -def cut_normalized(labels, rag, thresh=0.001, num_cuts=10): +def cut_normalized(labels, rag, thresh=0.001, num_cuts=10, in_place=True): """Perform Normalized Graph cut on the Region Adjacency Graph. Given an image's labels and its similarity RAG, recursively perform @@ -85,6 +91,9 @@ def cut_normalized(labels, rag, thresh=0.001, num_cuts=10): value of the N-cut exceeds `thresh`. num_cuts : int The number or N-cuts to perform before determining the optimal one. + in_place : bool + If set, modifies `rag` in place. For each node `n` the function will + set a new attribute ``rag.node[n]['ncut label]``. Returns ------- @@ -106,8 +115,15 @@ def cut_normalized(labels, rag, thresh=0.001, num_cuts=10): IEEE Transactions on , vol.22, no.8, pp.888,905, August 2000 """ - map_array = np.arange(labels.max() + 1) - _ncut_relabel(rag, thresh, num_cuts, map_array) + if not in_place: + rag = rag.copy() + + _ncut_relabel(rag, thresh, num_cuts) + + map_array = np.zeros(labels.max() + 1) + # Mapping from old labels to new + for n, d in rag.nodes_iter(data=True): + map_array[d['labels']] = d['ncut label'] return map_array[labels] @@ -128,6 +144,16 @@ def partition_by_cut(cut, rag): sub1, sub2 : RAG The two resulting subgraphs from the bi-partition. """ + # `cut` is derived from `D` and `W` matrices, which also follow the + # ordering returned by `rag.nodes()` because we use + # nx.to_scipy_sparce_matrix. + + # Example + # rag.nodes() = [3, 7, 9, 13] + # cut = [True, False, True, False] + # nodes1 = [3, 9] + # nodes2 = [7, 10] + nodes1 = [n for i, n in enumerate(rag.nodes()) if cut[i]] nodes2 = [n for i, n in enumerate(rag.nodes()) if not cut[i]] @@ -153,9 +179,10 @@ def get_min_ncut(ev, d, w, num_cuts): Returns ------- - threshold, mcut : float - The threshold which produced the minimum ncut, and the value of the - ncut itself. + mask : array + The array of booleans which denotes the bi-partition. + mcut : float + The value of the minimum ncut. """ mcut = np.inf @@ -165,13 +192,21 @@ def get_min_ncut(ev, d, w, num_cuts): mask = ev > t cost = _ncut.ncut_cost(mask, d, w) if cost < mcut: + min_mask = mask mcut = cost - threshold = t - return threshold, mcut + return min_mask, mcut -def _ncut_relabel(rag, thresh, num_cuts, map_array): +def _label_all(rag, attr_name): + node = rag.nodes()[0] + new_label = rag.node[node]['labels'][0] + for n, d in rag.nodes_iter(data=True): + for l in d['labels']: + d[attr_name] = new_label + + +def _ncut_relabel(rag, thresh, num_cuts): """Perform Normalized Graph cut on the Region Adjacency Graph. Recursively partition the graph into 2, until further subdivision @@ -195,46 +230,37 @@ def _ncut_relabel(rag, thresh, num_cuts, map_array): the function. """ d, w = _ncut.DW_matrices(rag) - stop = False m = w.shape[0] if m > 2: d2 = d.copy() # Since d is diagonal, we can directly operate on it's data - # the inverse - d2.data = 1.0 / d2.data - # the square root - d2.data = np.sqrt(d2.data) + # the inverse of the square root + d2.data = np.reciprocal(np.sqrt(d2.data, out=d2.data), out=d2.data) + # Refer Shi & Malik 2001, Equation 7, Page 891 vals, vectors = linalg.eigsh(d2 * (d - w) * d2, which='SM', k=min(100, m - 2)) - else: - stop = True - if not stop: - # Pick second smalles eigenvector. + # Pick second smallest eigenvector. # Refer Shi & Malik 2001, Section 3.2.3, Page 893 vals, vectors = np.real(vals), np.real(vectors) index2 = _ncut_cy.argmin2(vals) ev = _ncut.normalize(vectors[:, index2]) - threshold, mcut = get_min_ncut(ev, d, w, num_cuts) + cut_mask, mcut = get_min_ncut(ev, d, w, num_cuts) if (mcut < thresh): - cut_mask = ev > threshold # Sub divide and perform N-cut again # Refer Shi & Malik 2001, Section 3.2.5, Page 893 sub1, sub2 = partition_by_cut(cut_mask, rag) - _ncut_relabel(sub1, thresh, num_cuts, map_array) - _ncut_relabel(sub2, thresh, num_cuts, map_array) + _ncut_relabel(sub1, thresh, num_cuts) + _ncut_relabel(sub2, thresh, num_cuts) return # The N-cut wasn't small enough, or could not be computed. # The remaining graph is a region. # Assign `ncut label` by picking any label from the existing nodes, since # `labels` are unique, `new_label` is also unique. - node = rag.nodes()[0] - new_label = rag.node[node]['labels'][0] - for n, d in rag.nodes_iter(data=True): - for l in d['labels']: - map_array[l] = new_label + + _label_all(rag, 'ncut label') diff --git a/skimage/graph/rag.py b/skimage/graph/rag.py index f33936c2..1366a6d9 100644 --- a/skimage/graph/rag.py +++ b/skimage/graph/rag.py @@ -2,6 +2,7 @@ try: import networkx as nx except ImportError: msg = "Graph functions require networkx, which is not installed" + class nx: class Graph: def __init__(self, *args, **kwargs): @@ -119,7 +120,7 @@ def _add_edge_filter(values, graph): return 0 -def rag_mean_color(image, labels, connectivity=2, mode='dissimilarity', +def rag_mean_color(image, labels, connectivity=2, mode='distance', sigma=255.0): """Compute the Region Adjacency Graph using mean colors. @@ -142,15 +143,15 @@ def rag_mean_color(image, labels, connectivity=2, mode='dissimilarity', are considered adjacent. It can range from 1 to `labels.ndim`. Its behavior is the same as `connectivity` parameter in `scipy.ndimage.filters.generate_binary_structure`. - mode : str['similarity' | 'dissimilarity'] + mode : {'distance', 'similarity'}, optional The strategy to assign edge weights. - 'similarity' : The weight between two adjacent regions is the + 'distance' : The weight between two adjacent regions is the :math:`|c_1 - c_2|`, where :math:`c_1` and :math:`c_2` are the mean - colors of the two regions. It represents how different two regions - are. + colors of the two regions. It represents the Euclidian distance in + their average color. - 'dissimilarity' : The weight between two adjacent is + 'similarity' : The weight between two adjacent is :math:`e^{-d^2/sigma}` where :math:`d=|c_1 - c_2|`, where :math:`c_1` and :math:`c_2` are the mean colors of the two regions. It represents how similar two regions are. @@ -229,9 +230,9 @@ def rag_mean_color(image, labels, connectivity=2, mode='dissimilarity', diff = np.linalg.norm(diff) if mode == 'similarity': d['weight'] = math.e ** (-(diff ** 2) / sigma) - elif mode == 'dissimilarity': + elif mode == 'distance': d['weight'] = diff else: - raise ValueError("The mode '%s' is not recodnized" % mode) + raise ValueError("The mode '%s' is not recognised" % mode) return graph From ded138cbcb3040beea4947f800ac3bb306200a11 Mon Sep 17 00:00:00 2001 From: Vighnesh Birodkar Date: Sat, 9 Aug 2014 21:35:07 +0530 Subject: [PATCH 29/36] test case for in place --- skimage/graph/tests/test_rag.py | 20 +++++++++++++++----- 1 file changed, 15 insertions(+), 5 deletions(-) diff --git a/skimage/graph/tests/test_rag.py b/skimage/graph/tests/test_rag.py index 1e643fc4..0e62aa1a 100644 --- a/skimage/graph/tests/test_rag.py +++ b/skimage/graph/tests/test_rag.py @@ -64,11 +64,15 @@ def test_threshold_cut(): labels[50:, 50:] = 3 rag = graph.rag_mean_color(img, labels) - new_labels = graph.cut_threshold(labels, rag, 10) - + new_labels = graph.cut_threshold(labels, rag, 10, in_place=False) # Two labels assert new_labels.max() == 1 + new_labels = graph.cut_threshold(labels, rag, 10) + # Two labels + assert new_labels.max() == 1 + + @skipif(not is_installed('networkx')) def test_cut_normalized(): @@ -85,14 +89,20 @@ def test_cut_normalized(): labels[50:, 50:] = 3 rag = graph.rag_mean_color(img, labels, mode='similarity') - new_labels = graph.cut_normalized(labels, rag) - new_labels, _, _ = segmentation.relabel_sequential(new_labels) + new_labels = graph.cut_normalized(labels, rag, in_place=False) + new_labels, _, _ = segmentation.relabel_sequential(new_labels) # Two labels assert new_labels.max() == 1 + new_labels = graph.cut_normalized(labels, rag) + new_labels, _, _ = segmentation.relabel_sequential(new_labels) + assert new_labels.max() == 1 + + @skipif(not is_installed('networkx')) def test_rag_error(): img = np.zeros((10, 10, 3), dtype='uint8') labels = np.zeros((10, 10), dtype='uint8') - testing.assert_raises(ValueError, graph.rag_mean_color, img, labels, 2, 'non existant mode') + testing.assert_raises(ValueError, graph.rag_mean_color, img, labels, + 2, 'non existant mode') From 03d3873cc43cb2f9725da27fd85d8704625a64af Mon Sep 17 00:00:00 2001 From: Vighnesh Birodkar Date: Sat, 9 Aug 2014 22:08:27 +0530 Subject: [PATCH 30/36] docstring of _label_all --- skimage/graph/graph_cut.py | 12 +++++++++++- 1 file changed, 11 insertions(+), 1 deletion(-) diff --git a/skimage/graph/graph_cut.py b/skimage/graph/graph_cut.py index 1e3d31dc..f69df469 100644 --- a/skimage/graph/graph_cut.py +++ b/skimage/graph/graph_cut.py @@ -199,6 +199,17 @@ def get_min_ncut(ev, d, w, num_cuts): def _label_all(rag, attr_name): + """Assign a uique integer to the given attribute in the RAG. + + This function assumes that all labels in `rag` are unique. It + picks up a random label from them and assigns it to the `attr_name` + attribute of all the nodes. + + rag : RAG + The Region Adjacency Graph. + attr_name : string + The attribute to which a unique integer is assigned. + """ node = rag.nodes()[0] new_label = rag.node[node]['labels'][0] for n, d in rag.nodes_iter(data=True): @@ -262,5 +273,4 @@ def _ncut_relabel(rag, thresh, num_cuts): # The remaining graph is a region. # Assign `ncut label` by picking any label from the existing nodes, since # `labels` are unique, `new_label` is also unique. - _label_all(rag, 'ncut label') From 0c185629d3002ffcb9d5306c7f848243585f8341 Mon Sep 17 00:00:00 2001 From: Vighnesh Birodkar Date: Sat, 9 Aug 2014 22:09:12 +0530 Subject: [PATCH 31/36] corrected similarity --- skimage/graph/rag.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/skimage/graph/rag.py b/skimage/graph/rag.py index 1366a6d9..c1108aa2 100644 --- a/skimage/graph/rag.py +++ b/skimage/graph/rag.py @@ -156,7 +156,7 @@ def rag_mean_color(image, labels, connectivity=2, mode='distance', :math:`c_1` and :math:`c_2` are the mean colors of the two regions. It represents how similar two regions are. sigma : float, optional - Used for computation when `mode` is "dissimilarity". It governs how + Used for computation when `mode` is "similarity". It governs how close to each other two colors should be, for their corresponding edge weight to be significant. A very large value of `sigma` could make any two colors behave as though they were similar. From 4a231cfd3587e6e043aafd7aab59df28e185a217 Mon Sep 17 00:00:00 2001 From: Vighnesh Birodkar Date: Sun, 10 Aug 2014 20:07:34 +0530 Subject: [PATCH 32/36] Minor changes --- skimage/graph/graph_cut.py | 13 ++++++------- skimage/graph/rag.py | 2 +- 2 files changed, 7 insertions(+), 8 deletions(-) diff --git a/skimage/graph/graph_cut.py b/skimage/graph/graph_cut.py index f69df469..56939b7e 100644 --- a/skimage/graph/graph_cut.py +++ b/skimage/graph/graph_cut.py @@ -50,10 +50,10 @@ def cut_threshold(labels, rag, thresh, in_place=True): http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.11.5274 """ - # Because deleting edges while iterating through them produces an error. if not in_place: rag = rag.copy() + # Because deleting edges while iterating through them produces an error. to_remove = [(x, y) for x, y, d in rag.edges_iter(data=True) if d['weight'] >= thresh] rag.remove_edges_from(to_remove) @@ -93,7 +93,7 @@ def cut_normalized(labels, rag, thresh=0.001, num_cuts=10, in_place=True): The number or N-cuts to perform before determining the optimal one. in_place : bool If set, modifies `rag` in place. For each node `n` the function will - set a new attribute ``rag.node[n]['ncut label]``. + set a new attribute ``rag.node[n]['ncut label']``. Returns ------- @@ -146,7 +146,7 @@ def partition_by_cut(cut, rag): """ # `cut` is derived from `D` and `W` matrices, which also follow the # ordering returned by `rag.nodes()` because we use - # nx.to_scipy_sparce_matrix. + # nx.to_scipy_sparse_matrix. # Example # rag.nodes() = [3, 7, 9, 13] @@ -199,7 +199,7 @@ def get_min_ncut(ev, d, w, num_cuts): def _label_all(rag, attr_name): - """Assign a uique integer to the given attribute in the RAG. + """Assign a unique integer to the given attribute in the RAG. This function assumes that all labels in `rag` are unique. It picks up a random label from them and assigns it to the `attr_name` @@ -213,8 +213,7 @@ def _label_all(rag, attr_name): node = rag.nodes()[0] new_label = rag.node[node]['labels'][0] for n, d in rag.nodes_iter(data=True): - for l in d['labels']: - d[attr_name] = new_label + d[attr_name] = new_label def _ncut_relabel(rag, thresh, num_cuts): @@ -245,7 +244,7 @@ def _ncut_relabel(rag, thresh, num_cuts): if m > 2: d2 = d.copy() - # Since d is diagonal, we can directly operate on it's data + # Since d is diagonal, we can directly operate on its data # the inverse of the square root d2.data = np.reciprocal(np.sqrt(d2.data, out=d2.data), out=d2.data) diff --git a/skimage/graph/rag.py b/skimage/graph/rag.py index c1108aa2..5ef9a83b 100644 --- a/skimage/graph/rag.py +++ b/skimage/graph/rag.py @@ -148,7 +148,7 @@ def rag_mean_color(image, labels, connectivity=2, mode='distance', 'distance' : The weight between two adjacent regions is the :math:`|c_1 - c_2|`, where :math:`c_1` and :math:`c_2` are the mean - colors of the two regions. It represents the Euclidian distance in + colors of the two regions. It represents the Euclidean distance in their average color. 'similarity' : The weight between two adjacent is From ab60edcdf73f3c7a4e5e0b91866ec1e5a08e8295 Mon Sep 17 00:00:00 2001 From: Vighnesh Birodkar Date: Sun, 10 Aug 2014 20:31:22 +0530 Subject: [PATCH 33/36] Fix Reference format --- skimage/graph/graph_cut.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/skimage/graph/graph_cut.py b/skimage/graph/graph_cut.py index 56939b7e..97ca7507 100644 --- a/skimage/graph/graph_cut.py +++ b/skimage/graph/graph_cut.py @@ -112,7 +112,7 @@ def cut_normalized(labels, rag, thresh=0.001, num_cuts=10, in_place=True): ---------- .. [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, August 2000 + IEEE Transactions on, vol. 22, no. 8, pp. 888-905, August 2000. """ if not in_place: From 6ef01cb4dca2ddc84e99f5c0117cf05dce2341b0 Mon Sep 17 00:00:00 2001 From: Vighnesh Birodkar Date: Wed, 13 Aug 2014 02:12:51 +0530 Subject: [PATCH 34/36] change ref format in example --- doc/examples/plot_ncut.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/examples/plot_ncut.py b/doc/examples/plot_ncut.py index 3160f710..6b8a62ab 100644 --- a/doc/examples/plot_ncut.py +++ b/doc/examples/plot_ncut.py @@ -10,7 +10,7 @@ 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 + IEEE Transactions on, vol. 22, no. 8, pp. 888-905, August 2000. """ from skimage import graph, data, io, segmentation, color from matplotlib import pyplot as plt From b61763117d4bcf4ad30de6b2f67cbfb3d58e2ca1 Mon Sep 17 00:00:00 2001 From: Vighnesh Birodkar Date: Thu, 14 Aug 2014 23:45:30 +0530 Subject: [PATCH 35/36] Update _ncut_cy.pyx --- skimage/graph/_ncut_cy.pyx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/skimage/graph/_ncut_cy.pyx b/skimage/graph/_ncut_cy.pyx index de3b5c82..14f3a94a 100644 --- a/skimage/graph/_ncut_cy.pyx +++ b/skimage/graph/_ncut_cy.pyx @@ -11,7 +11,7 @@ def argmin2(cnp.double_t[:] array): Parameters ---------- - a : array + array : array The array to process. Returns From afa345f363151f72d42976090ad16bac008f4cce Mon Sep 17 00:00:00 2001 From: Vighnesh Birodkar Date: Thu, 14 Aug 2014 23:47:00 +0530 Subject: [PATCH 36/36] Update graph_cut.py --- skimage/graph/graph_cut.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/skimage/graph/graph_cut.py b/skimage/graph/graph_cut.py index 97ca7507..3f973e22 100644 --- a/skimage/graph/graph_cut.py +++ b/skimage/graph/graph_cut.py @@ -77,7 +77,7 @@ def cut_normalized(labels, rag, thresh=0.001, num_cuts=10, in_place=True): Given an image's labels and its similarity RAG, recursively perform a 2-way normalized cut on it. All nodes belonging to a subgraph - which cannot be cut further, are assigned a unique label in the + that cannot be cut further are assigned a unique label in the output. Parameters