Подгонка точек к плоскости алгоритмов, как ее интерпретировать результаты?
Обновить: Я изменил методы Optimize и Eigen and Solve, чтобы отразить изменения. Все теперь возвращают «тот же» вектор, учитывающий точность станка.Я все еще нахожусь в тупике на метод Эйгена. Конкретно как / почему я выбираю срез собственного вектора, не имеет смысла. Это было просто методом проб и ошибок, пока обычные не совпали с другими решениями. Если кто-то может исправить / объяснить, что я действительно должен делать, или почему то, что я сделал, работает, я был бы признателен..
Спасибо Александр Крамер за объяснение, почему я беру кусочек, позволил выбрать только один правильный ответ
У меня есть изображение глубины. Я хочу вычислить грубую нормаль поверхности для пикселя в изображении глубины. Я рассматриваю окружающие пиксели, в простейшем случае матрицу 3х3, и подгоняю плоскость к этим точкам, и вычисляю вектор нормали к этой плоскости.
Звучит просто, но лучше всего сначала проверить алгоритмы подгонки плоскости. В поисках 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)