triangulation.py 8.1 KB

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