Plan des moindres carrés 3D

Question

Quel est l'algorithme permettant de calculer un plan des moindres carrés dans l'espace (x, y, z), à partir d'un ensemble de points de données 3D? En d’autres termes, si j’avais plusieurs points comme (1, 2, 3), (4, 5, 6), (7, 8, 9), etc., comment pourrions-nous calculer le meilleur plan d’ajustement f (x, y) = ax + by + c? Quel est l'algorithme permettant d'extraire a, b et c d'un ensemble de points 3D?

Answer 1

Si vous avez n points de données (x [i], y [i], z [i]), calculez la matrice symétrique 3x3 A dont les entrées sont:

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

Calculez également le vecteur des 3 éléments b:

{sum_i x[i]*z[i],   sum_i y[i]*z[i],    sum_i z[i]}

Puis résolvez Ax = b pour A et b donnés. Les trois composantes du vecteur solution sont les coefficients du plan d’ajustement des moindres carrés {a, b, c}.

Notez qu'il s'agit des "moindres carrés ordinaires". fit, qui ne convient que lorsque z est censé être une fonction linéaire de x et y. Si vous recherchez plus généralement un "avion parfaitement adapté", sur 3-espace, vous voudrez peut-être en savoir plus sur " géométrique " moindres carrés.

Notez également que cela échouera si vos points sont alignés, comme le sont vos exemples de points.

Answer 2

À moins que quelqu'un ne me dise comment taper des équations ici, laissez-moi simplement écrire les calculs finaux que vous devez faire:

d’abord, en fonction des points r_i \ n \ R, i = 1..N, calculez le centre de gravité de tous les points:

r_G = \frac{\sum_{i=1}^N r_i}{N}

puis, calculez le vecteur normal n qui, avec le vecteur de base r_G, définit le plan en calculant la matrice 3x3 A comme

A = \sum_{i=1}^N (r_i - r_G)(r_i - r_G)^T

avec cette matrice, le vecteur normal n est maintenant donné par le vecteur propre de A correspondant à la valeur propre minimale de A.

Pour en savoir plus sur les paires vecteur propre / valeur propre, utilisez n'importe quelle bibliothèque d'algèbre linéaire de votre choix.

Cette solution est basée sur le théorème de Rayleight-Ritz pour la matrice hermitienne A.

Answer 3

Voir "Ajustement des données par les moindres carrés" de David Eberly pour connaître la méthode proposée pour minimiser l'ajustement géométrique (distance orthogonale des points au plan).

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;
}

Programmé à 37 micro secondes pour ajuster un avion à 1 000 points (programme Windows 7, i7, 32 bits)

Answer 4

Comme pour toute approche par les moindres carrés, procédez comme suit:

Avant de commencer à coder

1. Ecrivez une équation pour un plan dans une paramétrisation, par exemple 0 = ax + de + z + d dans les paramètres (a, b, d) .

2. Rechercher une expression D (\ vec {v}; a, b, d) pour la distance par rapport à un point arbitraire \ vec {v} .

3. Notez la somme S = \ sigma_i = 0, n D ^ 2 (\ vec {x} _i) et simplifiez-la jusqu'à ce qu'elle soit exprimée en simples sommes du composants de v comme \ sigma v_x , \ sigma v_y ^ 2 , \ sigma v_x * v_z ...

4. Notez les expressions de minimisation par paramètre dS / dx_0 = 0 , dS / dy_0 = 0 ... ce qui vous donne un ensemble de trois équations dans trois paramètres et les sommes de l'étape précédente.

5. Résolvez cet ensemble d'équations pour les paramètres.

(ou pour les cas simples, il suffit de regarder le formulaire). Utiliser un paquet algèbre symbolique (comme Mathematica) pourrait vous rendre la vie beaucoup plus facile.

Le codage

• Écrivez le code pour former les sommes nécessaires et trouver les paramètres du dernier jeu ci-dessus.

Alternatives

Notez que si vous n'aviez que trois points, vous feriez mieux de trouver l'avion qui les traverse.

En outre, si la solution analytique est irréalisable (ce qui n'est pas le cas pour un avion, mais possible en général), vous pouvez effectuer les étapes 1 et 2 et utiliser un minimiseur de Monte Carlo sur la somme de l'étape 3.

Answer 5

CGAL::linear_least_squares_fitting_3

  

La fonction linear_least_squares_fitting_3 calcule la 3D la mieux ajustée   ligne ou plan (au sens des moindres carrés) d'un ensemble d'objets 3D tels que   sous forme de points, de segments, de triangles, de sphères, de billes, de cuboïdes ou de tétraèdres.

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

Answer 6

L’équation d’un plan est la suivante: ax + de + c = z. Donc, configurez des matrices comme celle-ci avec toutes vos données:

    x_0   y_0   1  
A = x_1   y_1   1  
          ... 
    x_n   y_n   1  

Et

    a  
x = b  
    c

Et

    z_0   
B = z_1   
    ...   
    z_n

En d'autres termes: Ax = B. Maintenant, résolvez pour x quels sont vos coefficients. Mais puisque (je suppose) que vous avez plus de 3 points, le système est surdéterminé, vous devez donc utiliser le pseudo-inverse de gauche. Donc, la réponse est:

a 
b = (A^T A)^-1 A^T B
c

Et voici un code Python simple avec un exemple:

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

Cela réduit le problème Total des moindres carrés , qui peut être résolu à l'aide de décomposition SVD .

Code C ++ utilisant 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

Tout ce que vous avez à faire est de résoudre le système d'équations.

Si ce sont vos arguments: (1, 2, 3), (4, 5, 6), (7, 8, 9)

Cela vous donne les équations:

3=a*1 + b*2 + c
6=a*4 + b*5 + c
9=a*7 + b*8 + c

Votre question devrait donc être: comment résoudre un système d'équations?

Par conséquent, je vous recommande de lire la question sur les SO.

.

Si j'ai mal compris votre question, faites-le nous savoir.

MODIFIER :

Ignorez ma réponse car vous vouliez probablement dire autre chose.

Answer 9

On dirait que tout ce que vous voulez faire est une régression linéaire avec 2 régresseurs. La page wikipedia sur le sujet devrait vous indiquer tout ce que vous devez savoir et plus encore.