Solución de mínimos cuadrados para ecuaciones simultáneas

Question

Estoy tratando de ajustar una transformación de un conjunto de coordenadas a otro.

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

Existe una fórmula conocida 'a mano' para ajustar P, Q, R, S a un conjunto de puntos correspondientes. Pero necesito una estimación de error sobre el ajuste, por lo que necesito una solución de mínimos cuadrados.

Lea 'Recetas numéricas', pero tengo problemas para encontrar la forma de hacerlo para conjuntos de datos con xey en ellas.

¿Alguien puede apuntar a una muestra de ejemplo/tutorial/código de cómo hacer esto?
No me molesta demasiado el idioma.
Pero, solo use la función integrada de Matlab/Lapack/numpy/R probablemente no sea útil.

editar: Tengo un conjunto grande de viejo (x, y) nuevo (x, y) para encajar. El problema está sobredeterminado (más puntos de datos que incógnitas) por lo que la inversión de matriz simple no es suficiente, y como dije, realmente necesito el error en el ajuste.

Answer 1

El siguiente código debería ser el truco. He utilizado la siguiente fórmula para los residuos:

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

Y luego derivan las fórmulas de mínimos cuadrados sobre la base de la general procedure descrito en MathWorld de Wolfram.

Probé este algoritmo en Excel y funciona como se esperaba. Usé una colección de diez puntos aleatorios que luego se rotaron, tradujeron y escalaron mediante una matriz de transformación generada aleatoriamente.

sin ruido aleatorio aplicado a los datos de salida, este programa produce cuatro parámetros (P, Q, R, y S) que son idénticos a los parámetros de entrada, y un valor rSquared de cero.

A medida que se aplica más y más ruido aleatorio a los puntos de salida, las constantes comienzan a alejarse de los valores correctos, y el valor rSquared aumenta en consecuencia.

Aquí está el código:

// 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 es que las cosas numéricas como esta suelen ser complicadas. Incluso cuando los algoritmos son simples, a menudo aparecen problemas en el cálculo real.

Por esa razón, si hay un sistema que puede obtener fácilmente que tiene una función incorporada, podría ser mejor usarlo.

Answer 3

Para encontrar P, Q, R y S, puede usar los mínimos cuadrados. Creo que lo confuso es que esa descripción habitual de mínimos cuadrados usa xey, pero no coinciden con la xey en su problema. Solo necesita traducir su problema cuidadosamente en el marco de mínimos cuadrados. En su caso, las variables independientes son las coordenadas no transformadas xey, las variables dependientes son las coordenadas transformadas x 'e y', y los parámetros ajustables son P, Q, R y S. (Si esto no es lo suficientemente claro, hágamelo saber y publicaré más detalles.)

Una vez que haya encontrado P, Q, R y S, entonces scale = sqrt (P^2 + Q^2) y podrá encontrar la rotación de sin rotation = Q/scale y cos rotation = P/scale.

Answer 4

Defina la matriz T 3x3 (P, Q, R, S) tal que (x',y',1) = T (x,y,1). Luego calcule

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

y minimice A contra (P, Q, R, S).

La codificación de esto usted mismo es un proyecto de tamaño mediano a grande, a menos que pueda garantizar que los datos están bien condicionados, especialmente cuando desea obtener buenas estimaciones de error del procedimiento. Usted es probablemente mejor fuera utilizando un minimizador existente que es compatible con las estimaciones de error ..

de partículas de tipo física usaría minuit ya sea directamente desde CERNLIB (con la codificación más fácil de hacer en FORTRAN77), o desde ROOT (con la codificación en C++ , o debería ser accesible a través de los enlaces de python). Pero esa es una gran instalación si ya no tienes una de estas herramientas.

Estoy seguro de que otros pueden sugerir otros minimizadores.

Answer 5

Puede usar el programa levmar para calcular esto. Es probado e integrado en múltiples productos, incluido el mío.Está licenciado bajo la GPL, pero si se trata de un proyecto que no es de código abierto, cambiará la licencia por usted (con una tarifa)

Answer 6

, gracias eJames, eso es casi exaclty lo que tengo. Lo codifiqué de un viejo manual de levantamiento del ejército que estaba basado en una nota anterior de "Instrucciones para topógrafos" que debe tener 100 años de antigüedad. (Utiliza N y E para Norte y Este en lugar de x/y)

El parámetro bondad de ajuste será muy útil: puedo tirar de forma interactiva los puntos seleccionados si empeoran.

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