scikit-image

1
import numpy as np
2
from scipy import sparse
3
from scipy.sparse import csgraph
4
from ..morphology._util import _raveled_offsets_and_distances
5
from ..util._map_array import map_array
6

7

8
def _weighted_abs_diff(values0, values1, distances):
9
    """A default edge function for complete image graphs.
10

11
    A pixel graph on an image with no edge values and no mask is a very
12
    boring regular lattice, so we define a default edge weight to be the
13
    absolute difference between values *weighted* by the distance
14
    between them.
15

16
    Parameters
17
    ----------
18
    values0 : array
19
        The pixel values for each node.
20
    values1 : array
21
        The pixel values for each neighbor.
22
    distances : array
23
        The distance between each node and its neighbor.
24

25
    Returns
26
    -------
27
    edge_values : array of float
28
        The computed values: abs(values0 - values1) * distances.
29
    """
30
    return np.abs(values0 - values1) * distances
31

32

33
def pixel_graph(image, *, mask=None, edge_function=None, connectivity=1, spacing=None):
34
    """Create an adjacency graph of pixels in an image.
35

36
    Pixels where the mask is True are nodes in the returned graph, and they are
37
    connected by edges to their neighbors according to the connectivity
38
    parameter. By default, the *value* of an edge when a mask is given, or when
39
    the image is itself the mask, is the Euclidean distance between the pixels.
40

41
    However, if an int- or float-valued image is given with no mask, the value
42
    of the edges is the absolute difference in intensity between adjacent
43
    pixels, weighted by the Euclidean distance.
44

45
    Parameters
46
    ----------
47
    image : array
48
        The input image. If the image is of type bool, it will be used as the
49
        mask as well.
50
    mask : array of bool
51
        Which pixels to use. If None, the graph for the whole image is used.
52
    edge_function : callable
53
        A function taking an array of pixel values, and an array of neighbor
54
        pixel values, and an array of distances, and returning a value for the
55
        edge. If no function is given, the value of an edge is just the
56
        distance.
57
    connectivity : int
58
        The square connectivity of the pixel neighborhood: the number of
59
        orthogonal steps allowed to consider a pixel a neighbor. See
60
        `scipy.ndimage.generate_binary_structure` for details.
61
    spacing : tuple of float
62
        The spacing between pixels along each axis.
63

64
    Returns
65
    -------
66
    graph : scipy.sparse.csr_matrix
67
        A sparse adjacency matrix in which entry (i, j) is 1 if nodes i and j
68
        are neighbors, 0 otherwise.
69
    nodes : array of int
70
        The nodes of the graph. These correspond to the raveled indices of the
71
        nonzero pixels in the mask.
72
    """
73
    if mask is None:
74
        if image.dtype == bool:
75
            mask = image
76
        else:
77
            mask = np.ones_like(image, dtype=bool)
78

79
    if edge_function is None:
80
        if image.dtype == bool:
81

82
            def edge_function(x, y, distances):
83
                return distances
84

85
        else:
86
            edge_function = _weighted_abs_diff
87

88
    # Strategy: we are going to build the (i, j, data) arrays of a scipy
89
    # sparse COO matrix, then convert to CSR (which is fast).
90
    # - grab the raveled IDs of the foreground (mask == True) parts of the
91
    #   image **in the padded space**.
92
    # - broadcast them together with the raveled offsets to their neighbors.
93
    #   This gives us for each foreground pixel a list of neighbors (that
94
    #   may or may not be selected by the mask). (We also track the *distance*
95
    #   to each neighbor.)
96
    # - select "valid" entries in the neighbors and distance arrays by indexing
97
    #   into the mask, which we can do since these are raveled indices.
98
    # - use np.repeat() to repeat each source index according to the number
99
    #   of neighbors selected by the mask it has. Each of these repeated
100
    #   indices will be lined up with its neighbor, i.e. **this is the i
101
    #   array** of the COO format matrix.
102
    # - use the mask as a boolean index to get a 1D view of the selected
103
    #   neighbors. **This is the j array.**
104
    # - by default, the same boolean indexing can be applied to the distances
105
    #   to each neighbor, to give the **data array.** Optionally, a
106
    #   provided edge function can be computed on the pixel values and the
107
    #   distances to give a different value for the edges.
108
    # Note, we use map_array to map the raveled coordinates in the padded
109
    # image to the ones in the original image, and those are the returned
110
    # nodes.
111
    padded = np.pad(mask, 1, mode='constant', constant_values=False)
112
    nodes_padded = np.flatnonzero(padded)
113
    neighbor_offsets_padded, distances_padded = _raveled_offsets_and_distances(
114
        padded.shape, connectivity=connectivity, spacing=spacing
115
    )
116
    neighbors_padded = nodes_padded[:, np.newaxis] + neighbor_offsets_padded
117
    neighbor_distances_full = np.broadcast_to(distances_padded, neighbors_padded.shape)
118
    nodes = np.flatnonzero(mask)
119
    nodes_sequential = np.arange(nodes.size)
120
    # neighbors outside the mask get mapped to 0, which is a valid index,
121
    # BUT, they will be masked out in the next step.
122
    neighbors = map_array(neighbors_padded, nodes_padded, nodes)
123
    neighbors_mask = padded.reshape(-1)[neighbors_padded]
124
    num_neighbors = np.sum(neighbors_mask, axis=1)
125
    indices = np.repeat(nodes, num_neighbors)
126
    indices_sequential = np.repeat(nodes_sequential, num_neighbors)
127
    neighbor_indices = neighbors[neighbors_mask]
128
    neighbor_distances = neighbor_distances_full[neighbors_mask]
129
    neighbor_indices_sequential = map_array(neighbor_indices, nodes, nodes_sequential)
130

131
    image_r = image.reshape(-1)
132
    data = edge_function(
133
        image_r[indices], image_r[neighbor_indices], neighbor_distances
134
    )
135

136
    m = nodes_sequential.size
137
    mat = sparse.coo_matrix(
138
        (data, (indices_sequential, neighbor_indices_sequential)), shape=(m, m)
139
    )
140
    graph = mat.tocsr()
141
    return graph, nodes
142

143

144
def central_pixel(graph, nodes=None, shape=None, partition_size=100):
145
    """Find the pixel with the highest closeness centrality.
146

147
    Closeness centrality is the inverse of the total sum of shortest distances
148
    from a node to every other node.
149

150
    Parameters
151
    ----------
152
    graph : scipy.sparse.csr_matrix
153
        The sparse matrix representation of the graph.
154
    nodes : array of int
155
        The raveled index of each node in graph in the image. If not provided,
156
        the returned value will be the index in the input graph.
157
    shape : tuple of int
158
        The shape of the image in which the nodes are embedded. If provided,
159
        the returned coordinates are a NumPy multi-index of the same
160
        dimensionality as the input shape. Otherwise, the returned coordinate
161
        is the raveled index provided in `nodes`.
162
    partition_size : int
163
        This function computes the shortest path distance between every pair
164
        of nodes in the graph. This can result in a very large (N*N) matrix.
165
        As a simple performance tweak, the distance values are computed in
166
        lots of `partition_size`, resulting in a memory requirement of only
167
        partition_size*N.
168

169
    Returns
170
    -------
171
    position : int or tuple of int
172
        If shape is given, the coordinate of the central pixel in the image.
173
        Otherwise, the raveled index of that pixel.
174
    distances : array of float
175
        The total sum of distances from each node to each other reachable
176
        node.
177
    """
178
    if nodes is None:
179
        nodes = np.arange(graph.shape[0])
180
    if partition_size is None:
181
        num_splits = 1
182
    else:
183
        num_splits = max(2, graph.shape[0] // partition_size)
184
    idxs = np.arange(graph.shape[0])
185
    total_shortest_path_len_list = []
186
    for partition in np.array_split(idxs, num_splits):
187
        shortest_paths = csgraph.shortest_path(graph, directed=False, indices=partition)
188
        shortest_paths_no_inf = np.nan_to_num(shortest_paths)
189
        total_shortest_path_len_list.append(np.sum(shortest_paths_no_inf, axis=1))
190
    total_shortest_path_len = np.concatenate(total_shortest_path_len_list)
191
    nonzero = np.flatnonzero(total_shortest_path_len)
192
    min_sp = np.argmin(total_shortest_path_len[nonzero])
193
    raveled_index = nodes[nonzero[min_sp]]
194
    if shape is not None:
195
        central = np.unravel_index(raveled_index, shape)
196
    else:
197
        central = raveled_index
198
    return central, total_shortest_path_len
199