Skip to content

Commit

Permalink
Use qhull for delaunay triangulation
Browse files Browse the repository at this point in the history
  • Loading branch information
@ianthomas23
ianthomas23 committed Nov 30, 2013
1 parent b35ad23 commit a65066d
Show file tree
Hide file tree
Showing 13 changed files with 807 additions and 1,514 deletions.
44 changes: 24 additions & 20 deletions CHANGELOG
View file
Original file line number Diff line number Diff line change
@@ -1,6 +1,11 @@
2013-11-28 Added qhull extension module to perform Delaunay triangulation more
robustly than before. It is used by tri.Triangulation (and hence
all pyplot.tri* methods) and mlab.griddata. Deprecated
matplotlib.delaunay module. - IMT

2013-10-27 Added get_rlabel_position and set_rlabel_position methods to
PolarAxes to control angular position of radial tick labels.
PolarAxes to control angular position of radial tick labels.

2013-10-06 Add stride-based functions to mlab for easy creation of 2D arrays
with less memory.

Expand All @@ -27,16 +32,16 @@
and matplotlib.units.Registry.

2013-06-26 Refactored the axes module: the axes module is now a folder,
containing the following submodule:
- _subplots.py, containing all the subplots helper methods
- _base.py, containing several private methods and a new
_AxesBase class. This _AxesBase class contains all the methods
that are not directly linked to plots of the "old" Axes
- _axes.py contains the Axes class. This class now inherits from
_AxesBase: it contains all "plotting" methods and labelling
methods.
This refactoring should not affect the API. Only private methods
are not importable from the axes module anymore.
containing the following submodule:
- _subplots.py, containing all the subplots helper methods
- _base.py, containing several private methods and a new
_AxesBase class. This _AxesBase class contains all the methods
that are not directly linked to plots of the "old" Axes
- _axes.py contains the Axes class. This class now inherits from
_AxesBase: it contains all "plotting" methods and labelling
methods.
This refactoring should not affect the API. Only private methods
are not importable from the axes module anymore.

2013-05-18 Added support for arbitrary rasterization resolutions to the
SVG backend. Previously the resolution was hard coded to 72
Expand All @@ -52,22 +57,21 @@

2013-04-25 Changed all instances of:

from matplotlib import MatplotlibDeprecationWarning as mplDeprecation
to:

from cbook import mplDeprecation
from matplotlib import MatplotlibDeprecationWarning as mplDeprecation
to:

and removed the import into the matplotlib namespace in __init__.py
Thomas Caswell
from cbook import mplDeprecation

and removed the import into the matplotlib namespace in __init__.py
Thomas Caswell

2013-04-15 Added 'axes.xmargin' and 'axes.ymargin' to rpParams to set default
margins on auto-scaleing. - TAC
margins on auto-scaleing. - TAC

2013-04-16 Added patheffect support for Line2D objects. -JJL

2013-03-19 Added support for passing `linestyle` kwarg to `step` so all `plot`
kwargs are passed to the underlying `plot` call. -TAC
kwargs are passed to the underlying `plot` call. -TAC

2013-02-25 Added classes CubicTriInterpolator, UniformTriRefiner, TriAnalyzer
to matplotlib.tri module. - GBy
Expand Down
4 changes: 1 addition & 3 deletions MANIFEST.in
View file
Original file line number Diff line number Diff line change
Expand Up @@ -14,7 +14,5 @@ recursive-include LICENSE *
recursive-include examples *
recursive-include doc *
recursive-include src *.cpp *.c *.h *.m
recursive-include CXX *.cxx *.hxx *.c *.h
recursive-include agg24 *
recursive-include lib *
recursive-include ttconv *.cpp *.h
recursive-include extern *
38 changes: 12 additions & 26 deletions examples/pylab_examples/griddata_demo.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -6,36 +6,22 @@
#npts = int(raw_input('enter # of random points to plot:'))
seed(0)
npts = 200
x = uniform(-2,2,npts)
y = uniform(-2,2,npts)
z = x*np.exp(-x**2-y**2)
x = uniform(-2, 2, npts)
y = uniform(-2, 2, npts)
z = x*np.exp(-x**2 - y**2)
# define grid.
xi = np.linspace(-2.1,2.1,100)
yi = np.linspace(-2.1,2.1,200)
xi = np.linspace(-2.1, 2.1, 100)
yi = np.linspace(-2.1, 2.1, 200)
# grid the data.
zi = griddata(x,y,z,xi,yi,interp='linear')
zi = griddata(x, y, z, xi, yi, interp='linear')
# contour the gridded data, plotting dots at the nonuniform data points.
CS = plt.contour(xi,yi,zi,15,linewidths=0.5,colors='k')
CS = plt.contourf(xi,yi,zi,15,cmap=plt.cm.rainbow,
CS = plt.contour(xi, yi, zi, 15, linewidths=0.5, colors='k')
CS = plt.contourf(xi, yi, zi, 15, cmap=plt.cm.rainbow,
vmax=abs(zi).max(), vmin=-abs(zi).max())
plt.colorbar() # draw colorbar
plt.colorbar() # draw colorbar
# plot data points.
plt.scatter(x,y,marker='o',c='b',s=5,zorder=10)
plt.xlim(-2,2)
plt.ylim(-2,2)
plt.scatter(x, y, marker='o', c='b', s=5, zorder=10)
plt.xlim(-2, 2)
plt.ylim(-2, 2)
plt.title('griddata test (%d points)' % npts)
plt.show()

# test case that scikits.delaunay fails on, but natgrid passes..
#data = np.array([[-1, -1], [-1, 0], [-1, 1],
# [ 0, -1], [ 0, 0], [ 0, 1],
# [ 1, -1 - np.finfo(np.float_).eps], [ 1, 0], [ 1, 1],
# ])
#x = data[:,0]
#y = data[:,1]
#z = x*np.exp(-x**2-y**2)
## define grid.
#xi = np.linspace(-1.1,1.1,100)
#yi = np.linspace(-1.1,1.1,100)
## grid the data.
#zi = griddata(x,y,z,xi,yi)
5 changes: 5 additions & 0 deletions lib/matplotlib/delaunay/triangulate.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -11,6 +11,11 @@

from matplotlib._delaunay import delaunay
from .interpolate import LinearInterpolator, NNInterpolator
from matplotlib.cbook import warn_deprecated
warn_deprecated('1.4',
name='matplotlib.delaunay',
alternative='matplotlib.tri.Triangulation',
obj_type='module')

__all__ = ['Triangulation', 'DuplicatePointWarning']

Expand Down
239 changes: 123 additions & 116 deletions lib/matplotlib/mlab.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3356,124 +3356,131 @@ def newfunc(val, mask, mval):
if opened:
fh.close()

def griddata(x,y,z,xi,yi,interp='nn'):
"""
``zi = griddata(x,y,z,xi,yi)`` fits a surface of the form *z* =
*f*(*x*, *y*) to the data in the (usually) nonuniformly spaced
vectors (*x*, *y*, *z*). :func:`griddata` interpolates this
surface at the points specified by (*xi*, *yi*) to produce
*zi*. *xi* and *yi* must describe a regular grid, can be either 1D
or 2D, but must be monotonically increasing.
A masked array is returned if any grid points are outside convex
hull defined by input data (no extrapolation is done).
If interp keyword is set to '`nn`' (default),
uses natural neighbor interpolation based on Delaunay
triangulation. By default, this algorithm is provided by the
:mod:`matplotlib.delaunay` package, written by Robert Kern. The
triangulation algorithm in this package is known to fail on some
nearly pathological cases. For this reason, a separate toolkit
(:mod:`mpl_tookits.natgrid`) has been created that provides a more
robust algorithm fof triangulation and interpolation. This
toolkit is based on the NCAR natgrid library, which contains code
that is not redistributable under a BSD-compatible license. When
installed, this function will use the :mod:`mpl_toolkits.natgrid`
algorithm, otherwise it will use the built-in
:mod:`matplotlib.delaunay` package.
If the interp keyword is set to '`linear`', then linear interpolation
is used instead of natural neighbor. In this case, the output grid
is assumed to be regular with a constant grid spacing in both the x and
y directions. For regular grids with nonconstant grid spacing, you
must use natural neighbor interpolation. Linear interpolation is only valid if
:mod:`matplotlib.delaunay` package is used - :mod:`mpl_tookits.natgrid`
only provides natural neighbor interpolation.
The natgrid matplotlib toolkit can be downloaded from
http://sourceforge.net/project/showfiles.php?group_id=80706&package_id=142792
"""
try:
from mpl_toolkits.natgrid import _natgrid, __version__
_use_natgrid = True
except ImportError:
import matplotlib.delaunay as delaunay
from matplotlib.delaunay import __version__
_use_natgrid = False
if not griddata._reported:
if _use_natgrid:
verbose.report('using natgrid version %s' % __version__)
else:
verbose.report('using delaunay version %s' % __version__)
griddata._reported = True

def griddata(x, y, z, xi, yi, interp='nn'):
"""Interpolates from a nonuniformly spaced grid to some other
grid.
Fits a surface of the form z = f(`x`, `y`) to the data in the
(usually) nonuniformly spaced vectors (`x`, `y`, `z`), then
interpolates this surface at the points specified by
(`xi`, `yi`) to produce `zi`.
Parameters
----------
x, y, z : 1d array_like
Coordinates of grid points to interpolate from.
xi, yi : 1d or 2d array_like
Coordinates of grid points to interpolate to.
interp : string key from {'nn', 'linear'}
Interpolation algorithm, either 'nn' for natural neighbor, or
'linear' for linear interpolation.
Returns
-------
2d float array
Array of values interpolated at (`xi`, `yi`) points. Array
will be masked is any of (`xi`, `yi`) are outside the convex
hull of (`x`, `y`).
Notes
-----
If `interp` is 'nn' (the default), uses natural neighbor
interpolation based on Delaunay triangulation. This option is
only available if the mpl_toolkits.natgrid module is installed.
This can be downloaded from https://github.com/matplotlib/natgrid.
The (`xi`, `yi`) grid must be regular and monotonically increasing
in this case.
If `interp` is 'linear', linear interpolation is used via
matplotlib.tri.LinearTriInterpolator.
Instead of using `griddata`, more flexible functionality and other
interpolation options are available using a
matplotlib.tri.Triangulation and a matplotlib.tri.TriInterpolator.
"""
# Check input arguments.
x = np.asanyarray(x, dtype=np.float64)
y = np.asanyarray(y, dtype=np.float64)
z = np.asanyarray(z, dtype=np.float64)
if x.shape != y.shape or x.shape != z.shape or x.ndim != 1:
raise ValueError("x, y and z must be equal-length 1-D arrays")

xi = np.asanyarray(xi, dtype=np.float64)
yi = np.asanyarray(yi, dtype=np.float64)
if xi.ndim != yi.ndim:
raise TypeError("inputs xi and yi must have same number of dimensions (1 or 2)")
if xi.ndim != 1 and xi.ndim != 2:
raise TypeError("inputs xi and yi must be 1D or 2D.")
if not len(x)==len(y)==len(z):
raise TypeError("inputs x,y,z must all be 1D arrays of the same length")
# remove masked points.
if hasattr(z,'mask'):
# make sure mask is not a scalar boolean array.
if z.mask.ndim:
x = x.compress(z.mask == False)
y = y.compress(z.mask == False)
z = z.compressed()
if _use_natgrid: # use natgrid toolkit if available.
if interp != 'nn':
raise ValueError("only natural neighor interpolation"
" allowed when using natgrid toolkit in griddata.")
raise ValueError("xi and yi must be arrays with the same number of "
"dimensions (1 or 2)")
if xi.ndim == 2 and xi.shape != yi.shape:
raise ValueError("if xi and yi are 2D arrays, they must have the same "
"shape")
if xi.ndim == 1:
xi, yi = np.meshgrid(xi, yi)

if interp == 'nn':
use_nn_interpolation = True
elif interp == 'linear':
use_nn_interpolation = False
else:
raise ValueError("interp keyword must be one of 'linear' (for linear "
"interpolation) or 'nn' (for natural neighbor "
"interpolation). Default is 'nn'.")

# Remove masked points.
mask = np.ma.getmask(z)
if not (mask is np.ma.nomask):
x = x.compress(~mask)
y = y.compress(~mask)
z = z.compressed()

if use_nn_interpolation:
try:
from mpl_toolkits.natgrid import _natgrid
except ImportError:
raise RuntimeError("To use interp='nn' (Natural Neighbor "
"interpolation) in griddata, natgrid must be installed. "
"Either install it from http://sourceforge.net/projects/"
"matplotlib/files/matplotlib-toolkits, or use interp='linear' "
"instead.")

if xi.ndim == 2:
xi = xi[0,:]
yi = yi[:,0]
# override default natgrid internal parameters.
_natgrid.seti('ext',0)
_natgrid.setr('nul',np.nan)
# cast input arrays to doubles (this makes a copy)
x = x.astype(np.float)
y = y.astype(np.float)
z = z.astype(np.float)
xo = xi.astype(np.float)
yo = yi.astype(np.float)
if min(xo[1:]-xo[0:-1]) < 0 or min(yo[1:]-yo[0:-1]) < 0:
raise ValueError('output grid defined by xi,yi must be monotone increasing')
# allocate array for output (buffer will be overwritten by nagridd)
zo = np.empty((yo.shape[0],xo.shape[0]), np.float)
_natgrid.natgridd(x,y,z,xo,yo,zo)
else: # use Robert Kern's delaunay package from scikits (default)
if xi.ndim != yi.ndim:
raise TypeError("inputs xi and yi must have same number of dimensions (1 or 2)")
if xi.ndim != 1 and xi.ndim != 2:
raise TypeError("inputs xi and yi must be 1D or 2D.")
if xi.ndim == 1:
xi,yi = np.meshgrid(xi,yi)
# triangulate data
tri = delaunay.Triangulation(x,y)
# interpolate data
if interp == 'nn':
interp = tri.nn_interpolator(z)
zo = interp(xi,yi)
elif interp == 'linear':
# make sure grid has constant dx, dy
dx = xi[0,1:]-xi[0,0:-1]
dy = yi[1:,0]-yi[0:-1,0]
epsx = np.finfo(xi.dtype).resolution
epsy = np.finfo(yi.dtype).resolution
if dx.max()-dx.min() > epsx or dy.max()-dy.min() > epsy:
raise ValueError("output grid must have constant spacing"
" when using interp='linear'")
interp = tri.linear_interpolator(z)
zo = interp[yi.min():yi.max():complex(0,yi.shape[0]),
xi.min():xi.max():complex(0,xi.shape[1])]
else:
raise ValueError("interp keyword must be one of"
" 'linear' (for linear interpolation) or 'nn'"
" (for natural neighbor interpolation). Default is 'nn'.")
# mask points on grid outside convex hull of input data.
if np.any(np.isnan(zo)):
zo = np.ma.masked_where(np.isnan(zo),zo)
return zo
griddata._reported = False
# natgrid expects 1D xi and yi arrays.
xi = xi[0, :]
yi = yi[:, 0]

# Override default natgrid internal parameters.
_natgrid.seti('ext', 0)
_natgrid.setr('nul', np.nan)

if np.min(np.diff(xi)) < 0 or np.min(np.diff(yi)) < 0:
raise ValueError("Output grid defined by xi,yi must be monotone "
"increasing")

# Allocate array for output (buffer will be overwritten by natgridd)
zi = np.empty((yi.shape[0], xi.shape[0]), np.float64)

# Natgrid requires each array to be contiguous rather than e.g. a view
# that is a non-contiguous slice of another array. Use numpy.require
# to deal with this, which will copy if necessary.
x = np.require(x, requirements=['C'])
y = np.require(y, requirements=['C'])
z = np.require(z, requirements=['C'])
xi = np.require(xi, requirements=['C'])
yi = np.require(yi, requirements=['C'])
_natgrid.natgridd(x, y, z, xi, yi, zi)

# Mask points on grid outside convex hull of input data.
if np.any(np.isnan(zi)):
zi = np.ma.masked_where(np.isnan(zi), zi)
return zi
else:
# Linear interpolation performed using a matplotlib.tri.Triangulation
# and a matplotlib.tri.LinearTriInterpolator.
from .tri import Triangulation, LinearTriInterpolator
triang = Triangulation(x, y)
interpolator = LinearTriInterpolator(triang, z)
return interpolator(xi, yi)


##################################################
# Linear interpolation algorithms
Expand Down
Binary file removed lib/matplotlib/tests/baseline_images/test_triangulation/tripcolor1.pdf
View file
Binary file not shown.
Loading

0 comments on commit a65066d

Please sign in to comment.