Cubic interpolation for triangular grids

CHANGELOG

@@ -1,3 +1,6 @@
+2013-02-25 Added classes CubicTriInterpolator, UniformTriRefiner, TriAnalyzer
+           to matplotlib.tri module. - GBy
+
 2013-01-23 Add 'savefig.directory' to rcParams to remember and fill in the last
            directory saved to for figure save dialogs - Martin Spacek
 

doc/api/tri_api.rst

@@ -4,18 +4,32 @@ triangular grids
 
 :mod:`matplotlib.tri`
 =====================
+.. automodule:: matplotlib.tri
 
 .. autoclass:: matplotlib.tri.Triangulation
-   :members:
+    :members:
 
 .. autoclass:: matplotlib.tri.TriFinder
-   :members:
 
 .. autoclass:: matplotlib.tri.TrapezoidMapTriFinder
-   :members: __call__
+    :members: __call__
+    :show-inheritance:
 
 .. autoclass:: matplotlib.tri.TriInterpolator
-   :members:
-
+
 .. autoclass:: matplotlib.tri.LinearTriInterpolator
-   :members: __call__
+    :members: __call__, gradient
+    :show-inheritance:
+
+.. autoclass:: matplotlib.tri.CubicTriInterpolator
+    :members: __call__, gradient
+    :show-inheritance:
+
+.. autoclass:: matplotlib.tri.TriRefiner
+
+.. autoclass:: matplotlib.tri.UniformTriRefiner
+    :show-inheritance:
+    :members: 
+
+.. autoclass:: matplotlib.tri.TriAnalyzer
+    :members: 

doc/users/whats_new.rst

@@ -75,10 +75,17 @@ 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
+Geoffroy Billotey and Ian Thomas added classes to perform interpolation within
+triangular grids: (:class:`~matplotlib.tri.LinearTriInterpolator` and
+:class:`~matplotlib.tri.CubicTriInterpolator`) and a utility class to find
 the triangles in which points lie (
 :class:`~matplotlib.tri.TrapezoidMapTriFinder`).
+A helper class to perform mesh refinement and smooth contouring was also added
+(:class:`~matplotlib.tri.UniformTriRefiner`).
+Finally, a class implementing some basic tools for triangular mesh improvement
+was added (:class:`~matplotlib.tri.TriAnalyzer`).
+
+.. plot:: mpl_examples/pylab_examples/tricontour_smooth_user.py
 
 Left and right side axes titles
 -------------------------------

examples/pylab_examples/tricontour_smooth_delaunay.py

@@ -0,0 +1,133 @@
+"""
+Demonstrates high-resolution tricontouring of a random set of points ;
+a matplotlib.tri.TriAnalyzer is used to improve the plot quality.
+
+The initial data points and triangular grid for this demo are:
+    - a set of random points is instantiated, inside [-1, 1] x [-1, 1] square
+    - A Delaunay triangulation of these points is then computed, of which a
+    random subset of triangles is masked out by the user (based on
+    *init_mask_frac* parameter). This simulates invalidated data.
+
+The proposed generic procedure to obtain a high resolution contouring of such
+a data set is the following:
+    1) Compute an extended mask with a matplotlib.tri.TriAnalyzer, which will
+    exclude badly shaped (flat) triangles from the border of the
+    triangulation. Apply the mask to the triangulation (using set_mask).
+    2) Refine and interpolate the data using a
+    matplotlib.tri.UniformTriRefiner.
+    3) Plot the refined data with tricontour.
+
+"""
+from matplotlib.tri import Triangulation, TriAnalyzer, UniformTriRefiner
+import matplotlib.pyplot as plt
+import matplotlib.cm as cm
+import numpy as np
+
+
+#-----------------------------------------------------------------------------
+# Analytical test function
+#-----------------------------------------------------------------------------
+def experiment_res(x, y):
+    """ An analytic function representing experiment results """
+    x = 2.*x
+    r1 = np.sqrt((0.5-x)**2 + (0.5-y)**2)
+    theta1 = np.arctan2(0.5-x, 0.5-y)
+    r2 = np.sqrt((-x-0.2)**2 + (-y-0.2)**2)
+    theta2 = np.arctan2(-x-0.2, -y-0.2)
+    z = (4*(np.exp((r1/10)**2)-1)*30. * np.cos(3*theta1) +
+         (np.exp((r2/10)**2)-1)*30. * np.cos(5*theta2) +
+         2*(x**2 + y**2))
+    return (np.max(z)-z)/(np.max(z)-np.min(z))
+
+#-----------------------------------------------------------------------------
+# Generating the initial data test points and triangulation for the demo
+#-----------------------------------------------------------------------------
+# User parameters for data test points
+n_test = 200  # Number of test data points, tested from 3 to 5000 for subdiv=3
+
+subdiv = 3  # Number of recursive subdivisions of the initial mesh for smooth
+            # plots. Values >3 might result in a very high number of triangles
+            # for the refine mesh: new triangles numbering = (4**subdiv)*ntri
+
+init_mask_frac = 0.0    # Float > 0. adjusting the proportion of
+                        # (invalid) initial triangles which will be masked
+                        # out. Enter 0 for no mask.
+
+min_circle_ratio = .01  # Minimum circle ratio - border triangles with circle
+                        # ratio below this will be masked if they touch a
+                        # border. Suggested value 0.01 ; Use -1 to keep
+                        # all triangles.
+
+# Random points
+random_gen = np.random.mtrand.RandomState(seed=127260)
+x_test = random_gen.uniform(-1., 1., size=n_test)
+y_test = random_gen.uniform(-1., 1., size=n_test)
+z_test = experiment_res(x_test, y_test)
+
+# meshing with Delaunay triangulation
+tri = Triangulation(x_test, y_test)
+ntri = tri.triangles.shape[0]
+
+# Some invalid data are masked out
+mask_init = np.zeros(ntri, dtype=np.bool)
+masked_tri = random_gen.randint(0, ntri, int(ntri*init_mask_frac))
+mask_init[masked_tri] = True
+tri.set_mask(mask_init)
+
+
+#-----------------------------------------------------------------------------
+# Improving the triangulation before high-res plots: removing flat triangles
+#-----------------------------------------------------------------------------
+# masking badly shaped triangles at the border of the triangular mesh.
+mask = TriAnalyzer(tri).get_flat_tri_mask(min_circle_ratio)
+tri.set_mask(mask)
+
+# refining the data
+refiner = UniformTriRefiner(tri)
+tri_refi, z_test_refi = refiner.refine_field(z_test, subdiv=subdiv)
+
+# analytical 'results' for comparison
+z_expected = experiment_res(tri_refi.x, tri_refi.y)
+
+# for the demo: loading the 'flat' triangles for plot
+flat_tri = Triangulation(x_test, y_test)
+flat_tri.set_mask(~mask)
+
+
+#-----------------------------------------------------------------------------
+# Now the plots
+#-----------------------------------------------------------------------------
+# User options for plots
+plot_tri = True          # plot of the base triangulation
+plot_masked_tri = True   # plot of the excessively flat excluded triangles
+plot_refi_tri = False    # plot of the refined triangulation
+plot_expected = False    # plot of the analytical function values for comparison
+
+
+# Graphical options for tricontouring
+levels = np.arange(0., 1., 0.025)
+cmap = cm.get_cmap(name='Blues', lut=None)
+
+plt.figure()
+plt.gca().set_aspect('equal')
+plt.title("Filtering a Delaunay mesh\n" +
+          "(application to high-resolution tricontouring)")
+
+# 1) plot of the refined (computed) data countours:
+plt.tricontour(tri_refi, z_test_refi, levels=levels, cmap=cmap,
+               linewidths=[2.0, 0.5, 1.0, 0.5])
+# 2) plot of the expected (analytical) data countours (dashed):
+if plot_expected:
+    plt.tricontour(tri_refi, z_expected, levels=levels, cmap=cmap,
+                   linestyles='--')
+# 3) plot of the fine mesh on which interpolation was done:
+if plot_refi_tri:
+    plt.triplot(tri_refi, color='0.97')
+# 4) plot of the initial 'coarse' mesh:
+if plot_tri:
+    plt.triplot(tri, color='0.7')
+# 4) plot of the unvalidated triangles from naive Delaunay Triangulation:
+if plot_masked_tri:
+    plt.triplot(flat_tri, color='red')
+
+plt.show()

examples/pylab_examples/tricontour_smooth_user.py

@@ -0,0 +1,76 @@
+"""
+Demonstrates high-resolution tricontouring on user-defined triangular grids
+with matplotlib.tri.UniformTriRefiner
+"""
+from matplotlib.tri import Triangulation, UniformTriRefiner
+import matplotlib.pyplot as plt
+import matplotlib.cm as cm
+import numpy as np
+import math
+
+
+#-----------------------------------------------------------------------------
+# Analytical test function
+#-----------------------------------------------------------------------------
+def function_z(x, y):
+    """ A function of 2 variables """
+    r1 = np.sqrt((0.5-x)**2 + (0.5-y)**2)
+    theta1 = np.arctan2(0.5-x, 0.5-y)
+    r2 = np.sqrt((-x-0.2)**2 + (-y-0.2)**2)
+    theta2 = np.arctan2(-x-0.2, -y-0.2)
+    z = -(2*(np.exp((r1/10)**2)-1)*30. * np.cos(7.*theta1) +
+          (np.exp((r2/10)**2)-1)*30. * np.cos(11.*theta2) +
+          0.7*(x**2 + y**2))
+    return (np.max(z)-z)/(np.max(z)-np.min(z))
+
+#-----------------------------------------------------------------------------
+# Creating a Triangulation 
+#-----------------------------------------------------------------------------
+# First create the x and y coordinates of the points.
+n_angles = 20
+n_radii = 10
+min_radius = 0.15
+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()
+z = function_z(x, y)
+
+# Now create the Triangulation.
+# (Creating a Triangulation without specifying the triangles results in the
+# Delaunay triangulation of the points.)
+triang = Triangulation(x, y)
+
+# Mask off unwanted triangles.
+xmid = x[triang.triangles].mean(axis=1)
+ymid = y[triang.triangles].mean(axis=1)
+mask = np.where(xmid*xmid + ymid*ymid < min_radius*min_radius, 1, 0)
+triang.set_mask(mask)
+
+#-----------------------------------------------------------------------------
+# Refine data
+#-----------------------------------------------------------------------------
+refiner = UniformTriRefiner(triang)
+tri_refi, z_test_refi = refiner.refine_field(z, subdiv=3)
+
+#-----------------------------------------------------------------------------
+# Plot the triangulation and the high-res iso-contours
+#-----------------------------------------------------------------------------
+plt.figure()
+plt.gca().set_aspect('equal')
+plt.triplot(triang, lw=0.5, color='white')
+
+levels = np.arange(0., 1., 0.025)
+cmap = cm.get_cmap(name='terrain', lut=None)
+plt.tricontourf(tri_refi, z_test_refi, levels=levels, cmap=cmap)
+plt.tricontour(tri_refi, z_test_refi, levels=levels,
+               colors=['0.25', '0.5', '0.5', '0.5', '0.5'],
+               linewidths=[1.0, 0.5, 0.5, 0.5, 0.5])
+
+plt.title("High-resolution tricontouring")
+
+plt.show()

examples/pylab_examples/trigradient_demo.py

@@ -0,0 +1,81 @@
+"""
+Demonstrates computation of gradient with matplotlib.tri.CubicTriInterpolator.
+"""
+from matplotlib.tri import Triangulation, UniformTriRefiner,\
+    CubicTriInterpolator
+import matplotlib.pyplot as plt
+import matplotlib.cm as cm
+import numpy as np
+import math
+
+
+#-----------------------------------------------------------------------------
+# Electrical potential of a dipole
+#-----------------------------------------------------------------------------
+def dipole_potential(x, y):
+    """ The electric dipole potential V """
+    r_sq = x**2 + y**2
+    theta = np.arctan2(y, x)
+    z = np.cos(theta)/r_sq
+    return (np.max(z)-z) / (np.max(z)-np.min(z))
+
+
+#-----------------------------------------------------------------------------
+# Creating a Triangulation
+#-----------------------------------------------------------------------------
+# First create the x and y coordinates of the points.
+n_angles = 30
+n_radii = 10
+min_radius = 0.2
+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()
+V = dipole_potential(x, y)
+
+# Create the Triangulation; no triangles specified so Delaunay triangulation
+# created.
+triang = Triangulation(x, y)
+
+# Mask off unwanted triangles.
+xmid = x[triang.triangles].mean(axis=1)
+ymid = y[triang.triangles].mean(axis=1)
+mask = np.where(xmid*xmid + ymid*ymid < min_radius*min_radius, 1, 0)
+triang.set_mask(mask)
+
+#-----------------------------------------------------------------------------
+# Refine data - interpolates the electrical potential V
+#-----------------------------------------------------------------------------
+refiner = UniformTriRefiner(triang)
+tri_refi, z_test_refi = refiner.refine_field(V, subdiv=3)
+
+#-----------------------------------------------------------------------------
+# Computes the electrical field (Ex, Ey) as gradient of electrical potential
+#-----------------------------------------------------------------------------
+tci = CubicTriInterpolator(triang, -V)
+# Gradient requested here at the mesh nodes but could be anywhere else:
+(Ex, Ey) = tci.gradient(triang.x, triang.y)
+E_norm = np.sqrt(Ex**2 + Ey**2)
+
+#-----------------------------------------------------------------------------
+# Plot the triangulation, the potential iso-contours and the vector field
+#-----------------------------------------------------------------------------
+plt.figure()
+plt.gca().set_aspect('equal')
+plt.triplot(triang, color='0.8')
+
+levels = np.arange(0., 1., 0.01)
+cmap = cm.get_cmap(name='hot', lut=None)
+plt.tricontour(tri_refi, z_test_refi, levels=levels, cmap=cmap,
+               linewidths=[2.0, 1.0, 1.0, 1.0])
+# Plots direction of the electrical vector field
+plt.quiver(triang.x, triang.y, Ex/E_norm, Ey/E_norm,
+           units='xy', scale=10., zorder=3, color='blue',
+           width=0.007, headwidth=3., headlength=4.)
+
+plt.title('Gradient plot: an electrical dipole')
+plt.show()