This commit is contained in:
2020-01-13 23:34:26 -07:00
commit 51e97c74a4
3096 changed files with 887666 additions and 0 deletions

View File

@@ -0,0 +1,12 @@
"""
Unstructured triangular grid functions.
"""
from .triangulation import *
from .tricontour import *
from .tritools import *
from .trifinder import *
from .triinterpolate import *
from .trirefine import *
from .tripcolor import *
from .triplot import *

View File

@@ -0,0 +1,221 @@
import numpy as np
class Triangulation(object):
"""
An unstructured triangular grid consisting of npoints points and
ntri triangles. The triangles can either be specified by the user
or automatically generated using a Delaunay triangulation.
Parameters
----------
x, y : array-like of shape (npoints)
Coordinates of grid points.
triangles : integer array_like of shape (ntri, 3), optional
For each triangle, the indices of the three points that make
up the triangle, ordered in an anticlockwise manner. If not
specified, the Delaunay triangulation is calculated.
mask : boolean array-like of shape (ntri), optional
Which triangles are masked out.
Attributes
----------
edges : int array of shape (nedges, 2)
See `~.Triangulation.edges`
neighbors : int array of shape (ntri, 3)
See `~.Triangulation.neighbors`
mask : bool array of shape (ntri, 3)
Masked out triangles.
is_delaunay : bool
Whether the Triangulation is a calculated Delaunay
triangulation (where `triangles` was not specified) or not.
Notes
-----
For a Triangulation to be valid it must not have duplicate points,
triangles formed from colinear points, or overlapping triangles.
"""
def __init__(self, x, y, triangles=None, mask=None):
from matplotlib import _qhull
self.x = np.asarray(x, dtype=np.float64)
self.y = np.asarray(y, dtype=np.float64)
if self.x.shape != self.y.shape or self.x.ndim != 1:
raise ValueError("x and y must be equal-length 1-D arrays")
self.mask = None
self._edges = None
self._neighbors = None
self.is_delaunay = False
if triangles is None:
# No triangulation specified, so use matplotlib._qhull to obtain
# Delaunay triangulation.
self.triangles, self._neighbors = _qhull.delaunay(x, y)
self.is_delaunay = True
else:
# Triangulation specified. Copy, since we may correct triangle
# orientation.
self.triangles = np.array(triangles, dtype=np.int32, order='C')
if self.triangles.ndim != 2 or self.triangles.shape[1] != 3:
raise ValueError('triangles must be a (?,3) array')
if self.triangles.max() >= len(self.x):
raise ValueError('triangles max element is out of bounds')
if self.triangles.min() < 0:
raise ValueError('triangles min element is out of bounds')
if mask is not None:
self.mask = np.asarray(mask, dtype=bool)
if self.mask.shape != (self.triangles.shape[0],):
raise ValueError('mask array must have same length as '
'triangles array')
# Underlying C++ object is not created until first needed.
self._cpp_triangulation = None
# Default TriFinder not created until needed.
self._trifinder = None
def calculate_plane_coefficients(self, z):
"""
Calculate plane equation coefficients for all unmasked triangles from
the point (x, y) coordinates and specified z-array of shape (npoints).
The returned array has shape (npoints, 3) and allows z-value at (x, y)
position in triangle tri to be calculated using
``z = array[tri, 0] * x + array[tri, 1] * y + array[tri, 2]``.
"""
return self.get_cpp_triangulation().calculate_plane_coefficients(z)
@property
def edges(self):
"""
Return integer array of shape (nedges, 2) containing all edges of
non-masked triangles.
Each row defines an edge by it's start point index and end point
index. Each edge appears only once, i.e. for an edge between points
*i* and *j*, there will only be either *(i, j)* or *(j, i)*.
"""
if self._edges is None:
self._edges = self.get_cpp_triangulation().get_edges()
return self._edges
def get_cpp_triangulation(self):
"""
Return the underlying C++ Triangulation object, creating it
if necessary.
"""
from matplotlib import _tri
if self._cpp_triangulation is None:
self._cpp_triangulation = _tri.Triangulation(
self.x, self.y, self.triangles, self.mask, self._edges,
self._neighbors, not self.is_delaunay)
return self._cpp_triangulation
def get_masked_triangles(self):
"""
Return an array of triangles that are not masked.
"""
if self.mask is not None:
return self.triangles[~self.mask]
else:
return self.triangles
@staticmethod
def get_from_args_and_kwargs(*args, **kwargs):
"""
Return a Triangulation object from the args and kwargs, and
the remaining args and kwargs with the consumed values removed.
There are two alternatives: either the first argument is a
Triangulation object, in which case it is returned, or the args
and kwargs are sufficient to create a new Triangulation to
return. In the latter case, see Triangulation.__init__ for
the possible args and kwargs.
"""
if isinstance(args[0], Triangulation):
triangulation = args[0]
args = args[1:]
else:
x = args[0]
y = args[1]
args = args[2:] # Consumed first two args.
# Check triangles in kwargs then args.
triangles = kwargs.pop('triangles', None)
from_args = False
if triangles is None and args:
triangles = args[0]
from_args = True
if triangles is not None:
try:
triangles = np.asarray(triangles, dtype=np.int32)
except ValueError:
triangles = None
if triangles is not None and (triangles.ndim != 2 or
triangles.shape[1] != 3):
triangles = None
if triangles is not None and from_args:
args = args[1:] # Consumed first item in args.
# Check for mask in kwargs.
mask = kwargs.pop('mask', None)
triangulation = Triangulation(x, y, triangles, mask)
return triangulation, args, kwargs
def get_trifinder(self):
"""
Return the default :class:`matplotlib.tri.TriFinder` of this
triangulation, creating it if necessary. This allows the same
TriFinder object to be easily shared.
"""
if self._trifinder is None:
# Default TriFinder class.
from matplotlib.tri.trifinder import TrapezoidMapTriFinder
self._trifinder = TrapezoidMapTriFinder(self)
return self._trifinder
@property
def neighbors(self):
"""
Return integer array of shape (ntri, 3) containing neighbor
triangles.
For each triangle, the indices of the three triangles that
share the same edges, or -1 if there is no such neighboring
triangle. neighbors[i,j] is the triangle that is the neighbor
to the edge from point index triangles[i,j] to point index
triangles[i,(j+1)%3].
"""
if self._neighbors is None:
self._neighbors = self.get_cpp_triangulation().get_neighbors()
return self._neighbors
def set_mask(self, mask):
"""
Set or clear the mask array. This is either None, or a boolean
array of shape (ntri).
"""
if mask is None:
self.mask = None
else:
self.mask = np.asarray(mask, dtype=bool)
if self.mask.shape != (self.triangles.shape[0],):
raise ValueError('mask array must have same length as '
'triangles array')
# Set mask in C++ Triangulation.
if self._cpp_triangulation is not None:
self._cpp_triangulation.set_mask(self.mask)
# Clear derived fields so they are recalculated when needed.
self._edges = None
self._neighbors = None
# Recalculate TriFinder if it exists.
if self._trifinder is not None:
self._trifinder._initialize()

View File

@@ -0,0 +1,278 @@
import numpy as np
from matplotlib.contour import ContourSet
from matplotlib.tri.triangulation import Triangulation
class TriContourSet(ContourSet):
"""
Create and store a set of contour lines or filled regions for
a triangular grid.
User-callable method: clabel
Useful attributes:
ax:
the axes object in which the contours are drawn
collections:
a silent_list of LineCollections or PolyCollections
levels:
contour levels
layers:
same as levels for line contours; half-way between
levels for filled contours. See _process_colors method.
"""
def __init__(self, ax, *args, **kwargs):
"""
Draw triangular grid contour lines or filled regions,
depending on whether keyword arg 'filled' is False
(default) or True.
The first argument of the initializer must be an axes
object. The remaining arguments and keyword arguments
are described in the docstring of `tricontour`.
"""
ContourSet.__init__(self, ax, *args, **kwargs)
def _process_args(self, *args, **kwargs):
"""
Process args and kwargs.
"""
if isinstance(args[0], TriContourSet):
C = args[0].cppContourGenerator
if self.levels is None:
self.levels = args[0].levels
else:
from matplotlib import _tri
tri, z = self._contour_args(args, kwargs)
C = _tri.TriContourGenerator(tri.get_cpp_triangulation(), z)
self._mins = [tri.x.min(), tri.y.min()]
self._maxs = [tri.x.max(), tri.y.max()]
self.cppContourGenerator = C
return kwargs
def _get_allsegs_and_allkinds(self):
"""
Create and return allsegs and allkinds by calling underlying C code.
"""
allsegs = []
if self.filled:
lowers, uppers = self._get_lowers_and_uppers()
allkinds = []
for lower, upper in zip(lowers, uppers):
segs, kinds = self.cppContourGenerator.create_filled_contour(
lower, upper)
allsegs.append([segs])
allkinds.append([kinds])
else:
allkinds = None
for level in self.levels:
segs = self.cppContourGenerator.create_contour(level)
allsegs.append(segs)
return allsegs, allkinds
def _contour_args(self, args, kwargs):
if self.filled:
fn = 'contourf'
else:
fn = 'contour'
tri, args, kwargs = Triangulation.get_from_args_and_kwargs(*args,
**kwargs)
z = np.ma.asarray(args[0])
if z.shape != tri.x.shape:
raise ValueError('z array must have same length as triangulation x'
' and y arrays')
# z values must be finite, only need to check points that are included
# in the triangulation.
z_check = z[np.unique(tri.get_masked_triangles())]
if np.ma.is_masked(z_check):
raise ValueError('z must not contain masked points within the '
'triangulation')
if not np.isfinite(z_check).all():
raise ValueError('z array must not contain non-finite values '
'within the triangulation')
z = np.ma.masked_invalid(z, copy=False)
self.zmax = float(z_check.max())
self.zmin = float(z_check.min())
if self.logscale and self.zmin <= 0:
raise ValueError('Cannot %s log of negative values.' % fn)
self._contour_level_args(z, args[1:])
return (tri, z)
def tricontour(ax, *args, **kwargs):
"""
Draw contours on an unstructured triangular grid.
`.tricontour` and `.tricontourf` draw contour lines and filled contours,
respectively. Except as noted, function signatures and return values are
the same for both versions.
The triangulation can be specified in one of two ways; either ::
tricontour(triangulation, ...)
where *triangulation* is a `matplotlib.tri.Triangulation` object, or ::
tricontour(x, y, ...)
tricontour(x, y, triangles, ...)
tricontour(x, y, triangles=triangles, ...)
tricontour(x, y, mask=mask, ...)
tricontour(x, y, triangles, mask=mask, ...)
in which case a `.Triangulation` object will be created. See that class'
docstring for an explanation of these cases.
The remaining arguments may be::
tricontour(..., Z)
where *Z* is the array of values to contour, one per point in the
triangulation. The level values are chosen automatically.
::
tricontour(..., Z, N)
contour up to *N+1* automatically chosen contour levels (*N* intervals).
::
tricontour(..., Z, V)
draw contour lines at the values specified in sequence *V*,
which must be in increasing order.
::
tricontourf(..., Z, V)
fill the (len(*V*)-1) regions between the values in *V*,
which must be in increasing order.
::
tricontour(Z, **kwargs)
Use keyword args to control colors, linewidth, origin, cmap ... see
below for more details.
`.tricontour(...)` returns a `~matplotlib.contour.TriContourSet` object.
Optional keyword arguments:
*colors*: [ *None* | string | (mpl_colors) ]
If *None*, the colormap specified by cmap will be used.
If a string, like 'r' or 'red', all levels will be plotted in this
color.
If a tuple of matplotlib color args (string, float, rgb, etc),
different levels will be plotted in different colors in the order
specified.
*alpha*: float
The alpha blending value
*cmap*: [ *None* | Colormap ]
A cm :class:`~matplotlib.colors.Colormap` instance or
*None*. If *cmap* is *None* and *colors* is *None*, a
default Colormap is used.
*norm*: [ *None* | Normalize ]
A :class:`matplotlib.colors.Normalize` instance for
scaling data values to colors. If *norm* is *None* and
*colors* is *None*, the default linear scaling is used.
*levels* [level0, level1, ..., leveln]
A list of floating point numbers indicating the level
curves to draw, in increasing order; e.g., to draw just
the zero contour pass ``levels=[0]``
*origin*: [ *None* | 'upper' | 'lower' | 'image' ]
If *None*, the first value of *Z* will correspond to the
lower left corner, location (0,0). If 'image', the rc
value for ``image.origin`` will be used.
This keyword is not active if *X* and *Y* are specified in
the call to contour.
*extent*: [ *None* | (x0,x1,y0,y1) ]
If *origin* is not *None*, then *extent* is interpreted as
in :func:`matplotlib.pyplot.imshow`: it gives the outer
pixel boundaries. In this case, the position of Z[0,0]
is the center of the pixel, not a corner. If *origin* is
*None*, then (*x0*, *y0*) is the position of Z[0,0], and
(*x1*, *y1*) is the position of Z[-1,-1].
This keyword is not active if *X* and *Y* are specified in
the call to contour.
*locator*: [ *None* | ticker.Locator subclass ]
If *locator* is None, the default
:class:`~matplotlib.ticker.MaxNLocator` is used. The
locator is used to determine the contour levels if they
are not given explicitly via the *V* argument.
*extend*: [ 'neither' | 'both' | 'min' | 'max' ]
Unless this is 'neither', contour levels are automatically
added to one or both ends of the range so that all data
are included. These added ranges are then mapped to the
special colormap values which default to the ends of the
colormap range, but can be set via
:meth:`matplotlib.colors.Colormap.set_under` and
:meth:`matplotlib.colors.Colormap.set_over` methods.
*xunits*, *yunits*: [ *None* | registered units ]
Override axis units by specifying an instance of a
:class:`matplotlib.units.ConversionInterface`.
tricontour-only keyword arguments:
*linewidths*: [ *None* | number | tuple of numbers ]
If *linewidths* is *None*, defaults to rc:`lines.linewidth`.
If a number, all levels will be plotted with this linewidth.
If a tuple, different levels will be plotted with different
linewidths in the order specified
*linestyles*: [ *None* | 'solid' | 'dashed' | 'dashdot' | 'dotted' ]
If *linestyles* is *None*, the 'solid' is used.
*linestyles* can also be an iterable of the above strings
specifying a set of linestyles to be used. If this
iterable is shorter than the number of contour levels
it will be repeated as necessary.
If contour is using a monochrome colormap and the contour
level is less than 0, then the linestyle specified
in :rc:`contour.negative_linestyle` will be used.
tricontourf-only keyword arguments:
*antialiased*: bool
enable antialiasing
Note: `.tricontourf` fills intervals that are closed at the top; that is,
for boundaries *z1* and *z2*, the filled region is::
z1 < Z <= z2
except for the lowest interval, which is closed on both sides (i.e. it
includes the lowest value).
"""
kwargs['filled'] = False
return TriContourSet(ax, *args, **kwargs)
def tricontourf(ax, *args, **kwargs):
kwargs['filled'] = True
return TriContourSet(ax, *args, **kwargs)
tricontourf.__doc__ = tricontour.__doc__

View File

@@ -0,0 +1,92 @@
import numpy as np
from matplotlib.tri import Triangulation
class TriFinder(object):
"""
Abstract base class for classes used to find the triangles of a
Triangulation in which (x,y) points lie.
Rather than instantiate an object of a class derived from TriFinder, it is
usually better to use the function
:func:`matplotlib.tri.Triangulation.get_trifinder`.
Derived classes implement __call__(x,y) where x,y are array_like point
coordinates of the same shape.
"""
def __init__(self, triangulation):
if not isinstance(triangulation, Triangulation):
raise ValueError('Expected a Triangulation object')
self._triangulation = triangulation
class TrapezoidMapTriFinder(TriFinder):
"""
:class:`~matplotlib.tri.TriFinder` class implemented using the trapezoid
map algorithm from the book "Computational Geometry, Algorithms and
Applications", second edition, by M. de Berg, M. van Kreveld, M. Overmars
and O. Schwarzkopf.
The triangulation must be valid, i.e. it must not have duplicate points,
triangles formed from colinear points, or overlapping triangles. The
algorithm has some tolerance to triangles formed from colinear points, but
this should not be relied upon.
"""
def __init__(self, triangulation):
from matplotlib import _tri
TriFinder.__init__(self, triangulation)
self._cpp_trifinder = _tri.TrapezoidMapTriFinder(
triangulation.get_cpp_triangulation())
self._initialize()
def __call__(self, x, y):
"""
Return an array containing the indices of the triangles in which the
specified x,y points lie, or -1 for points that do not lie within a
triangle.
*x*, *y* are array_like x and y coordinates of the same shape and any
number of dimensions.
Returns integer array with the same shape and *x* and *y*.
"""
x = np.asarray(x, dtype=np.float64)
y = np.asarray(y, dtype=np.float64)
if x.shape != y.shape:
raise ValueError("x and y must be array-like with the same shape")
# C++ does the heavy lifting, and expects 1D arrays.
indices = (self._cpp_trifinder.find_many(x.ravel(), y.ravel())
.reshape(x.shape))
return indices
def _get_tree_stats(self):
"""
Return a python list containing the statistics about the node tree:
0: number of nodes (tree size)
1: number of unique nodes
2: number of trapezoids (tree leaf nodes)
3: number of unique trapezoids
4: maximum parent count (max number of times a node is repeated in
tree)
5: maximum depth of tree (one more than the maximum number of
comparisons needed to search through the tree)
6: mean of all trapezoid depths (one more than the average number
of comparisons needed to search through the tree)
"""
return self._cpp_trifinder.get_tree_stats()
def _initialize(self):
"""
Initialize the underlying C++ object. Can be called multiple times if,
for example, the triangulation is modified.
"""
self._cpp_trifinder.initialize()
def _print_tree(self):
"""
Print a text representation of the node tree, which is useful for
debugging purposes.
"""
self._cpp_trifinder.print_tree()

File diff suppressed because it is too large Load Diff

View File

@@ -0,0 +1,138 @@
import numpy as np
from matplotlib import cbook
from matplotlib.collections import PolyCollection, TriMesh
from matplotlib.colors import Normalize
from matplotlib.tri.triangulation import Triangulation
def tripcolor(ax, *args, alpha=1.0, norm=None, cmap=None, vmin=None,
vmax=None, shading='flat', facecolors=None, **kwargs):
"""
Create a pseudocolor plot of an unstructured triangular grid.
The triangulation can be specified in one of two ways; either::
tripcolor(triangulation, ...)
where triangulation is a :class:`matplotlib.tri.Triangulation`
object, or
::
tripcolor(x, y, ...)
tripcolor(x, y, triangles, ...)
tripcolor(x, y, triangles=triangles, ...)
tripcolor(x, y, mask=mask, ...)
tripcolor(x, y, triangles, mask=mask, ...)
in which case a Triangulation object will be created. See
:class:`~matplotlib.tri.Triangulation` for a explanation of these
possibilities.
The next argument must be *C*, the array of color values, either
one per point in the triangulation if color values are defined at
points, or one per triangle in the triangulation if color values
are defined at triangles. If there are the same number of points
and triangles in the triangulation it is assumed that color
values are defined at points; to force the use of color values at
triangles use the kwarg ``facecolors=C`` instead of just ``C``.
*shading* may be 'flat' (the default) or 'gouraud'. If *shading*
is 'flat' and C values are defined at points, the color values
used for each triangle are from the mean C of the triangle's
three points. If *shading* is 'gouraud' then color values must be
defined at points.
The remaining kwargs are the same as for
:meth:`~matplotlib.axes.Axes.pcolor`.
"""
cbook._check_in_list(['flat', 'gouraud'], shading=shading)
tri, args, kwargs = Triangulation.get_from_args_and_kwargs(*args, **kwargs)
# C is the colors array defined at either points or faces (i.e. triangles).
# If facecolors is None, C are defined at points.
# If facecolors is not None, C are defined at faces.
if facecolors is not None:
C = facecolors
else:
C = np.asarray(args[0])
# If there are a different number of points and triangles in the
# triangulation, can omit facecolors kwarg as it is obvious from
# length of C whether it refers to points or faces.
# Do not do this for gouraud shading.
if (facecolors is None and len(C) == len(tri.triangles) and
len(C) != len(tri.x) and shading != 'gouraud'):
facecolors = C
# Check length of C is OK.
if ((facecolors is None and len(C) != len(tri.x)) or
(facecolors is not None and len(C) != len(tri.triangles))):
raise ValueError('Length of color values array must be the same '
'as either the number of triangulation points '
'or triangles')
# Handling of linewidths, shading, edgecolors and antialiased as
# in Axes.pcolor
linewidths = (0.25,)
if 'linewidth' in kwargs:
kwargs['linewidths'] = kwargs.pop('linewidth')
kwargs.setdefault('linewidths', linewidths)
edgecolors = 'none'
if 'edgecolor' in kwargs:
kwargs['edgecolors'] = kwargs.pop('edgecolor')
ec = kwargs.setdefault('edgecolors', edgecolors)
if 'antialiased' in kwargs:
kwargs['antialiaseds'] = kwargs.pop('antialiased')
if 'antialiaseds' not in kwargs and ec.lower() == "none":
kwargs['antialiaseds'] = False
if shading == 'gouraud':
if facecolors is not None:
raise ValueError('Gouraud shading does not support the use '
'of facecolors kwarg')
if len(C) != len(tri.x):
raise ValueError('For gouraud shading, the length of color '
'values array must be the same as the '
'number of triangulation points')
collection = TriMesh(tri, **kwargs)
else:
# Vertices of triangles.
maskedTris = tri.get_masked_triangles()
verts = np.stack((tri.x[maskedTris], tri.y[maskedTris]), axis=-1)
# Color values.
if facecolors is None:
# One color per triangle, the mean of the 3 vertex color values.
C = C[maskedTris].mean(axis=1)
elif tri.mask is not None:
# Remove color values of masked triangles.
C = C[~tri.mask]
collection = PolyCollection(verts, **kwargs)
collection.set_alpha(alpha)
collection.set_array(C)
if norm is not None and not isinstance(norm, Normalize):
raise ValueError("'norm' must be an instance of 'Normalize'")
collection.set_cmap(cmap)
collection.set_norm(norm)
if vmin is not None or vmax is not None:
collection.set_clim(vmin, vmax)
else:
collection.autoscale_None()
ax.grid(False)
minx = tri.x.min()
maxx = tri.x.max()
miny = tri.y.min()
maxy = tri.y.max()
corners = (minx, miny), (maxx, maxy)
ax.update_datalim(corners)
ax.autoscale_view()
ax.add_collection(collection)
return collection

View File

@@ -0,0 +1,85 @@
import numpy as np
from matplotlib.tri.triangulation import Triangulation
def triplot(ax, *args, **kwargs):
"""
Draw a unstructured triangular grid as lines and/or markers.
The triangulation to plot can be specified in one of two ways;
either::
triplot(triangulation, ...)
where triangulation is a :class:`matplotlib.tri.Triangulation`
object, or
::
triplot(x, y, ...)
triplot(x, y, triangles, ...)
triplot(x, y, triangles=triangles, ...)
triplot(x, y, mask=mask, ...)
triplot(x, y, triangles, mask=mask, ...)
in which case a Triangulation object will be created. See
:class:`~matplotlib.tri.Triangulation` for a explanation of these
possibilities.
The remaining args and kwargs are the same as for
:meth:`~matplotlib.axes.Axes.plot`.
Return a list of 2 :class:`~matplotlib.lines.Line2D` containing
respectively:
- the lines plotted for triangles edges
- the markers plotted for triangles nodes
"""
import matplotlib.axes
tri, args, kwargs = Triangulation.get_from_args_and_kwargs(*args, **kwargs)
x, y, edges = (tri.x, tri.y, tri.edges)
# Decode plot format string, e.g., 'ro-'
fmt = args[0] if args else ""
linestyle, marker, color = matplotlib.axes._base._process_plot_format(fmt)
# Insert plot format string into a copy of kwargs (kwargs values prevail).
kw = kwargs.copy()
for key, val in zip(('linestyle', 'marker', 'color'),
(linestyle, marker, color)):
if val is not None:
kw[key] = kwargs.get(key, val)
# Draw lines without markers.
# Note 1: If we drew markers here, most markers would be drawn more than
# once as they belong to several edges.
# Note 2: We insert nan values in the flattened edges arrays rather than
# plotting directly (triang.x[edges].T, triang.y[edges].T)
# as it considerably speeds-up code execution.
linestyle = kw['linestyle']
kw_lines = {
**kw,
'marker': 'None', # No marker to draw.
'zorder': kw.get('zorder', 1), # Path default zorder is used.
}
if linestyle not in [None, 'None', '', ' ']:
tri_lines_x = np.insert(x[edges], 2, np.nan, axis=1)
tri_lines_y = np.insert(y[edges], 2, np.nan, axis=1)
tri_lines = ax.plot(tri_lines_x.ravel(), tri_lines_y.ravel(),
**kw_lines)
else:
tri_lines = ax.plot([], [], **kw_lines)
# Draw markers separately.
marker = kw['marker']
kw_markers = {
**kw,
'linestyle': 'None', # No line to draw.
}
if marker not in [None, 'None', '', ' ']:
tri_markers = ax.plot(x, y, **kw_markers)
else:
tri_markers = ax.plot([], [], **kw_markers)
return tri_lines + tri_markers

View File

@@ -0,0 +1,315 @@
"""
Mesh refinement for triangular grids.
"""
import numpy as np
from matplotlib.tri.triangulation import Triangulation
import matplotlib.tri.triinterpolate
class TriRefiner(object):
"""
Abstract base class for classes implementing mesh refinement.
A TriRefiner encapsulates a Triangulation object and provides tools for
mesh refinement and interpolation.
Derived classes must implements:
- ``refine_triangulation(return_tri_index=False, **kwargs)`` , where
the optional keyword arguments *kwargs* are defined in each
TriRefiner concrete implementation, and which returns:
- a refined triangulation
- optionally (depending on *return_tri_index*), for each
point of the refined triangulation: the index of
the initial triangulation triangle to which it belongs.
- ``refine_field(z, triinterpolator=None, **kwargs)`` , where:
- *z* array of field values (to refine) defined at the base
triangulation nodes
- *triinterpolator* is a
:class:`~matplotlib.tri.TriInterpolator` (optional)
- the other optional keyword arguments *kwargs* are defined in
each TriRefiner concrete implementation
and which returns (as a tuple) a refined triangular mesh and the
interpolated values of the field at the refined triangulation nodes.
"""
def __init__(self, triangulation):
if not isinstance(triangulation, Triangulation):
raise ValueError("Expected a Triangulation object")
self._triangulation = triangulation
class UniformTriRefiner(TriRefiner):
"""
Uniform mesh refinement by recursive subdivisions.
Parameters
----------
triangulation : :class:`~matplotlib.tri.Triangulation`
The encapsulated triangulation (to be refined)
"""
# See Also
# --------
# :class:`~matplotlib.tri.CubicTriInterpolator` and
# :class:`~matplotlib.tri.TriAnalyzer`.
# """
def __init__(self, triangulation):
TriRefiner.__init__(self, triangulation)
def refine_triangulation(self, return_tri_index=False, subdiv=3):
"""
Computes an uniformly refined triangulation *refi_triangulation* of
the encapsulated :attr:`triangulation`.
This function refines the encapsulated triangulation by splitting each
father triangle into 4 child sub-triangles built on the edges midside
nodes, recursively (level of recursion *subdiv*).
In the end, each triangle is hence divided into ``4**subdiv``
child triangles.
The default value for *subdiv* is 3 resulting in 64 refined
subtriangles for each triangle of the initial triangulation.
Parameters
----------
return_tri_index : boolean, optional
Boolean indicating whether an index table indicating the father
triangle index of each point will be returned. Default value
False.
subdiv : integer, optional
Recursion level for the subdivision. Defaults value 3.
Each triangle will be divided into ``4**subdiv`` child triangles.
Returns
-------
refi_triangulation : :class:`~matplotlib.tri.Triangulation`
The returned refined triangulation
found_index : array-like of integers
Index of the initial triangulation containing triangle, for each
point of *refi_triangulation*.
Returned only if *return_tri_index* is set to True.
"""
refi_triangulation = self._triangulation
ntri = refi_triangulation.triangles.shape[0]
# Computes the triangulation ancestors numbers in the reference
# triangulation.
ancestors = np.arange(ntri, dtype=np.int32)
for _ in range(subdiv):
refi_triangulation, ancestors = self._refine_triangulation_once(
refi_triangulation, ancestors)
refi_npts = refi_triangulation.x.shape[0]
refi_triangles = refi_triangulation.triangles
# Now we compute found_index table if needed
if return_tri_index:
# We have to initialize found_index with -1 because some nodes
# may very well belong to no triangle at all, e.g., in case of
# Delaunay Triangulation with DuplicatePointWarning.
found_index = np.full(refi_npts, -1, dtype=np.int32)
tri_mask = self._triangulation.mask
if tri_mask is None:
found_index[refi_triangles] = np.repeat(ancestors,
3).reshape(-1, 3)
else:
# There is a subtlety here: we want to avoid whenever possible
# that refined points container is a masked triangle (which
# would result in artifacts in plots).
# So we impose the numbering from masked ancestors first,
# then overwrite it with unmasked ancestor numbers.
ancestor_mask = tri_mask[ancestors]
found_index[refi_triangles[ancestor_mask, :]
] = np.repeat(ancestors[ancestor_mask],
3).reshape(-1, 3)
found_index[refi_triangles[~ancestor_mask, :]
] = np.repeat(ancestors[~ancestor_mask],
3).reshape(-1, 3)
return refi_triangulation, found_index
else:
return refi_triangulation
def refine_field(self, z, triinterpolator=None, subdiv=3):
"""
Refines a field defined on the encapsulated triangulation.
Returns *refi_tri* (refined triangulation), *refi_z* (interpolated
values of the field at the node of the refined triangulation).
Parameters
----------
z : 1d-array-like of length ``n_points``
Values of the field to refine, defined at the nodes of the
encapsulated triangulation. (``n_points`` is the number of points
in the initial triangulation)
triinterpolator : :class:`~matplotlib.tri.TriInterpolator`, optional
Interpolator used for field interpolation. If not specified,
a :class:`~matplotlib.tri.CubicTriInterpolator` will
be used.
subdiv : integer, optional
Recursion level for the subdivision. Defaults to 3.
Each triangle will be divided into ``4**subdiv`` child triangles.
Returns
-------
refi_tri : :class:`~matplotlib.tri.Triangulation` object
The returned refined triangulation
refi_z : 1d array of length: *refi_tri* node count.
The returned interpolated field (at *refi_tri* nodes)
"""
if triinterpolator is None:
interp = matplotlib.tri.CubicTriInterpolator(
self._triangulation, z)
else:
if not isinstance(triinterpolator,
matplotlib.tri.TriInterpolator):
raise ValueError("Expected a TriInterpolator object")
interp = triinterpolator
refi_tri, found_index = self.refine_triangulation(
subdiv=subdiv, return_tri_index=True)
refi_z = interp._interpolate_multikeys(
refi_tri.x, refi_tri.y, tri_index=found_index)[0]
return refi_tri, refi_z
@staticmethod
def _refine_triangulation_once(triangulation, ancestors=None):
"""
This function refines a matplotlib.tri *triangulation* by splitting
each triangle into 4 child-masked_triangles built on the edges midside
nodes.
The masked triangles, if present, are also split but their children
returned masked.
If *ancestors* is not provided, returns only a new triangulation:
child_triangulation.
If the array-like key table *ancestor* is given, it shall be of shape
(ntri,) where ntri is the number of *triangulation* masked_triangles.
In this case, the function returns
(child_triangulation, child_ancestors)
child_ancestors is defined so that the 4 child masked_triangles share
the same index as their father: child_ancestors.shape = (4 * ntri,).
"""
x = triangulation.x
y = triangulation.y
# According to tri.triangulation doc:
# neighbors[i,j] is the triangle that is the neighbor
# to the edge from point index masked_triangles[i,j] to point
# index masked_triangles[i,(j+1)%3].
neighbors = triangulation.neighbors
triangles = triangulation.triangles
npts = np.shape(x)[0]
ntri = np.shape(triangles)[0]
if ancestors is not None:
ancestors = np.asarray(ancestors)
if np.shape(ancestors) != (ntri,):
raise ValueError(
"Incompatible shapes provide for triangulation"
".masked_triangles and ancestors: {0} and {1}".format(
np.shape(triangles), np.shape(ancestors)))
# Initiating tables refi_x and refi_y of the refined triangulation
# points
# hint: each apex is shared by 2 masked_triangles except the borders.
borders = np.sum(neighbors == -1)
added_pts = (3*ntri + borders) // 2
refi_npts = npts + added_pts
refi_x = np.zeros(refi_npts)
refi_y = np.zeros(refi_npts)
# First part of refi_x, refi_y is just the initial points
refi_x[:npts] = x
refi_y[:npts] = y
# Second part contains the edge midside nodes.
# Each edge belongs to 1 triangle (if border edge) or is shared by 2
# masked_triangles (interior edge).
# We first build 2 * ntri arrays of edge starting nodes (edge_elems,
# edge_apexes); we then extract only the masters to avoid overlaps.
# The so-called 'master' is the triangle with biggest index
# The 'slave' is the triangle with lower index
# (can be -1 if border edge)
# For slave and master we will identify the apex pointing to the edge
# start
edge_elems = np.tile(np.arange(ntri, dtype=np.int32), 3)
edge_apexes = np.repeat(np.arange(3, dtype=np.int32), ntri)
edge_neighbors = neighbors[edge_elems, edge_apexes]
mask_masters = (edge_elems > edge_neighbors)
# Identifying the "masters" and adding to refi_x, refi_y vec
masters = edge_elems[mask_masters]
apex_masters = edge_apexes[mask_masters]
x_add = (x[triangles[masters, apex_masters]] +
x[triangles[masters, (apex_masters+1) % 3]]) * 0.5
y_add = (y[triangles[masters, apex_masters]] +
y[triangles[masters, (apex_masters+1) % 3]]) * 0.5
refi_x[npts:] = x_add
refi_y[npts:] = y_add
# Building the new masked_triangles; each old masked_triangles hosts
# 4 new masked_triangles
# there are 6 pts to identify per 'old' triangle, 3 new_pt_corner and
# 3 new_pt_midside
new_pt_corner = triangles
# What is the index in refi_x, refi_y of point at middle of apex iapex
# of elem ielem ?
# If ielem is the apex master: simple count, given the way refi_x was
# built.
# If ielem is the apex slave: yet we do not know; but we will soon
# using the neighbors table.
new_pt_midside = np.empty([ntri, 3], dtype=np.int32)
cum_sum = npts
for imid in range(3):
mask_st_loc = (imid == apex_masters)
n_masters_loc = np.sum(mask_st_loc)
elem_masters_loc = masters[mask_st_loc]
new_pt_midside[:, imid][elem_masters_loc] = np.arange(
n_masters_loc, dtype=np.int32) + cum_sum
cum_sum += n_masters_loc
# Now dealing with slave elems.
# for each slave element we identify the master and then the inode
# once slave_masters is identified, slave_masters_apex is such that:
# neighbors[slaves_masters, slave_masters_apex] == slaves
mask_slaves = np.logical_not(mask_masters)
slaves = edge_elems[mask_slaves]
slaves_masters = edge_neighbors[mask_slaves]
diff_table = np.abs(neighbors[slaves_masters, :] -
np.outer(slaves, np.ones(3, dtype=np.int32)))
slave_masters_apex = np.argmin(diff_table, axis=1)
slaves_apex = edge_apexes[mask_slaves]
new_pt_midside[slaves, slaves_apex] = new_pt_midside[
slaves_masters, slave_masters_apex]
# Builds the 4 child masked_triangles
child_triangles = np.empty([ntri*4, 3], dtype=np.int32)
child_triangles[0::4, :] = np.vstack([
new_pt_corner[:, 0], new_pt_midside[:, 0],
new_pt_midside[:, 2]]).T
child_triangles[1::4, :] = np.vstack([
new_pt_corner[:, 1], new_pt_midside[:, 1],
new_pt_midside[:, 0]]).T
child_triangles[2::4, :] = np.vstack([
new_pt_corner[:, 2], new_pt_midside[:, 2],
new_pt_midside[:, 1]]).T
child_triangles[3::4, :] = np.vstack([
new_pt_midside[:, 0], new_pt_midside[:, 1],
new_pt_midside[:, 2]]).T
child_triangulation = Triangulation(refi_x, refi_y, child_triangles)
# Builds the child mask
if triangulation.mask is not None:
child_triangulation.set_mask(np.repeat(triangulation.mask, 4))
if ancestors is None:
return child_triangulation
else:
return child_triangulation, np.repeat(ancestors, 4)

View File

@@ -0,0 +1,299 @@
"""
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)