插入由其角节点已知的 3D 表面并使用颜色图对其进行着色 [英] Interpolating a 3D surface known by its corner nodes and coloring it with a colormap

问题描述

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

I want to construct a 3D representation of experimental data to track the deformation of a membrane. Experimentally, only the corner nodes are known. However I want to plot the deformaiton of the overall structure and this why I want to interpolate the membrane to enable a nice colormap of it. By searching around, I came almost close to it with the following code:

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.

The blue polygon is the surface known by its corner nodes and that I want to colormap. The colormapped surface is my best result so far. However, there are the black polygons near the top of the surface that are troubling me. I think it might be due to the fact that the surface doesn't fit the meshgrid and so the fourth corner is here a Nan.

是否有避免这些黑色三角形的解决方法，或者更好的颜色映射仅由其角节点知道的表面的方法?

Is there a workaround to avoid these black triangles or even better a better way of colormapping a surface known only by its corner nodes?

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

Here is the figure with the triangulation solution given in my first comment by using the following command

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()))

推荐答案

问题归结为如何在 matplotlib 中对表面进行插值着色，即等效于 Matlab 的 shading('interp')功能.简短的回答是:你不能.它不受本机支持，因此最好的方法是手工完成，这是目前提出的解决方案的目标.

The question boils down to how to do interpolated shading of a surface in matplotlib, i.e., the equivalent of Matlab's shading('interp') feature. The short answer is: You can't. It's not supported natively, so the best one can hope for is to do it by hand, which is what the solutions presented so far are aiming at.

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

I went down this road a few years ago, when I was getting frustrated with Matlab's shading('interp') as well: It works by simply interpolating the 4 corner colors on each quadrilateral, which means that the direction of the color gradient can be different on neighboring quadrilaterals. What I wanted was that each color band would be exactly between two well defined values on the z axis, with no visual breaks between neighboring cells.

进行三角测量绝对是正确的想法.但不是简单地细化网格并希望达到相邻三角形的颜色在视觉上无法区分的点(没有达到首先出现伪影的点)，我的方法是计算三角剖分上的轮廓带，然后在 3D 中绘制它们.

Working on a triangulation is definitely the right idea. But instead of simply refining the grid and hope to reach a point where the colors of neighboring triangles get visually indistinguishable (without reaching the point where artifacts appear first), my approach was to calculate the contour bands on the triangulation and then plot them in 3D.

当我第一次实现这个时，matplotlib 不支持在三角剖分上绘制轮廓.现在它通过 _tri.TriContourGenerator 实现.如果这也提供了提取的多边形顶点的 z 值，我们就完成了.不幸的是，它们在 Python 级别上是不可访问的，因此我们需要尝试通过比较 create_filled_contours() 和 create_contours() 的输出来重建它们，这是在以下代码:

When I first implemented this, matplotlib didn't support contouring on a triangulation. Now it does via _tri.TriContourGenerator. If this was providing the z values of the extracted polygon vertices as well, we would be done. Unfortunately, they are not accessible on the Python level, so we need to try to reconstruct them by comparing the outputs of create_filled_contours() and create_contours(), which is done in the following code:

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 值的重建并不完全正确，因为我们还需要检查 x、y 顶点是否实际上是原始数据集的一部分，在这种情况下，必须采用其原始 z 值.但是，修改轮廓算法的 C++ 代码以跟踪 z 值会容易得多.这将是一个很小的变化，而试图覆盖 Python 中的所有情况简直就是一场噩梦.

Note that the reconstruction of the z values is not fully correct, since we would also need to check if a x, y vertex is in fact part of the original data set, in which case its original z value must be taken. However, it would be much easier to modify the C++ code of the contouring algorithm to keep track of the z values. This would be a small change, while trying to cover all cases in Python is nothing short of a nightmare.

关于效率，嗯，我们正在尝试在 Python 级别上完成显卡的工作，所以它会很糟糕.但这与所有 mplot3d 相同.如果需要性能实现，我推荐 VTK 中的 BandedContourFilter().这工作得非常快，也可以从 Python 中使用.

Regarding efficiency, well, we are trying to do the job of a graphics card on the Python level, so it's going to be horrible. But that's the same with all of mplot3d. If one needs a performance implementation, I recommend BandedContourFilter() from VTK. This works blazingly fast and can be used from Python as well.

这篇关于插入由其角节点已知的 3D 表面并使用颜色图对其进行着色的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！