拟合点到平面算法，如何迭代结果? [英] Fit points to a plane algorithms, how to iterpret results?

问题描述

更新:我修改了Optimize和Eigen and Solve方法以反映更改.现在，所有函数都返回相同"向量，以实现机器精度. 我仍然对本征方法感到困惑.具体来说，如何/为什么选择特征向量的切片没有意义.直到正常情况与其他解决方案匹配之前，这只是反复试验.如果有人可以纠正/解释我真正应该做的事情，或者为什么我做过的事情行得通，那么我将不胜感激..

谢谢 亚历山大·克雷默(Alexander Kramer)解释了我为什么要切片的原因，只允许选择一个正确的答案

我有一个深度图像.我想计算深度图像中像素的粗表面法线.我考虑周围的像素(在最简单的情况下为3x3矩阵)，然后将平面拟合到这些点，然后计算该平面的法向单位矢量.

I have a depth image. I want to calculate a crude surface normal for a pixel in the depth image. I consider the surrounding pixels, in the simplest case a 3x3 matrix, and fit a plane to these point, and calculate the normal unit vector to this plane.

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

Sounds easy, but thought best to verify the plane fitting algorithms first. Searching SO and various other sites I see methods using least squares, singlualar value decomposition, eigenvectors/values etc.

尽管我不完全理解数学，但我已经能够使各种片段/示例正常工作.我遇到的问题是，每种方法的答案都不同.我原以为各种答案将是相似的(不完全是)，但是它们似乎有很大的不同.也许有些方法不适合我的数据，但是不确定为什么我得到不同的结果.有什么想法吗?

Although I don't fully understand the maths I have been able to get the various fragments/example to work. The problem I am having, is that I am getting different answers for each method. I was expecting the various answers would be similar (not exact), but they seem significantly different. Perhaps some methods are not suited to my data, but not sure why I am getting different results. Any ideas why?

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

Here is the Updated output of the code:

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]

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

The following code implements five different methods to calculate the surface normal of a plane. The algorithms/code were sourced from various forums on the internet.

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)

推荐答案

优化

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

Optimize

The normal vector of a plane a*x + b*y +c*z = 0, equals (a,b,c)

优化方法找到a和b的值，使得a * x + b * y〜z(〜表示近似值)在计算中完全不使用c的值.我没有在这台计算机上安装numpy，但我希望将模型更改为(a * x + b * y)/c应该可以解决此问题.对于所有数据集，它不会给出相同的结果.此方法将始终假定一个平面穿过原点.

The optimize method finds a values for a and b such that a*x+b*y~z (~ denotes approximates) It omits to use the value of c in the calculation at all. I don't have numpy installed on this machine but I expect that changing the model to (a*x+b*y)/c should fix this method. It will not give the same result for all data-sets. This method will always assume a plane that goes through the origin.

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

produce the same results. (The difference is about the size of machine precision).

选择了错误的特征向量.就像在SVD和LTSQ中一样，对应于最大特征值(lambda = 1.50)的特征向量为x=[0, sqrt(2)/2, sqrt(2)/2].

The wrong eigenvector is chosen. The eigenvector corresponding to the greatest eigenvalue (lambda = 1.50) is x=[0, sqrt(2)/2, sqrt(2)/2] just as in the SVD and LTSQ.

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

I have no clue how this is supposed to work.

这篇关于拟合点到平面算法，如何迭代结果?的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！