根据径向分量减少 3D 插值尺寸

Question

我有以下代表 3D 网格中的点的数据：

x = array([ 0, 0.08885313,  0.05077321,  0.05077321,  0.03807991,  0.03807991,
        0.03807991,  0.02538661,  0.02538661,  0.0126933 ,  0.0126933 ,
       -0.        , -0.        , -0.0126933 , -0.0126933 , -0.02538661,
       -0.02538661, -0.02538661, -0.03807991, -0.05077321, -0.05077321,
       -0.05077321, -0.06346652, -0.07615982, -0.11423973, -0.12693304,
       -0.13962634])
y = array([ 0, -0.15231964, -0.08885313, -0.17770625, -0.08885313, -0.10154643,
       -0.12693304, -0.07615982, -0.08885313, -0.07615982, -0.08885313,
       -0.07615982, -0.10154643, -0.08885313, -0.10154643, -0.07615982,
       -0.08885313, -0.17770625, -0.10154643, -0.08885313, -0.11423973,
       -0.12693304, -0.10154643, -0.08885313, -0.17770625, -0.12693304,
       -0.11423973])
z = array([ 0, 1.21839241, 0.78673339, 1.21839241, 0.70318648, 0.82850684,
       0.96078945, 0.64748854, 0.71014872, 0.63356406, 0.71711096,
       0.65445078, 0.77977115, 0.73103545, 0.77977115, 0.68926199,
       0.73799769, 1.19750569, 0.85635581, 0.84243133, 0.96078945,
       1.02344963, 0.93294048, 0.97471393, 1.24624138, 1.20446793,
       1.22535466])

我正在尝试将曲面拟合到这些点。我想拟合曲面，使其仅在数据点和原点之间延伸。

我尝试过以下代码：

data = np.array([x, y, z]).T 
  
# Define mathematical function for curve fitting 
def func(xy, a, b, c, d, e, f): 
    x, y = xy 
    return a + b*x + c*y + d*x**2 + e*y**2 + f*x*y 
  
# Perform curve fitting 
popt, pcov = curve_fit(func, (x, y), z) 
  
# Print optimized parameters 
print(popt) 
  
# Create 3D plot of the data points and the fitted curve 
fig = plt.figure() 
ax = fig.add_subplot(111, projection='3d') 
ax.scatter(x, y, z, color='blue') 
x_range = np.linspace(-0.2, 0.2, 50) 
y_range = np.linspace(-0.2, 0.2, 50) 
X, Y = np.meshgrid(x_range, y_range) 
Z = func((X, Y), *popt) 
ax.plot_surface(X, Y, Z, color='red', alpha=0.5) 

ax.set_xlim([-0.2, 0.2])  
ax.set_ylim([-0.2, 0.2])
ax.set_zlim([0,1.2])

plt.show()

这段代码给了我以下情节。我用蓝线绘制了我希望在其之间进行插值的位置。我如何添加两个不同的角度来限制我在径向方向上的绘图，以便插值位于这些线之间？

Answer 1

限制线似乎位于零点 (0, 0, 0) 的左侧和右侧，穿过 x 轴。这意味着具有 x 的最小值和最大值的两个点可用于定义这两条线。您可以使用 numpy.linalg.lstsq 找到它们的系数，就像这里的答案一样。线的 XY 投影可用于在绘制的表面上设置限制条件（使用 numpy.where）：

import matplotlib.pyplot as plt
import numpy as np
from scipy.optimize import curve_fit
from numpy import ones, vstack
from numpy.linalg import lstsq

x = np.array([ 0, 0.08885313,  0.05077321,  0.05077321,  0.03807991,  0.03807991,
        0.03807991,  0.02538661,  0.02538661,  0.0126933 ,  0.0126933 ,
       -0.        , -0.        , -0.0126933 , -0.0126933 , -0.02538661,
       -0.02538661, -0.02538661, -0.03807991, -0.05077321, -0.05077321,
       -0.05077321, -0.06346652, -0.07615982, -0.11423973, -0.12693304,
       -0.13962634])
y = np.array([ 0, -0.15231964, -0.08885313, -0.17770625, -0.08885313, -0.10154643,
       -0.12693304, -0.07615982, -0.08885313, -0.07615982, -0.08885313,
       -0.07615982, -0.10154643, -0.08885313, -0.10154643, -0.07615982,
       -0.08885313, -0.17770625, -0.10154643, -0.08885313, -0.11423973,
       -0.12693304, -0.10154643, -0.08885313, -0.17770625, -0.12693304,
       -0.11423973])
z = np.array([ 0, 1.21839241, 0.78673339, 1.21839241, 0.70318648, 0.82850684,
       0.96078945, 0.64748854, 0.71014872, 0.63356406, 0.71711096,
       0.65445078, 0.77977115, 0.73103545, 0.77977115, 0.68926199,
       0.73799769, 1.19750569, 0.85635581, 0.84243133, 0.96078945,
       1.02344963, 0.93294048, 0.97471393, 1.24624138, 1.20446793,
       1.22535466])

# the data set seems to be within two boundaries located to left and right of the (0, 0, 0) point across the x-axis
# find the maximum and minimum x values
i_max, = np.where(np.isclose(x, np.max(x)))
i_min, = np.where(np.isclose(x, np.min(x)))
x_max = x[i_max][0]
x_min = x[i_min][0]

# find the corresponding y values (i.e. y values of those points where x is maximum)
y_max = y[i_max][0]
y_min = y[i_min][0]

# function finding a 2D-line coefficients, based on this answer https://stackoverflow.com/a/21566184/3715182
def line_coeffs(points):
    x_coords, y_coords = zip(*points)
    A = vstack([x_coords, ones(len(x_coords))]).T
    # y = a*x + b
    a, b = lstsq(A, y_coords, rcond=None)[0]
    return (a, b)

# find coefficients of the two lines "limiting" all points left and right across the x-axis in the XY-plane
k1_max, k2_max = line_coeffs([(x_max, y_max), (0, 0)])
k1_min, k2_min = line_coeffs([(x_min, y_min), (0, 0)])
  
# Define mathematical function for curve fitting 
def func(xy, a, b, c, d, e, f):
    x, y = xy
    return a + b*x + c*y + d*x**2 + e*y**2 + f*x*y
  
# Perform curve fitting 
popt, pcov = curve_fit(func, (x, y), z)
  
# Print optimized parameters 
print(popt)
  
# Create 3D plot of the data points and the fitted curve 
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.scatter(x, y, z, color='blue')
x_range = np.linspace(-0.2, 0.2, 1000)
y_range = np.linspace(-0.2, 0.2, 1000)
X, Y = np.meshgrid(x_range, y_range)
Z = func((X, Y), *popt)
# limit the surface with a condition, forcing its XY-projections to be within the area limited by two lines
ax.plot_surface(np.where(X >= (Y - k2_min)/k1_min, np.where(X <= (Y - k2_max)/k1_max, X, np.nan), np.nan), Y, Z, color='red', alpha=0.5)

ax.set_xlim([-0.2, 0.2])
ax.set_ylim([-0.2, 0.2])
ax.set_zlim([0,1.2])

ax.set_xlabel('x')  
ax.set_ylabel('y')
ax.set_zlabel('z')

plt.show()

结果：