如何从scipy.interpolate.Rbf创建的样条中计算任意值？

Question

我在3维空间(x, y, z)中有几个数据点，并已使用scipy.interpolate.Rbf对其进行插值。这给了我一个很好的代表我3D对象表面的样条线。我现在想确定几对x和y对，它们具有相同的任意z值。我想这样做是为了在任何给定的z值下计算3D对象的横截面。有人知道该怎么做吗？也许还有比使用scipy.interpolate.Rbf更好的方法了。

到目前为止，我已经通过使用matplotlib.pyplot绘制轮廓图并提取显示的线段来评估横截面。3D points and interpolated splinesegments extracted using a contour plot

Answer 1

我能够解决问题。我已经通过对x-y数据进行三角测量并用z平面切割三角形来计算了面积，我想计算出[z=z0）的横截面积。具体来说，我已经搜索了z值都在z0上方和下方的三角形。然后，我计算了这些三角形的边等于z0的边的x和y值。然后，我使用scipy.spatial.ConvexHull对相交的点进行排序。然后，使用鞋带公式即可确定面积。

我在此处附加了示例代码：

import numpy as np
from scipy import spatial
import matplotlib.pyplot as plt

# Generation of random test data
n = 500
x = np.random.random(n) 
y = np.random.random(n)
z = np.exp(-2*(x-.5)**2-4*(y-.5)**2)

z0 = .75

# Triangulation of the test data
triang= spatial.Delaunay(np.array([x, y]).T)


# Determine all triangles where not all points are above or below z0, i.e. the triangles that intersect z0
tri_inter=np.zeros_like(triang.simplices, dtype=np.int) # The triangles which intersect the plane at z0, filled below

i = 0
for tri in triang.simplices:    
    if ~np.all(z[tri] > z0) and ~np.all(z[tri] < z0):
        tri_inter[i,:] = tri
        i += 1
tri_inter = tri_inter[~np.all(tri_inter==0, axis=1)] # Remove all rows with only 0 

# The number of interpolated values for x and y has twice the length of the triangles
# Because each triangle intersects the plane at z0 twice
x_inter = np.zeros(tri_inter.shape[0]*2)
y_inter = np.zeros(tri_inter.shape[0]*2)

for j, tri in enumerate(tri_inter):

    # Determine which of the three points are above and which are below z0
    points_above = []
    points_below = []

    for i in tri:    
        if z[i] > z0:
            points_above.append(i)
        else:
            points_below.append(i)     

    # Calculate the intersections and put the values into x_inter and y_inter
    t = (z0-z[points_below[0]])/(z[points_above[0]]-z[points_below[0]])
    x_new = t * (x[points_above[0]]-x[points_below[0]]) + x[points_below[0]]
    y_new = t * (y[points_above[0]]-y[points_below[0]]) + y[points_below[0]]

    x_inter[j*2] = x_new
    y_inter[j*2] = y_new

    if len(points_above) > len(points_below): 
        t = (z0-z[points_below[0]])/(z[points_above[1]]-z[points_below[0]])
        x_new = t * (x[points_above[1]]-x[points_below[0]]) + x[points_below[0]]
        y_new = t * (y[points_above[1]]-y[points_below[0]]) + y[points_below[0]]

    else:
        t = (z0-z[points_below[1]])/(z[points_above[0]]-z[points_below[1]])
        x_new = t * (x[points_above[0]]-x[points_below[1]]) + x[points_below[1]]
        y_new = t * (y[points_above[0]]-y[points_below[1]]) + y[points_below[1]]

    x_inter[j*2+1] = x_new
    y_inter[j*2+1] = y_new

# sort points to calculate area
hull = spatial.ConvexHull(np.array([x_inter, y_inter]).T)

x_hull, y_hull = x_inter[hull.vertices], y_inter[hull.vertices]

# Calculation of are using the shoelace formula
area = 0.5*np.abs(np.dot(x_hull,np.roll(y_hull,1))-np.dot(y_hull,np.roll(x_hull,1)))

print('Area:', area)

plt.figure()
plt.plot(x_inter, y_inter, 'ro')
plt.plot(x_hull, y_hull, 'b--')
plt.triplot(x, y, triangles=tri_inter, color='k')
plt.show()