algorithm - libros - Plano de mínimos cuadrados 3D
libros en español pdf (9)
Al igual que con cualquier enfoque de mínimos cuadrados, procedemos así:
Antes de comenzar a codificar
Escriba una ecuación para un plano en alguna parametrización, digamos
0 = ax + by + z + den los parámetros(a, b, d).Encuentre una expresión
D(/vec{v};a, b, d)para la distancia desde un punto arbitrario/vec{v}.Escriba la suma
S = /sigma_i=0,n D^2(/vec{x}_i), y simplifique hasta que se exprese en términos de sumas simples de los componentes devcomo/sigma v_x,/sigma v_y^2,/sigma v_x*v_z...Escriba las expresiones de minimización por parámetro
dS/dx_0 = 0,dS/dy_0 = 0... que le da un conjunto de tres ecuaciones en tres parámetros y las sumas del paso anterior.Resuelve este conjunto de ecuaciones para los parámetros.
(o para casos simples, solo busque el formulario). Usar un paquete de álgebra simbólica (como Mathematica) podría hacerte la vida mucho más fácil.
La codificación
- Escriba el código para formar las sumas necesarias y encuentre los parámetros del último conjunto anterior.
Alternativas
Ten en cuenta que si tuvieras solo tres puntos, sería mejor que encuentres el avión que los atraviesa.
Además, si la solución analítica es inviable (no es el caso de un avión, pero es posible en general) puede hacer los pasos 1 y 2, y usar un minimizador de Monte Carlo en la suma del paso 3.
¿Cuál es el algoritmo para calcular un plano de mínimos cuadrados en el espacio (x, y, z), dado un conjunto de puntos de datos 3D? En otras palabras, si tuviera un montón de puntos como (1, 2, 3), (4, 5, 6), (7, 8, 9), etc., ¿cómo sería posible calcular el plano de mejor ajuste f (x, y) = ax + por + c? ¿Cuál es el algoritmo para obtener a, byc de un conjunto de puntos 3D?
Consulte "Adaptación de datos de los mínimos cuadrados" de David Eberly para ver cómo se me ocurrió este para minimizar el ajuste geométrico (distancia ortogonal de los puntos al plano).
bool Geom_utils::Fit_plane_direct(const arma::mat& pts_in, Plane& plane_out)
{
bool success(false);
int K(pts_in.n_cols);
if(pts_in.n_rows == 3 && K > 2) // check for bad sizing and indeterminate case
{
plane_out._p_3 = (1.0/static_cast<double>(K))*arma::sum(pts_in,1);
arma::mat A(pts_in);
A.each_col() -= plane_out._p_3; //[x1-p, x2-p, ..., xk-p]
arma::mat33 M(A*A.t());
arma::vec3 D;
arma::mat33 V;
if(arma::eig_sym(D,V,M))
{
// diagonalization succeeded
plane_out._n_3 = V.col(0); // in ascending order by default
if(plane_out._n_3(2) < 0)
{
plane_out._n_3 = -plane_out._n_3; // upward pointing
}
success = true;
}
}
return success;
}
Programado a 37 microsegundos para adaptar un avión a 1000 puntos (programa Windows 7, i7, 32 bits)
Esto se reduce al problema de Mínimos cuadrados totales , que se puede resolver usando la descomposición SVD .
Código C ++ usando OpenCV:
float fitPlaneToSetOfPoints(const std::vector<cv::Point3f> &pts, cv::Point3f &p0, cv::Vec3f &nml) {
const int SCALAR_TYPE = CV_32F;
typedef float ScalarType;
// Calculate centroid
p0 = cv::Point3f(0,0,0);
for (int i = 0; i < pts.size(); ++i)
p0 = p0 + conv<cv::Vec3f>(pts[i]);
p0 *= 1.0/pts.size();
// Compose data matrix subtracting the centroid from each point
cv::Mat Q(pts.size(), 3, SCALAR_TYPE);
for (int i = 0; i < pts.size(); ++i) {
Q.at<ScalarType>(i,0) = pts[i].x - p0.x;
Q.at<ScalarType>(i,1) = pts[i].y - p0.y;
Q.at<ScalarType>(i,2) = pts[i].z - p0.z;
}
// Compute SVD decomposition and the Total Least Squares solution, which is the eigenvector corresponding to the least eigenvalue
cv::SVD svd(Q, cv::SVD::MODIFY_A|cv::SVD::FULL_UV);
nml = svd.vt.row(2);
// Calculate the actual RMS error
float err = 0;
for (int i = 0; i < pts.size(); ++i)
err += powf(nml.dot(pts[i] - p0), 2);
err = sqrtf(err / pts.size());
return err;
}
La ecuación para un plano es: ax + by + c = z. Así que configure matrices como esta con todos sus datos:
x_0 y_0 1
A = x_1 y_1 1
...
x_n y_n 1
Y
a
x = b
c
Y
z_0
B = z_1
...
z_n
En otras palabras: Ax = B. Ahora resuelve para x cuáles son tus coeficientes. Pero como (supongo) que tiene más de 3 puntos, el sistema está sobredeterminado, por lo que debe usar el pseudo inverso inverso. Entonces la respuesta es:
a
b = (A^T A)^-1 A^T B
c
Y aquí hay un código Python simple con un ejemplo:
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import numpy as np
N_POINTS = 10
TARGET_X_SLOPE = 2
TARGET_y_SLOPE = 3
TARGET_OFFSET = 5
EXTENTS = 5
NOISE = 5
# create random 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)
fit = (A.T * A).I * A.T * b
errors = b - A * fit
residual = np.linalg.norm(errors)
print "solution:"
print "%f x + %f y + %f = z" % (fit[0], fit[1], fit[2])
print "errors:"
print errors
print "residual:"
print 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()
Parece que todo lo que quiere hacer es una regresión lineal con 2 regresores. La página de wikipedia sobre el tema debe decirle todo lo que necesita saber y más.
Si tiene n puntos de datos (x [i], y [i], z [i]), calcule la matriz simétrica A 3x3 cuyas entradas son:
sum_i x[i]*x[i], sum_i x[i]*y[i], sum_i x[i]
sum_i x[i]*y[i], sum_i y[i]*y[i], sum_i y[i]
sum_i x[i], sum_i y[i], n
Calcule también el vector de 3 elementos b:
{sum_i x[i]*z[i], sum_i y[i]*z[i], sum_i z[i]}
Luego resuelve Ax = b para el dado A y b. Los tres componentes del vector solución son los coeficientes del plano de ajuste de mínimos cuadrados {a, b, c}.
Tenga en cuenta que este es el ajuste de "mínimos cuadrados ordinarios", que es apropiado solo cuando se espera que z sea una función lineal de xey. Si está buscando más generalmente un "plano de mejor ajuste" en 3 espacios, es posible que desee aprender sobre mínimos cuadrados "geométricos".
Tenga en cuenta también que esto fallará si sus puntos están en una línea, como son sus puntos de ejemplo.
Todo lo que tendrás que hacer es resolver el sistema de ecuaciones.
Si esos son sus puntos: (1, 2, 3), (4, 5, 6), (7, 8, 9)
Eso te da las ecuaciones:
3=a*1 + b*2 + c
6=a*4 + b*5 + c
9=a*7 + b*8 + c
Entonces, tu pregunta debería ser: ¿Cómo resuelvo un sistema de ecuaciones?
Por lo tanto, recomiendo leer this pregunta SO.
Si he entendido mal tu pregunta, avísanos.
EDITAR :
Ignora mi respuesta ya que probablemente querías decir algo más.
a menos que alguien me diga cómo escribir ecuaciones aquí, permítame anotar los cálculos finales que tiene que hacer:
primero, dados los puntos r_i / n / R, i = 1..N, calcule el centro de masa de todos los puntos:
r_G = /frac{/sum_{i=1}^N r_i}{N}
luego, calcule el vector normal n, que junto con el vector base r_G define el plano calculando la matriz A 3x3 como
A = /sum_{i=1}^N (r_i - r_G)(r_i - r_G)^T
con esta matriz, el vector normal n ahora viene dado por el vector propio de A que corresponde al autovalor mínimo de A.
Para conocer los pares eigenvector / eigenvalue, use cualquier biblioteca de álgebra lineal de su elección.
Esta solución se basa en el teorema de Rayleight-Ritz para la matriz hermitiana A.
CGAL::linear_least_squares_fitting_3
Función linear_least_squares_fitting_3 calcula la línea o plano 3D que mejor se ajusta (en el sentido de mínimos cuadrados) de un conjunto de objetos tridimensionales como puntos, segmentos, triángulos, esferas, bolas, cuboides o tetraedros.