Soluzione dei minimi quadrati alle equazioni simultanee

Question

Sto cercando di adattare una trasformazione da un insieme di coordinate a un altro.

x' = R + Px + Qy 
y' = S - Qx + Py 
Where P,Q,R,S are constants, P = scale*cos(rotation). Q=scale*sin(rotation) 

C'è una ben nota formula "per mano" per il montaggio di P, Q, R, S su un insieme di punti corrispondenti. Ma ho bisogno di avere una stima dell'errore sulla misura - quindi ho bisogno di una soluzione dei minimi quadrati.

Leggi 'Numerical Recipes', ma sto avendo problemi a lavorare fuori come fare questo per i set di dati con xey in loro.

può puntare a chiunque di un campione di esempio/tutorial/codice di come fare questo?
Non troppo preoccupato per la lingua.
Ma - usa solo la funzione incorporata di Matlab/Lapack/numpy/R probabilmente non utile!

modifica: Ho un ampio set di vecchi (x, y) nuovi (x, y) per adattarli. Il problema è sovradeterminato (più punti di dati che incognite), quindi l'inversione di matrice non è sufficiente e, come ho detto, ho davvero bisogno dell'errore sull'adattamento.

Answer 1

Il seguente codice dovrebbe fare il trucco. Ho usato la seguente formula per i residui:

residual[i] = (computed_x[i] - actual_x[i])^2 
       + (computed_y[i] - actual_y[i])^2 

E poi derivato il formule minimi quadrati in base alla general procedure descritto in MathWorld di Wolfram.

Ho testato questo algoritmo in Excel e funziona come previsto. Ho usato una raccolta di dieci punti casuali che sono stati poi ruotati, tradotti e ridimensionati da una matrice di trasformazione generata in modo casuale.

Senza rumore casuale applicato ai dati di uscita, questo programma produce quattro parametri (P, Q, R, e S) che sono identici ai parametri di ingresso, e un valore rSquared zero.

Poiché sempre più rumore casuale viene applicato ai punti di uscita, le costanti iniziano ad allontanarsi dai valori corretti e il valore rSquared aumenta di conseguenza.

Ecco il codice: eJames

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

Answer 2

Un problema è che le cose numeriche come questa sono spesso complicate. Anche quando gli algoritmi sono semplici, spesso ci sono problemi che appaiono nel calcolo reale.

Per questo motivo, se c'è un sistema si può ottenere facilmente che ha un built-in funzione, potrebbe essere meglio usare quella.

Answer 3

Per trovare P, Q, R, S e, quindi è possibile utilizzare dei minimi quadrati. Penso che la cosa confusa sia che la solita descrizione dei minimi quadrati usa x e y, ma non corrispondono alla xey nel tuo problema. Hai solo bisogno di tradurre il tuo problema attentamente nel quadro dei minimi quadrati. Nel tuo caso le variabili indipendenti sono le coordinate non trasformate x e y, le variabili dipendenti sono le coordinate trasformate x 'e y', ei parametri regolabili sono P, Q, R e S. (Se questo non è abbastanza chiaro, fammi sapere e posterò maggiori dettagli.)

Una volta trovato P, Q, R e S, quindi ridimensiona = sqrt (P^2 + Q^2) e puoi trovare la rotazione dalla rotazione sin = Q/scala e cos rotazione = P/scala.

Answer 4

definire la matrice 3x3 T (P, Q, R, S) tale che (x',y',1) = T (x,y,1). Poi calcolare

A = \sum_i |(T (x_i,y_i,1)) - (x'_i,y'_i,1)|^2 

e minimizzare A contro (P, Q, R, S).

codifica da soli è un mezzo al progetto di grandi dimensioni a meno che non si può guarntee che i dati sono ben condizionati, soprattutto quando si vogliono ottenere buoni stime di errore fuori della procedura. Sei probabilmente meglio fuori utilizzando un Minimizer esistente che supporta le stime di errore ..

particelle tipo di fisica userebbe minuit sia direttamente dal CERNLIB (con la codifica più facilmente realizzabile in Fortran77), o da ROOT (con la codifica in C++ , o dovrebbe essere accessibile attraverso i collegamenti Python). Ma questa è una grande installazione se non si dispone già di uno di questi strumenti.

sono sicuro che gli altri possono suggerire altre minimi.

Answer 5

È possibile utilizzare il programma levmar per calcolare questo. È testato e integrato in più prodotti, incluso il mio.La sua licenza sotto licenza GPL, ma se questo è un progetto non-opensource, cambierà la licenza per te (a pagamento)

Answer 6

Grazie, questo è quasi exaclty quello che ho. L'ho codificato da un vecchio manuale di rilevamento dell'esercito basato su una precedente "Istruzioni per Geometri", nota che deve essere vecchia di 100 anni! (Usa N ed E per Nord e Est piuttosto che x/y)

Il parametro di bontà dell'adattamento sarà molto utile: posso eliminare i punti selezionati in modo interattivo se peggiorano la misura.

FindTransformation(vector<Point2D> known,vector<Point2D> unknown) { 
{ 
    // sums 
    for (unsigned int ii=0;ii<known.size();ii++) { 
     sum_e += unknown[ii].x; 
     sum_n += unknown[ii].y; 
     sum_E += known[ii].x; 
     sum_N += known[ii].y;        
     ++n;   
    } 

    // mean position 
    me = sum_e/(double)n; 
    mn = sum_n/(double)n; 
    mE = sum_E/(double)n; 
    mN = sum_N/(double)n; 

    // differences 
    for (unsigned int ii=0;ii<known.size();ii++) { 

     de = unknown[ii].x - me; 
     dn = unknown[ii].y - mn; 

     // for P 
     sum_deE += (de*known[ii].x); 
     sum_dnN += (dn*known[ii].y); 
     sum_dee += (de*unknown[ii].x); 
     sum_dnn += (dn*unknown[ii].y); 

     // for Q 
     sum_dnE += (dn*known[ii].x); 
     sum_deN += (de*known[ii].y);      
    } 

double P = (sum_deE + sum_dnN)/(sum_dee + sum_dnn); 
double Q = (sum_dnE - sum_deN)/(sum_dee + sum_dnn); 

double R = mE - (P*me) - (Q*mn); 
double S = mN + (Q*me) - (P*mn); 
}