Plano de mínimos cuadrados 3D

Question

¿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 se podría calcular el mejor plano de ajuste f (x, y) = ax + by + c? ¿Cuál es el algoritmo para obtener a, byc de un conjunto de puntos 3D?

Answer 1

Si tiene n puntos de datos (x [i], y [i], z [i]), calcule la matriz simétrica A de 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 las A y b dadas. 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 " mínimos cuadrados ordinarios " ajuste, que es apropiado solo cuando se espera que z sea una función lineal de x e y. Si está buscando de manera más general un " plano de mejor ajuste " en 3 espacios, es posible que desee aprender sobre " geométrico " mínimos cuadrados.

Tenga en cuenta también que esto fallará si sus puntos están en una línea, como lo están sus puntos de ejemplo.

Answer 2

a menos que alguien me diga cómo escribir ecuaciones aquí, permítame anotar los cálculos finales que tiene que hacer:

primero, los puntos dados 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 de 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 correspondiente al valor propio mínimo de A.

Para obtener información sobre 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 A de Hermitian.

Answer 3

Vea 'Ajuste de datos de mínimos cuadrados' de David Eberly para ver cómo se me ocurrió este para minimizar el ajuste geométrico (distancia ortogonal desde 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 micro segundos ajustando un avión a 1000 puntos (Windows 7, i7, programa de 32 bits)

Answer 4

Al igual que con cualquier enfoque de mínimos cuadrados, proceda así:

Antes de comenzar a codificar

1. Escriba una ecuación para un plano en alguna parametrización, digamos 0 = ax + by + z + d en los parámetros (a, b, d) .

2. Encuentre una expresión D (\ vec {v}; a, b, d) para la distancia desde un punto arbitrario \ vec {v} .

3. 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 componentes de v como \ sigma v_x , \ sigma v_y ^ 2 , \ sigma v_x * v_z ...

4. Escriba las expresiones de minimización por parámetro dS / dx_0 = 0 , dS / dy_0 = 0 ... que le proporciona un conjunto de tres ecuaciones en tres parámetros y las sumas del paso anterior.

5. Resuelva este conjunto de ecuaciones para los parámetros.

(o para casos simples, simplemente busque el formulario). El uso de 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

Tenga en cuenta que si realmente tuviera solo tres puntos, sería mejor simplemente encontrar 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 realizar los pasos 1 y 2, y utilizar un Minimizador de Monte Carlo en la suma del paso 3.

Answer 5

CGAL::linear_least_squares_fitting_3

  

La función linear_least_squares_fitting_3 calcula el mejor ajuste 3D   línea o plano (en el sentido de mínimos cuadrados) de un conjunto de objetos 3D como   como puntos, segmentos, triángulos, esferas, bolas, cuboides o tetraedros.

http://www.cgal.org/Manual/ latest / doc_html / cgal_manual / Principal_component_analysis_ref / Function_linear_least_squares_fitting_3.html

Answer 6

La ecuación para un plano es: ax + by + c = z. 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()

Answer 7

Esto se reduce al problema Total de mínimos cuadrados , que se puede resolver usando 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;
}

Answer 8

Todo lo que tendrá que hacer es resolver el sistema de ecuaciones.

Si esos son tus 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 su pregunta debería ser: ¿Cómo resuelvo un sistema de ecuaciones?

Por lo tanto, recomiendo leer esta pregunta SO.

Si no he entendido bien su pregunta, háganoslo saber.

EDITAR :

Ignora mi respuesta ya que probablemente quisiste decir otra cosa.

Answer 9

Parece que todo lo que quieres 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 algo más.