用奇异值分解对三维点进行平面拟合

6
亲爱的stackoverflow用户们,
我正在尝试计算由一组3D点定义的任意(但光滑)表面上的法向量。为此,我使用平面拟合算法,该算法基于我计算法向量的点的10个最近邻点找到局部最小二乘平面。
然而,它并不总是找到看起来最好的平面。因此,我想知道我的实现或算法是否存在缺陷。我使用奇异值分解,因为我在有关平面拟合的几个链接中发现了推荐。下面是一段可以在我的机器上重现这种行为的代码:
#library imports
import numpy as np
import math
import matplotlib.pyplot    as     plt
from   mpl_toolkits.mplot3d import Axes3D

#values used for best plane fit
xyz = np.array([[-1.04194694, -1.17965867,  1.09517722],
[-0.39947906, -1.37104542,  1.36019265],
[-1.0634807 , -1.35020616,  0.46773962],
[-0.48640524, -1.64476106,  0.2726187 ],
[-0.05720509, -1.6791781 ,  0.76964551],
[-1.27522669, -1.10240358,  0.33761405],
[-0.61274031, -1.52709874, -0.09945502],
[-1.402693  , -0.86807757,  0.88866091],
[-0.72520241, -0.86800727,  1.69729388]])

''' best plane fit'''
#1.calculate centroid of points and make points relative to it
centroid         = xyz.mean(axis = 0)
xyzT             = np.transpose(xyz)
xyzR             = xyz - centroid                         #points relative to centroid
xyzRT            = np.transpose(xyzR)                       

#2. calculate the singular value decomposition of the xyzT matrix and get the normal as the last column of u matrix
u, sigma, v       = np.linalg.svd(xyzRT)
normal            = u[2]                                 
normal            = normal / np.linalg.norm(normal)       #we want normal vectors normalized to unity

'''matplotlib display'''
#prepare normal vector for display
forGraphs = list()
forGraphs.append(np.array([centroid[0],centroid[1],centroid[2],normal[0],normal[1], normal[2]]))

#get d coefficient to plane for display
d = normal[0] * centroid[0] + normal[1] * centroid[1] + normal[2] * centroid[2]

# create x,y for display
minPlane = int(math.floor(min(min(xyzT[0]), min(xyzT[1]), min(xyzT[2]))))
maxPlane = int(math.ceil(max(max(xyzT[0]), max(xyzT[1]), max(xyzT[2]))))
xx, yy = np.meshgrid(range(minPlane,maxPlane), range(minPlane,maxPlane))

# calculate corresponding z for display
z = (-normal[0] * xx - normal[1] * yy + d) * 1. /normal[2]

#matplotlib display code
forGraphs = np.asarray(forGraphs)
X, Y, Z, U, V, W = zip(*forGraphs)
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot_surface(xx, yy, z, alpha=0.2)
ax.scatter(xyzT[0],xyzT[1],xyzT[2])
ax.quiver(X, Y, Z, U, V, W)
ax.set_xlim([min(xyzT[0])- 0.1, max(xyzT[0]) + 0.1])
ax.set_ylim([min(xyzT[1])- 0.1, max(xyzT[1]) + 0.1])
ax.set_zlim([min(xyzT[2])- 0.1, max(xyzT[2]) + 0.1])   
plt.show() 

结果是: enter image description here 我希望它更像这样: enter image description here (对于不好的草图,我很抱歉)
那么,这里有什么问题呢?可能是我的matplotlib代码中的显示错误吗?
祝一切顺利!
2个回答

3
维基百科文章中,可以了解到右奇异向量是最小化“正交”的向量。所以我想你不需要转置,使用v[2]代替u[2]即可;对我来说可以运行。请注意,使用第二个元素,即最后一个元素,依赖于numpy(LAPACK)以降序返回奇异值的事实。

谢谢!我在SVD中使用了非转置矩阵(xyzR),并将v [2] 作为法线,它按照我的期望工作。你的答案解决了我的问题! - Glxblt76

0

实际上,我认为主要问题在于您在获取法向量时使用了左奇异矩阵的第三而不是第三,这可能是您遇到问题的原因。在您的代码中使用normal=u[:,2]应该可以解决问题。除此之外,我认为您做得很好。

如果您对非转置点矩阵进行SVD,则确实需要查看右奇异向量(但是v[2]仍然是第三行)。您应该注意的是,在这种情况下,转置矩阵v.T与您原始代码中的u矩阵相同。


网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接