10

更新：我修改了 Optimize 和 Eigen 和 Solve 方法以反映变化。所有现在都返回“相同”的向量，允许机器精度。 我仍然对 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)
4

2 回答 2

6

优化

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

优化方法为 a 和 b 找到 a 值，使得 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 中的一样。

解决

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

于 2013-04-11T23:31:34.607 回答
5

特征解中平面的法向量是最小特征值的特征向量。一些 Eigen 实现对特征值和特征向量进行排序，而另一些则没有。因此，在某些实现中，将第一个（或最后一个）特征向量作为正常值就足够了。在其他实现中，您必须先对它们进行排序。另一方面，大多数 SVD 实现提供排序值，因此它是简单的第一个（或最后一个）向量。

于 2013-08-22T16:01:38.273 回答