Added LinearTriInterpolator class.

Author: ianthomas23

CHANGELOG

@@ -1,3 +1,7 @@
+2012-12-22 Added classes for interpolation within triangular grids
+           (LinearTriInterpolator) and to find the triangles in which points
+           lie (TrapezoidMapTriFinder) to matplotlib.tri module. - IMT
+
 2012-12-05 Added MatplotlibDeprecationWarning class for signaling deprecation.
            Matplotlib developers can use this class as follows:
 

doc/api/tri_api.rst

@@ -13,3 +13,9 @@ triangular grids
 
 .. autoclass:: matplotlib.tri.TrapezoidMapTriFinder
    :members: __call__
+
+.. autoclass:: matplotlib.tri.TriInterpolator
+   :members:
+
+.. autoclass:: matplotlib.tri.LinearTriInterpolator
+   :members: __call__

doc/users/whats_new.rst

@@ -49,6 +49,13 @@ using Inkscape for example, while preserving their intended position. For
 `svg` please note that you'll have to disable the default text-to-path
 conversion (`mpl.rc('svg', fonttype='none')`).
 
+Triangular grid interpolation
+-----------------------------
+Ian Thomas added classes to perform interpolation within triangular grids
+(:class:`~matplotlib.tri.LinearTriInterpolator`) and a utility class to find
+the triangles in which points lie (
+:class:`~matplotlib.tri.TrapezoidMapTriFinder`).
+
 .. _whats-new-1-2:
 
 new in matplotlib-1.2
@@ -298,7 +305,7 @@ to address the most common layout issues.
     fig, axes_list = plt.subplots(2, 1)
     for ax in axes_list.flat:
         ax.set(xlabel="x-label", ylabel="y-label", title="before tight_layout")
-	ax.locator_params(nbins=3)
+    ax.locator_params(nbins=3)
 
     plt.show()
 
@@ -308,7 +315,7 @@ to address the most common layout issues.
     fig, axes_list = plt.subplots(2, 1)
     for ax in axes_list.flat:
         ax.set(xlabel="x-label", ylabel="y-label", title="after tight_layout")
-	ax.locator_params(nbins=3)
+    ax.locator_params(nbins=3)
 
     plt.tight_layout()
     plt.show()

examples/pylab_examples/triinterp_demo.py

@@ -0,0 +1,36 @@
+"""
+Interpolation from triangular grid to quad grid.
+"""
+import matplotlib.pyplot as plt
+import matplotlib.tri as mtri
+import numpy as np
+import math
+
+
+# Create triangulation.
+x = np.asarray([0, 1, 2, 3, 0.5, 1.5, 2.5, 1, 2, 1.5])
+y = np.asarray([0, 0, 0, 0, 1,   1,   1,   2, 2, 3])
+triangles = [[0,1,4], [1,2,5], [2,3,6], [1,5,4], [2,6,5], [4,5,7], [5,6,8],
+             [5,8,7], [7,8,9]]
+triang = mtri.Triangulation(x, y, triangles)
+
+# Interpolate to regularly-spaced quad grid.
+z = np.cos(1.5*x)*np.cos(1.5*y)
+interp = mtri.LinearTriInterpolator(triang, z)
+xi, yi = np.meshgrid(np.linspace(0, 3, 20), np.linspace(0, 3, 20))
+zi = interp(xi, yi)
+
+# Plot the triangulation.
+plt.subplot(121)
+plt.tricontourf(triang, z)
+plt.triplot(triang, 'ko-')
+plt.title('Triangular grid')
+
+# Plot interpolation to quad grid.
+plt.subplot(122)
+plt.contourf(xi, yi, zi)
+plt.plot(xi, yi, 'k-', alpha=0.5)
+plt.plot(xi.T, yi.T, 'k-', alpha=0.5)
+plt.title('Linear interpolation')
+
+plt.show()

lib/matplotlib/tests/test_triangulation.py

@@ -219,3 +219,34 @@ def test_trifinder():
     assert_equal(trifinder, triang.get_trifinder())
     tris = trifinder(xs, ys)
     assert_array_equal(tris, [-1, -1, 1, -1])
+
+def test_triinterp():
+    # Test points within triangles of masked triangulation.
+    x,y = np.meshgrid(np.arange(4), np.arange(4))
+    x = x.ravel()
+    y = y.ravel()
+    z = 1.23*x - 4.79*y
+    triangles = [[0,1,4], [1,5,4], [1,2,5], [2,6,5], [2,3,6], [3,7,6], [4,5,8],
+                 [5,9,8], [5,6,9], [6,10,9], [6,7,10], [7,11,10], [8,9,12],
+                 [9,13,12], [9,10,13], [10,14,13], [10,11,14], [11,15,14]]
+    mask = np.zeros(len(triangles))
+    mask[8:10] = 1
+    triang = mtri.Triangulation(x, y, triangles, mask)
+    linear_interp = mtri.LinearTriInterpolator(triang, z)
+
+    xs = np.linspace(0.25, 2.75, 6)
+    ys = [0.25, 0.75, 2.25, 2.75]
+    xs,ys = np.meshgrid(xs,ys)
+    xs = xs.ravel()
+    ys = ys.ravel()
+    zs = linear_interp(xs, ys)
+    assert_array_almost_equal(zs, (1.23*xs - 4.79*ys))
+
+    # Test points outside triangulation.
+    xs = [-0.25, 1.25, 1.75, 3.25]
+    ys = xs
+    xs, ys = np.meshgrid(xs,ys)
+    xs = xs.ravel()
+    ys = ys.ravel()
+    zs = linear_interp(xs, ys)
+    assert_array_equal(zs.mask, [True]*16)

lib/matplotlib/tri/__init__.py

@@ -6,5 +6,6 @@
 from triangulation import *
 from tricontour import *
 from trifinder import *
+from triinterpolate import *
 from tripcolor import *
 from triplot import *

lib/matplotlib/tri/_tri.cpp

@@ -137,6 +137,11 @@ double XYZ::dot(const XYZ& other) const
     return x*other.x + y*other.y + z*other.z;
 }
 
+double XYZ::length_squared() const
+{
+    return x*x + y*y + z*z;
+}
+
 XYZ XYZ::operator-(const XYZ& other) const
 {
     return XYZ(x - other.x, y - other.y, z - other.z);
@@ -375,6 +380,95 @@ void Triangulation::calculate_neighbors()
     // boundary edges, but the boundaries are calculated separately elsewhere.
 }
 
+Py::Object Triangulation::calculate_plane_coefficients(const Py::Tuple &args)
+{
+    _VERBOSE("Triangulation::calculate_plane_coefficients");
+    args.verify_length(1);
+
+    PyArrayObject* z = (PyArrayObject*)PyArray_ContiguousFromObject(
+                           args[0].ptr(), PyArray_DOUBLE, 1, 1);
+    if (z == 0 || PyArray_DIM(z,0) != PyArray_DIM(_x,0)) {
+        Py_XDECREF(z);
+        throw Py::ValueError(
+            "z array must have same length as triangulation x and y arrays");
+    }
+    const double* zs = (const double*)PyArray_DATA(z);
+
+    npy_intp dims[2] = {_ntri, 3};
+    PyArrayObject* planes_array = (PyArrayObject*)PyArray_SimpleNew(
+                                                      2, dims, PyArray_DOUBLE);
+    double* planes = (double*)PyArray_DATA(planes_array);
+    const int* tris = get_triangles_ptr();
+    const double* xs = (const double*)PyArray_DATA(_x);
+    const double* ys = (const double*)PyArray_DATA(_y);
+    for (int tri = 0; tri < _ntri; ++tri)
+    {
+        if (is_masked(tri))
+        {
+            *planes++ = 0.0;
+            *planes++ = 0.0;
+            *planes++ = 0.0;
+            tris += 3;
+        }
+        else
+        {
+            // Equation of plane for all points r on plane is r.normal = p
+            // where normal is vector normal to the plane, and p is a constant.
+            // Rewrite as r_x*normal_x + r_y*normal_y + r_z*normal_z = p
+            // and rearrange to give
+            // r_z = (-normal_x/normal_z)*r_x + (-normal_y/normal_z)*r_y +
+            //       p/normal_z
+            XYZ point0(xs[*tris], ys[*tris], zs[*tris]);
+            tris++;
+            XYZ point1(xs[*tris], ys[*tris], zs[*tris]);
+            tris++;
+            XYZ point2(xs[*tris], ys[*tris], zs[*tris]);
+            tris++;
+
+            XYZ normal = (point1 - point0).cross(point2 - point0);
+
+            if (normal.z == 0.0)
+            {
+                // Normal is in x-y plane which means triangle consists of
+                // colinear points.  Try to do the best we can by taking plane
+                // through longest side of triangle.
+                double length_sqr_01 = (point1 - point0).length_squared();
+                double length_sqr_12 = (point2 - point1).length_squared();
+                double length_sqr_20 = (point0 - point2).length_squared();
+                if (length_sqr_01 > length_sqr_12)
+                {
+                    if (length_sqr_01 > length_sqr_20)
+                        normal = normal.cross(point1 - point0);
+                    else
+                        normal = normal.cross(point0 - point2);
+                }
+                else
+                {
+                    if (length_sqr_12 > length_sqr_20)
+                        normal = normal.cross(point2 - point1);
+                    else
+                        normal = normal.cross(point0 - point2);
+                }
+
+                if (normal.z == 0.0)
+                {
+                    // The 3 triangle points have identical x and y!  The best
+                    // we can do here is take normal = (0,0,1) and for the
+                    // constant p take the mean of the 3 points' z-values.
+                    normal = XYZ(0.0, 0.0, 1.0);
+                    point0.z = (point0.z + point1.z + point2.z) / 3.0;
+                }
+            }
+
+            *planes++ = -normal.x / normal.z;           // x
+            *planes++ = -normal.y / normal.z;           // y
+            *planes++ = normal.dot(point0) / normal.z;  // constant
+        }
+    }
+
+    return Py::asObject((PyObject*)planes_array);
+}
+
 void Triangulation::correct_triangles()
 {
     int* triangles_ptr = (int*)PyArray_DATA(_triangles);
@@ -506,6 +600,9 @@ void Triangulation::init_type()
     behaviors().name("Triangulation");
     behaviors().doc("Triangulation");
 
+    add_varargs_method("calculate_plane_coefficients",
+                       &Triangulation::calculate_plane_coefficients,
+                       "calculate_plane_coefficients(z)");
     add_noargs_method("get_edges", &Triangulation::get_edges,
                       "get_edges()");
     add_noargs_method("get_neighbors", &Triangulation::get_neighbors,

lib/matplotlib/tri/_tri.h

@@ -115,6 +115,7 @@ struct XYZ
     XYZ(const double& x_, const double& y_, const double& z_);
     XYZ cross(const XYZ& other) const;
     double dot(const XYZ& other) const;
+    double length_squared() const;
     XYZ operator-(const XYZ& other) const;
     friend std::ostream& operator<<(std::ostream& os, const XYZ& xyz);
 
@@ -189,6 +190,13 @@ class Triangulation : public Py::PythonExtension<Triangulation>
 
     virtual ~Triangulation();
 
+    /* Calculate plane equation coefficients for all unmasked triangles from
+     * the point (x,y) coordinates and point z-array of shape (npoints) passed
+     * in via the args.  Returned array has shape (npoints,3) and allows
+     * z-value at (x,y) coordinates in triangle tri to be calculated using
+     *      z = array[tri,0]*x + array[tri,1]*y + array[tri,2]. */
+    Py::Object calculate_plane_coefficients(const Py::Tuple &args);
+
     // Return the boundaries collection, creating it if necessary.
     const Boundaries& get_boundaries() const;
 

lib/matplotlib/tri/triangulation.py

@@ -88,6 +88,16 @@ def __init__(self, x, y, triangles=None, mask=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).
+        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):
         if self._edges is None:

lib/matplotlib/tri/triinterpolate.py

@@ -0,0 +1,101 @@
+from __future__ import print_function
+from matplotlib.tri import Triangulation
+from matplotlib.tri.trifinder import TriFinder
+import numpy as np
+
+
+class TriInterpolator(object):
+    """
+    Abstract base class for classes used to perform interpolation on
+    triangular grids.
+
+    Derived classes implement __call__(x,y) where x,y are array_like point
+    coordinates of the same shape, and that returns a masked array of the same
+    shape containing the interpolated z-values.
+    """
+    def __init__(self, triangulation, z, trifinder=None):
+        if not isinstance(triangulation, Triangulation):
+            raise ValueError('Expected a Triangulation object')
+        self._triangulation = triangulation
+
+        self._z = np.asarray(z)
+        if self._z.shape != self._triangulation.x.shape:
+            raise ValueError('z array must have same length as triangulation x'
+                             ' and y arrays')
+
+        if trifinder is not None and not isinstance(trifinder, TriFinder):
+            raise ValueError('Expected a TriFinder object')
+        self._trifinder = trifinder or self._triangulation.get_trifinder()
+
+
+class LinearTriInterpolator(TriInterpolator):
+    """
+    A LinearTriInterpolator performs linear interpolation on a triangular grid.
+
+    Each triangle is represented by a plane so that an interpolated value at
+    point (x,y) lies on the plane of the triangle containing (x,y).
+    Interpolated values are therefore continuous across the triangulation, but
+    their first derivatives are discontinuous at edges between triangles.
+    """
+    def __init__(self, triangulation, z, trifinder=None):
+        """
+        *triangulation*: the :class:`~matplotlib.tri.Triangulation` to
+        interpolate over.
+
+        *z*: array_like of shape (npoints).
+          Array of values, defined at grid points, to interpolate between.
+
+        *trifinder*: optional :class:`~matplotlib.tri.TriFinder` object.
+          If this is not specified, the Triangulation's default TriFinder will
+          be used by calling :func:`matplotlib.tri.Triangulation.get_trifinder`.
+        """
+        TriInterpolator.__init__(self, triangulation, z, trifinder)
+
+        # Store plane coefficients for fast interpolation calculations.
+        self._plane_coefficients = \
+            self._triangulation.calculate_plane_coefficients(self._z)
+
+        # Store vectorized interpolation function, so can pass in arbitrarily
+        # shape arrays of x, y and tri and the _single_interp function is
+        # called in turn with scalar x, y and tri.
+        self._multi_interp = np.vectorize(self._single_interp,
+                                          otypes=[np.float])
+
+    def __call__(self, x, y):
+        """
+        Return a masked array containing linearly interpolated values at the
+        specified x,y points.
+
+        *x*, *y* are array_like x and y coordinates of the same shape and any
+        number of dimensions.
+
+        Returned masked array has the same shape as *x* and *y*; values
+        corresponding to (x,y) points outside of the triangulation are masked
+        out.
+        """
+        # Check arguments.
+        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 equal-shaped arrays")
+
+        # Indices of triangles containing x, y points, or -1 for no triangles.
+        tris = self._trifinder(x, y)
+
+        # Perform interpolation.
+        z = self._multi_interp(x, y, tris)
+
+        # Return masked array.
+        return np.ma.masked_invalid(z, copy=False)
+
+    def _single_interp(self, x, y, tri):
+        """
+        Return single interpolated value at specified (*x*, *y*) coordinates
+        within triangle with index *tri*.  Returns np.nan if tri == -1.
+        """
+        if tri == -1:
+            return np.nan
+        else:
+            return (self._plane_coefficients[tri,0] * x +
+                    self._plane_coefficients[tri,1] * y +
+                    self._plane_coefficients[tri,2])