Shortcuts on this page
r m x toggle line displays
j k next/prev highlighted chunk
0 (zero) top of page
1 (one) first highlighted chunk
1"""
2Tools for triangular grids.
3"""
5import numpy as np
7from matplotlib import _api
8from matplotlib.tri import Triangulation
11class TriAnalyzer:
12 """
13 Define basic tools for triangular mesh analysis and improvement.
15 A TriAnalyzer encapsulates a `.Triangulation` object and provides basic
16 tools for mesh analysis and mesh improvement.
18 Attributes
19 ----------
20 scale_factors
22 Parameters
23 ----------
24 triangulation : `~matplotlib.tri.Triangulation`
25 The encapsulated triangulation to analyze.
26 """
28 def __init__(self, triangulation):
29 _api.check_isinstance(Triangulation, triangulation=triangulation)
30 self._triangulation = triangulation
32 @property
33 def scale_factors(self):
34 """
35 Factors to rescale the triangulation into a unit square.
37 Returns
38 -------
39 (float, float)
40 Scaling factors (kx, ky) so that the triangulation
41 ``[triangulation.x * kx, triangulation.y * ky]``
42 fits exactly inside a unit square.
43 """
44 compressed_triangles = self._triangulation.get_masked_triangles()
45 node_used = (np.bincount(np.ravel(compressed_triangles),
46 minlength=self._triangulation.x.size) != 0)
47 return (1 / np.ptp(self._triangulation.x[node_used]),
48 1 / np.ptp(self._triangulation.y[node_used]))
50 def circle_ratios(self, rescale=True):
51 """
52 Return a measure of the triangulation triangles flatness.
54 The ratio of the incircle radius over the circumcircle radius is a
55 widely used indicator of a triangle flatness.
56 It is always ``<= 0.5`` and ``== 0.5`` only for equilateral
57 triangles. Circle ratios below 0.01 denote very flat triangles.
59 To avoid unduly low values due to a difference of scale between the 2
60 axis, the triangular mesh can first be rescaled to fit inside a unit
61 square with `scale_factors` (Only if *rescale* is True, which is
62 its default value).
64 Parameters
65 ----------
66 rescale : bool, default: True
67 If True, internally rescale (based on `scale_factors`), so that the
68 (unmasked) triangles fit exactly inside a unit square mesh.
70 Returns
71 -------
72 masked array
73 Ratio of the incircle radius over the circumcircle radius, for
74 each 'rescaled' triangle of the encapsulated triangulation.
75 Values corresponding to masked triangles are masked out.
77 """
78 # Coords rescaling
79 if rescale:
80 (kx, ky) = self.scale_factors
81 else:
82 (kx, ky) = (1.0, 1.0)
83 pts = np.vstack([self._triangulation.x*kx,
84 self._triangulation.y*ky]).T
85 tri_pts = pts[self._triangulation.triangles]
86 # Computes the 3 side lengths
87 a = tri_pts[:, 1, :] - tri_pts[:, 0, :]
88 b = tri_pts[:, 2, :] - tri_pts[:, 1, :]
89 c = tri_pts[:, 0, :] - tri_pts[:, 2, :]
90 a = np.hypot(a[:, 0], a[:, 1])
91 b = np.hypot(b[:, 0], b[:, 1])
92 c = np.hypot(c[:, 0], c[:, 1])
93 # circumcircle and incircle radii
94 s = (a+b+c)*0.5
95 prod = s*(a+b-s)*(a+c-s)*(b+c-s)
96 # We have to deal with flat triangles with infinite circum_radius
97 bool_flat = (prod == 0.)
98 if np.any(bool_flat):
99 # Pathologic flow
100 ntri = tri_pts.shape[0]
101 circum_radius = np.empty(ntri, dtype=np.float64)
102 circum_radius[bool_flat] = np.inf
103 abc = a*b*c
104 circum_radius[~bool_flat] = abc[~bool_flat] / (
105 4.0*np.sqrt(prod[~bool_flat]))
106 else:
107 # Normal optimized flow
108 circum_radius = (a*b*c) / (4.0*np.sqrt(prod))
109 in_radius = (a*b*c) / (4.0*circum_radius*s)
110 circle_ratio = in_radius/circum_radius
111 mask = self._triangulation.mask
112 if mask is None:
113 return circle_ratio
114 else:
115 return np.ma.array(circle_ratio, mask=mask)
117 def get_flat_tri_mask(self, min_circle_ratio=0.01, rescale=True):
118 """
119 Eliminate excessively flat border triangles from the triangulation.
121 Returns a mask *new_mask* which allows to clean the encapsulated
122 triangulation from its border-located flat triangles
123 (according to their :meth:`circle_ratios`).
124 This mask is meant to be subsequently applied to the triangulation
125 using `.Triangulation.set_mask`.
126 *new_mask* is an extension of the initial triangulation mask
127 in the sense that an initially masked triangle will remain masked.
129 The *new_mask* array is computed recursively; at each step flat
130 triangles are removed only if they share a side with the current mesh
131 border. Thus, no new holes in the triangulated domain will be created.
133 Parameters
134 ----------
135 min_circle_ratio : float, default: 0.01
136 Border triangles with incircle/circumcircle radii ratio r/R will
137 be removed if r/R < *min_circle_ratio*.
138 rescale : bool, default: True
139 If True, first, internally rescale (based on `scale_factors`) so
140 that the (unmasked) triangles fit exactly inside a unit square
141 mesh. This rescaling accounts for the difference of scale which
142 might exist between the 2 axis.
144 Returns
145 -------
146 array of bool
147 Mask to apply to encapsulated triangulation.
148 All the initially masked triangles remain masked in the
149 *new_mask*.
151 Notes
152 -----
153 The rationale behind this function is that a Delaunay
154 triangulation - of an unstructured set of points - sometimes contains
155 almost flat triangles at its border, leading to artifacts in plots
156 (especially for high-resolution contouring).
157 Masked with computed *new_mask*, the encapsulated
158 triangulation would contain no more unmasked border triangles
159 with a circle ratio below *min_circle_ratio*, thus improving the
160 mesh quality for subsequent plots or interpolation.
161 """
162 # Recursively computes the mask_current_borders, true if a triangle is
163 # at the border of the mesh OR touching the border through a chain of
164 # invalid aspect ratio masked_triangles.
165 ntri = self._triangulation.triangles.shape[0]
166 mask_bad_ratio = self.circle_ratios(rescale) < min_circle_ratio
168 current_mask = self._triangulation.mask
169 if current_mask is None:
170 current_mask = np.zeros(ntri, dtype=bool)
171 valid_neighbors = np.copy(self._triangulation.neighbors)
172 renum_neighbors = np.arange(ntri, dtype=np.int32)
173 nadd = -1
174 while nadd != 0:
175 # The active wavefront is the triangles from the border (unmasked
176 # but with a least 1 neighbor equal to -1
177 wavefront = (np.min(valid_neighbors, axis=1) == -1) & ~current_mask
178 # The element from the active wavefront will be masked if their
179 # circle ratio is bad.
180 added_mask = wavefront & mask_bad_ratio
181 current_mask = added_mask | current_mask
182 nadd = np.sum(added_mask)
184 # now we have to update the tables valid_neighbors
185 valid_neighbors[added_mask, :] = -1
186 renum_neighbors[added_mask] = -1
187 valid_neighbors = np.where(valid_neighbors == -1, -1,
188 renum_neighbors[valid_neighbors])
190 return np.ma.filled(current_mask, True)
192 def _get_compressed_triangulation(self):
193 """
194 Compress (if masked) the encapsulated triangulation.
196 Returns minimal-length triangles array (*compressed_triangles*) and
197 coordinates arrays (*compressed_x*, *compressed_y*) that can still
198 describe the unmasked triangles of the encapsulated triangulation.
200 Returns
201 -------
202 compressed_triangles : array-like
203 the returned compressed triangulation triangles
204 compressed_x : array-like
205 the returned compressed triangulation 1st coordinate
206 compressed_y : array-like
207 the returned compressed triangulation 2nd coordinate
208 tri_renum : int array
209 renumbering table to translate the triangle numbers from the
210 encapsulated triangulation into the new (compressed) renumbering.
211 -1 for masked triangles (deleted from *compressed_triangles*).
212 node_renum : int array
213 renumbering table to translate the point numbers from the
214 encapsulated triangulation into the new (compressed) renumbering.
215 -1 for unused points (i.e. those deleted from *compressed_x* and
216 *compressed_y*).
218 """
219 # Valid triangles and renumbering
220 tri_mask = self._triangulation.mask
221 compressed_triangles = self._triangulation.get_masked_triangles()
222 ntri = self._triangulation.triangles.shape[0]
223 if tri_mask is not None:
224 tri_renum = self._total_to_compress_renum(~tri_mask)
225 else:
226 tri_renum = np.arange(ntri, dtype=np.int32)
228 # Valid nodes and renumbering
229 valid_node = (np.bincount(np.ravel(compressed_triangles),
230 minlength=self._triangulation.x.size) != 0)
231 compressed_x = self._triangulation.x[valid_node]
232 compressed_y = self._triangulation.y[valid_node]
233 node_renum = self._total_to_compress_renum(valid_node)
235 # Now renumbering the valid triangles nodes
236 compressed_triangles = node_renum[compressed_triangles]
238 return (compressed_triangles, compressed_x, compressed_y, tri_renum,
239 node_renum)
241 @staticmethod
242 def _total_to_compress_renum(valid):
243 """
244 Parameters
245 ----------
246 valid : 1D bool array
247 Validity mask.
249 Returns
250 -------
251 int array
252 Array so that (`valid_array` being a compressed array
253 based on a `masked_array` with mask ~*valid*):
255 - For all i with valid[i] = True:
256 valid_array[renum[i]] = masked_array[i]
257 - For all i with valid[i] = False:
258 renum[i] = -1 (invalid value)
259 """
260 renum = np.full(np.size(valid), -1, dtype=np.int32)
261 n_valid = np.sum(valid)
262 renum[valid] = np.arange(n_valid, dtype=np.int32)
263 return renum