Der folgende Code sollte den Trick tun. I verwendet, um die folgende Formel für das Residuen:
residual[i] = (computed_x[i] - actual_x[i])^2
+ (computed_y[i] - actual_y[i])^2
und abgeleitete dann die Least-Squares-Formeln auf der Grundlage der bei general procedure Wolframs MathWorld beschrieben.
Ich habe diesen Algorithmus in Excel getestet und es funktioniert wie erwartet. Ich verwendete eine Sammlung von zehn zufälligen Punkten, die dann rotiert, übersetzt und skaliert wurden durch eine zufällig erzeugte Transformationsmatrix.
Ohne zufälligem Rauschen auf die Ausgangsdaten angewendet wird, erzeugt dieses Programm vier Parameter (P, Q, R und S), die mit den Eingangsparametern identisch sind, und einen Wert von Null rSquared.
Da mehr und mehr zufälliges Rauschen auf die Ausgangspunkte angewendet wird, beginnen die Konstanten von den korrekten Werten abzuweichen, und der Wert rSquared erhöht sich entsprechend. Hier
ist der Code:
// test data
const int N = 1000;
float oldPoints_x[N] = { ... };
float oldPoints_y[N] = { ... };
float newPoints_x[N] = { ... };
float newPoints_y[N] = { ... };
// compute various sums and sums of products
// across the entire set of test data
float Ex = Sum(oldPoints_x, N);
float Ey = Sum(oldPoints_y, N);
float Exn = Sum(newPoints_x, N);
float Eyn = Sum(newPoints_y, N);
float Ex2 = SumProduct(oldPoints_x, oldPoints_x, N);
float Ey2 = SumProduct(oldPoints_y, oldPoints_y, N);
float Exxn = SumProduct(oldPoints_x, newPoints_x, N);
float Exyn = SumProduct(oldPoints_x, newPoints_y, N);
float Eyxn = SumProduct(oldPoints_y, newPoints_x, N);
float Eyyn = SumProduct(oldPoints_y, newPoints_y, N);
// compute the transformation constants
// using least-squares regression
float divisor = Ex*Ex + Ey*Ey - N*(Ex2 + Ey2);
float P = (Exn*Ex + Eyn*Ey - N*(Exxn + Eyyn))/divisor;
float Q = (Exn*Ey + Eyn*Ex + N*(Exyn - Eyxn))/divisor;
float R = (Exn - P*Ex - Q*Ey)/N;
float S = (Eyn - P*Ey + Q*Ex)/N;
// compute the rSquared error value
// low values represent a good fit
float rSquared = 0;
float x;
float y;
for (int i = 0; i < N; i++)
{
x = R + P*oldPoints_x[i] + Q*oldPoints_y[i];
y = S - Q*oldPoints_x[i] + P*oldPoints_y[i];
rSquared += (x - newPoints_x[i])^2;
rSquared += (y - newPoints_y[i])^2;
}
Haben Sie eine Menge von (x_i, y_i, x'_i, y'_i) s oder x +/- dx, y +/- dy ... oder was? Wenn Sie jeweils genau eine von x, y, x ', y' haben, können Sie * nur * eine exakte Lösung erstellen, und es gibt keine Möglichkeit, eine Fehlerschätzung zu extrahieren ... – dmckee
Keine der Funktionen vom Minimierungstyp LMA/Gauss- Newton gibt einen direkten Fehler. Ich nehme an, ich könnte die beste Anpassung berechnen und dann den Fehler von jedem Punkt aus berechnen. Ich dachte, es wäre viel einfacher als das (dh. Ein einfacher Modus eines linerar LSquers) und ich war nur dumm –
Interessantes Problem! Ich postete einen Code, der den Trick machen sollte, aber der spaßige Teil bestand darin, meine rostigen alten mathematischen Fähigkeiten wieder auszubrechen :) –