Skip to content

Commit cd45a80

Browse files
GBilloteyGBillotey
committed
Issue 1521: Cubic interpolation for triangular grids
1 parent a78d7e8 commit cd45a80

13 files changed

+3035
-122
lines changed

‎CHANGELOG

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1,3 +1,6 @@
1+
2013-02-25 Added classes CubicTriInterpolator, UniformTriRefiner, TriAnalyzer
2+
to matplotlib.tri module. - GBy
3+
14
2013-01-23 Add 'savefig.directory' to rcParams to remember and fill in the last
25
directory saved to for figure save dialogs - Martin Spacek
36

‎doc/api/tri_api.rst

Lines changed: 20 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -4,18 +4,32 @@ triangular grids
44

55
:mod:`matplotlib.tri`
66
=====================
7+
.. automodule:: matplotlib.tri
78

89
.. autoclass:: matplotlib.tri.Triangulation
9-
:members:
10+
:members:
1011

1112
.. autoclass:: matplotlib.tri.TriFinder
12-
:members:
1313

1414
.. autoclass:: matplotlib.tri.TrapezoidMapTriFinder
15-
:members: __call__
15+
:members: __call__
16+
:show-inheritance:
1617

1718
.. autoclass:: matplotlib.tri.TriInterpolator
18-
:members:
19-
19+
2020
.. autoclass:: matplotlib.tri.LinearTriInterpolator
21-
:members: __call__
21+
:members: __call__, gradient
22+
:show-inheritance:
23+
24+
.. autoclass:: matplotlib.tri.CubicTriInterpolator
25+
:members: __call__, gradient
26+
:show-inheritance:
27+
28+
.. autoclass:: matplotlib.tri.TriRefiner
29+
30+
.. autoclass:: matplotlib.tri.UniformTriRefiner
31+
:show-inheritance:
32+
:members:
33+
34+
.. autoclass:: matplotlib.tri.TriAnalyzer
35+
:members:

‎doc/users/whats_new.rst

Lines changed: 9 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -75,10 +75,17 @@ conversion (`mpl.rc('svg', fonttype='none')`).
7575

7676
Triangular grid interpolation
7777
-----------------------------
78-
Ian Thomas added classes to perform interpolation within triangular grids
79-
(:class:`~matplotlib.tri.LinearTriInterpolator`) and a utility class to find
78+
Geoffroy Billotey and Ian Thomas added classes to perform interpolation within
79+
triangular grids: (:class:`~matplotlib.tri.LinearTriInterpolator` and
80+
:class:`~matplotlib.tri.CubicTriInterpolator`) and a utility class to find
8081
the triangles in which points lie (
8182
:class:`~matplotlib.tri.TrapezoidMapTriFinder`).
83+
A helper class to perform mesh refinement and smooth contouring was also added
84+
(:class:`~matplotlib.tri.UniformTriRefiner`).
85+
Finally, a class implementing some basic tools for triangular mesh improvement
86+
was added (:class:`~matplotlib.tri.TriAnalyzer`).
87+
88+
.. plot:: mpl_examples/pylab_examples/tricontour_smooth_user.py
8289

8390
Left and right side axes titles
8491
-------------------------------

‎examples/pylab_examples/tricontour_smooth_delaunay.py

Lines changed: 133 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,133 @@
1+
"""
2+
Demonstrates high-resolution tricontouring of a random set of points ;
3+
a matplotlib.tri.TriAnalyzer is used to improve the plot quality.
4+
5+
The initial data points and triangular grid for this demo are:
6+
- a set of random points is instantiated, inside [-1, 1] x [-1, 1] square
7+
- A Delaunay triangulation of these points is then computed, of which a
8+
random subset of triangles is masked out by the user (based on
9+
*init_mask_frac* parameter). This simulates invalidated data.
10+
11+
The proposed generic procedure to obtain a high resolution contouring of such
12+
a data set is the following:
13+
1) Compute an extended mask with a matplotlib.tri.TriAnalyzer, which will
14+
exclude badly shaped (flat) triangles from the border of the
15+
triangulation. Apply the mask to the triangulation (using set_mask).
16+
2) Refine and interpolate the data using a
17+
matplotlib.tri.UniformTriRefiner.
18+
3) Plot the refined data with tricontour.
19+
20+
"""
21+
from matplotlib.tri import Triangulation, TriAnalyzer, UniformTriRefiner
22+
import matplotlib.pyplot as plt
23+
import matplotlib.cm as cm
24+
import numpy as np
25+
26+
27+
#-----------------------------------------------------------------------------
28+
# Analytical test function
29+
#-----------------------------------------------------------------------------
30+
def experiment_res(x, y):
31+
""" An analytic function representing experiment results """
32+
x = 2.*x
33+
r1 = np.sqrt((0.5-x)**2 + (0.5-y)**2)
34+
theta1 = np.arctan2(0.5-x, 0.5-y)
35+
r2 = np.sqrt((-x-0.2)**2 + (-y-0.2)**2)
36+
theta2 = np.arctan2(-x-0.2, -y-0.2)
37+
z = (4*(np.exp((r1/10)**2)-1)*30. * np.cos(3*theta1) +
38+
(np.exp((r2/10)**2)-1)*30. * np.cos(5*theta2) +
39+
2*(x**2 + y**2))
40+
return (np.max(z)-z)/(np.max(z)-np.min(z))
41+
42+
#-----------------------------------------------------------------------------
43+
# Generating the initial data test points and triangulation for the demo
44+
#-----------------------------------------------------------------------------
45+
# User parameters for data test points
46+
n_test = 200 # Number of test data points, tested from 3 to 5000 for subdiv=3
47+
48+
subdiv = 3 # Number of recursive subdivisions of the initial mesh for smooth
49+
# plots. Values >3 might result in a very high number of triangles
50+
# for the refine mesh: new triangles numbering = (4**subdiv)*ntri
51+
52+
init_mask_frac = 0.0 # Float > 0. adjusting the proportion of
53+
# (invalid) initial triangles which will be masked
54+
# out. Enter 0 for no mask.
55+
56+
min_circle_ratio = .01 # Minimum circle ratio - border triangles with circle
57+
# ratio below this will be masked if they touch a
58+
# border. Suggested value 0.01 ; Use -1 to keep
59+
# all triangles.
60+
61+
# Random points
62+
random_gen = np.random.mtrand.RandomState(seed=127260)
63+
x_test = random_gen.uniform(-1., 1., size=n_test)
64+
y_test = random_gen.uniform(-1., 1., size=n_test)
65+
z_test = experiment_res(x_test, y_test)
66+
67+
# meshing with Delaunay triangulation
68+
tri = Triangulation(x_test, y_test)
69+
ntri = tri.triangles.shape[0]
70+
71+
# Some invalid data are masked out
72+
mask_init = np.zeros(ntri, dtype=np.bool)
73+
masked_tri = random_gen.randint(0, ntri, int(ntri*init_mask_frac))
74+
mask_init[masked_tri] = True
75+
tri.set_mask(mask_init)
76+
77+
78+
#-----------------------------------------------------------------------------
79+
# Improving the triangulation before high-res plots: removing flat triangles
80+
#-----------------------------------------------------------------------------
81+
# masking badly shaped triangles at the border of the triangular mesh.
82+
mask = TriAnalyzer(tri).get_flat_tri_mask(min_circle_ratio)
83+
tri.set_mask(mask)
84+
85+
# refining the data
86+
refiner = UniformTriRefiner(tri)
87+
tri_refi, z_test_refi = refiner.refine_field(z_test, subdiv=subdiv)
88+
89+
# analytical 'results' for comparison
90+
z_expected = experiment_res(tri_refi.x, tri_refi.y)
91+
92+
# for the demo: loading the 'flat' triangles for plot
93+
flat_tri = Triangulation(x_test, y_test)
94+
flat_tri.set_mask(~mask)
95+
96+
97+
#-----------------------------------------------------------------------------
98+
# Now the plots
99+
#-----------------------------------------------------------------------------
100+
# User options for plots
101+
plot_tri = True # plot of the base triangulation
102+
plot_masked_tri = True # plot of the excessively flat excluded triangles
103+
plot_refi_tri = False # plot of the refined triangulation
104+
plot_expected = False # plot of the analytical function values for comparison
105+
106+
107+
# Graphical options for tricontouring
108+
levels = np.arange(0., 1., 0.025)
109+
cmap = cm.get_cmap(name='Blues', lut=None)
110+
111+
plt.figure()
112+
plt.gca().set_aspect('equal')
113+
plt.title("Filtering a Delaunay mesh\n" +
114+
"(application to high-resolution tricontouring)")
115+
116+
# 1) plot of the refined (computed) data countours:
117+
plt.tricontour(tri_refi, z_test_refi, levels=levels, cmap=cmap,
118+
linewidths=[2.0, 0.5, 1.0, 0.5])
119+
# 2) plot of the expected (analytical) data countours (dashed):
120+
if plot_expected:
121+
plt.tricontour(tri_refi, z_expected, levels=levels, cmap=cmap,
122+
linestyles='--')
123+
# 3) plot of the fine mesh on which interpolation was done:
124+
if plot_refi_tri:
125+
plt.triplot(tri_refi, color='0.97')
126+
# 4) plot of the initial 'coarse' mesh:
127+
if plot_tri:
128+
plt.triplot(tri, color='0.7')
129+
# 4) plot of the unvalidated triangles from naive Delaunay Triangulation:
130+
if plot_masked_tri:
131+
plt.triplot(flat_tri, color='red')
132+
133+
plt.show()

‎examples/pylab_examples/tricontour_smooth_user.py

Lines changed: 76 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,76 @@
1+
"""
2+
Demonstrates high-resolution tricontouring on user-defined triangular grids
3+
with matplotlib.tri.UniformTriRefiner
4+
"""
5+
from matplotlib.tri import Triangulation, UniformTriRefiner
6+
import matplotlib.pyplot as plt
7+
import matplotlib.cm as cm
8+
import numpy as np
9+
import math
10+
11+
12+
#-----------------------------------------------------------------------------
13+
# Analytical test function
14+
#-----------------------------------------------------------------------------
15+
def function_z(x, y):
16+
""" A function of 2 variables """
17+
r1 = np.sqrt((0.5-x)**2 + (0.5-y)**2)
18+
theta1 = np.arctan2(0.5-x, 0.5-y)
19+
r2 = np.sqrt((-x-0.2)**2 + (-y-0.2)**2)
20+
theta2 = np.arctan2(-x-0.2, -y-0.2)
21+
z = -(2*(np.exp((r1/10)**2)-1)*30. * np.cos(7.*theta1) +
22+
(np.exp((r2/10)**2)-1)*30. * np.cos(11.*theta2) +
23+
0.7*(x**2 + y**2))
24+
return (np.max(z)-z)/(np.max(z)-np.min(z))
25+
26+
#-----------------------------------------------------------------------------
27+
# Creating a Triangulation
28+
#-----------------------------------------------------------------------------
29+
# First create the x and y coordinates of the points.
30+
n_angles = 20
31+
n_radii = 10
32+
min_radius = 0.15
33+
radii = np.linspace(min_radius, 0.95, n_radii)
34+
35+
angles = np.linspace(0, 2*math.pi, n_angles, endpoint=False)
36+
angles = np.repeat(angles[..., np.newaxis], n_radii, axis=1)
37+
angles[:, 1::2] += math.pi/n_angles
38+
39+
x = (radii*np.cos(angles)).flatten()
40+
y = (radii*np.sin(angles)).flatten()
41+
z = function_z(x, y)
42+
43+
# Now create the Triangulation.
44+
# (Creating a Triangulation without specifying the triangles results in the
45+
# Delaunay triangulation of the points.)
46+
triang = Triangulation(x, y)
47+
48+
# Mask off unwanted triangles.
49+
xmid = x[triang.triangles].mean(axis=1)
50+
ymid = y[triang.triangles].mean(axis=1)
51+
mask = np.where(xmid*xmid + ymid*ymid < min_radius*min_radius, 1, 0)
52+
triang.set_mask(mask)
53+
54+
#-----------------------------------------------------------------------------
55+
# Refine data
56+
#-----------------------------------------------------------------------------
57+
refiner = UniformTriRefiner(triang)
58+
tri_refi, z_test_refi = refiner.refine_field(z, subdiv=3)
59+
60+
#-----------------------------------------------------------------------------
61+
# Plot the triangulation and the high-res iso-contours
62+
#-----------------------------------------------------------------------------
63+
plt.figure()
64+
plt.gca().set_aspect('equal')
65+
plt.triplot(triang, lw=0.5, color='white')
66+
67+
levels = np.arange(0., 1., 0.025)
68+
cmap = cm.get_cmap(name='terrain', lut=None)
69+
plt.tricontourf(tri_refi, z_test_refi, levels=levels, cmap=cmap)
70+
plt.tricontour(tri_refi, z_test_refi, levels=levels,
71+
colors=['0.25', '0.5', '0.5', '0.5', '0.5'],
72+
linewidths=[1.0, 0.5, 0.5, 0.5, 0.5])
73+
74+
plt.title("High-resolution tricontouring")
75+
76+
plt.show()

‎examples/pylab_examples/trigradient_demo.py

Lines changed: 81 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,81 @@
1+
"""
2+
Demonstrates computation of gradient with matplotlib.tri.CubicTriInterpolator.
3+
"""
4+
from matplotlib.tri import Triangulation, UniformTriRefiner,\
5+
CubicTriInterpolator
6+
import matplotlib.pyplot as plt
7+
import matplotlib.cm as cm
8+
import numpy as np
9+
import math
10+
11+
12+
#-----------------------------------------------------------------------------
13+
# Electrical potential of a dipole
14+
#-----------------------------------------------------------------------------
15+
def dipole_potential(x, y):
16+
""" The electric dipole potential V """
17+
r_sq = x**2 + y**2
18+
theta = np.arctan2(y, x)
19+
z = np.cos(theta)/r_sq
20+
return (np.max(z)-z) / (np.max(z)-np.min(z))
21+
22+
23+
#-----------------------------------------------------------------------------
24+
# Creating a Triangulation
25+
#-----------------------------------------------------------------------------
26+
# First create the x and y coordinates of the points.
27+
n_angles = 30
28+
n_radii = 10
29+
min_radius = 0.2
30+
radii = np.linspace(min_radius, 0.95, n_radii)
31+
32+
angles = np.linspace(0, 2*math.pi, n_angles, endpoint=False)
33+
angles = np.repeat(angles[..., np.newaxis], n_radii, axis=1)
34+
angles[:, 1::2] += math.pi/n_angles
35+
36+
x = (radii*np.cos(angles)).flatten()
37+
y = (radii*np.sin(angles)).flatten()
38+
V = dipole_potential(x, y)
39+
40+
# Create the Triangulation; no triangles specified so Delaunay triangulation
41+
# created.
42+
triang = Triangulation(x, y)
43+
44+
# Mask off unwanted triangles.
45+
xmid = x[triang.triangles].mean(axis=1)
46+
ymid = y[triang.triangles].mean(axis=1)
47+
mask = np.where(xmid*xmid + ymid*ymid < min_radius*min_radius, 1, 0)
48+
triang.set_mask(mask)
49+
50+
#-----------------------------------------------------------------------------
51+
# Refine data - interpolates the electrical potential V
52+
#-----------------------------------------------------------------------------
53+
refiner = UniformTriRefiner(triang)
54+
tri_refi, z_test_refi = refiner.refine_field(V, subdiv=3)
55+
56+
#-----------------------------------------------------------------------------
57+
# Computes the electrical field (Ex, Ey) as gradient of electrical potential
58+
#-----------------------------------------------------------------------------
59+
tci = CubicTriInterpolator(triang, -V)
60+
# Gradient requested here at the mesh nodes but could be anywhere else:
61+
(Ex, Ey) = tci.gradient(triang.x, triang.y)
62+
E_norm = np.sqrt(Ex**2 + Ey**2)
63+
64+
#-----------------------------------------------------------------------------
65+
# Plot the triangulation, the potential iso-contours and the vector field
66+
#-----------------------------------------------------------------------------
67+
plt.figure()
68+
plt.gca().set_aspect('equal')
69+
plt.triplot(triang, color='0.8')
70+
71+
levels = np.arange(0., 1., 0.01)
72+
cmap = cm.get_cmap(name='hot', lut=None)
73+
plt.tricontour(tri_refi, z_test_refi, levels=levels, cmap=cmap,
74+
linewidths=[2.0, 1.0, 1.0, 1.0])
75+
# Plots direction of the electrical vector field
76+
plt.quiver(triang.x, triang.y, Ex/E_norm, Ey/E_norm,
77+
units='xy', scale=10., zorder=3, color='blue',
78+
width=0.007, headwidth=3., headlength=4.)
79+
80+
plt.title('Gradient plot: an electrical dipole')
81+
plt.show()

0 commit comments

Comments
 (0)