Решение методом наименьших квадратов к одновременным уравнениям
Ваш материал будет оттянут в точном порядке, в котором Вы вызываете glBegin/glEnd функции. Можно получить буферизующее глубину использование z-буфера, и если 2-е объекты имеют различные значения z, можно получить эффект, Вы хотите тот путь. Единственным путем Вы видите поведение, которое Вы описываете на Mac, то, если программа тянет материал в наоборот порядке вручную или использует z-буфер для выполнения этого. OpenGL иначе не имеет никакой функциональности автоматически, как Вы описываете.
6 ответов
Следующие код должен сделать свое дело. Я использовал следующую формулу для остатков:
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 .
As все больше и больше случайного шума применяется к выходным точкам, константы начинают отклоняться от правильных значений, и значение 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;
}
Вы можете использовать программу levmar , чтобы вычислить это. Он протестирован и интегрирован во многие продукты, включая мой. Он под лицензией GPL, но если это проект без открытого исходного кода, он изменит лицензию для вас (за определенную плату)
Чтобы найти 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 / шкала и вращение cos = P / шкала.
Определите матрицу 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
и минимизируйте A против (P, Q, R, S).
Самостоятельное кодирование - это проект среднего и большого размера, если вы не можете гарантировать, что данные хорошо согласованы, особенно когда вам нужны хорошие оценки ошибок вне процедуры. Вам, вероятно, лучше всего использовать существующий минимизатор, который поддерживает оценки ошибок.
Тип физики элементарных частиц будет использовать minuit либо непосредственно из CERNLIB (кодирование проще всего выполнить в fortran77) ) или из ROOT (с кодировкой на c ++, или он должен быть доступен через привязки python). Но это большая установка, если у вас еще нет одного из этих инструментов.
Я уверен, что другие могут предложить другие минимизаторы.
Спасибо, eJames, это почти то, что у меня есть. Я закодировал его из старого армейского руководства по геодезии, основанного на более ранней записке из «Инструкций геодезистам», которой должно быть 100 лет! (Он использует N и E для севера и востока, а не x / y)
Параметр согласия будет очень полезен - я могу интерактивно выкинуть выбранные точки, если они ухудшают соответствие.
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);
}
Одна из проблем заключается в том, что подобные числовые данные часто являются сложными. Даже когда алгоритмы просты, есть