Решение наименьших квадратов для уравнений

Question

Я пытаюсь подогнать преобразование из одного набора координат в другой.

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

Существует хорошо известная формула «от руки» для подгонки P, Q, R, S к набору соответствующих точек. Но мне нужно оценить погрешность подбора, поэтому мне нужно решение наименьших квадратов.

Прочитайте 'Числовые рецепты', но у меня возникли проблемы при разработке, как это сделать для наборов данных с x и y в них.

Может кто-нибудь указать на пример / учебник / пример кода, как это сделать?
Не слишком беспокоюсь о языке.
Но - просто использовать встроенную функцию Matlab / Lapack / numpy / R, вероятно, не поможет!

редактирование: У меня есть большой набор старых (х, у) новых (х, у), чтобы соответствовать. Проблема переопределена (больше точек данных, чем неизвестных), поэтому простой инверсии матрицы недостаточно - и, как я сказал, мне действительно нужна ошибка при подборе.

Answer 1

Чтобы найти P, Q, R и S, вы можете использовать наименьшие квадраты. Я думаю, что сбивает с толку то, что в обычном описании наименьших квадратов используются x и y, но они не соответствуют x и y в вашей задаче. Вам просто нужно аккуратно перевести свою задачу в структуру наименьших квадратов. В вашем случае независимыми переменными являются нетрансформированные координаты x и y, зависимыми переменными являются преобразованные координаты x 'и y', а настраиваемыми параметрами являются P, Q, R и S. (Если это не достаточно ясно, дайте мне знать, и я опубликую более подробно.)

Как только вы нашли P, Q, R и S, тогда scale = sqrt (P ^ 2 + Q ^ 2), и затем вы можете найти вращение от sin вращение = Q / scale и cos вращение = P / шкала.

Answer 2

Следующий код должен помочь. Я использовал следующую формулу для остатков:

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

А затем вывели формулы наименьших квадратов на основе общей процедуры , описанной в Wolfram's MathWorld.

Я проверил этот алгоритм в Excel, и он работает как положено. Я использовал набор из десяти случайных точек, которые затем поворачивались, переводились и масштабировались с помощью случайно сгенерированной матрицы преобразования.

Без случайного шума, применяемого к выходным данным, эта программа создает четыре параметра (P, Q, R и S), которые идентичны входным параметрам, и значение rSquared нуля.

Поскольку к выходным точкам применяется все больше и больше случайных шумов, константы начинают отклоняться от правильных значений, и значение rSquared соответственно увеличивается.

Вот код:

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

Спасибо, eJames, это почти точно, что у меня есть. Я закодировал его из старого армейского руководства по геодезии, основанного на более ранней заметке «Инструкции для геодезистов», которой должно быть 100 лет! (Он использует N и E для севера и востока, а не х / у)

Параметр качества подгонки будет очень полезен - я могу в интерактивном режиме выбрасывать выбранные точки, если они ухудшают подгонку.

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

Answer 4

Вы можете использовать программу levmar для вычисления этого. Он протестирован и интегрирован в несколько продуктов, включая мой. Он лицензирован по лицензии GPL, но если это не проект с открытым исходным кодом, он заменит лицензию для вас (за плату)

Answer 5

Определите матрицу 3x3 T (P, Q, R, S) так, чтобы (x',y',1) = T (x,y,1). Затем вычислите

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

и минимизировать А против (P, Q, R, S).

Самостоятельное кодирование - это проект среднего и крупного размера, если только вы не можете гарантировать, что данные хорошо подготовлены, особенно если вы хотите получить хорошие оценки ошибок вне процедуры. Вам, вероятно, лучше всего использовать существующий минимизатор, который поддерживает оценки ошибок ..

Тип физики элементарных частиц будет использовать minuit либо непосредственно из CERNLIB (с кодировкой, проще всего сделать в fortran77), либо из ROOT (с кодировкой в c ++, или он должен быть доступен через привязки Python). Но это большая установка, если у вас еще нет одного из этих инструментов.

Я уверен, что другие могут предложить другие минимизаторы.

Answer 6

Одна из проблем заключается в том, что такие числовые вещи часто бывают хитрыми. Даже когда алгоритмы просты, в реальных вычислениях часто возникают проблемы.

По этой причине, если есть система, которую вы легко можете получить и которая имеет встроенную функцию, лучше использовать ее.