不同算法拟合平面结果差异的原因与解决方法

在计算机视觉中，经常需要估计表面的法线向量。一种常见的方法是将数据点拟合到一个平面上，并使用该平面的法线向量作为表面法线向量。然而，不同的拟合算法可能会产生不同的结果。

2、解决方案

2.1 理解算法差异

• LTSQ（最小二乘法）： 旨在最小化点到平面的距离之和。
• SVD（奇异值分解）： 使用奇异值分解来计算平面参数。
• Eigen（特征值分解）： 使用特征值分解来计算平面参数。
• Solve（求解线性方程组）： 使用求解线性方程组的方法来计算平面参数。
• Optimize（优化）： 使用优化方法来计算平面参数。

2.2 分析示例代码

import numpy as np
import scipy.optimize

def fitPLaneLTSQ(XYZ):
    # Fits a plane to a point cloud, 
    # Where Z = aX + bY + c        ----Eqn #1
    # Rearranging 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)

• LTSQ 和 SVD 产生相同的结果。 (差异约为机器精度的数量级)。
• Eigen 选择了错误的特征向量。 最大的特征值 (lambda = 1.50) 对应的特征向量是 x=[0, sqrt(2)/2, sqrt(2)/2]，与 SVD 和 LTSQ 的结果相同。
• Solve： 我不知道这是如何工作的。

2.3 总结

• 不同算法对平面拟合的结果可能不同，这是由算法本身的原理和数据本身的性质决定的。
• 在选择算法时，需要考虑算法的准确度、鲁棒性和计算复杂度等因素。
• 在使用算法时，需要对算法的原理和参数设置有充分的了解，以确保算法能够正确地工作。