Merge pull request matplotlib#1582 from ianthomas23/1521_linear_tri_interpolator

Linear tri interpolator

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

@@ -5,6 +5,17 @@ triangular grids
 :mod:`matplotlib.tri`
 =====================
 
-.. automodule:: matplotlib.tri
-   :members: Triangulation
+.. autoclass:: matplotlib.tri.Triangulation
+   :members:
 
+.. autoclass:: matplotlib.tri.TriFinder
+   :members:
+
+.. 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/event_handling/trifinder_event_demo.py

@@ -0,0 +1,59 @@
+"""
+Example showing the use of a TriFinder object.  As the mouse is moved over the
+triangulation, the triangle under the cursor is highlighted and the index of
+the triangle is displayed in the plot title.
+"""
+import matplotlib.pyplot as plt
+from matplotlib.tri import Triangulation
+from matplotlib.patches import Polygon
+import numpy as np
+import math
+
+
+def update_polygon(tri):
+    if tri == -1:
+        points = [0, 0, 0]
+    else:
+        points = triangulation.triangles[tri]
+    xs = triangulation.x[points]
+    ys = triangulation.y[points]
+    polygon.set_xy(zip(xs, ys))
+
+
+def motion_notify(event):
+    if event.inaxes is None:
+        tri = -1
+    else:
+        tri = trifinder(event.xdata, event.ydata)
+    update_polygon(tri)
+    plt.title('In triangle %i' % tri)
+    event.canvas.draw()
+
+
+# Create a Triangulation.
+n_angles = 16
+n_radii = 5
+min_radius = 0.25
+radii = np.linspace(min_radius, 0.95, n_radii)
+angles = np.linspace(0, 2*math.pi, n_angles, endpoint=False)
+angles = np.repeat(angles[..., np.newaxis], n_radii, axis=1)
+angles[:, 1::2] += math.pi / n_angles
+x = (radii*np.cos(angles)).flatten()
+y = (radii*np.sin(angles)).flatten()
+triangulation = Triangulation(x, y)
+xmid = x[triangulation.triangles].mean(axis=1)
+ymid = y[triangulation.triangles].mean(axis=1)
+mask = np.where(xmid*xmid + ymid*ymid < min_radius*min_radius, 1, 0)
+triangulation.set_mask(mask)
+
+# Use the triangulation's default TriFinder object.
+trifinder = triangulation.get_trifinder()
+
+# Setup plot and callbacks.
+plt.subplot(111, aspect='equal')
+plt.triplot(triangulation, 'bo-')
+polygon = Polygon([[0, 0], [0, 0]], facecolor='y')  # dummy data for xs,ys
+update_polygon(-1)
+plt.gca().add_patch(polygon)
+plt.gcf().canvas.mpl_connect('motion_notify_event', motion_notify)
+plt.show()

examples/pylab_examples/triinterp_demo.py

@@ -0,0 +1,35 @@
+"""
+Interpolation from triangular grid to quad grid.
+"""
+import matplotlib.pyplot as plt
+import matplotlib.tri as mtri
+import numpy as np
+
+
+# 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

@@ -131,3 +131,127 @@ def test_no_modify():
     tri = mtri.Triangulation(points[:,0], points[:,1], triangles)
     edges = tri.edges
     assert_array_equal(old_triangles, triangles)
+
+def test_trifinder():
+    # Test points within triangles of masked triangulation.
+    x, y = np.meshgrid(np.arange(4), np.arange(4))
+    x = x.ravel()
+    y = y.ravel()
+    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)
+    trifinder = triang.get_trifinder()
+
+    xs = [0.25, 1.25, 2.25, 3.25]
+    ys = [0.25, 1.25, 2.25, 3.25]
+    xs, ys = np.meshgrid(xs, ys)
+    xs = xs.ravel()
+    ys = ys.ravel()
+    tris = trifinder(xs, ys)
+    assert_array_equal(tris, [0, 2, 4, -1, 6, -1, 10, -1,
+                              12, 14, 16, -1, -1, -1, -1, -1])
+    tris = trifinder(xs-0.5, ys-0.5)
+    assert_array_equal(tris, [-1, -1, -1, -1, -1, 1, 3, 5,
+                              -1, 7, -1, 11, -1, 13, 15, 17])
+
+    # Test points exactly on boundary edges of masked triangulation.
+    xs = [0.5, 1.5, 2.5, 0.5, 1.5, 2.5, 1.5, 1.5, 0.0, 1.0, 2.0, 3.0]
+    ys = [0.0, 0.0, 0.0, 3.0, 3.0, 3.0, 1.0, 2.0, 1.5, 1.5, 1.5, 1.5]
+    tris = trifinder(xs, ys)
+    assert_array_equal(tris, [0, 2, 4, 13, 15, 17, 3, 14, 6, 7, 10, 11])
+
+    # Test points exactly on boundary corners of masked triangulation.
+    xs = [0.0, 3.0]
+    ys = [0.0, 3.0]
+    tris = trifinder(xs, ys)
+    assert_array_equal(tris, [0, 17])
+
+    # Test triangles with horizontal colinear points.  These are not valid
+    # triangulations, but we try to deal with the simplest violations.
+    delta = 0.0  # If +ve, triangulation is OK, if -ve triangulation invalid,
+                 # if zero have colinear points but should pass tests anyway.
+    x = [1.5, 0,  1,  2, 3, 1.5,   1.5]
+    y = [-1,  0,  0,  0, 0, delta, 1]
+    triangles = [[0, 2, 1], [0, 3, 2], [0, 4, 3], [1, 2, 5], [2, 3, 5],
+                 [3, 4, 5], [1, 5, 6], [4, 6, 5]]
+    triang = mtri.Triangulation(x, y, triangles)
+    trifinder = triang.get_trifinder()
+
+    xs = [-0.1, 0.4, 0.9, 1.4, 1.9, 2.4, 2.9]
+    ys = [-0.1, 0.1]
+    xs, ys = np.meshgrid(xs, ys)
+    tris = trifinder(xs, ys)
+    assert_array_equal(tris, [[-1, 0, 0, 1, 1, 2, -1],
+                              [-1, 6, 6, 6, 7, 7, -1]])
+
+    # Test triangles with vertical colinear points.  These are not valid
+    # triangulations, but we try to deal with the simplest violations.
+    delta = 0.0  # If +ve, triangulation is OK, if -ve triangulation invalid,
+                 # if zero have colinear points but should pass tests anyway.
+    x = [-1, -delta, 0,  0,  0, 0, 1]
+    y = [1.5, 1.5,   0,  1,  2, 3, 1.5]
+    triangles = [[0, 1, 2], [0, 1, 5], [1, 2, 3], [1, 3, 4], [1, 4, 5],
+                 [2, 6, 3], [3, 6, 4], [4, 6, 5]]
+    triang = mtri.Triangulation(x, y, triangles)
+    trifinder = triang.get_trifinder()
+
+    xs = [-0.1, 0.1]
+    ys = [-0.1, 0.4, 0.9, 1.4, 1.9, 2.4, 2.9]
+    xs, ys = np.meshgrid(xs, ys)
+    tris = trifinder(xs, ys)
+    assert_array_equal(tris, [[-1, -1], [0, 5], [0, 5], [0, 6], [1, 6], [1, 7],
+                              [-1, -1]])
+
+    # Test that changing triangulation by setting a mask causes the trifinder
+    # to be reinitialised.
+    x = [0, 1, 0, 1]
+    y = [0, 0, 1, 1]
+    triangles = [[0, 1, 2], [1, 3, 2]]
+    triang = mtri.Triangulation(x, y, triangles)
+    trifinder = triang.get_trifinder()
+
+    xs = [-0.2, 0.2, 0.8, 1.2]
+    ys = [0.5, 0.5, 0.5, 0.5]
+    tris = trifinder(xs, ys)
+    assert_array_equal(tris, [-1, 0, 1, -1])
+
+    triang.set_mask([1, 0])
+    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

@@ -5,5 +5,7 @@
 from __future__ import print_function
 from triangulation import *
 from tricontour import *
+from trifinder import *
+from triinterpolate import *
 from tripcolor import *
 from triplot import *