python - 如何在 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. 计算 pt 处的斜率(正切)向量 m = [ m_x, m_y ]
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建议一种方法可以为 1D 情况快速产生精度非常高的结果，但它适用于单一维度，对我来说可能很难转换为两个维度。

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^2(x,y) (z 只是 的函数code>x 和 y) 在两点之间，D(x,y)，因为它消除了平方根。为了找到 D(x,y) 的临界点，我们取它的偏导数 w.r.t 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])