Plano de mejor ajuste dado un conjunto de puntos

Question

No hay nada más que explicar. Es que no sé cómo encontrar el mejor plano de ajuste dado un conjunto de $N$ puntos en un $3D$ espacio. Luego tengo que escribir el algoritmo correspondiente. Gracias ;)

Answer 1

Resta el centroide, forma un $3\times N$ matriz $\mathbf X$ de las coordenadas resultantes y calcular su descomposición del valor singular . El vector normal del plano mejor ajustado es el vector singular izquierdo correspondiente al menor valor singular. Véase esta respuesta para explicar por qué esto es numéricamente preferible al cálculo del vector propio de $\mathbf X\mathbf X^\top$ correspondiente al menor valor propio.

Aquí hay una implementación en Python, como se pidió:

import numpy as np

# generate some random test points 
m = 20 # number of points
delta = 0.01 # size of random displacement
origin = np.random.rand(3, 1) # random origin for the plane
basis = np.random.rand(3, 2) # random basis vectors for the plane
coefficients = np.random.rand(2, m) # random coefficients for points on the plane

# generate random points on the plane and add random displacement
points = basis @ coefficients \
         + np.tile(origin, (1, m)) \
         + delta * np.random.rand(3, m)

# now find the best-fitting plane for the test points

# subtract out the centroid and take the SVD
svd = np.linalg.svd(points - np.mean(points, axis=1, keepdims=True))

# Extract the left singular vectors
left = svd[0]

1 2

# the corresponding left singular vector is the normal vector of the best-fitting plane

left[:, -1]

2

# its dot product with the basis vectors of the plane is approximately zero

left[:, -1] @ basis

2

Answer 2

Una simple solución de mínimos cuadrados debería ser suficiente.

La ecuación de un plano es: $ax + by + c = z$ . Así que establece matrices como esta con todos tus datos:

$$ \begin{bmatrix} x_0 & y_0 & 1 \\ x_1 & y_1 & 1 \\ &... \\ x_n & y_n & 1 \\ \end{bmatrix} \begin{bmatrix} a \\ b \\ c \end{bmatrix} = \begin{bmatrix} z_0 \\ z_1 \\ ... \\ z_n \end{bmatrix} $$ En otras palabras:

$$Ax = B$$

Ahora resuelve para $x$ que son sus coeficientes. Pero como (supongo) tienes más de 3 puntos, el sistema está sobredeterminado, así que necesitas usar la pseudoinversa de la izquierda: $A^+ = (A^T A)^{-1} A^T$ . Así que la respuesta es: $$ \begin{bmatrix} a \\ b \\ c \end{bmatrix} = (A^T A)^{-1} A^T B $$

Y aquí hay un sencillo código Python con un ejemplo:

import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import numpy as np

# These constants are to create random data for the sake of this example
N_POINTS = 10
TARGET_X_SLOPE = 2
TARGET_y_SLOPE = 3
TARGET_OFFSET  = 5
EXTENTS = 5
NOISE = 5

# Create random data.
# In your solution, you would provide your own xs, ys, and zs data.
xs = [np.random.uniform(2*EXTENTS)-EXTENTS for i in range(N_POINTS)]
ys = [np.random.uniform(2*EXTENTS)-EXTENTS for i in range(N_POINTS)]
zs = []
for i in range(N_POINTS):
    zs.append(xs[i]*TARGET_X_SLOPE + \
              ys[i]*TARGET_y_SLOPE + \
              TARGET_OFFSET + np.random.normal(scale=NOISE))

# plot raw data
plt.figure()
ax = plt.subplot(111, projection='3d')
ax.scatter(xs, ys, zs, color='b')

# do fit
tmp_A = []
tmp_b = []
for i in range(len(xs)):
    tmp_A.append([xs[i], ys[i], 1])
    tmp_b.append(zs[i])
b = np.matrix(tmp_b).T
A = np.matrix(tmp_A)

# Manual solution
fit = (A.T * A).I * A.T * b
errors = b - A * fit
residual = np.linalg.norm(errors)

# Or use Scipy
# from scipy.linalg import lstsq
# fit, residual, rnk, s = lstsq(A, b)

print("solution: %f x + %f y + %f = z" % (fit[0], fit[1], fit[2]))
print("errors: \n", errors)
print("residual:", residual)

# plot plane
xlim = ax.get_xlim()
ylim = ax.get_ylim()
X,Y = np.meshgrid(np.arange(xlim[0], xlim[1]),
                  np.arange(ylim[0], ylim[1]))
Z = np.zeros(X.shape)
for r in range(X.shape[0]):
    for c in range(X.shape[1]):
        Z[r,c] = fit[0] * X[r,c] + fit[1] * Y[r,c] + fit[2]
ax.plot_wireframe(X,Y,Z, color='k')

ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('z')
plt.show()

Answer 3

Considerando un plano de ecuación $Ax+By+Cz=0$ y un punto de coordenadas $(x_i,y_i,z_i)$ la distancia viene dada por $$d_i=\pm\frac{Ax_i+By_i+Cz_i}{\sqrt{A^2+B^2+C^2}}$$ y supongo que quieres minimizar $$F=\sum_{i=1}^n d_i^2=\sum_{i=1}^n\frac{(Ax_i+By_i+Cz_i)^2}{{A^2+B^2+C^2}}$$ Configurar $C=1$ entonces tenemos que minimizar con respecto a $A,B$ $$F=\sum_{i=1}^n\frac{(Ax_i+By_i+z_i)^2}{{A^2+B^2+1}}$$ Tomar derivados $$F'_A=\sum _{i=1}^n \left(\frac{2 x_i (A x_i+B y_i+z_i)}{A^2+B^2+1}-\frac{2 A (A x_i+B y_i+z_i)^2}{\left(A^2+B^2+1\right)^2}\right)$$ $$F'_B=\sum _{i=1}^n \left(\frac{2 y_i (A x_i+B y_i+z_i)}{A^2+B^2+1}-\frac{2 B (A x_i+B y_i+z_i)^2}{\left(A^2+B^2+1\right)^2}\right)$$ Dado que fijaremos estas derivadas iguales a $0$ las ecuaciones pueden simplificarse a $$\sum _{i=1}^n \left({ x_i (A x_i+B y_i+z_i)}-\frac{ A (A x_i+B y_i+z_i)^2}{\left(A^2+B^2+1\right)}\right)=0$$ $$\sum _{i=1}^n \left({ y_i (A x_i+B y_i+z_i)}-\frac{ B (A x_i+B y_i+z_i)^2}{\left(A^2+B^2+1\right)}\right)=0$$ que son no lineales con respecto a los parámetros $A,B$ ; entonces, se requieren buenas estimaciones ya que probablemente se utilizará Newton-Raphson para pulir las soluciones.

Se pueden obtener haciendo primero una regresión multilineal (sin intercepción en su caso) $$z=\alpha x+\beta y$$ y utilizar $A=-\alpha$ y $B=-\beta$ para iniciar el proceso iterativo. Los valores vienen dados por $$A=-\frac{\text{Sxy} \,\text{Syz}-\text{Sxz}\, \text{Syy}}{\text{Sxy}^2-\text{Sxx}\, \text{Syy}}\qquad B=-\frac{\text{Sxy}\, \text{Sxz}-\text{Sxx}\, \text{Syz}}{\text{Sxy}^2-\text{Sxx}\, \text{Syy}}$$

A modo de ejemplo, consideraré los siguientes datos $$\left( \begin{array}{ccc} x & y & z \\ 1 & 1 & 9 \\ 1 & 2 & 14 \\ 1 & 3 & 20 \\ 2 & 1 & 11 \\ 2 & 2 & 17 \\ 2 & 3 & 23 \\ 3 & 1 & 15 \\ 3 & 2 & 20 \\ 3 & 3 & 26 \end{array} \right)$$

El paso previo da $z=2.97436 x+5.64103 y$ y resolviendo las ecuaciones rigurosas se llega a $A=-2.97075$ , $B=-5.64702$ que se acercan bastante a las estimaciones (debido a los pequeños errores).

Answer 4

Para completar la respuesta de Claude Leibovici :

Con el ejemplo numérico propuesto por Claude Leibovici que calculó los parámetros de un plano ajustado $\quad z=Ax+By\quad$ el ajuste del plano $\quad Z=\alpha X+\beta Y+\gamma\quad$ puede llevarse a cabo gracias al método de componentes principales (como sugiere joriki).

La teoría se puede encontrar en muchos libros. En las páginas 24-25 del documento se ofrece una sinopsis: https://fr.scribd.com/doc/31477970/Regressions-et-trajectoires-3D

Los símbolos utilizados a continuación corresponden a los de las fórmulas del documento de referencia.