определение дисторсии по точкам

[mrgloom_]

есть набор точек типа такого.


есть теоретическое преобразование которое его описывает, на самом деле может быть и не совсем так.
$u= (m_{11}*x + m_{12}*y + m_{13})/(m_{31}*x + m_{32}*y + m_{33})+c_{1}+x+k_{1}*(x*x+y*y)*x; v= (m_{21}*x + m_{22}*y + m_{23})/(m_{31}*x + m_{32}*y + m_{33})+c_{2}+y+k_{2}*(x*x+y*y)*y;$

так вот надо по точкам определить параметры, это вроде регрессия называется.
но получается у меня есть только точки (u,v), а (х,у) нету и вообще уравнения заданы в параметрическом виде, не очень понятно как с этим работать.
впринципе меня устроила бы аппроксимация полиномом каким либо, но я не уверен, какой он должен быть.

хотелось бы понять как это можно решить и какую библиотеку можно использовать для с++.

[ewert]

mrgloom_ в сообщении #598133 писал(а):

так вот надо по точкам определить параметры, это вроде регрессия называется.

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

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

Модель должна быть, безусловно, полярной, т.е. пересчёту подлежит только расстояние до центра снимка типа $r'=f(r)$, где $r$ -- расстояние до центра на оригинале и $r'$ -- расстояние на пересчитанной фотографии, или наоборот, это не особо важно (надо выбрать вариант, который проще для реализации). Простейшая модель -- кубическая: $f(r)=r(1-\alpha\frac{r^2}{h^2})$, где $h$ -- размер самого снимка. Модель -- однопараметрическая (с параметром $\alpha$; коэффициент при первой степени не имеет значения, т.к. влияет лишь на масштаб и тем самым никак не сказывается на углах). Поэтому минимизацию можно проводить сколь угодно тупым способом, ну хоть делением на три -- это на эффективности в диалоговом режиме практически не скажется.

Вульгарно, конечно, но результат должно дать (хотя честно признаюсь, что сам лично не пробовал; хотя желание такое лет пять назад было, но потом заглохло).

[mrgloom_]

Простейшая модель -- кубическая:
у меня преобразование задано как
$u = c_{1} + x + k_{1}*(x*x + y*y)*x v = c_{2} + y + k_{2}*(x*x + y*y)*y$
и я не понимаю как его свести к вашему уравнению, опять же непонятно какое уравнение использовать для регрессии (у меня есть только набор точек (u,v)).


попытался решить эти уравнения в коэффициентах c1,c2,k1,k2 ,как бы при известных 2-х пар точек (x,y) и (u,v) до преобразования и после.
но почему то при переводе на с++ и последующей прогонке получилось деление на 0 и даже координаты центра отрицательные, чего быть не может.

Цитата:

Solve[v1 == c2 + y1 (1 + k2 (x1^2 + y1^2)) &&
v2 == c2 + y2 (1 + k2 (x2^2 + y2^2)) &&
u1 == c1 + x1 (1 + k1 (x1^2 + y1^2)) &&
u2 == c1 + x2 (1 + k2 (x2^2 + y2^2)), {c1, c2, k1, k2}]

ответ:

Цитата:

{{c1 -> -((-v1 x2^3 + v2 x2^3 + u2 x1^2 y1 - x1^2 x2 y1 + x2^3 y1 +
u2 y1^3 - x2 y1^3 - u2 x2^2 y2 - v1 x2 y2^2 + v2 x2 y2^2 +
x2 y1 y2^2 - u2 y2^3)/(-x1^2 y1 - y1^3 + x2^2 y2 + y2^3)),
c2 -> -((-y1 (x1^2 + y1^2) (v2 - y2) + (v1 - y1) y2 (x2^2 + y2^2))/(
y1 (x1^2 + y1^2) - y2 (x2^2 + y2^2))),
k1 -> -((-u1 + u2 + x1 - x2)/(
x1 (x1^2 + y1^2))) + ((-v1 + y1) (-x2^3 - x2 y2^2))/(
x1 y1 (x1^2 + y1^2)^2) + ((x2^3 +
x2 y2^2) (-y1 (x1^2 + y1^2) (v2 - y2) + (v1 - y1) y2 (x2^2 +
y2^2)))/(
x1 y1 (x1^2 + y1^2)^2 (y1 (x1^2 + y1^2) - y2 (x2^2 + y2^2))),
k2 -> -((-v1 + v2 + y1 - y2)/(x1^2 y1 + y1^3 - x2^2 y2 - y2^3))}}

возможно надо было решать с параметром Reals? там ответ отличается.

[mrgloom_]

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

запихиваем в математику (пробовал 4 разных метода Solve,NSolve и с Reals и без)
так же в математике можно переводить тоже на С, с помощью CForm, только как то не до конца адекватно, но возможно я не до конца разобрался.

Цитата:

CForm[c1 /.
Solve[v1 == c2 + y1 (1 + k2 (x1^2 + y1^2)) &&
v2 == c2 + y2 (1 + k2 (x2^2 + y2^2)) &&
u1 == c1 + x1 (1 + k1 (x1^2 + y1^2)) &&
u2 == c1 + x2 (1 + k2 (x2^2 + y2^2)), {c1, c2, k1, k2}]]
CForm[c2 /.
Solve[v1 == c2 + y1 (1 + k2 (x1^2 + y1^2)) &&
v2 == c2 + y2 (1 + k2 (x2^2 + y2^2)) &&
u1 == c1 + x1 (1 + k1 (x1^2 + y1^2)) &&
u2 == c1 + x2 (1 + k2 (x2^2 + y2^2)), {c1, c2, k1, k2}]]
CForm[k1 /.
Solve[v1 == c2 + y1 (1 + k2 (x1^2 + y1^2)) &&
v2 == c2 + y2 (1 + k2 (x2^2 + y2^2)) &&
u1 == c1 + x1 (1 + k1 (x1^2 + y1^2)) &&
u2 == c1 + x2 (1 + k2 (x2^2 + y2^2)), {c1, c2, k1, k2}]]
CForm[k2 /.
Solve[v1 == c2 + y1 (1 + k2 (x1^2 + y1^2)) &&
v2 == c2 + y2 (1 + k2 (x2^2 + y2^2)) &&
u1 == c1 + x1 (1 + k1 (x1^2 + y1^2)) &&
u2 == c1 + x2 (1 + k2 (x2^2 + y2^2)), {c1, c2, k1, k2}]]

решение в реалз

CForm[c1 /.
Solve[v1 == c2 + y1 (1 + k2 (x1^2 + y1^2)) &&
v2 == c2 + y2 (1 + k2 (x2^2 + y2^2)) &&
u1 == c1 + x1 (1 + k1 (x1^2 + y1^2)) &&
u2 == c1 + x2 (1 + k2 (x2^2 + y2^2)), {c1, c2, k1, k2},Reals]]
CForm[c2 /.
Solve[v1 == c2 + y1 (1 + k2 (x1^2 + y1^2)) &&
v2 == c2 + y2 (1 + k2 (x2^2 + y2^2)) &&
u1 == c1 + x1 (1 + k1 (x1^2 + y1^2)) &&
u2 == c1 + x2 (1 + k2 (x2^2 + y2^2)), {c1, c2, k1, k2},Reals]]
CForm[k1 /.
Solve[v1 == c2 + y1 (1 + k2 (x1^2 + y1^2)) &&
v2 == c2 + y2 (1 + k2 (x2^2 + y2^2)) &&
u1 == c1 + x1 (1 + k1 (x1^2 + y1^2)) &&
u2 == c1 + x2 (1 + k2 (x2^2 + y2^2)), {c1, c2, k1, k2},Reals]]
CForm[k2 /.
Solve[v1 == c2 + y1 (1 + k2 (x1^2 + y1^2)) &&
v2 == c2 + y2 (1 + k2 (x2^2 + y2^2)) &&
u1 == c1 + x1 (1 + k1 (x1^2 + y1^2)) &&
u2 == c1 + x2 (1 + k2 (x2^2 + y2^2)), {c1, c2, k1, k2},Reals]]

Еще и численный метод решения:
CForm[c1 /.
NSolve[v1 == c2 + y1 (1 + k2 (x1^2 + y1^2)) &&
v2 == c2 + y2 (1 + k2 (x2^2 + y2^2)) &&
u1 == c1 + x1 (1 + k1 (x1^2 + y1^2)) &&
u2 == c1 + x2 (1 + k2 (x2^2 + y2^2)), {c1, c2, k1, k2}]]
CForm[c2 /.
NSolve[v1 == c2 + y1 (1 + k2 (x1^2 + y1^2)) &&
v2 == c2 + y2 (1 + k2 (x2^2 + y2^2)) &&
u1 == c1 + x1 (1 + k1 (x1^2 + y1^2)) &&
u2 == c1 + x2 (1 + k2 (x2^2 + y2^2)), {c1, c2, k1, k2}]]
CForm[k1 /.
NSolve[v1 == c2 + y1 (1 + k2 (x1^2 + y1^2)) &&
v2 == c2 + y2 (1 + k2 (x2^2 + y2^2)) &&
u1 == c1 + x1 (1 + k1 (x1^2 + y1^2)) &&
u2 == c1 + x2 (1 + k2 (x2^2 + y2^2)), {c1, c2, k1, k2}]]
CForm[k2 /.
NSolve[v1 == c2 + y1 (1 + k2 (x1^2 + y1^2)) &&
v2 == c2 + y2 (1 + k2 (x2^2 + y2^2)) &&
u1 == c1 + x1 (1 + k1 (x1^2 + y1^2)) &&
u2 == c1 + x2 (1 + k2 (x2^2 + y2^2)), {c1, c2, k1, k2}]]

тоже самое только на реалз
CForm[c1 /.
NSolve[v1 == c2 + y1 (1 + k2 (x1^2 + y1^2)) &&
v2 == c2 + y2 (1 + k2 (x2^2 + y2^2)) &&
u1 == c1 + x1 (1 + k1 (x1^2 + y1^2)) &&
u2 == c1 + x2 (1 + k2 (x2^2 + y2^2)), {c1, c2, k1, k2},Reals]]
CForm[c2 /.
NSolve[v1 == c2 + y1 (1 + k2 (x1^2 + y1^2)) &&
v2 == c2 + y2 (1 + k2 (x2^2 + y2^2)) &&
u1 == c1 + x1 (1 + k1 (x1^2 + y1^2)) &&
u2 == c1 + x2 (1 + k2 (x2^2 + y2^2)), {c1, c2, k1, k2},Reals]]
CForm[k1 /.
NSolve[v1 == c2 + y1 (1 + k2 (x1^2 + y1^2)) &&
v2 == c2 + y2 (1 + k2 (x2^2 + y2^2)) &&
u1 == c1 + x1 (1 + k1 (x1^2 + y1^2)) &&
u2 == c1 + x2 (1 + k2 (x2^2 + y2^2)), {c1, c2, k1, k2},Reals]]
CForm[k2 /.
NSolve[v1 == c2 + y1 (1 + k2 (x1^2 + y1^2)) &&
v2 == c2 + y2 (1 + k2 (x2^2 + y2^2)) &&
u1 == c1 + x1 (1 + k1 (x1^2 + y1^2)) &&
u2 == c1 + x2 (1 + k2 (x2^2 + y2^2)), {c1, c2, k1, k2},Reals]]

далее собственно пробуем восстановить параметры по формуле экспортировав ответы на С.

void barrel_parameters(double x1,double y1,double x2,double y2,double u1,double v1,double u2,double v2)
{
   //Solve
   double c1= -((-(v1*pow(x2,3)) + v2*pow(x2,3) + u2*pow(x1,2)*y1 - pow(x1,2)*x2*y1 + pow(x2,3)*y1 + u2*pow(y1,3) - x2*pow(y1,3) - u2*pow(x2,2)*y2 - v1*x2*pow(y2,2) + 
        v2*x2*pow(y2,2) + x2*y1*pow(y2,2) - u2*pow(y2,3))/(-(pow(x1,2)*y1) - pow(y1,3) + pow(x2,2)*y2 + pow(y2,3)));
   double c2= -((-(y1*(pow(x1,2) + pow(y1,2))*(v2 - y2)) + (v1 - y1)*y2*(pow(x2,2) + pow(y2,2)))/(y1*(pow(x1,2) + pow(y1,2)) - y2*(pow(x2,2) + pow(y2,2))));
   double k1= -((-u1 + u2 + x1 - x2)/(x1*(pow(x1,2) + pow(y1,2)))) + ((-v1 + y1)*(-pow(x2,3) - x2*pow(y2,2)))/(x1*y1*pow(pow(x1,2) + pow(y1,2),2)) + 
      ((pow(x2,3) + x2*pow(y2,2))*(-(y1*(pow(x1,2) + pow(y1,2))*(v2 - y2)) + (v1 - y1)*y2*(pow(x2,2) + pow(y2,2))))/
       (x1*y1*pow(pow(x1,2) + pow(y1,2),2)*(y1*(pow(x1,2) + pow(y1,2)) - y2*(pow(x2,2) + pow(y2,2))));
   double k2= -((-v1 + v2 + y1 - y2)/(pow(x1,2)*y1 + pow(y1,3) - pow(x2,2)*y2 - pow(y2,3)));

   //Solve reals
   double с1_r= u1 - x1 - (pow(x1,3)*(-u1 + u2 + x1 - x2 - (pow(x2,3)*(-v1 + v2 + y1 - y2))/(-(pow(x1,2)*y1) - pow(y1,3) + pow(x2,2)*y2 + pow(y2,3)) - 
          (x2*(-v1 + v2 + y1 - y2)*pow(y2,2))/(-(pow(x1,2)*y1) - pow(y1,3) + pow(x2,2)*y2 + pow(y2,3))))/(-pow(x1,3) - x1*pow(y1,2)) - 
     (x1*pow(y1,2)*(-u1 + u2 + x1 - x2 - (pow(x2,3)*(-v1 + v2 + y1 - y2))/(-(pow(x1,2)*y1) - pow(y1,3) + pow(x2,2)*y2 + pow(y2,3)) - 
          (x2*(-v1 + v2 + y1 - y2)*pow(y2,2))/(-(pow(x1,2)*y1) - pow(y1,3) + pow(x2,2)*y2 + pow(y2,3))))/(-pow(x1,3) - x1*pow(y1,2));
   double c2_r= v1 - y1 - (pow(x1,2)*y1*(-v1 + v2 + y1 - y2))/(-(pow(x1,2)*y1) - pow(y1,3) + pow(x2,2)*y2 + pow(y2,3)) - 
       (pow(y1,3)*(-v1 + v2 + y1 - y2))/(-(pow(x1,2)*y1) - pow(y1,3) + pow(x2,2)*y2 + pow(y2,3));
   double k1_r= (-u1 + u2 + x1 - x2 - (pow(x2,3)*(-v1 + v2 + y1 - y2))/(-(pow(x1,2)*y1) - pow(y1,3) + pow(x2,2)*y2 + pow(y2,3)) - 
         (x2*(-v1 + v2 + y1 - y2)*pow(y2,2))/(-(pow(x1,2)*y1) - pow(y1,3) + pow(x2,2)*y2 + pow(y2,3)))/(-pow(x1,3) - x1*pow(y1,2));
   double k2_r= (-v1 + v2 + y1 - y2)/(-(pow(x1,2)*y1) - pow(y1,3) + pow(x2,2)*y2 + pow(y2,3));

   //NSolve
   double c1_n= 1.*u2 - 1.*x2 + (1.*(1.*(v1 - 1.*y1) - 1.*(v2 - 1.*y2))*(1.*pow(x2,3) + 1.*x2*pow(y2,2)))/(-1.*y1*(pow(x1,2) + pow(y1,2)) + 1.*y2*(pow(x2,2) + pow(y2,2)));
   double c2_n= 1.*v1 - 1.*y1 + (1.*(1.*pow(x1,2)*y1 + 1.*pow(y1,3))*(1.*(v1 - 1.*y1) - 1.*(v2 - 1.*y2)))/(-1.*y1*(pow(x1,2) + pow(y1,2)) + 1.*y2*(pow(x2,2) + pow(y2,2)));
   double k1_n= -((-1.*u1 + 1.*u2 + 1.*x1 - 1.*x2)/(x1*(pow(x1,2) + pow(y1,2)))) + 
      (1.*(1.*(v1 - 1.*y1) - 1.*(v2 - 1.*y2))*(-1.*pow(x2,3) - 1.*x2*pow(y2,2)))/
       (x1*(pow(x1,2) + pow(y1,2))*(-1.*y1*(pow(x1,2) + pow(y1,2)) + 1.*y2*(pow(x2,2) + pow(y2,2))));
   double k2_n= -((1.*(v1 - 1.*y1) - 1.*(v2 - 1.*y2))/(-1.*y1*(pow(x1,2) + pow(y1,2)) + 1.*y2*(pow(x2,2) + pow(y2,2))));

   //NSolve reals
   double c1_nr= u1 - 1.*x1 - (1.*pow(x1,3)*(-1.*u1 + u2 + x1 - 1.*x2 - 
          (1.*pow(x2,3)*(-1.*v1 + v2 + y1 - 1.*y2))/(-1.*pow(x1,2)*y1 - 1.*pow(y1,3) + pow(x2,2)*y2 + pow(y2,3)) - 
          (1.*x2*(-1.*v1 + v2 + y1 - 1.*y2)*pow(y2,2))/(-1.*pow(x1,2)*y1 - 1.*pow(y1,3) + pow(x2,2)*y2 + pow(y2,3))))/(-1.*pow(x1,3) - 1.*x1*pow(y1,2)) - 
     (1.*x1*pow(y1,2)*(-1.*u1 + u2 + x1 - 1.*x2 - (1.*pow(x2,3)*(-1.*v1 + v2 + y1 - 1.*y2))/(-1.*pow(x1,2)*y1 - 1.*pow(y1,3) + pow(x2,2)*y2 + pow(y2,3)) - 
          (1.*x2*(-1.*v1 + v2 + y1 - 1.*y2)*pow(y2,2))/(-1.*pow(x1,2)*y1 - 1.*pow(y1,3) + pow(x2,2)*y2 + pow(y2,3))))/(-1.*pow(x1,3) - 1.*x1*pow(y1,2));
   double c2_nr= v1 - 1.*y1 - (1.*pow(x1,2)*y1*(-1.*v1 + v2 + y1 - 1.*y2))/(-1.*pow(x1,2)*y1 - 1.*pow(y1,3) + pow(x2,2)*y2 + pow(y2,3)) - 
       (1.*pow(y1,3)*(-1.*v1 + v2 + y1 - 1.*y2))/(-1.*pow(x1,2)*y1 - 1.*pow(y1,3) + pow(x2,2)*y2 + pow(y2,3));
   double k1_nr= (-1.*u1 + u2 + x1 - 1.*x2 - (1.*pow(x2,3)*(-1.*v1 + v2 + y1 - 1.*y2))/(-1.*pow(x1,2)*y1 - 1.*pow(y1,3) + pow(x2,2)*y2 + pow(y2,3)) - 
         (1.*x2*(-1.*v1 + v2 + y1 - 1.*y2)*pow(y2,2))/(-1.*pow(x1,2)*y1 - 1.*pow(y1,3) + pow(x2,2)*y2 + pow(y2,3)))/(-1.*pow(x1,3) - 1.*x1*pow(y1,2));
   double k2_nr= (-1.*v1 + v2 + y1 - 1.*y2)/(-1.*pow(x1,2)*y1 - 1.*pow(y1,3) + pow(x2,2)*y2 + pow(y2,3));


//тестируем
   IplImage* src=cvLoadImage("сетка.PNG");//чистая сетка
   cvSaveImage("test_barrel.png",barrel_pincusion_dist(src, c1,c2,k1,k2));
}

для теста использовал картинки



вызываю

barrel_parameters(3,2,284,283,26.4,25.2,259,259);

как моделирую само искажение

IplImage* barrel_pincusion_dist(IplImage* img, int Cx,int Cy,double kx,double ky)
{
   IplImage* mapx = cvCreateImage( cvGetSize(img), IPL_DEPTH_32F, 1 );
   IplImage* mapy = cvCreateImage( cvGetSize(img), IPL_DEPTH_32F, 1 );

   int w= img->width;
   int h= img->height;

   float* pbuf = (float*)mapx->imageData;
   for (int y = 0; y < h; y++)
   {
      int ty= y-Cy;
      for (int x = 0; x < w; x++)
      {
         int tx= x-Cx;
         int rt= tx*tx+ty*ty;
                        
         *pbuf = (float)(tx*(1+kx*rt)+Cx);
         ++pbuf;
      }
   }

   pbuf = (float*)mapy->imageData;
   for (int y = 0;y < h; y++)
   {
      int ty= y-Cy;
      for (int x = 0; x < w; x++) 
      {
         int tx= x-Cx;
         int rt= tx*tx+ty*ty;

         *pbuf = (float)(ty*(1+ky*rt)+Cy);
         ++pbuf;
      }
   }

   IplImage* temp = cvCloneImage(img);
    cvRemap( temp, img, mapx, mapy ); 
   cvReleaseImage(&temp);
   cvReleaseImage(&mapx);
   cvReleaseImage(&mapy);

   return img;
}

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

[mrgloom_]

кстати я понял как получить параметры искажения имея только 1 дугу, учитывая, что эта дуга была до искажения прямой.
для этого потребуется 6 точек лежащих на дуге.

u1= c1+x1(1+k1(x1^2+y1^2))
v1= c2+y1(1+k2(x1^2+y1^2))
y1=k*x1+b

повторяем эти уравнения 6 раз для разных точек и получаем 18 уравнений из которых можем получить x,y точки и параметры дисторсии.
вот только как решить потом 18 уравнений?

а если мы знаем где находится центр, обычно в центре изображения, т.е. знаем с1,с2, то всё еще проще.

u1= c1+x1(1+k1(x1^2+y1^2))
v1= c2+y1(1+k2(x1^2+y1^2))
u2= c1+x2(1+k1(x2^2+y2^2))
v2= c2+y2(1+k2(x2^2+y2^2))
x1=u1
y2=v2