python - 如何从使用 scipy.interpolate.Rbf 创建的样条曲线计算任意值？

Question

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

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

Answer 1

我能够解决这个问题。我通过对 x-y 数据进行三角测量并用 z 平面切割三角形来计算面积，我想计算 (z=z0) 的横截面积。具体来说，我搜索了 z 值既高于又低于 z0 的三角形。然后我计算了这些三角形边的 x 和 y 值，其中边等于 z0。然后我使用 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()