3D-Least-Squares-Flugzeug

Question

Was ist der Algorithmus zur Berechnung einer Least-Squares-Ebene in (x, y, z) Raum, da ein Satz von 3D-Datenpunkten? Mit anderen Worten, wenn hatte ich eine Reihe von Punkten, wie (1, 2, 3), (4, 5, 6), (7, 8, 9) usw., wie würde ein über die beste Passform f die Berechnung Ebene (x, y) = ax + by + c? Was ist der Algorithmus a, b für immer, und c aus einem Satz von 3D-Punkten?

Answer 1

Wenn Sie n Datenpunkte (x [i], Y [i], z [i]), berechne die 3x3 symmetrische Matrix A, deren Einträge sind:

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

Auch berechnen die 3-Element-Vektor b:

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

Dann Ax = b für die gegebenen A und b lösen. Die drei Komponenten des Lösungsvektors sind die Koeffizienten der Least-Square-Fit-Ebene {a, b, c}.

Beachten Sie, dass dies die „ordinary least squares“ geeignet ist, die nur geeignet ist, wenn z erwartet wird, eine lineare Funktion von x und y zu sein. Wenn Sie ganz allgemein für eine „best fit Ebene“ in 3-Raum suchen, möchten Sie vielleicht über „geometrische“ kleinsten Quadrate lernen.

Beachten Sie auch, dass dies fehlschlagen, wenn Sie Ihre Punkte in einer Linie sind, als Beispiel Punkte sind.

Answer 2

es sei denn, jemand mir sagt, wie Gleichungen hier zu geben, lassen Sie mich aufschreiben nur die letzten Berechnungen was Sie tun müssen:

ersten, gegebenen Punkte r_i \ n \ r, i = 1..N, die Berechnung der Massenmittelpunkt aller Punkte:

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

Dann berechnet die Normalvektor n, das zusammen mit dem Basisvektor r_G definiert die Ebene, die durch die 3x3-Matrix A als Berechnung

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

mit dieser Matrix, der Normalenvektor n nun von dem Eigenvektor von A entsprechend den minimalen Eigenwert von A gegeben.

Um herauszufinden, über den Eigenvektor / Eigenwert Paare sind, kann jede lineare Algebra Bibliothek Ihrer Wahl.

Diese Lösung basiert auf dem Rayleight-Ritz Satz für die hermitesche Matrix A.

Answer 3

Siehe ‚Least Squares Fitting von Daten‘ von David Eberly, wie ich mit diesem kam die geometrische Anpassung (orthogonalen Abstand von Punkten auf die Ebene) zu minimieren.

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

Timed bei 37 Mikrosekunden Anpassen einer Ebene auf 1000 Punkte (Windows 7, i7, 32-Bit-Programm)

Answer 4

Wie bei jedem kleinsten Quadrate Ansatz, gehen Sie wie folgt aus:

Bevor Sie beginnen Codierung

1. Schreiben Sie eine Gleichung für eine Ebene in einem gewissen Parametrisierung nach unten, sagen 0 = ax + by + z + d in dir Parameter (a, b, d).

2. ein Ausdruck D(\vec{v};a, b, d) für den Abstand von einem beliebigen Punkt \vec{v} finden.

3. Schreiben Sie die Summe S = \sigma_i=0,n D^2(\vec{x}_i) nach unten, und vereinfachen, bis es im Hinblick auf die einfache Summe der Komponenten von v wie \sigma v_x, \sigma v_y^2 ausgedrückt wird, \sigma v_x*v_z ...

4. Notieren Sie sich den pro Parameter Minimierung Ausdrücke dS/dx_0 = 0, dS/dy_0 = 0 ... die Ihnen einen Satz von drei Gleichungen in drei Parametern und die Summen aus dem vorherigen Schritt.

gibt

5. Lösen Sie diesen Satz von Gleichungen für die Parameter.

(oder für einfache Fälle, schauen Sie einfach das Formular oben). ein symbolisches Algebra-Paket (wie Mathematica) verwendet, könnte man das Leben viel einfacher machen.

Die Kodierung

• Code über die benötigten Summen zu bilden, und die Parameter aus dem letzten Satz oben finden.

Alternativen

Beachten Sie, dass, wenn Sie tatsächlich nur drei Punkte haben, man wäre besser, nur das Flugzeug zu finden, die durch sie geht.

Auch wenn die analytische Lösung in undurchführbar (nicht der Fall für ein Flugzeug, aber generell möglich) Sie tun können, die Schritte 1 und 2, und verwenden Sie ein Monte Carlo minimizer in Schritt auf die Summe 3.

Answer 5

CGAL::linear_least_squares_fitting_3

  

Funktion linear_least_squares_fitting_3 berechnet die am besten passenden 3D   Linie oder Ebene (im Sinne der kleinsten Quadrate) von einem Satz von 3D-Objekten, wie beispiel   als Punkte, Segmente, Dreiecke, Kugeln, Kugeln, Quader oder Tetraedern.

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

Answer 6

Die Gleichung für eine Ebene ist: ax + by + c = z. So gründen Matrizen wie dies mit allen Daten:

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

Und

    a  
x = b  
    c

Und

    z_0   
B = z_1   
    ...   
    z_n

Mit anderen Worten: Ax = B. Jetzt für x lösen, die Ihre Koeffizienten sind. Aber da (ich nehme an) Sie mehr als 3 Punkte haben, ist das System überbestimmt, so dass Sie die linke Pseudo inverse verwenden müssen. So lautet die Antwort:

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

Und hier ist ein einfacher Python-Code mit einem Beispiel:

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

Dies reduziert das Total Least Squares Problem, das unter Verwendung gelöst werden kann SVD Zersetzung.

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

Alles, was Sie tun müssen, ist das System von Gleichungen zu lösen.

Wenn das sind Ihre Punkte: (1, 2, 3), (4, 5, 6), (7, 8, 9)

Das gibt Ihnen die Gleichungen:

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

So Ihre Frage eigentlich sein sollte: Wie löse ich ein Gleichungssystem

Deshalb empfehle ich diese SO Frage.

Wenn ich das falsch verstanden habe Ihre Frage lassen Sie uns wissen.

Bearbeiten :

Ignorieren meiner Antwort, wie Sie wahrscheinlich etwas anderes gemeint.

Answer 9

Es klingt wie alles, was Sie tun wollen lineare Regression mit 2 Regressoren ist. Die Wikipedia-Seite zu diesem Thema sollten Sie alles, was Sie wissen müssen, sagen, und dann einige.