python - 插值由其角节点已知的 3D 表面并使用颜色图对其着色

Question

我想构建实验数据的 3D 表示来跟踪膜的变形。实验上，只有角节点是已知的。但是我想绘制整体结构的变形，这就是为什么我想插入膜以启用它的漂亮颜色图。通过四处搜索，我使用以下代码几乎接近它:

import numpy
from mpl_toolkits.mplot3d import Axes3D
from mpl_toolkits.mplot3d.art3d import Poly3DCollection
import matplotlib.pyplot as plt
from matplotlib import cm
from scipy.interpolate import griddata

x=numpy.array([0, 0, 1, 1])
y=numpy.array([0.5, 0.75, 1, 0.5])
z=numpy.array([0, 0.5, 1,0])

fig = plt.figure()
ax = Axes3D(fig)
verts = [zip(x, y, z)]
PC = Poly3DCollection(verts)
ax.add_collection3d(PC)

xi = numpy.linspace(x.min(),x.max(),20)
yi = numpy.linspace(y.min(),y.max(),20)
zi = griddata((x,y),z, (xi[None,:], yi[:,None]), method='linear')
xig, yig = numpy.meshgrid(xi, -yi)
ax.plot_surface(xig, yig, zi, rstride=1, cstride=1,  linewidth=0,cmap=plt.cm.jet,norm=plt.Normalize(vmax=abs(yi).max(), vmin=-abs(yi).max()))
plt.show()

得到如下图:

蓝色多边形是其角节点已知的表面，我想对其进行颜色映射。彩色贴图表面是我迄今为止最好的结果。然而，表面顶部附近的黑色多边形让我很困扰。我认为这可能是因为表面不适合网格，所以第四个角在这里是 Nan。

是否有避免这些黑色三角形的解决方法，或者更好的方法来对仅由角节点识别的表面进行颜色映射？

编辑:这是我在第一条评论中使用以下命令给出的三角测量解决方案的图

triang = tri.Triangulation(x, y)
ax.plot_trisurf(x, y, z, triangles=triang.triangles, cmap=cm.jet,norm=plt.Normalize(vmax=abs(yi).max(), vmin=-abs(yi).max()))

Answer 1

的确，plot_trisurf 似乎非常适合这项任务!此外，您可以利用 tri.UniformTriRefiner 获得具有较小三角形的 Triangulation:

import numpy
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
from matplotlib import tri, cm

x = numpy.array([0, 0, 1, 1])
y = numpy.array([0.5, 0.75, 1, 0.5])
z = numpy.array([0, 0.5, 1, 0])

triang = tri.Triangulation(x, y)
refiner = tri.UniformTriRefiner(triang)
new, new_z = refiner.refine_field(z, subdiv=4)

norm = plt.Normalize(vmax=abs(y).max(), vmin=-abs(y).max())
kwargs = dict(triangles=new.triangles, cmap=cm.jet, norm=norm, linewidth=0.2)

fig = plt.figure()
ax = Axes3D(fig)
pt = ax.plot_trisurf(new.x, new.y, new_z, **kwargs)
plt.show()

结果如下图:

三角网格细化最近才添加到 matplotlib 中，因此您需要 1.3 版才能使用它。虽然如果你会被 1.2 版卡住，你也应该能够使用 source from Github直接，如果你注释掉 import matplotlib.tri.triinterpolate 行和所有 refine_field 方法。然后，您需要使用 refine_triangulation 方法并使用 griddata 对新的对应 Z 值进行插值。

---

编辑:上面的代码使用三次插值来确定新三角形的 Z 值，但对于线性插值，您可以替换/添加这些行:

interpolator = tri.LinearTriInterpolator(triang, z)
new, new_z = refiner.refine_field(z, interpolator, subdiv=4)

或者，使用 scipy.interpolate.griddata 进行插值:

from scipy.interpolate import griddata

new = refiner.refine_triangulation(subdiv = 4)
new_z = griddata((x,y),z, (new.x, new.y), method='linear')