_graph_cut.py 10 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319
  1. import networkx as nx
  2. import numpy as np
  3. from scipy.sparse import linalg
  4. from skimage._shared.compat import SCIPY_GE_1_17_0_DEV0
  5. from . import _ncut, _ncut_cy
  6. def cut_threshold(labels, rag, thresh, in_place=True):
  7. """Combine regions separated by weight less than threshold.
  8. Given an image's labels and its RAG, output new labels by
  9. combining regions whose nodes are separated by a weight less
  10. than the given threshold.
  11. Parameters
  12. ----------
  13. labels : ndarray
  14. The array of labels.
  15. rag : RAG
  16. The region adjacency graph.
  17. thresh : float
  18. The threshold. Regions connected by edges with smaller weights are
  19. combined.
  20. in_place : bool
  21. If set, modifies `rag` in place. The function will remove the edges
  22. with weights less that `thresh`. If set to `False` the function
  23. makes a copy of `rag` before proceeding.
  24. Returns
  25. -------
  26. out : ndarray
  27. The new labelled array.
  28. Examples
  29. --------
  30. >>> from skimage import data, segmentation, graph
  31. >>> img = data.astronaut()
  32. >>> labels = segmentation.slic(img)
  33. >>> rag = graph.rag_mean_color(img, labels)
  34. >>> new_labels = graph.cut_threshold(labels, rag, 10)
  35. References
  36. ----------
  37. .. [1] Alain Tremeau and Philippe Colantoni
  38. "Regions Adjacency Graph Applied To Color Image Segmentation"
  39. :DOI:`10.1109/83.841950`
  40. """
  41. if not in_place:
  42. rag = rag.copy()
  43. # Because deleting edges while iterating through them produces an error.
  44. to_remove = [(x, y) for x, y, d in rag.edges(data=True) if d['weight'] >= thresh]
  45. rag.remove_edges_from(to_remove)
  46. comps = nx.connected_components(rag)
  47. # We construct an array which can map old labels to the new ones.
  48. # All the labels within a connected component are assigned to a single
  49. # label in the output.
  50. map_array = np.arange(labels.max() + 1, dtype=labels.dtype)
  51. for i, nodes in enumerate(comps):
  52. for node in nodes:
  53. for label in rag.nodes[node]['labels']:
  54. map_array[label] = i
  55. return map_array[labels]
  56. def cut_normalized(
  57. labels,
  58. rag,
  59. thresh=0.001,
  60. num_cuts=10,
  61. in_place=True,
  62. max_edge=1.0,
  63. *,
  64. rng=None,
  65. ):
  66. """Perform Normalized Graph cut on the Region Adjacency Graph.
  67. Given an image's labels and its similarity RAG, recursively perform
  68. a 2-way normalized cut on it. All nodes belonging to a subgraph
  69. that cannot be cut further are assigned a unique label in the
  70. output.
  71. Parameters
  72. ----------
  73. labels : ndarray
  74. The array of labels.
  75. rag : RAG
  76. The region adjacency graph.
  77. thresh : float
  78. The threshold. A subgraph won't be further subdivided if the
  79. value of the N-cut exceeds `thresh`.
  80. num_cuts : int
  81. The number or N-cuts to perform before determining the optimal one.
  82. in_place : bool
  83. If set, modifies `rag` in place. For each node `n` the function will
  84. set a new attribute ``rag.nodes[n]['ncut label']``.
  85. max_edge : float, optional
  86. The maximum possible value of an edge in the RAG. This corresponds to
  87. an edge between identical regions. This is used to put self
  88. edges in the RAG.
  89. rng : {`numpy.random.Generator`, int}, optional
  90. Pseudo-random number generator.
  91. By default, a PCG64 generator is used (see :func:`numpy.random.default_rng`).
  92. If `rng` is an int, it is used to seed the generator.
  93. The `rng` is used to determine the starting point
  94. of `scipy.sparse.linalg.eigsh`.
  95. Returns
  96. -------
  97. out : ndarray
  98. The new labeled array.
  99. Examples
  100. --------
  101. >>> from skimage import data, segmentation, graph
  102. >>> img = data.astronaut()
  103. >>> labels = segmentation.slic(img)
  104. >>> rag = graph.rag_mean_color(img, labels, mode='similarity')
  105. >>> new_labels = graph.cut_normalized(labels, rag)
  106. References
  107. ----------
  108. .. [1] Shi, J.; Malik, J., "Normalized cuts and image segmentation",
  109. Pattern Analysis and Machine Intelligence,
  110. IEEE Transactions on, vol. 22, no. 8, pp. 888-905, August 2000.
  111. """
  112. rng = np.random.default_rng(rng)
  113. if not in_place:
  114. rag = rag.copy()
  115. for node in rag.nodes():
  116. rag.add_edge(node, node, weight=max_edge)
  117. _ncut_relabel(rag, thresh, num_cuts, rng)
  118. map_array = np.zeros(labels.max() + 1, dtype=labels.dtype)
  119. # Mapping from old labels to new
  120. for n, d in rag.nodes(data=True):
  121. map_array[d['labels']] = d['ncut label']
  122. return map_array[labels]
  123. def partition_by_cut(cut, rag):
  124. """Compute resulting subgraphs from given bi-partition.
  125. Parameters
  126. ----------
  127. cut : array
  128. A array of booleans. Elements set to `True` belong to one
  129. set.
  130. rag : RAG
  131. The Region Adjacency Graph.
  132. Returns
  133. -------
  134. sub1, sub2 : RAG
  135. The two resulting subgraphs from the bi-partition.
  136. """
  137. # `cut` is derived from `D` and `W` matrices, which also follow the
  138. # ordering returned by `rag.nodes()` because we use
  139. # nx.to_scipy_sparse_array.
  140. # Example
  141. # rag.nodes() = [3, 7, 9, 13]
  142. # cut = [True, False, True, False]
  143. # nodes1 = [3, 9]
  144. # nodes2 = [7, 10]
  145. nodes1 = [n for i, n in enumerate(rag.nodes()) if cut[i]]
  146. nodes2 = [n for i, n in enumerate(rag.nodes()) if not cut[i]]
  147. sub1 = rag.subgraph(nodes1)
  148. sub2 = rag.subgraph(nodes2)
  149. return sub1, sub2
  150. def get_min_ncut(ev, d, w, num_cuts):
  151. """Threshold an eigenvector evenly, to determine minimum ncut.
  152. Parameters
  153. ----------
  154. ev : array
  155. The eigenvector to threshold.
  156. d : ndarray
  157. The diagonal matrix of the graph.
  158. w : ndarray
  159. The weight matrix of the graph.
  160. num_cuts : int
  161. The number of evenly spaced thresholds to check for.
  162. Returns
  163. -------
  164. mask : array
  165. The array of booleans which denotes the bi-partition.
  166. mcut : float
  167. The value of the minimum ncut.
  168. """
  169. mcut = np.inf
  170. mn = ev.min()
  171. mx = ev.max()
  172. # If all values in `ev` are equal, it implies that the graph can't be
  173. # further sub-divided. In this case the bi-partition is the the graph
  174. # itself and an empty set.
  175. min_mask = np.zeros_like(ev, dtype=bool)
  176. if np.allclose(mn, mx):
  177. return min_mask, mcut
  178. # Refer Shi & Malik 2001, Section 3.1.3, Page 892
  179. # Perform evenly spaced n-cuts and determine the optimal one.
  180. for t in np.linspace(mn, mx, num_cuts, endpoint=False):
  181. mask = ev > t
  182. cost = _ncut.ncut_cost(mask, d, w)
  183. if cost < mcut:
  184. min_mask = mask
  185. mcut = cost
  186. return min_mask, mcut
  187. def _label_all(rag, attr_name):
  188. """Assign a unique integer to the given attribute in the RAG.
  189. This function assumes that all labels in `rag` are unique. It
  190. picks up a random label from them and assigns it to the `attr_name`
  191. attribute of all the nodes.
  192. rag : RAG
  193. The Region Adjacency Graph.
  194. attr_name : string
  195. The attribute to which a unique integer is assigned.
  196. """
  197. node = min(rag.nodes())
  198. new_label = rag.nodes[node]['labels'][0]
  199. for n, d in rag.nodes(data=True):
  200. d[attr_name] = new_label
  201. def _ncut_relabel(rag, thresh, num_cuts, random_generator):
  202. """Perform Normalized Graph cut on the Region Adjacency Graph.
  203. Recursively partition the graph into 2, until further subdivision
  204. yields a cut greater than `thresh` or such a cut cannot be computed.
  205. For such a subgraph, indices to labels of all its nodes map to a single
  206. unique value.
  207. Parameters
  208. ----------
  209. rag : RAG
  210. The region adjacency graph.
  211. thresh : float
  212. The threshold. A subgraph won't be further subdivided if the
  213. value of the N-cut exceeds `thresh`.
  214. num_cuts : int
  215. The number or N-cuts to perform before determining the optimal one.
  216. random_generator : `numpy.random.Generator`
  217. Provides initial values for eigenvalue solver.
  218. """
  219. d, w = _ncut.DW_matrices(rag)
  220. m = w.shape[0]
  221. if (m > 2) and (d != w).nnz > 0:
  222. # This avoids further segmenting a graph that is too small,
  223. # and the degenerate case (d == w), which typically occurs
  224. # when only three single pixels remain.
  225. #
  226. # We're not sure exactly why this latter case arises. For
  227. # SciPy <= 0.14, SciPy continued to compute an eigenvector,
  228. # but newer versions (correctly) won't. We refuse to guess,
  229. # and stop further segmentation.
  230. #
  231. # It may make sense to a warning here; on the other hand segmentations
  232. # are not a ground truth, so this level of "noise" should be acceptable.
  233. d2 = d.copy()
  234. # Since d is diagonal, we can directly operate on its data
  235. # the inverse of the square root
  236. d2.data = np.reciprocal(np.sqrt(d2.data, out=d2.data), out=d2.data)
  237. # Refer Shi & Malik 2001, Equation 7, Page 891
  238. A = d2 @ (d - w) @ d2
  239. # Initialize the vector to ensure reproducibility.
  240. v0 = random_generator.random(A.shape[0])
  241. # SciPy 1.17.0.dev0 adds the new `rng` keyword, allowing `eigsh` to
  242. # become deterministic
  243. rng_kw = {"rng": random_generator} if SCIPY_GE_1_17_0_DEV0 else {}
  244. vals, vectors = linalg.eigsh(A, which='SM', v0=v0, k=min(100, m - 2), **rng_kw)
  245. # Pick second smallest eigenvector.
  246. # Refer Shi & Malik 2001, Section 3.2.3, Page 893
  247. vals, vectors = np.real(vals), np.real(vectors)
  248. index2 = _ncut_cy.argmin2(vals)
  249. ev = vectors[:, index2]
  250. cut_mask, mcut = get_min_ncut(ev, d, w, num_cuts)
  251. if mcut < thresh:
  252. # Sub divide and perform N-cut again
  253. # Refer Shi & Malik 2001, Section 3.2.5, Page 893
  254. sub1, sub2 = partition_by_cut(cut_mask, rag)
  255. _ncut_relabel(sub1, thresh, num_cuts, random_generator)
  256. _ncut_relabel(sub2, thresh, num_cuts, random_generator)
  257. return
  258. # The N-cut wasn't small enough, or could not be computed.
  259. # The remaining graph is a region.
  260. # Assign `ncut label` by picking any label from the existing nodes, since
  261. # `labels` are unique, `new_label` is also unique.
  262. _label_all(rag, 'ncut label')