_find_contours.py 9.4 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218
  1. import numpy as np
  2. from ._find_contours_cy import _get_contour_segments
  3. from collections import deque
  4. _param_options = ('high', 'low')
  5. def find_contours(
  6. image, level=None, fully_connected='low', positive_orientation='low', *, mask=None
  7. ):
  8. """Find iso-valued contours in a 2D array for a given level value.
  9. Uses the "marching squares" method to compute the iso-valued contours of
  10. the input 2D array for a particular level value. Array values are linearly
  11. interpolated to provide better precision for the output contours.
  12. Parameters
  13. ----------
  14. image : (M, N) ndarray of double
  15. Input image in which to find contours.
  16. level : float, optional
  17. Value along which to find contours in the array. By default, the level
  18. is set to (max(image) + min(image)) / 2
  19. .. versionchanged:: 0.18
  20. This parameter is now optional.
  21. fully_connected : str, {'low', 'high'}
  22. Indicates whether array elements below the given level value are to be
  23. considered fully-connected (and hence elements above the value will
  24. only be face connected), or vice-versa. (See notes below for details.)
  25. positive_orientation : str, {'low', 'high'}
  26. Indicates whether the output contours will produce positively-oriented
  27. polygons around islands of low- or high-valued elements. If 'low' then
  28. contours will wind counter-clockwise around elements below the
  29. iso-value. Alternately, this means that low-valued elements are always
  30. on the left of the contour. (See below for details.)
  31. mask : (M, N) ndarray of bool or None
  32. A boolean mask, True where we want to draw contours.
  33. Note that NaN values are always excluded from the considered region
  34. (``mask`` is set to ``False`` wherever ``array`` is ``NaN``).
  35. Returns
  36. -------
  37. contours : list of (K, 2) ndarrays
  38. Each contour is a ndarray of ``(row, column)`` coordinates along the contour.
  39. See Also
  40. --------
  41. skimage.measure.marching_cubes
  42. Notes
  43. -----
  44. The marching squares algorithm is a special case of the marching cubes
  45. algorithm [1]_. A simple explanation is available here:
  46. https://users.polytech.unice.fr/~lingrand/MarchingCubes/algo.html
  47. There is a single ambiguous case in the marching squares algorithm: when
  48. a given ``2 x 2``-element square has two high-valued and two low-valued
  49. elements, each pair diagonally adjacent. (Where high- and low-valued is
  50. with respect to the contour value sought.) In this case, either the
  51. high-valued elements can be 'connected together' via a thin isthmus that
  52. separates the low-valued elements, or vice-versa. When elements are
  53. connected together across a diagonal, they are considered 'fully
  54. connected' (also known as 'face+vertex-connected' or '8-connected'). Only
  55. high-valued or low-valued elements can be fully-connected, the other set
  56. will be considered as 'face-connected' or '4-connected'. By default,
  57. low-valued elements are considered fully-connected; this can be altered
  58. with the 'fully_connected' parameter.
  59. Output contours are not guaranteed to be closed: contours which intersect
  60. the array edge or a masked-off region (either where mask is False or where
  61. array is NaN) will be left open. All other contours will be closed. (The
  62. closed-ness of a contours can be tested by checking whether the beginning
  63. point is the same as the end point.)
  64. Contours are oriented. By default, array values lower than the contour
  65. value are to the left of the contour and values greater than the contour
  66. value are to the right. This means that contours will wind
  67. counter-clockwise (i.e. in 'positive orientation') around islands of
  68. low-valued pixels. This behavior can be altered with the
  69. 'positive_orientation' parameter.
  70. The order of the contours in the output list is determined by the position
  71. of the smallest ``x,y`` (in lexicographical order) coordinate in the
  72. contour. This is a side effect of how the input array is traversed, but
  73. can be relied upon.
  74. .. warning::
  75. Array coordinates/values are assumed to refer to the *center* of the
  76. array element. Take a simple example input: ``[0, 1]``. The interpolated
  77. position of 0.5 in this array is midway between the 0-element (at
  78. ``x=0``) and the 1-element (at ``x=1``), and thus would fall at
  79. ``x=0.5``.
  80. This means that to find reasonable contours, it is best to find contours
  81. midway between the expected "light" and "dark" values. In particular,
  82. given a binarized array, *do not* choose to find contours at the low or
  83. high value of the array. This will often yield degenerate contours,
  84. especially around structures that are a single array element wide. Instead,
  85. choose a middle value, as above.
  86. References
  87. ----------
  88. .. [1] Lorensen, William and Harvey E. Cline. Marching Cubes: A High
  89. Resolution 3D Surface Construction Algorithm. Computer Graphics
  90. (SIGGRAPH 87 Proceedings) 21(4) July 1987, p. 163-170).
  91. :DOI:`10.1145/37401.37422`
  92. Examples
  93. --------
  94. >>> a = np.zeros((3, 3))
  95. >>> a[0, 0] = 1
  96. >>> a
  97. array([[1., 0., 0.],
  98. [0., 0., 0.],
  99. [0., 0., 0.]])
  100. >>> find_contours(a, 0.5)
  101. [array([[0. , 0.5],
  102. [0.5, 0. ]])]
  103. """
  104. if fully_connected not in _param_options:
  105. raise ValueError(
  106. 'Parameters "fully_connected" must be either ' '"high" or "low".'
  107. )
  108. if positive_orientation not in _param_options:
  109. raise ValueError(
  110. 'Parameters "positive_orientation" must be either ' '"high" or "low".'
  111. )
  112. if image.shape[0] < 2 or image.shape[1] < 2:
  113. raise ValueError("Input array must be at least 2x2.")
  114. if image.ndim != 2:
  115. raise ValueError('Only 2D arrays are supported.')
  116. if mask is not None:
  117. if mask.shape != image.shape:
  118. raise ValueError('Parameters "array" and "mask"' ' must have same shape.')
  119. if not np.can_cast(mask.dtype, bool, casting='safe'):
  120. raise TypeError('Parameter "mask" must be a binary array.')
  121. mask = mask.astype(np.uint8, copy=False)
  122. if level is None:
  123. level = (np.nanmin(image) + np.nanmax(image)) / 2.0
  124. segments = _get_contour_segments(
  125. image.astype(np.float64), float(level), fully_connected == 'high', mask=mask
  126. )
  127. contours = _assemble_contours(segments)
  128. if positive_orientation == 'high':
  129. contours = [c[::-1] for c in contours]
  130. return contours
  131. def _assemble_contours(segments):
  132. current_index = 0
  133. contours = {}
  134. starts = {}
  135. ends = {}
  136. for from_point, to_point in segments:
  137. # Ignore degenerate segments.
  138. # This happens when (and only when) one vertex of the square is
  139. # exactly the contour level, and the rest are above or below.
  140. # This degenerate vertex will be picked up later by neighboring
  141. # squares.
  142. if from_point == to_point:
  143. continue
  144. tail, tail_num = starts.pop(to_point, (None, None))
  145. head, head_num = ends.pop(from_point, (None, None))
  146. if tail is not None and head is not None:
  147. # We need to connect these two contours.
  148. if tail is head:
  149. # We need to closed a contour: add the end point
  150. head.append(to_point)
  151. else: # tail is not head
  152. # We need to join two distinct contours.
  153. # We want to keep the first contour segment created, so that
  154. # the final contours are ordered left->right, top->bottom.
  155. if tail_num > head_num:
  156. # tail was created second. Append tail to head.
  157. head.extend(tail)
  158. # Remove tail from the detected contours
  159. contours.pop(tail_num, None)
  160. # Update starts and ends
  161. starts[head[0]] = (head, head_num)
  162. ends[head[-1]] = (head, head_num)
  163. else: # tail_num <= head_num
  164. # head was created second. Prepend head to tail.
  165. tail.extendleft(reversed(head))
  166. # Remove head from the detected contours
  167. starts.pop(head[0], None) # head[0] can be == to_point!
  168. contours.pop(head_num, None)
  169. # Update starts and ends
  170. starts[tail[0]] = (tail, tail_num)
  171. ends[tail[-1]] = (tail, tail_num)
  172. elif tail is None and head is None:
  173. # We need to add a new contour
  174. new_contour = deque((from_point, to_point))
  175. contours[current_index] = new_contour
  176. starts[from_point] = (new_contour, current_index)
  177. ends[to_point] = (new_contour, current_index)
  178. current_index += 1
  179. elif head is None: # tail is not None
  180. # tail first element is to_point: the new segment should be
  181. # prepended.
  182. tail.appendleft(from_point)
  183. # Update starts
  184. starts[from_point] = (tail, tail_num)
  185. else: # tail is None and head is not None:
  186. # head last element is from_point: the new segment should be
  187. # appended
  188. head.append(to_point)
  189. # Update ends
  190. ends[to_point] = (head, head_num)
  191. return [np.array(contour) for _, contour in sorted(contours.items())]