小编典典

插值其角节点已知的3D曲面并用色图着色

python

我想构建实验数据的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()))

在此处输入图片说明


阅读 209

收藏
2021-01-16

共1个答案

小编典典

问题归结为如何在matplotlib中对曲面进行插值着色,即与Matlab的shading('interp')功能等效。简短的答案是:您不能。它本身不受支持,因此最好的办法是手工完成,这是迄今为止提供的解决方案的目标。

几年前,当我对Matlab感到沮丧时,我沿着这条路走了shading('interp'):它通过简单地在每个四边形上插入4个角颜色来工作,这意味着相邻的四边形上颜色梯度的方向可以不同。我想要的是每个色带都恰好在z轴上两个明确定义的值之间,相邻单元格之间没有视觉中断。

进行三角剖分绝对是正确的想法。但是,我不仅要精简网格并希望达到相邻三角形的颜色在视觉上无法区分的点(而不会达到首先出现伪像的点),我的方法是在三角剖分上计算轮廓带,然后以3D方式绘制它们。

当我第一次实现此功能时,matplotlib不支持在三角剖分上绘制轮廓。现在通过_tri.TriContourGenerator。如果这也提供了提取的多边形顶点的z值,那么我们就可以完成。不幸的是,他们不是在Python的水平接近,所以我们需要尝试通过比较器的输出来重建他们create_filled_contours()create_contours(),这在下面的代码完成:

import numpy as np
from mpl_toolkits.mplot3d import Axes3D
from mpl_toolkits.mplot3d.art3d import Poly3DCollection
import matplotlib.pyplot as plt
from matplotlib import _tri, tri, cm

def contour_bands_3d(x, y, z, nbands=20):
    # obtain the contouring engine on a triangulation
    TRI = tri.Triangulation(x, y)
    C = _tri.TriContourGenerator(TRI.get_cpp_triangulation(), z)

    # define the band breaks
    brks = np.linspace(z.min(), z.max(), nbands+1)

    # the contour lines
    lines = [C.create_contour(b) for b in brks]

    # the contour bands
    bands = [C.create_filled_contour(brks[i], brks[i+1]) for i in xrange(nbands)]

    # compare the x, y vertices of each band with the x, y vertices of the upper
    # contour line; if matching, z = z1, otherwise z = z0 (see text for caveats)
    eps = 1e-6
    verts = []
    for i in xrange(nbands):
        b = bands[i][0]
        l = lines[i+1][0]
        z0, z1 = brks[i:i+2]
        zi = np.array([z1 if (np.abs(bb - l) < eps).all(1).any() else z0 for bb in b])
        verts.append(np.c_[b, zi[:,None]])
    return brks, verts

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

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

# calculate the 3d contour bands
brks, verts = contour_bands_3d(x, -y, z)

cmap = cm.get_cmap('jet')
norm = plt.Normalize(vmax=abs(y).max(), vmin=-abs(y).max())

PC = Poly3DCollection(verts, cmap=cmap, norm=norm, edgecolors='none')
PC.set_array(brks[:-1])
ax.add_collection(PC)
ax.set_ylim((-1, 1))
plt.show()

结果如下:

带状轮廓膜

请注意,z值的重构并不完全正确,因为我们还需要检查ax,y顶点是否实际上是原始数据集的一部分,在这种情况下,必须采用其原始z值。但是,修改轮廓算法的C
++代码以跟踪z值会容易得多。这将是一个很小的变化,而试图涵盖Python中的所有情况无疑是一场噩梦。

关于效率,我们正在尝试在Python级别完成图形卡的工作,因此这将是可怕的。但是,所有都一样mplot3d。如果一个人需要一个性能的实现,我建议BandedContourFilter()VTK。这非常快,并且可以在Python中使用。

2021-01-16