如何在 Python 中快速估计点和双三次样条曲面之间的距离？

Question

如何在 Python 中快速估计点与双三次样条曲面之间的距离？是否有我可以在 SciPy、NumPy 或其他包中利用的现有解决方案？

我有一个由双三次插值定义的表面：

import numpy as np
import scipy.interpolate

# Define regular grid surface
xmin,xmax,ymin,ymax = 25, 125, -50, 50
x = np.linspace(xmin,xmax, 201)
y = np.linspace(ymin,ymax, 201)
xx, yy = np.meshgrid(x, y)
z_ideal = ( xx**2 + yy**2 ) / 400
z_ideal += z_ideal + np.random.uniform(-0.5, 0.5, z_ideal.shape)
s_ideal = scipy.interpolate.interp2d(x, y, z_ideal, kind='cubic')   

我已经得到了该表面的一些测量点：

# Fake some measured points on the surface
z_measured = z_ideal + np.random.uniform(-0.1, 0.1, z_ideal.shape)
s_measured = scipy.interpolate.interp2d(x, y, z_measured, kind='cubic')
p_x = np.random.uniform(xmin,xmax,10000)
p_y = np.random.uniform(ymin,ymax,10000)
p_z = s_measured( p_x, p_y )

我想找到表面s_ideal 上与p 中每个点最近的点。一般情况下，对于变化很大的样条曲线可能有多种解决方案，因此我将问题限制在已知在点沿 z 的投影附近只有一个解决方案的表面上。 这不是极少数的测量或表面定义点，所以我想优化速度，即使牺牲精度到可能1E-5。

想到的方法是使用梯度下降法，对每个测量点做类似p:

1. 使用pt = [p_x, p_y, p_z]作为初始测试点，其中p_z = s_ideal(pt)
2. 计算斜率（切线）矢量m = [ m_x, m_y ]pt
3. 计算从pt到p的向量r：r = p - pt
4. 如果r 和m 之间的角度theta 在90 度的某个阈值内，那么pt 就是最后一点。
5. 否则，将pt 更新为：

r_len = numpy.linalg.norm(r)
dx = r_len * m_x
dy = r_len * m_y
if theta > 90:
    pt = [ p_x + dx, p_y + dy ]
else:
    pt = [ p_x - dx, p_y - dy ]

我发现this 建议一种方法可以快速产生一维情况下的高精度结果，但它是单一维度的，我可能很难转换为二维。

Answer 1

该问题旨在最小化三维表面S(x,y,z) 和另一个点x0,y0,z0 之间的欧几里得距离。该表面定义在一个矩形(x,y) 网格上，其中z(x,y) = f(x,y) + random_noise(x,y)。将噪声引入“理想”曲面会大大增加问题的复杂性，因为它需要使用二维三阶样条对曲面进行插值。

不明白为什么在理想表面引入噪声实际上是必要的。如果理想表面确实是理想的，那么应该充分理解x 和y 中的真正多项式拟合可以确定，如果不是分析方法，至少是经验方法。如果随机噪声要模拟实际测量，则只需记录测量足够多次，直到噪声平均为零。同样，使用信号过滤可以帮助消除噪声并揭示信号的真实行为。

要找到表面上离另一个点最近的点，必须使用距离方程及其导数。如果曲面真的只能使用样条基来描述，那么必须reconstruct 样条表示并找到它的导数，这是不平凡的。或者，可以使用精细网格评估表面，但在这里，很快就会遇到内存问题，这就是首先使用插值的原因。

但是，如果我们同意可以使用x 和y 中的简单表达式来定义曲面，那么最小化就变得微不足道了：

为了最小化，看两点D(x,y)之间距离的平方更方便d^2(x,y)（z只是x和y的函数），因为它消除平方根。为了找到D(x,y) 的临界点，我们对其与x 和y 的偏导数，并通过设置= 0 找到它们的根：d/dx D(x,y) = f1(x,y) = 0 和d/dy D(x,y) = f2(x,y)=0。这是一个非线性方程组，我们可以使用scipy.optimize.root 求解。我们只需要传递root 一个猜测（感兴趣的pt 在表面上的投影）和方程组的Jacobian。

import numpy as np
import scipy.interpolate
import scipy.optimize

# Define regular grid surface
xmin,xmax,ymin,ymax = 25, 125, -50, 50
x = np.linspace(xmin,xmax, 201)
y = np.linspace(ymin,ymax, 201)
xx, yy = np.meshgrid(x, y)
z_ideal = ( xx**2 + yy**2 ) / 400

# Fake some measured points on the surface
z_measured = z_ideal + np.random.uniform(-0.1, 0.1, z_ideal.shape)
s_measured = scipy.interpolate.interp2d(x, y, z_measured, kind='cubic')
p_x = np.random.uniform(xmin,xmax,10000)
p_y = np.random.uniform(ymin,ymax,10000)

# z_ideal function
def z(x):
    return (x[0] ** 2 + x[1] ** 2) / 400

# returns the system of equations
def f(x,pt):
    x0,y0,z0 = pt
    f1 = 2*(x[0] - x0) + (z(x)-z0)*x[0]/100
    f2 = 2*(x[1] - y0) + (z(x)-z0)*x[1]/100
    return [f1,f2]

# returns Jacobian of the system of equations
def jac(x, pt):
    x0,y0,z0 = pt
    return [[2*x[0]+1/100*(1/400*(z(x)+2*x[0]**2))-z0, x[0]*x[1]/2e4],
    [2*x[1]+1/100*(1/400*(z(x)+2*x[1]**2))-z0, x[0]*x[1]/2e4]]

def minimize_distance(pt):
    guess = [pt[0],pt[1]]
    return scipy.optimize.root(f,guess,jac=jac, args=pt)

# select a random point from the measured data
x0,y0 = p_x[30], p_y[30]
z0 = float(s_measured(x0,y0))

minimize_distance([x0,y0,z0])

输出：

    fjac: array([[-0.99419141, -0.1076264 ],
       [ 0.1076264 , -0.99419141]])
     fun: array([ -1.05033229e-08,  -2.63163477e-07])
 message: 'The solution converged.'
    nfev: 19
    njev: 2
     qtf: array([  2.80642738e-07,   2.13792093e-06])
       r: array([-2.63044477, -0.48260582, -2.33011149])
  status: 1
 success: True
       x: array([ 110.6726472 ,   39.28642206])

Answer 2

是的！将 K-Means 与聚类一起使用就可以做到这一点。所以s_ideal 将成为目标，然后你在p_z 上进行训练。最终，您将获得质心，它将为您提供表面上最接近的点 s_ideal 到 p 中的每个点。

Here 是一个例子，它非常接近你想要的。