更新：我修改了优化、特征值和求解方法以反映更改。现在所有方法都返回“相同”的向量，允许机器精度。但是我仍然被特征值方法难住了。具体来说，我选择特征向量的切片的如何/为什么没有意义。直到正常匹配其他解决方案之前，这只是试错。如果有人可以纠正/解释我真正应该做什么，或者为什么我所做的有效，我将不胜感激。
谢谢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)
XYZ是一个二维的z值地图，例如z(x,y) = XYZ[x,y]，还是一个点列表，形式为xi, yi, zi = XYZ[i]？根据您的代码示例，我猜测应该是点列表，但在您的介绍中提到了图像，这表明它是z值地图。 - askewchanB = v[3,:]更改为B = np.array(v[3,:])[0]。第一个返回矩阵而不是向量。 - PanchaGil