300 lines
12 KiB
Python
300 lines
12 KiB
Python
|
"""
|
||
|
|
Tools for triangular grids.
|
||
|
|
"""
|
||
|
|
|
||
|
|
import numpy as np
|
||
|
|
|
||
|
|
from matplotlib.tri import Triangulation
|
||
|
|
|
||
|
|
|
||
|
|
class TriAnalyzer(object):
|
||
|
|
"""
|
||
|
|
Define basic tools for triangular mesh analysis and improvement.
|
||
|
|
|
||
|
|
A TriAnalyzer encapsulates a :class:`~matplotlib.tri.Triangulation`
|
||
|
|
object and provides basic tools for mesh analysis and mesh improvement.
|
||
|
|
|
||
|
|
Parameters
|
||
|
|
----------
|
||
|
|
triangulation : :class:`~matplotlib.tri.Triangulation` object
|
||
|
|
The encapsulated triangulation to analyze.
|
||
|
|
|
||
|
|
Attributes
|
||
|
|
----------
|
||
|
|
`scale_factors`
|
||
|
|
|
||
|
|
"""
|
||
|
|
def __init__(self, triangulation):
|
||
|
|
if not isinstance(triangulation, Triangulation):
|
||
|
|
raise ValueError("Expected a Triangulation object")
|
||
|
|
self._triangulation = triangulation
|
||
|
|
|
||
|
|
@property
|
||
|
|
def scale_factors(self):
|
||
|
|
"""
|
||
|
|
Factors to rescale the triangulation into a unit square.
|
||
|
|
|
||
|
|
Returns *k*, tuple of 2 scale factors.
|
||
|
|
|
||
|
|
Returns
|
||
|
|
-------
|
||
|
|
k : tuple of 2 floats (kx, ky)
|
||
|
|
Tuple of floats that would rescale the triangulation :
|
||
|
|
``[triangulation.x * kx, triangulation.y * ky]``
|
||
|
|
fits exactly inside a unit square.
|
||
|
|
|
||
|
|
"""
|
||
|
|
compressed_triangles = self._triangulation.get_masked_triangles()
|
||
|
|
node_used = (np.bincount(np.ravel(compressed_triangles),
|
||
|
|
minlength=self._triangulation.x.size) != 0)
|
||
|
|
return (1 / np.ptp(self._triangulation.x[node_used]),
|
||
|
|
1 / np.ptp(self._triangulation.y[node_used]))
|
||
|
|
|
||
|
|
def circle_ratios(self, rescale=True):
|
||
|
|
"""
|
||
|
|
Returns a measure of the triangulation triangles flatness.
|
||
|
|
|
||
|
|
The ratio of the incircle radius over the circumcircle radius is a
|
||
|
|
widely used indicator of a triangle flatness.
|
||
|
|
It is always ``<= 0.5`` and ``== 0.5`` only for equilateral
|
||
|
|
triangles. Circle ratios below 0.01 denote very flat triangles.
|
||
|
|
|
||
|
|
To avoid unduly low values due to a difference of scale between the 2
|
||
|
|
axis, the triangular mesh can first be rescaled to fit inside a unit
|
||
|
|
square with :attr:`scale_factors` (Only if *rescale* is True, which is
|
||
|
|
its default value).
|
||
|
|
|
||
|
|
Parameters
|
||
|
|
----------
|
||
|
|
rescale : boolean, optional
|
||
|
|
If True, a rescaling will be internally performed (based on
|
||
|
|
:attr:`scale_factors`, so that the (unmasked) triangles fit
|
||
|
|
exactly inside a unit square mesh. Default is True.
|
||
|
|
|
||
|
|
Returns
|
||
|
|
-------
|
||
|
|
circle_ratios : masked array
|
||
|
|
Ratio of the incircle radius over the
|
||
|
|
circumcircle radius, for each 'rescaled' triangle of the
|
||
|
|
encapsulated triangulation.
|
||
|
|
Values corresponding to masked triangles are masked out.
|
||
|
|
|
||
|
|
"""
|
||
|
|
# Coords rescaling
|
||
|
|
if rescale:
|
||
|
|
(kx, ky) = self.scale_factors
|
||
|
|
else:
|
||
|
|
(kx, ky) = (1.0, 1.0)
|
||
|
|
pts = np.vstack([self._triangulation.x*kx,
|
||
|
|
self._triangulation.y*ky]).T
|
||
|
|
tri_pts = pts[self._triangulation.triangles]
|
||
|
|
# Computes the 3 side lengths
|
||
|
|
a = tri_pts[:, 1, :] - tri_pts[:, 0, :]
|
||
|
|
b = tri_pts[:, 2, :] - tri_pts[:, 1, :]
|
||
|
|
c = tri_pts[:, 0, :] - tri_pts[:, 2, :]
|
||
|
|
a = np.hypot(a[:, 0], a[:, 1])
|
||
|
|
b = np.hypot(b[:, 0], b[:, 1])
|
||
|
|
c = np.hypot(c[:, 0], c[:, 1])
|
||
|
|
# circumcircle and incircle radii
|
||
|
|
s = (a+b+c)*0.5
|
||
|
|
prod = s*(a+b-s)*(a+c-s)*(b+c-s)
|
||
|
|
# We have to deal with flat triangles with infinite circum_radius
|
||
|
|
bool_flat = (prod == 0.)
|
||
|
|
if np.any(bool_flat):
|
||
|
|
# Pathologic flow
|
||
|
|
ntri = tri_pts.shape[0]
|
||
|
|
circum_radius = np.empty(ntri, dtype=np.float64)
|
||
|
|
circum_radius[bool_flat] = np.inf
|
||
|
|
abc = a*b*c
|
||
|
|
circum_radius[~bool_flat] = abc[~bool_flat] / (
|
||
|
|
4.0*np.sqrt(prod[~bool_flat]))
|
||
|
|
else:
|
||
|
|
# Normal optimized flow
|
||
|
|
circum_radius = (a*b*c) / (4.0*np.sqrt(prod))
|
||
|
|
in_radius = (a*b*c) / (4.0*circum_radius*s)
|
||
|
|
circle_ratio = in_radius/circum_radius
|
||
|
|
mask = self._triangulation.mask
|
||
|
|
if mask is None:
|
||
|
|
return circle_ratio
|
||
|
|
else:
|
||
|
|
return np.ma.array(circle_ratio, mask=mask)
|
||
|
|
|
||
|
|
def get_flat_tri_mask(self, min_circle_ratio=0.01, rescale=True):
|
||
|
|
"""
|
||
|
|
Eliminates excessively flat border triangles from the triangulation.
|
||
|
|
|
||
|
|
Returns a mask *new_mask* which allows to clean the encapsulated
|
||
|
|
triangulation from its border-located flat triangles
|
||
|
|
(according to their :meth:`circle_ratios`).
|
||
|
|
This mask is meant to be subsequently applied to the triangulation
|
||
|
|
using :func:`matplotlib.tri.Triangulation.set_mask`.
|
||
|
|
*new_mask* is an extension of the initial triangulation mask
|
||
|
|
in the sense that an initially masked triangle will remain masked.
|
||
|
|
|
||
|
|
The *new_mask* array is computed recursively; at each step flat
|
||
|
|
triangles are removed only if they share a side with the current mesh
|
||
|
|
border. Thus no new holes in the triangulated domain will be created.
|
||
|
|
|
||
|
|
Parameters
|
||
|
|
----------
|
||
|
|
min_circle_ratio : float, optional
|
||
|
|
Border triangles with incircle/circumcircle radii ratio r/R will
|
||
|
|
be removed if r/R < *min_circle_ratio*. Default value: 0.01
|
||
|
|
rescale : boolean, optional
|
||
|
|
If True, a rescaling will first be internally performed (based on
|
||
|
|
:attr:`scale_factors` ), so that the (unmasked) triangles fit
|
||
|
|
exactly inside a unit square mesh. This rescaling accounts for the
|
||
|
|
difference of scale which might exist between the 2 axis. Default
|
||
|
|
(and recommended) value is True.
|
||
|
|
|
||
|
|
Returns
|
||
|
|
-------
|
||
|
|
new_mask : array-like of booleans
|
||
|
|
Mask to apply to encapsulated triangulation.
|
||
|
|
All the initially masked triangles remain masked in the
|
||
|
|
*new_mask*.
|
||
|
|
|
||
|
|
Notes
|
||
|
|
-----
|
||
|
|
The rationale behind this function is that a Delaunay
|
||
|
|
triangulation - of an unstructured set of points - sometimes contains
|
||
|
|
almost flat triangles at its border, leading to artifacts in plots
|
||
|
|
(especially for high-resolution contouring).
|
||
|
|
Masked with computed *new_mask*, the encapsulated
|
||
|
|
triangulation would contain no more unmasked border triangles
|
||
|
|
with a circle ratio below *min_circle_ratio*, thus improving the
|
||
|
|
mesh quality for subsequent plots or interpolation.
|
||
|
|
"""
|
||
|
|
# Recursively computes the mask_current_borders, true if a triangle is
|
||
|
|
# at the border of the mesh OR touching the border through a chain of
|
||
|
|
# invalid aspect ratio masked_triangles.
|
||
|
|
ntri = self._triangulation.triangles.shape[0]
|
||
|
|
mask_bad_ratio = self.circle_ratios(rescale) < min_circle_ratio
|
||
|
|
|
||
|
|
current_mask = self._triangulation.mask
|
||
|
|
if current_mask is None:
|
||
|
|
current_mask = np.zeros(ntri, dtype=bool)
|
||
|
|
valid_neighbors = np.copy(self._triangulation.neighbors)
|
||
|
|
renum_neighbors = np.arange(ntri, dtype=np.int32)
|
||
|
|
nadd = -1
|
||
|
|
while nadd != 0:
|
||
|
|
# The active wavefront is the triangles from the border (unmasked
|
||
|
|
# but with a least 1 neighbor equal to -1
|
||
|
|
wavefront = (np.min(valid_neighbors, axis=1) == -1) & ~current_mask
|
||
|
|
# The element from the active wavefront will be masked if their
|
||
|
|
# circle ratio is bad.
|
||
|
|
added_mask = wavefront & mask_bad_ratio
|
||
|
|
current_mask = added_mask | current_mask
|
||
|
|
nadd = np.sum(added_mask)
|
||
|
|
|
||
|
|
# now we have to update the tables valid_neighbors
|
||
|
|
valid_neighbors[added_mask, :] = -1
|
||
|
|
renum_neighbors[added_mask] = -1
|
||
|
|
valid_neighbors = np.where(valid_neighbors == -1, -1,
|
||
|
|
renum_neighbors[valid_neighbors])
|
||
|
|
|
||
|
|
return np.ma.filled(current_mask, True)
|
||
|
|
|
||
|
|
def _get_compressed_triangulation(self, return_tri_renum=False,
|
||
|
|
return_node_renum=False):
|
||
|
|
"""
|
||
|
|
Compress (if masked) the encapsulated triangulation.
|
||
|
|
|
||
|
|
Returns minimal-length triangles array (*compressed_triangles*) and
|
||
|
|
coordinates arrays (*compressed_x*, *compressed_y*) that can still
|
||
|
|
describe the unmasked triangles of the encapsulated triangulation.
|
||
|
|
|
||
|
|
Parameters
|
||
|
|
----------
|
||
|
|
return_tri_renum : boolean, optional
|
||
|
|
Indicates whether a renumbering table to translate the triangle
|
||
|
|
numbers from the encapsulated triangulation numbering into the
|
||
|
|
new (compressed) renumbering will be returned.
|
||
|
|
return_node_renum : boolean, optional
|
||
|
|
Indicates whether a renumbering table to translate the nodes
|
||
|
|
numbers from the encapsulated triangulation numbering into the
|
||
|
|
new (compressed) renumbering will be returned.
|
||
|
|
|
||
|
|
Returns
|
||
|
|
-------
|
||
|
|
compressed_triangles : array-like
|
||
|
|
the returned compressed triangulation triangles
|
||
|
|
compressed_x : array-like
|
||
|
|
the returned compressed triangulation 1st coordinate
|
||
|
|
compressed_y : array-like
|
||
|
|
the returned compressed triangulation 2nd coordinate
|
||
|
|
tri_renum : array-like of integers
|
||
|
|
renumbering table to translate the triangle numbers from the
|
||
|
|
encapsulated triangulation into the new (compressed) renumbering.
|
||
|
|
-1 for masked triangles (deleted from *compressed_triangles*).
|
||
|
|
Returned only if *return_tri_renum* is True.
|
||
|
|
node_renum : array-like of integers
|
||
|
|
renumbering table to translate the point numbers from the
|
||
|
|
encapsulated triangulation into the new (compressed) renumbering.
|
||
|
|
-1 for unused points (i.e. those deleted from *compressed_x* and
|
||
|
|
*compressed_y*). Returned only if *return_node_renum* is True.
|
||
|
|
|
||
|
|
"""
|
||
|
|
# Valid triangles and renumbering
|
||
|
|
tri_mask = self._triangulation.mask
|
||
|
|
compressed_triangles = self._triangulation.get_masked_triangles()
|
||
|
|
ntri = self._triangulation.triangles.shape[0]
|
||
|
|
tri_renum = self._total_to_compress_renum(tri_mask, ntri)
|
||
|
|
|
||
|
|
# Valid nodes and renumbering
|
||
|
|
node_mask = (np.bincount(np.ravel(compressed_triangles),
|
||
|
|
minlength=self._triangulation.x.size) == 0)
|
||
|
|
compressed_x = self._triangulation.x[~node_mask]
|
||
|
|
compressed_y = self._triangulation.y[~node_mask]
|
||
|
|
node_renum = self._total_to_compress_renum(node_mask)
|
||
|
|
|
||
|
|
# Now renumbering the valid triangles nodes
|
||
|
|
compressed_triangles = node_renum[compressed_triangles]
|
||
|
|
|
||
|
|
# 4 cases possible for return
|
||
|
|
if not return_tri_renum:
|
||
|
|
if not return_node_renum:
|
||
|
|
return compressed_triangles, compressed_x, compressed_y
|
||
|
|
else:
|
||
|
|
return (compressed_triangles, compressed_x, compressed_y,
|
||
|
|
node_renum)
|
||
|
|
else:
|
||
|
|
if not return_node_renum:
|
||
|
|
return (compressed_triangles, compressed_x, compressed_y,
|
||
|
|
tri_renum)
|
||
|
|
else:
|
||
|
|
return (compressed_triangles, compressed_x, compressed_y,
|
||
|
|
tri_renum, node_renum)
|
||
|
|
|
||
|
|
@staticmethod
|
||
|
|
def _total_to_compress_renum(mask, n=None):
|
||
|
|
"""
|
||
|
|
Parameters
|
||
|
|
----------
|
||
|
|
mask : 1d boolean array or None
|
||
|
|
mask
|
||
|
|
n : integer
|
||
|
|
length of the mask. Useful only id mask can be None
|
||
|
|
|
||
|
|
Returns
|
||
|
|
-------
|
||
|
|
renum : integer array
|
||
|
|
array so that (`valid_array` being a compressed array
|
||
|
|
based on a `masked_array` with mask *mask*) :
|
||
|
|
|
||
|
|
- For all i such as mask[i] = False:
|
||
|
|
valid_array[renum[i]] = masked_array[i]
|
||
|
|
- For all i such as mask[i] = True:
|
||
|
|
renum[i] = -1 (invalid value)
|
||
|
|
|
||
|
|
"""
|
||
|
|
if n is None:
|
||
|
|
n = np.size(mask)
|
||
|
|
if mask is not None:
|
||
|
|
renum = np.full(n, -1, dtype=np.int32) # Default num is -1
|
||
|
|
valid = np.arange(n, dtype=np.int32)[~mask]
|
||
|
|
renum[valid] = np.arange(np.size(valid, 0), dtype=np.int32)
|
||
|
|
return renum
|
||
|
|
else:
|
||
|
|
return np.arange(n, dtype=np.int32)
|