_graph.py 9.0 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220
  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. from ..segmentation.random_walker_segmentation import _safe_downcast_indices
  7. def _weighted_abs_diff(values0, values1, distances):
  8. """A default edge function for complete image graphs.
  9. A pixel graph on an image with no edge values and no mask is a very
  10. boring regular lattice, so we define a default edge weight to be the
  11. absolute difference between values *weighted* by the distance
  12. between them.
  13. Parameters
  14. ----------
  15. values0 : array
  16. The pixel values for each node.
  17. values1 : array
  18. The pixel values for each neighbor.
  19. distances : array
  20. The distance between each node and its neighbor.
  21. Returns
  22. -------
  23. edge_values : array of float
  24. The computed values: abs(values0 - values1) * distances.
  25. """
  26. return np.abs(values0 - values1) * distances
  27. def pixel_graph(
  28. image,
  29. *,
  30. mask=None,
  31. edge_function=None,
  32. connectivity=1,
  33. spacing=None,
  34. sparse_type="matrix",
  35. ):
  36. """Create an adjacency graph of pixels in an image.
  37. Pixels where the mask is True are nodes in the returned graph, and they are
  38. connected by edges to their neighbors according to the connectivity
  39. parameter. By default, the *value* of an edge when a mask is given, or when
  40. the image is itself the mask, is the Euclidean distance between the pixels.
  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. Parameters
  45. ----------
  46. image : array
  47. The input image. If the image is of type bool, it will be used as the
  48. mask as well.
  49. mask : array of bool
  50. Which pixels to use. If None, the graph for the whole image is used.
  51. edge_function : callable
  52. A function taking an array of pixel values, and an array of neighbor
  53. pixel values, and an array of distances, and returning a value for the
  54. edge. If no function is given, the value of an edge is just the
  55. distance.
  56. connectivity : int
  57. The square connectivity of the pixel neighborhood: the number of
  58. orthogonal steps allowed to consider a pixel a neighbor. See
  59. `scipy.ndimage.generate_binary_structure` for details.
  60. spacing : tuple of float
  61. The spacing between pixels along each axis.
  62. sparse_type : {"matrix", "array"}, optional
  63. The return type of `graph`, either `scipy.sparse.csr_array` or
  64. `scipy.sparse.csr_matrix` (default).
  65. Returns
  66. -------
  67. graph : scipy.sparse.csr_matrix or scipy.sparse.csr_array
  68. A sparse adjacency matrix in which entry (i, j) is 1 if nodes i and j
  69. are neighbors, 0 otherwise. Depending on `sparse_type`, this can be
  70. returned as a `scipy.sparse.csr_array`.
  71. nodes : array of int
  72. The nodes of the graph. These correspond to the raveled indices of the
  73. nonzero pixels in the mask.
  74. """
  75. if mask is None:
  76. if image.dtype == bool:
  77. mask = image
  78. else:
  79. mask = np.ones_like(image, dtype=bool)
  80. if edge_function is None:
  81. if image.dtype == bool:
  82. def edge_function(x, y, distances):
  83. return distances
  84. else:
  85. edge_function = _weighted_abs_diff
  86. # Strategy: we are going to build the (i, j, data) arrays of a scipy
  87. # sparse CSR matrix.
  88. # - grab the raveled IDs of the foreground (mask == True) parts of the
  89. # image **in the padded space**.
  90. # - broadcast them together with the raveled offsets to their neighbors.
  91. # This gives us for each foreground pixel a list of neighbors (that
  92. # may or may not be selected by the mask). (We also track the *distance*
  93. # to each neighbor.)
  94. # - select "valid" entries in the neighbors and distance arrays by indexing
  95. # into the mask, which we can do since these are raveled indices.
  96. # - use np.repeat() to repeat each source index according to the number
  97. # of neighbors selected by the mask it has. Each of these repeated
  98. # indices will be lined up with its neighbor, i.e. **this is the row_ind
  99. # array** of the CSR format matrix.
  100. # - use the mask as a boolean index to get a 1D view of the selected
  101. # neighbors. **This is the col_ind array.**
  102. # - by default, the same boolean indexing can be applied to the distances
  103. # to each neighbor, to give the **data array.** Optionally, a
  104. # provided edge function can be computed on the pixel values and the
  105. # distances to give a different value for the edges.
  106. # Note, we use map_array to map the raveled coordinates in the padded
  107. # image to the ones in the original image, and those are the returned
  108. # nodes.
  109. padded = np.pad(mask, 1, mode='constant', constant_values=False)
  110. nodes_padded = np.flatnonzero(padded)
  111. neighbor_offsets_padded, distances_padded = _raveled_offsets_and_distances(
  112. padded.shape, connectivity=connectivity, spacing=spacing
  113. )
  114. neighbors_padded = nodes_padded[:, np.newaxis] + neighbor_offsets_padded
  115. neighbor_distances_full = np.broadcast_to(distances_padded, neighbors_padded.shape)
  116. nodes = np.flatnonzero(mask)
  117. nodes_sequential = np.arange(nodes.size)
  118. # neighbors outside the mask get mapped to 0, which is a valid index,
  119. # BUT, they will be masked out in the next step.
  120. neighbors = map_array(neighbors_padded, nodes_padded, nodes)
  121. neighbors_mask = padded.reshape(-1)[neighbors_padded]
  122. num_neighbors = np.sum(neighbors_mask, axis=1)
  123. indices = np.repeat(nodes, num_neighbors)
  124. indices_sequential = np.repeat(nodes_sequential, num_neighbors)
  125. neighbor_indices = neighbors[neighbors_mask]
  126. neighbor_distances = neighbor_distances_full[neighbors_mask]
  127. neighbor_indices_sequential = map_array(neighbor_indices, nodes, nodes_sequential)
  128. image_r = image.reshape(-1)
  129. data = edge_function(
  130. image_r[indices], image_r[neighbor_indices], neighbor_distances
  131. )
  132. m = nodes_sequential.size
  133. graph = sparse.csr_array(
  134. (data, (indices_sequential, neighbor_indices_sequential)), shape=(m, m)
  135. )
  136. if sparse_type == "matrix":
  137. graph = sparse.csr_matrix(graph)
  138. elif sparse_type != "array":
  139. msg = f"`sparse_type` must be 'array' or 'matrix', got {sparse_type}"
  140. raise ValueError(msg)
  141. return graph, nodes
  142. def central_pixel(graph, nodes=None, shape=None, partition_size=100):
  143. """Find the pixel with the highest closeness centrality.
  144. Closeness centrality is the inverse of the total sum of shortest distances
  145. from a node to every other node.
  146. Parameters
  147. ----------
  148. graph : scipy.sparse.csr_array or scipy.sparse.csr_matrix
  149. The sparse representation of the graph.
  150. nodes : array of int
  151. The raveled index of each node in graph in the image. If not provided,
  152. the returned value will be the index in the input graph.
  153. shape : tuple of int
  154. The shape of the image in which the nodes are embedded. If provided,
  155. the returned coordinates are a NumPy multi-index of the same
  156. dimensionality as the input shape. Otherwise, the returned coordinate
  157. is the raveled index provided in `nodes`.
  158. partition_size : int
  159. This function computes the shortest path distance between every pair
  160. of nodes in the graph. This can result in a very large (N*N) matrix.
  161. As a simple performance tweak, the distance values are computed in
  162. lots of `partition_size`, resulting in a memory requirement of only
  163. partition_size*N.
  164. Returns
  165. -------
  166. position : int or tuple of int
  167. If shape is given, the coordinate of the central pixel in the image.
  168. Otherwise, the raveled index of that pixel.
  169. distances : array of float
  170. The total sum of distances from each node to each other reachable
  171. node.
  172. """
  173. if nodes is None:
  174. nodes = np.arange(graph.shape[0])
  175. if partition_size is None:
  176. num_splits = 1
  177. else:
  178. num_splits = max(2, graph.shape[0] // partition_size)
  179. graph.indices, graph.indptr = _safe_downcast_indices(
  180. graph, np.int32, 'index values too large for csgraph'
  181. )
  182. idxs = np.arange(graph.shape[0])
  183. total_shortest_path_len_list = []
  184. for partition in np.array_split(idxs, num_splits):
  185. shortest_paths = csgraph.shortest_path(graph, directed=False, indices=partition)
  186. shortest_paths_no_inf = np.nan_to_num(shortest_paths)
  187. total_shortest_path_len_list.append(np.sum(shortest_paths_no_inf, axis=1))
  188. total_shortest_path_len = np.concatenate(total_shortest_path_len_list)
  189. nonzero = np.flatnonzero(total_shortest_path_len)
  190. min_sp = np.argmin(total_shortest_path_len[nonzero])
  191. raveled_index = nodes[nonzero[min_sp]]
  192. if shape is not None:
  193. central = np.unravel_index(raveled_index, shape)
  194. else:
  195. central = raveled_index
  196. return central, total_shortest_path_len