Ajustar los puntos a los algoritmos de un plano, ¿cómo interpretar los resultados?

Actualización : He modificado los métodos Optimize, Eigen y Solve para reflejar los cambios. Todos ahora devuelven el “mismo” vector permitiendo la precisión de la máquina. Todavía estoy perplejo con el método Eigen. Específicamente, cómo / por qué selecciono una porción del vector propio no tiene sentido. Fue solo prueba y error hasta que lo normal coincidió con las otras soluciones. Si alguien puede corregir / explicar lo que realmente debería hacer, o por qué lo que he hecho funciona, lo apreciaría. .

Gracias Alexander Kramer, por explicar por qué tomo una porción, solo se me permitió seleccionar una respuesta correcta

Tengo una imagen de profundidad. Quiero calcular una superficie bruta normal para un píxel en la imagen de profundidad. Considero que los píxeles circundantes, en el caso más simple, son una matriz de 3×3, y ajuste un plano a estos puntos, y calculo el vector de unidad normal para este plano.

Suena fácil, pero se pensó mejor para verificar primero los algoritmos de ajuste de planos. Buscando SO y varios otros sitios, veo métodos que utilizan mínimos cuadrados, descomposición de valores singlualar, vectores propios / valores, etc.

Aunque no comprendo completamente las matemáticas, he podido hacer funcionar los diversos fragmentos / ejemplos. El problema que tengo es que obtengo respuestas diferentes para cada método. Esperaba que las distintas respuestas fueran similares (no exactas), pero parecen significativamente diferentes. Quizás algunos métodos no sean adecuados para mis datos, pero no estoy seguro de por qué estoy obteniendo resultados diferentes. ¿Alguna idea de por qué?

    Aquí está la salida actualizada del 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] 

    El siguiente código implementa cinco métodos diferentes para calcular la normal de superficie de un plano. Los algoritmos / códigos se obtuvieron de varios foros en 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,BT) 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) 

    Optimizar

    El vector normal de un plano a * x + b * y + c * z = 0, es igual a (a, b, c)

    El método de optimización encuentra valores para a y b tales que a * x + b * y ~ z (~ denota aproximados) Omite el uso del valor de c en el cálculo. No tengo el número instalado en esta máquina, pero espero que cambiar el modelo a (a * x + b * y) / c debería solucionar este método. No dará el mismo resultado para todos los conjuntos de datos. Este método siempre asumirá un plano que pase por el origen.

    SVD y LTSQ

    Producen los mismos resultados. (La diferencia es sobre el tamaño de la precisión de la máquina).

    Eigen

    Se elige el vector propio equivocado. El vector propio correspondiente al mayor valor propio ( lambda = 1.50 ) es x=[0, sqrt(2)/2, sqrt(2)/2] tal como en la SVD y LTSQ.

    Resolver

    No tengo ni idea de cómo se supone que funciona esto.

    El vector normal del plano en la solución Eigen es el vector propio para el valor propio más pequeño. Algunas implementaciones de Eigen clasifican los valores propios y los vectores propios que otros no. Entonces, en algunas implementaciones, es suficiente tomar el primer (o último) vector propio para normal. En otras implementaciones tienes que ordenarlas primero. Por otro lado, la mayoría de las implementaciones SVD proporcionan valores ordenados, por lo que es el primer vector (o el último) simple.