将 "quadratic" 个表面拟合到 3D 中的数据点

Fitting "quadratic" surface to data points in 3D

我正在尝试将二次平面拟合到 python 中的数据点云。我的平面函数的形式是

   f(x,y,z) = a*x**2 + b*y**2 + c*x*y + d*x + e*y + f - z

目前,我的数据点没有与之相关的错误,但是,如果有必要,可以假设一些错误。遵循 , I work out the vertical distance from a point p0=(x0,y0,z0) (which correspond to my data points) to a point on the plane p=(x,y,z) following this method 的建议。然后我得到

def vertical_distance(params,p0):
    *** snip ***
    nominator = f + a*x**2 + b*y**2 + c*x*y - x0*(2*a*x-c*y-d) - y0*(2*b*y-c*x-e) + z0
    denominator = sqrt((2*a*x+c*y+d)**2 + (2*b*y+c*x+e)**2 + (-1)**2)
    return nominator/denominator

最终,我认为我需要最小化 vertical_distance 函数。我可以愉快地在二维中提供起始参数列表 (params) 和数据点数组,但是,我不确定如何在 3D 中实现这一点。 ODR 包似乎只允许包含 x、y 或二维的数据。此外,如何将平面 (p) 上的点实现到最小化例程中?我猜想在拟合操作期间,点会根据参数优化而变化,因此在那个时刻平面的精确方程。

我想 «quadratic surface» 比 «plane» 更准确。

而问题是拟合 z = ax^2 + by^2 + cxy + dx + ey + f 到给定的一组点 P.

要通过优化做到这一点,您需要制定残差函数(例如,垂直距离)。

对于来自 P 残差的每个 3D 点 p 是

|p_2 – ap_0^2 + bp_1^2 + c*p_0*p_1 + dp_0 + ep_1 + f|

您需要最小化所有残差,即它们的平方和,可变参数 a…f。

下面的代码在技术上应该可以解决上面的问题。但是解决问题是 multi-extremal 并且如果没有良好的起点或搜索全球化,这样的例程可能无法找到正确的参数集。

import numpy
import scipy.optimize

P = numpy.random.rand(3,10) # given point set

def quadratic(x,y, a, b, c, d, e, f): 
  #fit quadratic surface
  return a*x**2 + b*y**2 + c*x*y + d*x + e*y + f

def residual(params, points):
  #total residual
  residuals = [
    p[2] - quadratic(p[0], p[1],
    params[0], params[1], params[2], params[3], params[4], params[5]) for p in points]

  return numpy.linalg.norm(residuals)

result = scipy.optimize.minimize(residual, 
                                 (1, 1, 0, 0, 0, 0),#starting point
                                 args=P)