profile.py 6.6 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190
  1. import numpy as np
  2. from scipy import ndimage as ndi
  3. from .._shared.utils import _validate_interpolation_order, _fix_ndimage_mode
  4. def profile_line(
  5. image,
  6. src,
  7. dst,
  8. linewidth=1,
  9. order=None,
  10. mode='reflect',
  11. cval=0.0,
  12. *,
  13. reduce_func=np.mean,
  14. ):
  15. """Return the intensity profile of an image measured along a scan line.
  16. Parameters
  17. ----------
  18. image : ndarray, shape (M, N[, C])
  19. The image, either grayscale (2D array) or multichannel
  20. (3D array, where the final axis contains the channel
  21. information).
  22. src : array_like, shape (2,)
  23. The coordinates of the start point of the scan line.
  24. dst : array_like, shape (2,)
  25. The coordinates of the end point of the scan
  26. line. The destination point is *included* in the profile, in
  27. contrast to standard numpy indexing.
  28. linewidth : int, optional
  29. Width of the scan, perpendicular to the line
  30. order : int in {0, 1, 2, 3, 4, 5}, optional
  31. The order of the spline interpolation, default is 0 if
  32. image.dtype is bool and 1 otherwise. The order has to be in
  33. the range 0-5. See `skimage.transform.warp` for detail.
  34. mode : {'constant', 'nearest', 'reflect', 'mirror', 'wrap'}, optional
  35. How to compute any values falling outside of the image.
  36. cval : float, optional
  37. If `mode` is 'constant', what constant value to use outside the image.
  38. reduce_func : callable, optional
  39. Function used to calculate the aggregation of pixel values
  40. perpendicular to the profile_line direction when `linewidth` > 1.
  41. If set to None the unreduced array will be returned.
  42. Returns
  43. -------
  44. return_value : array
  45. The intensity profile along the scan line. The length of the profile
  46. is the ceil of the computed length of the scan line.
  47. Examples
  48. --------
  49. >>> x = np.array([[1, 1, 1, 2, 2, 2]])
  50. >>> img = np.vstack([np.zeros_like(x), x, x, x, np.zeros_like(x)])
  51. >>> img
  52. array([[0, 0, 0, 0, 0, 0],
  53. [1, 1, 1, 2, 2, 2],
  54. [1, 1, 1, 2, 2, 2],
  55. [1, 1, 1, 2, 2, 2],
  56. [0, 0, 0, 0, 0, 0]])
  57. >>> profile_line(img, (2, 1), (2, 4))
  58. array([1., 1., 2., 2.])
  59. >>> profile_line(img, (1, 0), (1, 6), cval=4)
  60. array([1., 1., 1., 2., 2., 2., 2.])
  61. The destination point is included in the profile, in contrast to
  62. standard numpy indexing.
  63. For example:
  64. >>> profile_line(img, (1, 0), (1, 6)) # The final point is out of bounds
  65. array([1., 1., 1., 2., 2., 2., 2.])
  66. >>> profile_line(img, (1, 0), (1, 5)) # This accesses the full first row
  67. array([1., 1., 1., 2., 2., 2.])
  68. For different reduce_func inputs:
  69. >>> profile_line(img, (1, 0), (1, 3), linewidth=3, reduce_func=np.mean)
  70. array([0.66666667, 0.66666667, 0.66666667, 1.33333333])
  71. >>> profile_line(img, (1, 0), (1, 3), linewidth=3, reduce_func=np.max)
  72. array([1, 1, 1, 2])
  73. >>> profile_line(img, (1, 0), (1, 3), linewidth=3, reduce_func=np.sum)
  74. array([2, 2, 2, 4])
  75. The unreduced array will be returned when `reduce_func` is None or when
  76. `reduce_func` acts on each pixel value individually.
  77. >>> profile_line(img, (1, 2), (4, 2), linewidth=3, order=0,
  78. ... reduce_func=None)
  79. array([[1, 1, 2],
  80. [1, 1, 2],
  81. [1, 1, 2],
  82. [0, 0, 0]])
  83. >>> profile_line(img, (1, 0), (1, 3), linewidth=3, reduce_func=np.sqrt)
  84. array([[1. , 1. , 0. ],
  85. [1. , 1. , 0. ],
  86. [1. , 1. , 0. ],
  87. [1.41421356, 1.41421356, 0. ]])
  88. """
  89. order = _validate_interpolation_order(image.dtype, order)
  90. mode = _fix_ndimage_mode(mode)
  91. perp_lines = _line_profile_coordinates(src, dst, linewidth=linewidth)
  92. if image.ndim == 3:
  93. pixels = [
  94. ndi.map_coordinates(
  95. image[..., i],
  96. perp_lines,
  97. prefilter=order > 1,
  98. order=order,
  99. mode=mode,
  100. cval=cval,
  101. )
  102. for i in range(image.shape[2])
  103. ]
  104. pixels = np.transpose(np.asarray(pixels), (1, 2, 0))
  105. else:
  106. pixels = ndi.map_coordinates(
  107. image, perp_lines, prefilter=order > 1, order=order, mode=mode, cval=cval
  108. )
  109. # The outputted array with reduce_func=None gives an array where the
  110. # row values (axis=1) are flipped. Here, we make this consistent.
  111. pixels = np.flip(pixels, axis=1)
  112. if reduce_func is None:
  113. intensities = pixels
  114. else:
  115. try:
  116. intensities = reduce_func(pixels, axis=1)
  117. except TypeError: # function doesn't allow axis kwarg
  118. intensities = np.apply_along_axis(reduce_func, arr=pixels, axis=1)
  119. return intensities
  120. def _line_profile_coordinates(src, dst, linewidth=1):
  121. """Return the coordinates of the profile of an image along a scan line.
  122. Parameters
  123. ----------
  124. src : 2-tuple of numeric scalar (float or int)
  125. The start point of the scan line.
  126. dst : 2-tuple of numeric scalar (float or int)
  127. The end point of the scan line.
  128. linewidth : int, optional
  129. Width of the scan, perpendicular to the line
  130. Returns
  131. -------
  132. coords : array, shape (2, N, C), float
  133. The coordinates of the profile along the scan line. The length of the
  134. profile is the ceil of the computed length of the scan line.
  135. Notes
  136. -----
  137. This is a utility method meant to be used internally by skimage functions.
  138. The destination point is included in the profile, in contrast to
  139. standard numpy indexing.
  140. """
  141. src_row, src_col = src = np.asarray(src, dtype=float)
  142. dst_row, dst_col = dst = np.asarray(dst, dtype=float)
  143. d_row, d_col = dst - src
  144. theta = np.arctan2(d_row, d_col)
  145. length = int(np.ceil(np.hypot(d_row, d_col) + 1))
  146. # we add one above because we include the last point in the profile
  147. # (in contrast to standard numpy indexing)
  148. line_col = np.linspace(src_col, dst_col, length)
  149. line_row = np.linspace(src_row, dst_row, length)
  150. # we subtract 1 from linewidth to change from pixel-counting
  151. # (make this line 3 pixels wide) to point distances (the
  152. # distance between pixel centers)
  153. col_width = (linewidth - 1) * np.sin(-theta) / 2
  154. row_width = (linewidth - 1) * np.cos(theta) / 2
  155. perp_rows = np.stack(
  156. [
  157. np.linspace(row_i - row_width, row_i + row_width, linewidth)
  158. for row_i in line_row
  159. ]
  160. )
  161. perp_cols = np.stack(
  162. [
  163. np.linspace(col_i - col_width, col_i + col_width, linewidth)
  164. for col_i in line_col
  165. ]
  166. )
  167. return np.stack([perp_rows, perp_cols])