scikit-image

1
import networkx as nx
2
import numpy as np
3
from scipy.sparse import linalg
4

5
from . import _ncut, _ncut_cy
6

7

8
def cut_threshold(labels, rag, thresh, in_place=True):
9
    """Combine regions separated by weight less than threshold.
10

11
    Given an image's labels and its RAG, output new labels by
12
    combining regions whose nodes are separated by a weight less
13
    than the given threshold.
14

15
    Parameters
16
    ----------
17
    labels : ndarray
18
        The array of labels.
19
    rag : RAG
20
        The region adjacency graph.
21
    thresh : float
22
        The threshold. Regions connected by edges with smaller weights are
23
        combined.
24
    in_place : bool
25
        If set, modifies `rag` in place. The function will remove the edges
26
        with weights less that `thresh`. If set to `False` the function
27
        makes a copy of `rag` before proceeding.
28

29
    Returns
30
    -------
31
    out : ndarray
32
        The new labelled array.
33

34
    Examples
35
    --------
36
    >>> from skimage import data, segmentation, graph
37
    >>> img = data.astronaut()
38
    >>> labels = segmentation.slic(img)
39
    >>> rag = graph.rag_mean_color(img, labels)
40
    >>> new_labels = graph.cut_threshold(labels, rag, 10)
41

42
    References
43
    ----------
44
    .. [1] Alain Tremeau and Philippe Colantoni
45
           "Regions Adjacency Graph Applied To Color Image Segmentation"
46
           :DOI:`10.1109/83.841950`
47

48
    """
49
    if not in_place:
50
        rag = rag.copy()
51

52
    # Because deleting edges while iterating through them produces an error.
53
    to_remove = [(x, y) for x, y, d in rag.edges(data=True) if d['weight'] >= thresh]
54
    rag.remove_edges_from(to_remove)
55

56
    comps = nx.connected_components(rag)
57

58
    # We construct an array which can map old labels to the new ones.
59
    # All the labels within a connected component are assigned to a single
60
    # label in the output.
61
    map_array = np.arange(labels.max() + 1, dtype=labels.dtype)
62
    for i, nodes in enumerate(comps):
63
        for node in nodes:
64
            for label in rag.nodes[node]['labels']:
65
                map_array[label] = i
66

67
    return map_array[labels]
68

69

70
def cut_normalized(
71
    labels,
72
    rag,
73
    thresh=0.001,
74
    num_cuts=10,
75
    in_place=True,
76
    max_edge=1.0,
77
    *,
78
    rng=None,
79
):
80
    """Perform Normalized Graph cut on the Region Adjacency Graph.
81

82
    Given an image's labels and its similarity RAG, recursively perform
83
    a 2-way normalized cut on it. All nodes belonging to a subgraph
84
    that cannot be cut further are assigned a unique label in the
85
    output.
86

87
    Parameters
88
    ----------
89
    labels : ndarray
90
        The array of labels.
91
    rag : RAG
92
        The region adjacency graph.
93
    thresh : float
94
        The threshold. A subgraph won't be further subdivided if the
95
        value of the N-cut exceeds `thresh`.
96
    num_cuts : int
97
        The number or N-cuts to perform before determining the optimal one.
98
    in_place : bool
99
        If set, modifies `rag` in place. For each node `n` the function will
100
        set a new attribute ``rag.nodes[n]['ncut label']``.
101
    max_edge : float, optional
102
        The maximum possible value of an edge in the RAG. This corresponds to
103
        an edge between identical regions. This is used to put self
104
        edges in the RAG.
105
    rng : {`numpy.random.Generator`, int}, optional
106
        Pseudo-random number generator.
107
        By default, a PCG64 generator is used (see :func:`numpy.random.default_rng`).
108
        If `rng` is an int, it is used to seed the generator.
109

110
        The `rng` is used to determine the starting point
111
        of `scipy.sparse.linalg.eigsh`.
112

113
    Returns
114
    -------
115
    out : ndarray
116
        The new labeled array.
117

118
    Examples
119
    --------
120
    >>> from skimage import data, segmentation, graph
121
    >>> img = data.astronaut()
122
    >>> labels = segmentation.slic(img)
123
    >>> rag = graph.rag_mean_color(img, labels, mode='similarity')
124
    >>> new_labels = graph.cut_normalized(labels, rag)
125

126
    References
127
    ----------
128
    .. [1] Shi, J.; Malik, J., "Normalized cuts and image segmentation",
129
           Pattern Analysis and Machine Intelligence,
130
           IEEE Transactions on, vol. 22, no. 8, pp. 888-905, August 2000.
131

132
    """
133
    rng = np.random.default_rng(rng)
134
    if not in_place:
135
        rag = rag.copy()
136

137
    for node in rag.nodes():
138
        rag.add_edge(node, node, weight=max_edge)
139

140
    _ncut_relabel(rag, thresh, num_cuts, rng)
141

142
    map_array = np.zeros(labels.max() + 1, dtype=labels.dtype)
143
    # Mapping from old labels to new
144
    for n, d in rag.nodes(data=True):
145
        map_array[d['labels']] = d['ncut label']
146

147
    return map_array[labels]
148

149

150
def partition_by_cut(cut, rag):
151
    """Compute resulting subgraphs from given bi-partition.
152

153
    Parameters
154
    ----------
155
    cut : array
156
        A array of booleans. Elements set to `True` belong to one
157
        set.
158
    rag : RAG
159
        The Region Adjacency Graph.
160

161
    Returns
162
    -------
163
    sub1, sub2 : RAG
164
        The two resulting subgraphs from the bi-partition.
165
    """
166
    # `cut` is derived from `D` and `W` matrices, which also follow the
167
    # ordering returned by `rag.nodes()` because we use
168
    # nx.to_scipy_sparse_matrix.
169

170
    # Example
171
    # rag.nodes() = [3, 7, 9, 13]
172
    # cut = [True, False, True, False]
173
    # nodes1 = [3, 9]
174
    # nodes2 = [7, 10]
175

176
    nodes1 = [n for i, n in enumerate(rag.nodes()) if cut[i]]
177
    nodes2 = [n for i, n in enumerate(rag.nodes()) if not cut[i]]
178

179
    sub1 = rag.subgraph(nodes1)
180
    sub2 = rag.subgraph(nodes2)
181

182
    return sub1, sub2
183

184

185
def get_min_ncut(ev, d, w, num_cuts):
186
    """Threshold an eigenvector evenly, to determine minimum ncut.
187

188
    Parameters
189
    ----------
190
    ev : array
191
        The eigenvector to threshold.
192
    d : ndarray
193
        The diagonal matrix of the graph.
194
    w : ndarray
195
        The weight matrix of the graph.
196
    num_cuts : int
197
        The number of evenly spaced thresholds to check for.
198

199
    Returns
200
    -------
201
    mask : array
202
        The array of booleans which denotes the bi-partition.
203
    mcut : float
204
        The value of the minimum ncut.
205
    """
206
    mcut = np.inf
207
    mn = ev.min()
208
    mx = ev.max()
209

210
    # If all values in `ev` are equal, it implies that the graph can't be
211
    # further sub-divided. In this case the bi-partition is the the graph
212
    # itself and an empty set.
213
    min_mask = np.zeros_like(ev, dtype=bool)
214
    if np.allclose(mn, mx):
215
        return min_mask, mcut
216

217
    # Refer Shi & Malik 2001, Section 3.1.3, Page 892
218
    # Perform evenly spaced n-cuts and determine the optimal one.
219
    for t in np.linspace(mn, mx, num_cuts, endpoint=False):
220
        mask = ev > t
221
        cost = _ncut.ncut_cost(mask, d, w)
222
        if cost < mcut:
223
            min_mask = mask
224
            mcut = cost
225

226
    return min_mask, mcut
227

228

229
def _label_all(rag, attr_name):
230
    """Assign a unique integer to the given attribute in the RAG.
231

232
    This function assumes that all labels in `rag` are unique. It
233
    picks up a random label from them and assigns it to the `attr_name`
234
    attribute of all the nodes.
235

236
    rag : RAG
237
        The Region Adjacency Graph.
238
    attr_name : string
239
        The attribute to which a unique integer is assigned.
240
    """
241
    node = min(rag.nodes())
242
    new_label = rag.nodes[node]['labels'][0]
243
    for n, d in rag.nodes(data=True):
244
        d[attr_name] = new_label
245

246

247
def _ncut_relabel(rag, thresh, num_cuts, random_generator):
248
    """Perform Normalized Graph cut on the Region Adjacency Graph.
249

250
    Recursively partition the graph into 2, until further subdivision
251
    yields a cut greater than `thresh` or such a cut cannot be computed.
252
    For such a subgraph, indices to labels of all its nodes map to a single
253
    unique value.
254

255
    Parameters
256
    ----------
257
    rag : RAG
258
        The region adjacency graph.
259
    thresh : float
260
        The threshold. A subgraph won't be further subdivided if the
261
        value of the N-cut exceeds `thresh`.
262
    num_cuts : int
263
        The number or N-cuts to perform before determining the optimal one.
264
    random_generator : `numpy.random.Generator`
265
        Provides initial values for eigenvalue solver.
266
    """
267
    d, w = _ncut.DW_matrices(rag)
268
    m = w.shape[0]
269

270
    if m > 2:
271
        d2 = d.copy()
272
        # Since d is diagonal, we can directly operate on its data
273
        # the inverse of the square root
274
        d2.data = np.reciprocal(np.sqrt(d2.data, out=d2.data), out=d2.data)
275

276
        # Refer Shi & Malik 2001, Equation 7, Page 891
277
        A = d2 * (d - w) * d2
278
        # Initialize the vector to ensure reproducibility.
279
        v0 = random_generator.random(A.shape[0])
280
        vals, vectors = linalg.eigsh(A, which='SM', v0=v0, k=min(100, m - 2))
281

282
        # Pick second smallest eigenvector.
283
        # Refer Shi & Malik 2001, Section 3.2.3, Page 893
284
        vals, vectors = np.real(vals), np.real(vectors)
285
        index2 = _ncut_cy.argmin2(vals)
286
        ev = vectors[:, index2]
287

288
        cut_mask, mcut = get_min_ncut(ev, d, w, num_cuts)
289
        if mcut < thresh:
290
            # Sub divide and perform N-cut again
291
            # Refer Shi & Malik 2001, Section 3.2.5, Page 893
292
            sub1, sub2 = partition_by_cut(cut_mask, rag)
293

294
            _ncut_relabel(sub1, thresh, num_cuts, random_generator)
295
            _ncut_relabel(sub2, thresh, num_cuts, random_generator)
296
            return
297

298
    # The N-cut wasn't small enough, or could not be computed.
299
    # The remaining graph is a region.
300
    # Assign `ncut label` by picking any label from the existing nodes, since
301
    # `labels` are unique, `new_label` is also unique.
302
    _label_all(rag, 'ncut label')
303