Interpolating a 3D surface known by its corner nod

问题:

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

and get the following plot:

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?

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

回答1:

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.

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.

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.

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

This is the result:

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.

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.

回答2:

Indeed it seems plot_trisurf should be perfect for this task! Additionally, you can make use of tri.UniformTriRefiner to get a Triangulation with smaller triangles:

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

Resulting in the following image:

Triangular grid refinement was only recently added to matplotlib so you will need version 1.3 to use it. Though if you would be stuck with version 1.2 you should also be able to use the source from Github directly, if you comment out the line import matplotlib.tri.triinterpolate and all of the refine_field method. Then you need to use the refine_triangulation method and use griddata to interpolate the new corresponding Z-values.

---

Edit: The above code uses cubic interpolation to determine the Z-values for the new triangles, but for linear interpolation you could substitute / add these lines:

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

Alternatively, to do the interpolation with 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')