Подгонка точек к плоскости алгоритмов, как ее интерпретировать результаты?

Обновить: Я изменил методы 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)

Ответы на вопрос(2)

Ваш ответ на вопрос