| 123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319 |
- import networkx as nx
- import numpy as np
- from scipy.sparse import linalg
- from skimage._shared.compat import SCIPY_GE_1_17_0_DEV0
- from . import _ncut, _ncut_cy
- def cut_threshold(labels, rag, thresh, in_place=True):
- """Combine regions separated by weight less than threshold.
- Given an image's labels and its RAG, output new labels by
- combining regions whose nodes are separated by a weight less
- than the given threshold.
- Parameters
- ----------
- labels : ndarray
- The array of labels.
- rag : RAG
- The region adjacency graph.
- 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
- -------
- out : ndarray
- The new labelled array.
- Examples
- --------
- >>> from skimage import data, segmentation, graph
- >>> img = data.astronaut()
- >>> labels = segmentation.slic(img)
- >>> rag = graph.rag_mean_color(img, labels)
- >>> new_labels = graph.cut_threshold(labels, rag, 10)
- References
- ----------
- .. [1] Alain Tremeau and Philippe Colantoni
- "Regions Adjacency Graph Applied To Color Image Segmentation"
- :DOI:`10.1109/83.841950`
- """
- 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(data=True) if d['weight'] >= thresh]
- rag.remove_edges_from(to_remove)
- comps = nx.connected_components(rag)
- # We construct an array which can map old labels to the new ones.
- # All the labels within a connected component are assigned to a single
- # label in the output.
- map_array = np.arange(labels.max() + 1, dtype=labels.dtype)
- for i, nodes in enumerate(comps):
- for node in nodes:
- for label in rag.nodes[node]['labels']:
- map_array[label] = i
- return map_array[labels]
- def cut_normalized(
- labels,
- rag,
- thresh=0.001,
- num_cuts=10,
- in_place=True,
- max_edge=1.0,
- *,
- rng=None,
- ):
- """Perform Normalized Graph cut on the Region Adjacency Graph.
- Given an image's labels and its similarity RAG, recursively perform
- a 2-way normalized cut on it. All nodes belonging to a subgraph
- that 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`.
- 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.nodes[n]['ncut label']``.
- max_edge : float, optional
- The maximum possible value of an edge in the RAG. This corresponds to
- an edge between identical regions. This is used to put self
- edges in the RAG.
- rng : {`numpy.random.Generator`, int}, optional
- Pseudo-random number generator.
- By default, a PCG64 generator is used (see :func:`numpy.random.default_rng`).
- If `rng` is an int, it is used to seed the generator.
- The `rng` is used to determine the starting point
- of `scipy.sparse.linalg.eigsh`.
- Returns
- -------
- out : ndarray
- The new labeled array.
- Examples
- --------
- >>> from skimage import data, segmentation, graph
- >>> img = data.astronaut()
- >>> labels = segmentation.slic(img)
- >>> rag = graph.rag_mean_color(img, labels, mode='similarity')
- >>> new_labels = graph.cut_normalized(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, August 2000.
- """
- rng = np.random.default_rng(rng)
- if not in_place:
- rag = rag.copy()
- for node in rag.nodes():
- rag.add_edge(node, node, weight=max_edge)
- _ncut_relabel(rag, thresh, num_cuts, rng)
- map_array = np.zeros(labels.max() + 1, dtype=labels.dtype)
- # Mapping from old labels to new
- for n, d in rag.nodes(data=True):
- map_array[d['labels']] = d['ncut label']
- return map_array[labels]
- def partition_by_cut(cut, rag):
- """Compute resulting subgraphs from given bi-partition.
- 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.
- """
- # `cut` is derived from `D` and `W` matrices, which also follow the
- # ordering returned by `rag.nodes()` because we use
- # nx.to_scipy_sparse_array.
- # 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]]
- sub1 = rag.subgraph(nodes1)
- sub2 = rag.subgraph(nodes2)
- return sub1, sub2
- def get_min_ncut(ev, d, w, num_cuts):
- """Threshold an eigenvector 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
- -------
- mask : array
- The array of booleans which denotes the bi-partition.
- mcut : float
- The value of the minimum ncut.
- """
- mcut = np.inf
- mn = ev.min()
- mx = ev.max()
- # If all values in `ev` are equal, it implies that the graph can't be
- # further sub-divided. In this case the bi-partition is the the graph
- # itself and an empty set.
- min_mask = np.zeros_like(ev, dtype=bool)
- if np.allclose(mn, mx):
- return min_mask, mcut
- # 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(mn, mx, num_cuts, endpoint=False):
- mask = ev > t
- cost = _ncut.ncut_cost(mask, d, w)
- if cost < mcut:
- min_mask = mask
- mcut = cost
- return min_mask, mcut
- def _label_all(rag, attr_name):
- """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`
- attribute of all the nodes.
- rag : RAG
- The Region Adjacency Graph.
- attr_name : string
- The attribute to which a unique integer is assigned.
- """
- node = min(rag.nodes())
- new_label = rag.nodes[node]['labels'][0]
- for n, d in rag.nodes(data=True):
- d[attr_name] = new_label
- def _ncut_relabel(rag, thresh, num_cuts, random_generator):
- """Perform Normalized Graph cut on the Region Adjacency Graph.
- Recursively partition the graph into 2, until further subdivision
- yields a cut greater 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.
- Parameters
- ----------
- 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`.
- num_cuts : int
- The number or N-cuts to perform before determining the optimal one.
- random_generator : `numpy.random.Generator`
- Provides initial values for eigenvalue solver.
- """
- d, w = _ncut.DW_matrices(rag)
- m = w.shape[0]
- if (m > 2) and (d != w).nnz > 0:
- # This avoids further segmenting a graph that is too small,
- # and the degenerate case (d == w), which typically occurs
- # when only three single pixels remain.
- #
- # We're not sure exactly why this latter case arises. For
- # SciPy <= 0.14, SciPy continued to compute an eigenvector,
- # but newer versions (correctly) won't. We refuse to guess,
- # and stop further segmentation.
- #
- # It may make sense to a warning here; on the other hand segmentations
- # are not a ground truth, so this level of "noise" should be acceptable.
- d2 = d.copy()
- # 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)
- # Refer Shi & Malik 2001, Equation 7, Page 891
- A = d2 @ (d - w) @ d2
- # Initialize the vector to ensure reproducibility.
- v0 = random_generator.random(A.shape[0])
- # SciPy 1.17.0.dev0 adds the new `rng` keyword, allowing `eigsh` to
- # become deterministic
- rng_kw = {"rng": random_generator} if SCIPY_GE_1_17_0_DEV0 else {}
- vals, vectors = linalg.eigsh(A, which='SM', v0=v0, k=min(100, m - 2), **rng_kw)
- # 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 = vectors[:, index2]
- cut_mask, mcut = get_min_ncut(ev, d, w, num_cuts)
- if mcut < thresh:
- # 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, random_generator)
- _ncut_relabel(sub2, thresh, num_cuts, random_generator)
- 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.
- _label_all(rag, 'ncut label')
|