1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113
| from matplotlib.tri import Triangulation, TriAnalyzer, UniformTriRefiner import matplotlib.pyplot as plt import matplotlib.cm as cm import numpy as np
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))
n_test = 200
subdiv = 3
init_mask_frac = 0.0
min_circle_ratio = .01
random_gen = np.random.RandomState(seed=19680801) 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)
tri = Triangulation(x_test, y_test) ntri = tri.triangles.shape[0]
mask_init = np.zeros(ntri, dtype=bool) masked_tri = random_gen.randint(0, ntri, int(ntri * init_mask_frac)) mask_init[masked_tri] = True tri.set_mask(mask_init)
mask = TriAnalyzer(tri).get_flat_tri_mask(min_circle_ratio) tri.set_mask(mask)
refiner = UniformTriRefiner(tri) tri_refi, z_test_refi = refiner.refine_field(z_test, subdiv=subdiv)
z_expected = experiment_res(tri_refi.x, tri_refi.y)
flat_tri = Triangulation(x_test, y_test) flat_tri.set_mask(~mask)
plot_tri = True plot_masked_tri = True plot_refi_tri = False plot_expected = False
levels = np.arange(0., 1., 0.025) cmap = cm.get_cmap(name='Blues', lut=None)
fig, ax = plt.subplots() ax.set_aspect('equal') ax.set_title("Filtering a Delaunay mesh\n" + "(application to high-resolution tricontouring)")
ax.tricontour(tri_refi, z_test_refi, levels=levels, cmap=cmap, linewidths=[2.0, 0.5, 1.0, 0.5])
if plot_expected: ax.tricontour(tri_refi, z_expected, levels=levels, cmap=cmap, linestyles='--')
if plot_refi_tri: ax.triplot(tri_refi, color='0.97')
if plot_tri: ax.triplot(tri, color='0.7')
if plot_masked_tri: ax.triplot(flat_tri, color='red')
plt.show()
|