python - 拟合点到平面算法，如何解释结果？

Question

更新:我修改了 Optimize 和 Eigen and Solve 方法以反射(reflect)变化。所有现在都返回允许机器精度的“相同”向量。 我仍然对 Eigen 方法感到困惑。特别是我如何/为什么选择特征向量的切片没有意义。在正常匹配其他解决方案之前，这只是反复试验。如果有人可以纠正/解释我真正应该做什么，或者为什么我所做的工作有效，我将不胜感激。。

谢谢 Alexander Kramer，解释了为什么我要分一杯羹，只允许选择一个正确答案

我有一张深度图。我想计算深度图像中像素的粗糙表面法线。我考虑周围的像素，在最简单的情况下是一个 3x3 矩阵，并为这些点拟合一个平面，并计算该平面的法线单位向量。

听起来很简单，但我认为最好先验证平面拟合算法。搜索 SO 和其他各种网站，我看到了使用最小二乘法、奇异值分解、特征向量/值等的方法。

虽然我不完全理解数学，但我已经能够让各种片段/示例发挥作用。我遇到的问题是，每种方法都得到不同的答案。我期待各种答案会相似(不准确)，但它们似乎有很大不同。也许有些方法不适合我的数据，但不确定为什么我会得到不同的结果。有什么想法吗？

这是代码的更新后的输出:

LTSQ:   [ -8.10792259e-17   7.07106781e-01  -7.07106781e-01]
SVD:    [ 0.                0.70710678      -0.70710678]
Eigen:  [ 0.                0.70710678      -0.70710678]
Solve:  [ 0.                0.70710678       0.70710678]
Optim:  [ -1.56069661e-09   7.07106781e-01   7.07106782e-01]

以下代码实现了五种不同的方法来计算平面的表面法线。算法/代码来自互联网上的各种论坛。

import numpy as np
import scipy.optimize

def fitPLaneLTSQ(XYZ):
    # Fits a plane to a point cloud, 
    # Where Z = aX + bY + c        ----Eqn #1
    # Rearanging Eqn1: aX + bY -Z +c =0
    # Gives normal (a,b,-1)
    # Normal = (a,b,-1)
    [rows,cols] = XYZ.shape
    G = np.ones((rows,3))
    G[:,0] = XYZ[:,0]  #X
    G[:,1] = XYZ[:,1]  #Y
    Z = XYZ[:,2]
    (a,b,c),resid,rank,s = np.linalg.lstsq(G,Z) 
    normal = (a,b,-1)
    nn = np.linalg.norm(normal)
    normal = normal / nn
    return normal


def fitPlaneSVD(XYZ):
    [rows,cols] = XYZ.shape
    # Set up constraint equations of the form  AB = 0,
    # where B is a column vector of the plane coefficients
    # in the form b(1)*X + b(2)*Y +b(3)*Z + b(4) = 0.
    p = (np.ones((rows,1)))
    AB = np.hstack([XYZ,p])
    [u, d, v] = np.linalg.svd(AB,0)        
    B = v[3,:];                    # Solution is last column of v.
    nn = np.linalg.norm(B[0:3])
    B = B / nn
    return B[0:3]


def fitPlaneEigen(XYZ):
    # Works, in this case but don't understand!
    average=sum(XYZ)/XYZ.shape[0]
    covariant=np.cov(XYZ - average)
    eigenvalues,eigenvectors = np.linalg.eig(covariant)
    want_max = eigenvectors[:,eigenvalues.argmax()]
    (c,a,b) = want_max[3:6]    # Do not understand! Why 3:6? Why (c,a,b)?
    normal = np.array([a,b,c])
    nn = np.linalg.norm(normal)
    return normal / nn  

def fitPlaneSolve(XYZ):
    X = XYZ[:,0]
    Y = XYZ[:,1]
    Z = XYZ[:,2] 
    npts = len(X)
    A = np.array([ [sum(X*X), sum(X*Y), sum(X)],
                   [sum(X*Y), sum(Y*Y), sum(Y)],
                   [sum(X),   sum(Y), npts] ])
    B = np.array([ [sum(X*Z), sum(Y*Z), sum(Z)] ])
    normal = np.linalg.solve(A,B.T)
    nn = np.linalg.norm(normal)
    normal = normal / nn
    return normal.ravel()

def fitPlaneOptimize(XYZ):
    def residiuals(parameter,f,x,y):
        return [(f[i] - model(parameter,x[i],y[i])) for i in range(len(f))]


    def model(parameter, x, y):
        a, b, c = parameter
        return a*x + b*y + c

    X = XYZ[:,0]
    Y = XYZ[:,1]
    Z = XYZ[:,2]
    p0 = [1., 1.,1.] # initial guess
    result = scipy.optimize.leastsq(residiuals, p0, args=(Z,X,Y))[0]
    normal = result[0:3]
    nn = np.linalg.norm(normal)
    normal = normal / nn
    return normal


if __name__=="__main__":
    XYZ = np.array([
        [0,0,1],
        [0,1,2],
        [0,2,3],
        [1,0,1],
        [1,1,2],
        [1,2,3],
        [2,0,1],
        [2,1,2],
        [2,2,3]
        ])
    print "Solve: ", fitPlaneSolve(XYZ)
    print "Optim: ",fitPlaneOptimize(XYZ)
    print "SVD:   ",fitPlaneSVD(XYZ)
    print "LTSQ:  ",fitPLaneLTSQ(XYZ)
    print "Eigen: ",fitPlaneEigen(XYZ)

Answer 1

优化

一个平面的法向量a*x + b*y +c*z = 0, 等于(a,b,c)

optimize 方法找到 a 和 b 的值，使得 a*x+b*y~z(~ 表示近似值)它在计算中完全省略了 c 的值。我没有在这台机器上安装 numpy，但我希望将模型更改为 (a*x+b*y)/c 应该可以修复此方法。它不会为所有数据集提供相同的结果。此方法将始终假设一个平面穿过原点。

SVD 和 LTSQ

产生相同的结果。 (区别在于机器精度的大小)。

本征

选择了错误的特征向量。与最大特征值 (lambda = 1.50) 对应的特征向量是 x=[0, sqrt(2)/2, sqrt(2)/2] 就像在SVD 和 LTSQ。

解决

我不知道这应该如何工作。