Rozwiązanie najmniejszych kwadratów do równoczesnych równań

Question

Próbuję dopasować transformację z jednego zestawu współrzędnych do drugiego.

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

Istnieje dobrze znana formuła "odręczna" do dopasowania P, Q, R, S do zestawu odpowiednich punktów. Ale muszę mieć oszacowanie błędu na dopasowanie - więc potrzebuję rozwiązania najmniejszych kwadratów.

Przeczytałem "Receptury numeryczne", ale mam problem z ustaleniem, jak to zrobić dla zestawów danych z X i Y w nich.

Czy ktoś może wskazać przykład/samouczek/próbkę kodu, jak to zrobić?
Nie przejmowałem się tym językiem.
Ale - po prostu użyj wbudowanej funkcji Matlab/Lapack/numpy/R prawdopodobnie nie jest pomocna!

edytuj: Mam duży zestaw starych (x, y) nowych (x, y) do dopasowania. Problem jest przesadnie określony (więcej punktów danych niż niewiadomych), więc prosta inwersja matrycy nie jest wystarczająca - i jak powiedziałem, naprawdę potrzebuję błędu w dopasowaniu.

Answer 1

Poniższy kod powinien wystarczyć. Że stosuje się następujący wzór dla reszt:

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

i wyprowadzany wzorach najmniejszych kwadratów na podstawie general procedure opisanego w Wolframa MathWorld.

Przetestowałem ten algorytm w programie Excel i działa on zgodnie z oczekiwaniami. Użyłem zbioru dziesięciu losowych punktów, które następnie zostały obrócone, przetłumaczone i przeskalowane przez losowo wygenerowaną macierz transformacji.

Bez szumu losowego stosowane do danych wyjściowych, program ten produkuje cztery parametry (P, Q, R i S), które są identyczne do parametrów wejściowych i wartości rSquared zero.

Gdy coraz więcej szumów losowych jest stosowanych do punktów wyjściowych, stałe zaczynają odbiegać od prawidłowych wartości, a wartość rSquared odpowiednio wzrasta.

Oto kod:

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

Jedną z kwestii jest to, że takie numeryczne elementy są często trudne. Nawet gdy algorytmy są proste, często pojawiają się problemy, które pojawiają się w rzeczywistych obliczeniach.

Z tego powodu, jeśli istnieje system, który można łatwo uzyskać z wbudowaną funkcją, najlepiej go użyć.

Answer 3

Aby znaleźć P, Q, R i S, możesz użyć najmniejszych kwadratów. Myślę, że mylące jest to, że zwykły opis najmniejszych kwadratów używa x i y, ale nie pasują one do x i y twojego problemu. Trzeba tylko dokładnie przetłumaczyć swój problem na strukturę najmniejszych kwadratów. W twoim przypadku zmiennymi niezależnymi są nieprzetworzone współrzędne x i y, zmienne zależne są przekształconymi współrzędnymi x 'i y', a parametrami zmiennymi są P, Q, R i S. (Jeśli nie jest to wystarczająco jasne, daj mi znać, a ja opublikuję więcej szczegółów.)

Po znalezieniu P, Q, R i S, a następnie skali = sqrt (P^2 + Q^2), a następnie można znaleźć obrót z obrotu sin = Q/skali i obrotu cos = P/skali.

Answer 4

Zdefiniuj macierz 3x3 T (P, Q, R, S) tak, aby (x',y',1) = T (x,y,1). Następnie obliczyć i zminimalizować A w stosunku do (P, Q, R, S).

Samo kodowanie to projekt od średniego do dużego, chyba że można zagwarantować, że dane są dobrze sprawdzone, szczególnie gdy potrzebne są dobre oceny błędów z procedury. Jesteś prawdopodobnie najlepiej wyłączyć za pomocą istniejącego Minimizer który obsługuje szacunków błędach ..

fizyka cząstek typu użyłby minuit bezpośrednio od CERNLIB (z kodowaniem najłatwiej zrobić w fortran77) lub z ROOT (z kodowania w C++ lub powinien być dostępny przez powiązania Pythona). Ale to jest duża instalacja, jeśli nie masz już jednego z tych narzędzi.

Jestem pewien, że inni mogą sugerować inne minimalizatory.

Answer 5

Możesz użyć programu levmar, aby to obliczyć. Został przetestowany i zintegrowany z wieloma produktami, w tym moimi.Jest licencjonowany na licencji GPL, ale jeśli jest to projekt inny niż open source, zmieni licencję dla ciebie (za opłatą)

Answer 6

Dzięki eJames, to jest to prawie exaclty co mam. Zakodowałem ją na podstawie starego podręcznika pomiarowego, który opierał się na wcześniejszej notce "Instrukcji dla geodetów", która musi mieć 100 lat! (Używa N i E dla północy i wschodu zamiast x/y)

Dobroć parametru dopasowania będzie bardzo przydatna - mogę interaktywnie wyrzucać wybrane punkty, jeśli poprawią dopasowanie.

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