_triangulation.py 9.6 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247
  1. import sys
  2. import numpy as np
  3. from matplotlib import _api
  4. class Triangulation:
  5. """
  6. An unstructured triangular grid consisting of npoints points and
  7. ntri triangles. The triangles can either be specified by the user
  8. or automatically generated using a Delaunay triangulation.
  9. Parameters
  10. ----------
  11. x, y : (npoints,) array-like
  12. Coordinates of grid points.
  13. triangles : (ntri, 3) array-like of int, optional
  14. For each triangle, the indices of the three points that make
  15. up the triangle, ordered in an anticlockwise manner. If not
  16. specified, the Delaunay triangulation is calculated.
  17. mask : (ntri,) array-like of bool, optional
  18. Which triangles are masked out.
  19. Attributes
  20. ----------
  21. triangles : (ntri, 3) array of int
  22. For each triangle, the indices of the three points that make
  23. up the triangle, ordered in an anticlockwise manner. If you want to
  24. take the *mask* into account, use `get_masked_triangles` instead.
  25. mask : (ntri, 3) array of bool or None
  26. Masked out triangles.
  27. is_delaunay : bool
  28. Whether the Triangulation is a calculated Delaunay
  29. triangulation (where *triangles* was not specified) or not.
  30. Notes
  31. -----
  32. For a Triangulation to be valid it must not have duplicate points,
  33. triangles formed from colinear points, or overlapping triangles.
  34. """
  35. def __init__(self, x, y, triangles=None, mask=None):
  36. from matplotlib import _qhull
  37. self.x = np.asarray(x, dtype=np.float64)
  38. self.y = np.asarray(y, dtype=np.float64)
  39. if self.x.shape != self.y.shape or self.x.ndim != 1:
  40. raise ValueError("x and y must be equal-length 1D arrays, but "
  41. f"found shapes {self.x.shape!r} and "
  42. f"{self.y.shape!r}")
  43. self.mask = None
  44. self._edges = None
  45. self._neighbors = None
  46. self.is_delaunay = False
  47. if triangles is None:
  48. # No triangulation specified, so use matplotlib._qhull to obtain
  49. # Delaunay triangulation.
  50. self.triangles, self._neighbors = _qhull.delaunay(x, y, sys.flags.verbose)
  51. self.is_delaunay = True
  52. else:
  53. # Triangulation specified. Copy, since we may correct triangle
  54. # orientation.
  55. try:
  56. self.triangles = np.array(triangles, dtype=np.int32, order='C')
  57. except ValueError as e:
  58. raise ValueError('triangles must be a (N, 3) int array, not '
  59. f'{triangles!r}') from e
  60. if self.triangles.ndim != 2 or self.triangles.shape[1] != 3:
  61. raise ValueError(
  62. 'triangles must be a (N, 3) int array, but found shape '
  63. f'{self.triangles.shape!r}')
  64. if self.triangles.max() >= len(self.x):
  65. raise ValueError(
  66. 'triangles are indices into the points and must be in the '
  67. f'range 0 <= i < {len(self.x)} but found value '
  68. f'{self.triangles.max()}')
  69. if self.triangles.min() < 0:
  70. raise ValueError(
  71. 'triangles are indices into the points and must be in the '
  72. f'range 0 <= i < {len(self.x)} but found value '
  73. f'{self.triangles.min()}')
  74. # Underlying C++ object is not created until first needed.
  75. self._cpp_triangulation = None
  76. # Default TriFinder not created until needed.
  77. self._trifinder = None
  78. self.set_mask(mask)
  79. def calculate_plane_coefficients(self, z):
  80. """
  81. Calculate plane equation coefficients for all unmasked triangles from
  82. the point (x, y) coordinates and specified z-array of shape (npoints).
  83. The returned array has shape (npoints, 3) and allows z-value at (x, y)
  84. position in triangle tri to be calculated using
  85. ``z = array[tri, 0] * x + array[tri, 1] * y + array[tri, 2]``.
  86. """
  87. return self.get_cpp_triangulation().calculate_plane_coefficients(z)
  88. @property
  89. def edges(self):
  90. """
  91. Return integer array of shape (nedges, 2) containing all edges of
  92. non-masked triangles.
  93. Each row defines an edge by its start point index and end point
  94. index. Each edge appears only once, i.e. for an edge between points
  95. *i* and *j*, there will only be either *(i, j)* or *(j, i)*.
  96. """
  97. if self._edges is None:
  98. self._edges = self.get_cpp_triangulation().get_edges()
  99. return self._edges
  100. def get_cpp_triangulation(self):
  101. """
  102. Return the underlying C++ Triangulation object, creating it
  103. if necessary.
  104. """
  105. from matplotlib import _tri
  106. if self._cpp_triangulation is None:
  107. self._cpp_triangulation = _tri.Triangulation(
  108. # For unset arrays use empty tuple which has size of zero.
  109. self.x, self.y, self.triangles,
  110. self.mask if self.mask is not None else (),
  111. self._edges if self._edges is not None else (),
  112. self._neighbors if self._neighbors is not None else (),
  113. not self.is_delaunay)
  114. return self._cpp_triangulation
  115. def get_masked_triangles(self):
  116. """
  117. Return an array of triangles taking the mask into account.
  118. """
  119. if self.mask is not None:
  120. return self.triangles[~self.mask]
  121. else:
  122. return self.triangles
  123. @staticmethod
  124. def get_from_args_and_kwargs(*args, **kwargs):
  125. """
  126. Return a Triangulation object from the args and kwargs, and
  127. the remaining args and kwargs with the consumed values removed.
  128. There are two alternatives: either the first argument is a
  129. Triangulation object, in which case it is returned, or the args
  130. and kwargs are sufficient to create a new Triangulation to
  131. return. In the latter case, see Triangulation.__init__ for
  132. the possible args and kwargs.
  133. """
  134. if isinstance(args[0], Triangulation):
  135. triangulation, *args = args
  136. if 'triangles' in kwargs:
  137. _api.warn_external(
  138. "Passing the keyword 'triangles' has no effect when also "
  139. "passing a Triangulation")
  140. if 'mask' in kwargs:
  141. _api.warn_external(
  142. "Passing the keyword 'mask' has no effect when also "
  143. "passing a Triangulation")
  144. else:
  145. x, y, triangles, mask, args, kwargs = \
  146. Triangulation._extract_triangulation_params(args, kwargs)
  147. triangulation = Triangulation(x, y, triangles, mask)
  148. return triangulation, args, kwargs
  149. @staticmethod
  150. def _extract_triangulation_params(args, kwargs):
  151. x, y, *args = args
  152. # Check triangles in kwargs then args.
  153. triangles = kwargs.pop('triangles', None)
  154. from_args = False
  155. if triangles is None and args:
  156. triangles = args[0]
  157. from_args = True
  158. if triangles is not None:
  159. try:
  160. triangles = np.asarray(triangles, dtype=np.int32)
  161. except ValueError:
  162. triangles = None
  163. if triangles is not None and (triangles.ndim != 2 or
  164. triangles.shape[1] != 3):
  165. triangles = None
  166. if triangles is not None and from_args:
  167. args = args[1:] # Consumed first item in args.
  168. # Check for mask in kwargs.
  169. mask = kwargs.pop('mask', None)
  170. return x, y, triangles, mask, args, kwargs
  171. def get_trifinder(self):
  172. """
  173. Return the default `matplotlib.tri.TriFinder` of this
  174. triangulation, creating it if necessary. This allows the same
  175. TriFinder object to be easily shared.
  176. """
  177. if self._trifinder is None:
  178. # Default TriFinder class.
  179. from matplotlib.tri._trifinder import TrapezoidMapTriFinder
  180. self._trifinder = TrapezoidMapTriFinder(self)
  181. return self._trifinder
  182. @property
  183. def neighbors(self):
  184. """
  185. Return integer array of shape (ntri, 3) containing neighbor triangles.
  186. For each triangle, the indices of the three triangles that
  187. share the same edges, or -1 if there is no such neighboring
  188. triangle. ``neighbors[i, j]`` is the triangle that is the neighbor
  189. to the edge from point index ``triangles[i, j]`` to point index
  190. ``triangles[i, (j+1)%3]``.
  191. """
  192. if self._neighbors is None:
  193. self._neighbors = self.get_cpp_triangulation().get_neighbors()
  194. return self._neighbors
  195. def set_mask(self, mask):
  196. """
  197. Set or clear the mask array.
  198. Parameters
  199. ----------
  200. mask : None or bool array of length ntri
  201. """
  202. if mask is None:
  203. self.mask = None
  204. else:
  205. self.mask = np.asarray(mask, dtype=bool)
  206. if self.mask.shape != (self.triangles.shape[0],):
  207. raise ValueError('mask array must have same length as '
  208. 'triangles array')
  209. # Set mask in C++ Triangulation.
  210. if self._cpp_triangulation is not None:
  211. self._cpp_triangulation.set_mask(
  212. self.mask if self.mask is not None else ())
  213. # Clear derived fields so they are recalculated when needed.
  214. self._edges = None
  215. self._neighbors = None
  216. # Recalculate TriFinder if it exists.
  217. if self._trifinder is not None:
  218. self._trifinder._initialize()