peak.py 14 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420
  1. from warnings import warn
  2. import numpy as np
  3. import scipy.ndimage as ndi
  4. from .. import measure
  5. from .._shared.coord import ensure_spacing
  6. def _get_high_intensity_peaks(image, mask, num_peaks, min_distance, p_norm):
  7. """
  8. Return the highest intensity peak coordinates.
  9. """
  10. # get coordinates of peaks
  11. coord = np.nonzero(mask)
  12. intensities = image[coord]
  13. # Highest peak first
  14. idx_maxsort = np.argsort(-intensities, kind="stable")
  15. coord = np.transpose(coord)[idx_maxsort]
  16. if np.isfinite(num_peaks):
  17. max_out = int(num_peaks)
  18. else:
  19. max_out = None
  20. if min_distance > 1:
  21. coord = ensure_spacing(
  22. coord, spacing=min_distance, p_norm=p_norm, max_out=max_out
  23. )
  24. if len(coord) > num_peaks:
  25. coord = coord[:num_peaks]
  26. return coord
  27. def _get_peak_mask(image, footprint, threshold, mask=None):
  28. """
  29. Return the mask containing all peak candidates above thresholds.
  30. """
  31. if footprint.size == 1 or image.size == 1:
  32. return image > threshold
  33. image_max = ndi.maximum_filter(image, footprint=footprint, mode='nearest')
  34. out = image == image_max
  35. # no peak for a trivial image
  36. image_is_trivial = np.all(out) if mask is None else np.all(out[mask])
  37. if image_is_trivial:
  38. out[:] = False
  39. if mask is not None:
  40. # isolated pixels in masked area are returned as peaks
  41. isolated_px = np.logical_xor(mask, ndi.binary_opening(mask))
  42. out[isolated_px] = True
  43. out &= image > threshold
  44. return out
  45. def _exclude_border(label, border_width):
  46. """Set label border values to 0."""
  47. # zero out label borders
  48. for i, width in enumerate(border_width):
  49. if width == 0:
  50. continue
  51. label[(slice(None),) * i + (slice(None, width),)] = 0
  52. label[(slice(None),) * i + (slice(-width, None),)] = 0
  53. return label
  54. def _get_threshold(image, threshold_abs, threshold_rel):
  55. """Return the threshold value according to an absolute and a relative
  56. value.
  57. """
  58. threshold = threshold_abs if threshold_abs is not None else image.min()
  59. if threshold_rel is not None:
  60. threshold = max(threshold, threshold_rel * image.max())
  61. return threshold
  62. def _get_excluded_border_width(image, min_distance, exclude_border):
  63. """Return border_width values relative to a min_distance if requested."""
  64. if isinstance(exclude_border, bool):
  65. border_width = (min_distance if exclude_border else 0,) * image.ndim
  66. elif isinstance(exclude_border, int):
  67. if exclude_border < 0:
  68. raise ValueError("`exclude_border` cannot be a negative value")
  69. border_width = (exclude_border,) * image.ndim
  70. elif isinstance(exclude_border, tuple):
  71. if len(exclude_border) != image.ndim:
  72. raise ValueError(
  73. "`exclude_border` should have the same length as the "
  74. "dimensionality of the image."
  75. )
  76. for exclude in exclude_border:
  77. if not isinstance(exclude, int):
  78. raise ValueError(
  79. "`exclude_border`, when expressed as a tuple, must only "
  80. "contain ints."
  81. )
  82. if exclude < 0:
  83. raise ValueError("`exclude_border` can not be a negative value")
  84. border_width = exclude_border
  85. else:
  86. raise TypeError(
  87. "`exclude_border` must be bool, int, or tuple with the same "
  88. "length as the dimensionality of the image."
  89. )
  90. return border_width
  91. def peak_local_max(
  92. image,
  93. min_distance=1,
  94. threshold_abs=None,
  95. threshold_rel=None,
  96. exclude_border=True,
  97. num_peaks=np.inf,
  98. footprint=None,
  99. labels=None,
  100. num_peaks_per_label=np.inf,
  101. p_norm=np.inf,
  102. ):
  103. """Find peaks in an image as coordinate list.
  104. Peaks are the local maxima in a region of `2 * min_distance + 1`
  105. (i.e. peaks are separated by at least `min_distance`).
  106. If both `threshold_abs` and `threshold_rel` are provided, the maximum
  107. of the two is chosen as the minimum intensity threshold of peaks.
  108. .. versionchanged:: 0.18
  109. Prior to version 0.18, peaks of the same height within a radius of
  110. `min_distance` were all returned, but this could cause unexpected
  111. behaviour. From 0.18 onwards, an arbitrary peak within the region is
  112. returned. See issue gh-2592.
  113. Parameters
  114. ----------
  115. image : ndarray
  116. Input image.
  117. min_distance : int, optional
  118. The minimal allowed distance separating peaks. To find the
  119. maximum number of peaks, use `min_distance=1`.
  120. threshold_abs : float or None, optional
  121. Minimum intensity of peaks. By default, the absolute threshold is
  122. the minimum intensity of the image.
  123. threshold_rel : float or None, optional
  124. Minimum intensity of peaks, calculated as
  125. ``max(image) * threshold_rel``.
  126. exclude_border : int, tuple of ints, or bool, optional
  127. If positive integer, `exclude_border` excludes peaks from within
  128. `exclude_border`-pixels of the border of the image.
  129. If tuple of non-negative ints, the length of the tuple must match the
  130. input array's dimensionality. Each element of the tuple will exclude
  131. peaks from within `exclude_border`-pixels of the border of the image
  132. along that dimension.
  133. If True, takes the `min_distance` parameter as value.
  134. If zero or False, peaks are identified regardless of their distance
  135. from the border.
  136. num_peaks : int, optional
  137. Maximum number of peaks. When the number of peaks exceeds `num_peaks`,
  138. return `num_peaks` peaks based on highest peak intensity.
  139. footprint : ndarray of bools, optional
  140. If provided, `footprint == 1` represents the local region within which
  141. to search for peaks at every point in `image`.
  142. labels : ndarray of ints, optional
  143. If provided, each unique region `labels == value` represents a unique
  144. region to search for peaks. Zero is reserved for background.
  145. num_peaks_per_label : int, optional
  146. Maximum number of peaks for each label.
  147. p_norm : float
  148. Which Minkowski p-norm to use. Should be in the range [1, inf].
  149. A finite large p may cause a ValueError if overflow can occur.
  150. ``inf`` corresponds to the Chebyshev distance and 2 to the
  151. Euclidean distance.
  152. Returns
  153. -------
  154. output : ndarray
  155. The coordinates of the peaks.
  156. Notes
  157. -----
  158. The peak local maximum function returns the coordinates of local peaks
  159. (maxima) in an image. Internally, a maximum filter is used for finding
  160. local maxima. This operation dilates the original image. After comparison
  161. of the dilated and original images, this function returns the coordinates
  162. of the peaks where the dilated image equals the original image.
  163. See also
  164. --------
  165. skimage.feature.corner_peaks
  166. Examples
  167. --------
  168. >>> img1 = np.zeros((7, 7))
  169. >>> img1[3, 4] = 1
  170. >>> img1[3, 2] = 1.5
  171. >>> img1
  172. array([[0. , 0. , 0. , 0. , 0. , 0. , 0. ],
  173. [0. , 0. , 0. , 0. , 0. , 0. , 0. ],
  174. [0. , 0. , 0. , 0. , 0. , 0. , 0. ],
  175. [0. , 0. , 1.5, 0. , 1. , 0. , 0. ],
  176. [0. , 0. , 0. , 0. , 0. , 0. , 0. ],
  177. [0. , 0. , 0. , 0. , 0. , 0. , 0. ],
  178. [0. , 0. , 0. , 0. , 0. , 0. , 0. ]])
  179. >>> peak_local_max(img1, min_distance=1)
  180. array([[3, 2],
  181. [3, 4]])
  182. >>> peak_local_max(img1, min_distance=2)
  183. array([[3, 2]])
  184. >>> img2 = np.zeros((20, 20, 20))
  185. >>> img2[10, 10, 10] = 1
  186. >>> img2[15, 15, 15] = 1
  187. >>> peak_idx = peak_local_max(img2, exclude_border=0)
  188. >>> peak_idx
  189. array([[10, 10, 10],
  190. [15, 15, 15]])
  191. >>> peak_mask = np.zeros_like(img2, dtype=bool)
  192. >>> peak_mask[tuple(peak_idx.T)] = True
  193. >>> np.argwhere(peak_mask)
  194. array([[10, 10, 10],
  195. [15, 15, 15]])
  196. """
  197. if (footprint is None or footprint.size == 1) and min_distance < 1:
  198. warn(
  199. "When min_distance < 1, peak_local_max acts as finding "
  200. "image > max(threshold_abs, threshold_rel * max(image)).",
  201. RuntimeWarning,
  202. stacklevel=2,
  203. )
  204. border_width = _get_excluded_border_width(image, min_distance, exclude_border)
  205. threshold = _get_threshold(image, threshold_abs, threshold_rel)
  206. if footprint is None:
  207. size = 2 * min_distance + 1
  208. footprint = np.ones((size,) * image.ndim, dtype=bool)
  209. else:
  210. footprint = np.asarray(footprint)
  211. if labels is None:
  212. # Non maximum filter
  213. mask = _get_peak_mask(image, footprint, threshold)
  214. mask = _exclude_border(mask, border_width)
  215. # Select highest intensities (num_peaks)
  216. coordinates = _get_high_intensity_peaks(
  217. image, mask, num_peaks, min_distance, p_norm
  218. )
  219. else:
  220. _labels = _exclude_border(labels.astype(int, casting="safe"), border_width)
  221. if np.issubdtype(image.dtype, np.floating):
  222. bg_val = np.finfo(image.dtype).min
  223. else:
  224. bg_val = np.iinfo(image.dtype).min
  225. # For each label, extract a smaller image enclosing the object of
  226. # interest, identify num_peaks_per_label peaks
  227. labels_peak_coord = []
  228. for label_idx, roi in enumerate(ndi.find_objects(_labels)):
  229. if roi is None:
  230. continue
  231. # Get roi mask
  232. label_mask = labels[roi] == label_idx + 1
  233. # Extract image roi
  234. img_object = image[roi].copy()
  235. # Ensure masked values don't affect roi's local peaks
  236. img_object[np.logical_not(label_mask)] = bg_val
  237. mask = _get_peak_mask(img_object, footprint, threshold, label_mask)
  238. coordinates = _get_high_intensity_peaks(
  239. img_object, mask, num_peaks_per_label, min_distance, p_norm
  240. )
  241. # transform coordinates in global image indices space
  242. for idx, s in enumerate(roi):
  243. coordinates[:, idx] += s.start
  244. labels_peak_coord.append(coordinates)
  245. if labels_peak_coord:
  246. coordinates = np.vstack(labels_peak_coord)
  247. else:
  248. coordinates = np.empty((0, 2), dtype=int)
  249. if len(coordinates) > num_peaks:
  250. out = np.zeros_like(image, dtype=bool)
  251. out[tuple(coordinates.T)] = True
  252. coordinates = _get_high_intensity_peaks(
  253. image, out, num_peaks, min_distance, p_norm
  254. )
  255. return coordinates
  256. def _prominent_peaks(
  257. image, min_xdistance=1, min_ydistance=1, threshold=None, num_peaks=np.inf
  258. ):
  259. """Return peaks with non-maximum suppression.
  260. Identifies most prominent features separated by certain distances.
  261. Non-maximum suppression with different sizes is applied separately
  262. in the first and second dimension of the image to identify peaks.
  263. Parameters
  264. ----------
  265. image : (M, N) ndarray
  266. Input image.
  267. min_xdistance : int
  268. Minimum distance separating features in the x dimension.
  269. min_ydistance : int
  270. Minimum distance separating features in the y dimension.
  271. threshold : float
  272. Minimum intensity of peaks. Default is `0.5 * max(image)`.
  273. num_peaks : int
  274. Maximum number of peaks. When the number of peaks exceeds `num_peaks`,
  275. return `num_peaks` coordinates based on peak intensity.
  276. Returns
  277. -------
  278. intensity, xcoords, ycoords : tuple of array
  279. Peak intensity values, x and y indices.
  280. """
  281. img = image.copy()
  282. rows, cols = img.shape
  283. if threshold is None:
  284. threshold = 0.5 * np.max(img)
  285. ycoords_size = 2 * min_ydistance + 1
  286. xcoords_size = 2 * min_xdistance + 1
  287. img_max = ndi.maximum_filter1d(
  288. img, size=ycoords_size, axis=0, mode='constant', cval=0
  289. )
  290. img_max = ndi.maximum_filter1d(
  291. img_max, size=xcoords_size, axis=1, mode='constant', cval=0
  292. )
  293. mask = img == img_max
  294. img *= mask
  295. img_t = img > threshold
  296. label_img = measure.label(img_t)
  297. props = measure.regionprops(label_img, img_max)
  298. # Sort the list of peaks by intensity, not left-right, so larger peaks
  299. # in Hough space cannot be arbitrarily suppressed by smaller neighbors
  300. props = sorted(props, key=lambda x: x.intensity_max)[::-1]
  301. coords = np.array([np.round(p.centroid) for p in props], dtype=int)
  302. img_peaks = []
  303. ycoords_peaks = []
  304. xcoords_peaks = []
  305. # relative coordinate grid for local neighborhood suppression
  306. ycoords_ext, xcoords_ext = np.mgrid[
  307. -min_ydistance : min_ydistance + 1, -min_xdistance : min_xdistance + 1
  308. ]
  309. for ycoords_idx, xcoords_idx in coords:
  310. accum = img_max[ycoords_idx, xcoords_idx]
  311. if accum > threshold:
  312. # absolute coordinate grid for local neighborhood suppression
  313. ycoords_nh = ycoords_idx + ycoords_ext
  314. xcoords_nh = xcoords_idx + xcoords_ext
  315. # no reflection for distance neighborhood
  316. ycoords_in = np.logical_and(ycoords_nh > 0, ycoords_nh < rows)
  317. ycoords_nh = ycoords_nh[ycoords_in]
  318. xcoords_nh = xcoords_nh[ycoords_in]
  319. # reflect xcoords and assume xcoords are continuous,
  320. # e.g. for angles:
  321. # (..., 88, 89, -90, -89, ..., 89, -90, -89, ...)
  322. xcoords_low = xcoords_nh < 0
  323. ycoords_nh[xcoords_low] = rows - ycoords_nh[xcoords_low]
  324. xcoords_nh[xcoords_low] += cols
  325. xcoords_high = xcoords_nh >= cols
  326. ycoords_nh[xcoords_high] = rows - ycoords_nh[xcoords_high]
  327. xcoords_nh[xcoords_high] -= cols
  328. # suppress neighborhood
  329. img_max[ycoords_nh, xcoords_nh] = 0
  330. # add current feature to peaks
  331. img_peaks.append(accum)
  332. ycoords_peaks.append(ycoords_idx)
  333. xcoords_peaks.append(xcoords_idx)
  334. img_peaks = np.array(img_peaks)
  335. ycoords_peaks = np.array(ycoords_peaks)
  336. xcoords_peaks = np.array(xcoords_peaks)
  337. if num_peaks < len(img_peaks):
  338. idx_maxsort = np.argsort(img_peaks)[::-1][:num_peaks]
  339. img_peaks = img_peaks[idx_maxsort]
  340. ycoords_peaks = ycoords_peaks[idx_maxsort]
  341. xcoords_peaks = xcoords_peaks[idx_maxsort]
  342. return img_peaks, xcoords_peaks, ycoords_peaks