In [1]:
>>> nx, ny = (10, 10)
>>> x = np.linspace(0, 1, nx)
>>> y = np.linspace(0, 1, ny)
>>> xv, yv = np.meshgrid(x, y)
In [2]:
points = np.array((xv.flatten(), yv.flatten())).transpose()
In [3]:
>>> from scipy.spatial import Delaunay
>>> tri = Delaunay(points)
In [4]:
>>> import matplotlib.pyplot as plt
>>> plt.triplot(points[:,0], points[:,1], tri.vertices.copy())
>>> plt.plot(points[:,0], points[:,1], 'o')
>>> plt.show()
In [5]:
indices, indptr = tri.vertex_neighbor_vertices
k=0
indptr[indices[k]:indices[k+1]]
Out[5]:
In [6]:
edges = []
for k in xrange(len(points)):
for p in indptr[indices[k]:indices[k+1]]:
if p > k:
edges.append([k, p])
print edges
In [9]:
tri.points
Out[9]:
In [4]:
from scipy.spatial import Delaunay
import numpy as np
def refine_contour_plot(points, values):
"""
Take a set of `points` and `values` used to plot a contour plot and returns new points to calculate in order of how the value difference between connecting points is reduced.
>>> refine_contour_plot([[0, 0], [1, 0], [0, 1]], [0., 1., 10.])
[[0., 0.5], [0.5, 0.5], [0.5, 0.]]
:Parameters:
- `points`: Set of (x, y) points.
- `values`: Values at the points
:Returns:
- `new_points`: set of points to recalculate and improve the contour plot
"""
tri = Delaunay(points)
indices, indptr = tri.vertex_neighbor_vertices
edges = []
for i in xrange(len(points)):
for j in indptr[indices[i]:indices[i + 1]]:
if j > i:
edges.append([i, j])
edges = np.array(edges)
vertexValues = np.array(values)[edges]
edgeValue = abs(vertexValues[:,0] - vertexValues[:,1])
return ((tri.points[edges[:,0],:] + tri.points[edges[:,1],:]) / 2)[np.argsort(edgeValue)][::-1]
In [42]:
refine_contour_plot([[0, 0], [1, 0], [0, 1]], [0., 1., 10.])
Out[42]:
In [1]:
def func(x, y):
return x**2 * y**2
x = np.linspace(0, 1, 1000)
y = np.linspace(0, 1, 1000)
xv, yv = np.meshgrid(x, y)
answer = func(xv, yv)
plt.contourf(x, y, answer)
Out[1]:
In [2]:
from matplotlib.mlab import griddata
def interpolate_values(points, values):
xi = np.linspace(0, 1, 1000)
yi = np.linspace(0, 1, 1000)
x = points[:,0]
y = points[:,1]
return xi, yi, griddata(x, y, values, xi, yi)
In [17]:
points = np.array([[0., 0.], [1., 0.], [0., 1.], [1., 1.]])
values = func(points[:,0], points[:,1])
N = 100
for refine_count in xrange(50):
xi, yi, valuesi = interpolate_values(points, values)
#diff = ((answer - valuesi) / np.maximum(1e-10, answer)).flatten()
diff = (answer - valuesi).flatten()
error1 = max(abs(diff))
error2 = np.sqrt(np.sum(diff**2))
print 'refinement level:',refine_count
print 'norm inf:',error1
print 'norm 2:',error2
tri = Delaunay(points)
fig = plt.figure(figsize=(10,3))
ax1 = fig.add_subplot(1, 3, 1)
ax1.contourf(xi, yi, valuesi)
ax2 = fig.add_subplot(1, 3, 3)
ax2.triplot(points[:,0], points[:,1], tri.vertices.copy())
ax3 = fig.add_subplot(1, 3, 2)
ax3.contourf(xi, yi, answer)
new_points = refine_contour_plot(points, values)
if len(new_points) >= N:
new_points = new_points[:N]
points = np.concatenate((points, new_points))
values = func(points[:,0], points[:,1])
plt.show()
In [ ]:
In [ ]:
In [ ]: