Atualizar: Eu modifiquei os métodos Optimize e Eigen and Solve para refletir as alterações. Todos agora retornam o "mesmo" vetor, permitindo a precisão da máquina.Eu ainda estou perplexo no método Eigen. Especificamente Como / Por que eu seleciono fatia do vetor próprio não faz sentido. Foi apenas tentativa e erro até que o normal correspondesse às outras soluções. Se alguém puder corrigir / explicar o que eu realmente deveria fazer, ou por que o que tenho feito, eu agradeceria..

obrigado Alexander Kramer, por explicar por que eu tomo uma fatia, só alowed para selecionar uma resposta correta

Eu tenho uma imagem de profundidade. Eu quero calcular uma superfície bruta normal para um pixel na imagem de profundidade. Eu considero os pixels circundantes, no caso mais simples, uma matriz 3x3, e encaixo um plano nesse ponto, e calculo o vetor unitário normal para esse plano.

Parece fácil, mas achei melhor verificar primeiro os algoritmos de ajuste do plano. Pesquisando SO e vários outros sites eu vejo métodos usando mínimos quadrados, decomposição do valor singlualar, autovetores / valores etc.

Embora eu não entenda completamente as contas, consegui que os vários fragmentos / exemplos funcionassem. O problema que estou tendo é que estou recebendo respostas diferentes para cada método. Eu estava esperando que as várias respostas seriam semelhantes (não exatas), mas elas parecem significativamente diferentes. Talvez alguns métodos não sejam adequados aos meus dados, mas não sei por que estou obtendo resultados diferentes. Alguma idéia por quê?

Aqui está oSaída atualizada do código:

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]

O código a seguir implementa cinco métodos diferentes para calcular a superfície normal de um plano. Os algoritmos / código foram obtidos de vários fóruns na internet.

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)