最适合 3D 数据的平面：两种不同方法的不同结果

Question

(虽然有很多关于如何使平面最适合 SO 上的某些 3D 数据的问题，但我找不到这个问题的答案。)

给定 N (x, y, z) 点，我需要最合适的平面

a*x + b*y + c*z + d = 0

通过 a, b, c, d 系数定义，最小化从点到平面的正交距离的平均值。点平面正交距离（对于给定的 (x0, y0, z0) 点）是 defined as:

d = |a*x0 + b*y0 + c*z0 + d|/sqrt(a^2 + b^2 + c^2)

我设置了两种方法（代码如下）：

Singular-Value Decomposition (source)
Basin-Hopping 平均正交距离的最小化

据我了解，SVD 方法应该通过最小化正交距离分析来产生精确最佳拟合平面。相反，我发现 BH 方法给出的结果 比所谓的精确 SVD 方法更好，即使 BH 运行次数很少也是如此。

"better" 我的意思是 BH 方法的最终平均正交距离值比 SVD 方法小。

我在这里错过了什么？

import numpy as np
import scipy.optimize as optimize


def perp_error(params, xyz):
    """
    Mean of the absolute values for the perpendicular distance of the
    'xyz' points, to the plane defined by the coefficients 'a,b,c,d' in
    'params'.
    """
    a, b, c, d = params
    x, y, z = xyz
    length = np.sqrt(a**2 + b**2 + c**2)
    return (np.abs(a * x + b * y + c * z + d) / length).mean()


def minPerpDist(x, y, z, N_min):
    """
    Basin-Hopping method, minimize mean absolute values of the
    orthogonal distances.
    """
    def unit_length(params):
        """
        Constrain the vector perpendicular to the plane to be of unit length.
        """
        a, b, c, d = params
        return a**2 + b**2 + c**2 - 1

    # Random initial guess for the a,b,c,d plane coefficients.
    initial_guess = np.random.uniform(-10., 10., 4)

    # Constrain the vector perpendicular to the plane to be of unit length.
    cons = ({'type': 'eq', 'fun': unit_length})
    min_kwargs = {"constraints": cons, "args": [x, y, z]}

    # Use Basin-Hopping to obtain the best fit coefficients.
    sol = optimize.basinhopping(
        perp_error, initial_guess, minimizer_kwargs=min_kwargs, niter=N_min)
    abcd = list(sol.x)

    return abcd


def SVD(X):
    """
    Singular value decomposition method.
    Source: https://gist.github.com/lambdalisue/7201028
    """
    # Find the average of points (centroid) along the columns
    C = np.average(X, axis=0)

    # Create CX vector (centroid to point) matrix
    CX = X - C
    # Singular value decomposition
    U, S, V = np.linalg.svd(CX)
    # The last row of V matrix indicate the eigenvectors of
    # smallest eigenvalues (singular values).
    N = V[-1]

    # Extract a, b, c, d coefficients.
    x0, y0, z0 = C
    a, b, c = N
    d = -(a * x0 + b * y0 + c * z0)

    return a, b, c, d


# Generate a random plane.
seed = np.random.randint(100000)
print("Seed: {}".format(seed))
np.random.seed(seed)
a, b, c, d = np.random.uniform(-10., 10., 4)
print("Orig abc(d=1): {:.3f} {:.3f} {:.3f}\n".format(a / d, b / d, c / d))

# Generate random (x, y, z) points.
N = 200
x, y = np.random.uniform(-5., 5., (2, N))
z = -(a * x + b * y + d) / c
# Add scatter in z.
z = z + np.random.uniform(-.2, .2, N)

# Solve using SVD method.
a, b, c, d = SVD(np.array([x, y, z]).T)
print("SVD  abc(d=1): {:.3f} {:.3f} {:.3f}".format(a / d, b / d, c / d))
# Orthogonal mean distance
print("Perp err: {:.5f}\n".format(perp_error((a, b, c, d), (x, y, z))))

# Solve using Basin-Hopping.
abcd = minPerpDist(x, y, z, 500)
a, b, c, d = abcd
print("BH   abc(d=1): {:.3f} {:.3f} {:.3f}".format(a / d, b / d, c / d))
print("Perp err: {:.5f}".format(perp_error(abcd, (x, y, z))))

Answer 1

我相信我找到了差异的原因。

当我使用 Basin-Hopping 最小化点到平面的垂直距离时，我使用的是 绝对值 点平面距离：

d_abs = |a*x0 + b*y0 + c*z0 + d| / sqrt(a^2 + b^2 + c^2)

另一方面，SVD 方法显然最小化了平方点平面距离：

d_sqr = (a*x0 + b*y0 + c*z0 + d)^2 / (a^2 + b^2 + c^2)

如果在问题共享的代码中，我在 perp_error() 函数中使用平方距离而不是绝对值距离，则两种方法给出完全相同的答案。

最适合 3D 数据的平面：两种不同方法的不同结果

Best fit plane to 3D data: different results for two different methods

python

3d

scipy

svd