Перейти к содержимому
Compvision.ru
mrgloom

определение дисторсии и преобразования по точкам

Recommended Posts

До этого был метод простых итераций

Здесь приведен Метод Ньютона:

Спасибо за ссылку сам пришел к этим методом. Буду изучать.

Поделиться сообщением


Ссылка на сообщение
Поделиться на других сайтах

а почему в одном случае работает, а в другом нет?

кстати если поставить

double Cx=100.0;

double Cy=120.0;

double x= 100;

double y= 120;

то выдаст

x=-1.#IND

Поделиться сообщением


Ссылка на сообщение
Поделиться на других сайтах

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

на деле же получается, что u,v = x,y с картинки?

пробовал так

IplImage* un_barrel_pincusion_dist(IplImage* img, int c1,int c2,double k1,double k2)

{

	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_x = (float*)mapx->imageData;

	float* pbuf_y = (float*)mapy->imageData;

	for (int y = 0; y < h; y++)

	{

		for (int x = 0; x < w; x++)

		{

			/*double u= x+k1*((x-c1)*(x-c1)+(y-c2)*(y-c2))*(x-c1);

			double v= y+k2*((x-c1)*(x-c1)+(y-c2)*(y-c2))*(y-c2);*/

			double u= x;

			double v= y;

			double x_prev=0,y_prev=0;

			double x_curr=x,y_curr=y;

			double delta=DBL_MAX;

			int iter=0;

			for(;delta>DBL_EPSILON;)

			{

				double t1 = x_curr - c1;

				double t2 = t1 * t1;

				double t3 = y_curr - c1;

				double t4 = t3 * t3;

				double deltax = (k1 * (t2 + t4) * t1 + x_curr - u) / k1 / t3 / t1 / 0.2e1;


				t2 = pow(x_curr - c1, 0.2e1);

				t3 = y_curr - c2;

				t4 = t3 * t3;

				double t6 = k2 * (t2 + t4);

				double deltay = (t6 * t3 + y_curr - v) / (0.2e1 * k2 * t3 * t3 + t6 + 0.1e1);


				x_curr=x_curr-deltax;

				y_curr=y_curr-deltay;


				//cout << deltax << endl;


				delta=sqrt(pow((deltax),2)+pow((deltay),2));

				iter++;

				//cout << delta << endl;

			}

			*pbuf_x = (float)(x_curr);

			*pbuf_y = (float)(y_curr);

			++pbuf_x;

			++pbuf_y;

		}

	}


	IplImage* temp = cvCloneImage(img);

    cvRemap( temp, img, mapx, mapy ); 

	cvReleaseImage(&temp);

	cvReleaseImage(&mapx);

	cvReleaseImage(&mapy);


	return img;

}

Поделиться сообщением


Ссылка на сообщение
Поделиться на других сайтах

И что получается?

Вроде должно работать, только при такой индексации

*pbuf_x = (float)(x_curr);
*pbuf_y = (float)(y_curr);
++pbuf_x;
++pbuf_y;[/code]

не учитывается выравнивание и все может ехать.

Поделиться сообщением


Ссылка на сообщение
Поделиться на других сайтах

черный экран с полоской выдает.

не учитывается выравнивание и все может ехать.

не понял, например в эту сторону норм работает.

IplImage* barrel_pincusion_dist(IplImage* img, double Cx,double 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_x = (float*)mapx->imageData;

	float* pbuf_y = (float*)mapy->imageData;

	for (int y = 0; y < h; y++)

	{

		for (int x = 0; x < w; x++)

		{         

			float u= Cx+(x-Cx)*(1+kx*((x-Cx)*(x-Cx)+(y-Cy)*(y-Cy)));

			float v= Cy+(y-Cy)*(1+ky*((x-Cx)*(x-Cx)+(y-Cy)*(y-Cy)));

			*pbuf_x = u;

			*pbuf_y = v;

			++pbuf_x;

			++pbuf_y;

		}

	}

IplImage* temp = cvCloneImage(img);

    cvRemap( temp, img, mapx, mapy ); 

	cvReleaseImage(&temp);

	cvReleaseImage(&mapx);

	cvReleaseImage(&mapy);


	return img;

}

Поделиться сообщением


Ссылка на сообщение
Поделиться на других сайтах

Тут опечатка:

double t1 = x_curr - c1;
double t2 = t1 * t1;
double t3 = y_curr - c1;[/code]

Зарядил себе этот код. Метод часто выдает в качестве результата #INF, при этом дельты увеличиваются.

Значит не выполняется какое-то условие сходимости.

Уменьшал шаг (умножал дельты на 0.0001), результат лучше, но все равно часть точек не считает.

Поделиться сообщением


Ссылка на сообщение
Поделиться на других сайтах

Вот держите (работает для положительных k1 и k2):

post-1-0-94538900-1343909181_thumb.jpg


#include <iostream>
#include <vector>
#include <stdio.h>
#include <stdarg.h>
#include "opencv2/core/core.hpp"
#include "opencv2/highgui/highgui.hpp"
#include "opencv2/imgproc/imgproc.hpp"
#include "fstream"
#include "iostream"
using namespace std;
using namespace cv;
Point uwarp1Point(double u, double v,double c1,double c2,double k1,double k2)
{
double x=u;
double y=v;

double t1 ;
double t4 ;
double t5 ;
double t6 ;
double t7 ;
double t10;
double t11;
double t14;
double t21;
double t23;
double t28;
double t32;
double t54;
double t56;
double t62;
double t65;
double t67;
double t72;
//x+k1*((x-c1)*(x-c1)+(y-c2)*(y-c2))*(x-c1)-u=0;
//y+k2*((x-c1)*(x-c1)+(y-c2)*(y-c2))*(y-c2)-v=0;
double x_prev=x,y_prev=y,delta=DBL_MAX;
for(;delta>1e-5;)
{

if(k1>0)
{
t1 = 1.0 / k1;
t4 = sqrt(0.3e1);
t5 = ( k1 * k1);
t6 = ( y_prev * y_prev);
t7 = ( t6 * t6);
t10 = ( c2 * c2);
t11 = ( t10 * t10);
t14 = ( k1 * t5);
t21 = ( k1 * t6);
t23 = ( k1 * t10);
t28 = ( y_prev * t6);
t32 = ( c2 * t10);
t54 = ( k1 * y_prev * c2);
t56 = ( c1 * c1);
t62 = ( u * u);
t65 = 4 + 12 * t5 * t7 + 12 * t5 * t11 + 4 * t14 * t6 * t7 + 4 * t14 * t10 * t11 + 12 * t21 + 12 * t23 + 72 * t5 * t6 * t10 - 48 * t5 * t28 * c2 - 48 * t5 * t32 * y_prev + 60 * t14 * t7 * t10 - 24 * t14 * y_prev * t7 * c2 + 60 * t14 * t6 * t11 - 80 * t14 * t28 * t32 - 24 * t14 * c2 * t11 * y_prev - 24 * t54 + 27 * k1 * t56 - 54 * k1 * c1 * u + 27 * k1 * t62;
t67 = sqrt( (t65 * t1));
t72 = pow((- (108 * c1) + (108 * u) + 0.12e2 * t4 * t67) * t5, 0.1e1 / 0.3e1);
x = t1 * t72 / 0.6e1 - 0.2e1 * (t21 + t23 + 1 - 2 * t54) / t72 + c1;
}

if(k2>0)
{
t1 = 1.0 / k2;
t4 = sqrt(0.3e1);
t5 = ( k2 * k2);
t6 = ( x_prev * x_prev);
t7 = ( t6 * t6);
t10 = ( c1 * c1);
t11 = ( t10 * t10);
t14 = ( k2 * t5);
t21 = ( k2 * t6);
t23 = ( k2 * t10);
t28 = ( x_prev * t6);
t32 = ( c1 * t10);
t54 = ( k2 * x_prev * c1);
t56 = ( c2 * c2);
t62 = ( v * v);
t65 = 4 + 12 * t5 * t7 + 12 * t5 * t11 + 4 * t14 * t6 * t7 + 4 * t14 * t10 * t11 + 12 * t21 + 12 * t23 + 72 * t5 * t6 * t10 - 48 * t5 * t28 * c1 - 48 * t5 * t32 * x_prev + 60 * t14 * t7 * t10 - 24 * t14 * x_prev * t7 * c1 + 60 * t14 * t6 * t11 - 80 * t14 * t28 * t32 - 24 * t14 * c1 * t11 * x_prev - 24 * t54 + 27 * k2 * t56 - 54 * k2 * c2 * v + 27 * k2 * t62;
t67 = sqrt( (t65 * t1));
t72 = pow((- (108 * c2) + (108 * v) + 0.12e2 * t4 * t67) * t5, 0.1e1 / 0.3e1);
y = t1 * t72 / 0.6e1 - 0.2e1 * (t21 + 1 + t23 - 2 * t54) / t72 + c2;
}
delta=sqrt(pow((x-x_prev),2)+pow((y-y_prev),2));

x_prev=x;
y_prev=y;
}



// cout << x << endl;
// cout << y << endl;
return Point(x,y);
}


Mat un_barrel_pincusion_dist(Mat& img, int c1,int c2,double k1,double k2)
{
Mat mapx=Mat::zeros(img.rows,img.cols,CV_32FC1);
Mat mapy=Mat::zeros(img.rows,img.cols,CV_32FC1);

int w= img.cols;
int h= img.rows;
for (int y = 0; y < h; y++)
{
for (int x = 0; x < w; x++)
{
Point p=uwarp1Point(x, y,c1,c2, k1,k2);
mapy.at<float>(y,x)=p.y;
mapx.at<float>(y,x)=p.x;
}
}
Mat IM;
remap(img,IM,mapx,mapy,1);
return IM;
}

//-----------------------------------------------------------------------------------------------------
//
//-----------------------------------------------------------------------------------------------------
int main( int argc, char** argv )
{
char filename[]="C:\\ImagesForTest\\cat.jpg";

namedWindow("Исправленное");

namedWindow("Исходное");

Mat IM=imread(filename,0);
Mat IMG;
IMG=un_barrel_pincusion_dist(IM,100.0,120.0,2e-4,3e-4);

imshow("Исходное",IM);

imshow("Исправленное",IMG);

waitKey(0);
cout << "OK" << endl;
return 0;
}
[/code]

Делал так: решил первое ур-ние относительно x, второе относительно y , взял действительные решения. И подставил в цикл.

  • Like 1

Поделиться сообщением


Ссылка на сообщение
Поделиться на других сайтах

спасибо, попробую.

кстати я понял как получить параметры искажения имея только 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.

known_center.png

то можно вообще по 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

Поделиться сообщением


Ссылка на сообщение
Поделиться на других сайтах

Ищите "numerical recipes in c" :) там есть готовые решалки на с.

ЗЫ: нормали кривых пересекаются в центре кривизны, так что и c1,c2 можно тоже найти, правда нужно еще найти эти нормали, а для этого еще 2 точки нужны.

Поделиться сообщением


Ссылка на сообщение
Поделиться на других сайтах
numerical recipes in c

так может есть какая нибудь готовая библиотека? просто не очень понятно по каким словам гуглить.

Поделиться сообщением


Ссылка на сообщение
Поделиться на других сайтах

Библиотек полно. Например: http://devernay.free.fr/hacks/cminpack/index.html

Гуглить можно по: nonlinear systems solver, gradient descent, Newtow-Raphson, nonlinear optimization c++ library.

Желательно, правда перед этим прочитать что есть Якобиан (имеется ввиду матрица Якоби), многие методы его хотят.

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

Обычно в качестве оного используется квадрат ошибки или какая-нибудь другая ошибка, которая в точке решения принимает экстремальное (мин или макс) значение.

У Вас квадрат ошибки может быть получен например так:

Выражаем уравнение в виде: f(X)=0;

Тогда минимизируемым функционалом будет: E(X)=(f(X))^2; видно, что в точке решения f(X)=0, E(X) имеет минимум.

Поделиться сообщением


Ссылка на сообщение
Поделиться на других сайтах

решил по 4-м точкам получить параметры перспективного искажения

у меня получилось

http://mathematica.stackexchange.com/questions/9244/solve-system-of-equations

математика зависает на обдумывании.

а в опенцв какая то магия

я так понял, что тут решается A*X=B

т.е. все таки в зд система получается линейная а если записать через формулу в 2д то нелинейная.

cv::Mat cv::getPerspectiveTransform( const Point2f src[], const Point2f dst[] )

{

    Mat M(3, 3, CV_64F), X(8, 1, CV_64F, M.data);

    double a[8][8], b[8];

    Mat A(8, 8, CV_64F, a), B(8, 1, CV_64F, ;


    for( int i = 0; i < 4; ++i )

    {

        a[i][0] = a[i+4][3] = src[i].x;

        a[i][1] = a[i+4][4] = src[i].y;

        a[i][2] = a[i+4][5] = 1;

        a[i][3] = a[i][4] = a[i][5] =

        a[i+4][0] = a[i+4][1] = a[i+4][2] = 0;

        a[i][6] = -src[i].x*dst[i].x;

        a[i][7] = -src[i].y*dst[i].x;

        a[i+4][6] = -src[i].x*dst[i].y;

        a[i+4][7] = -src[i].y*dst[i].y;

        b[i] = dst[i].x;

        b[i+4] = dst[i].y;

    }


    solve( A, B, X, DECOMP_SVD );

    ((double*)M.data)[8] = 1.;


    return M;

}[/code]

Поделиться сообщением


Ссылка на сообщение
Поделиться на других сайтах

все понял, на самом деле получается линейная система от mij всё правильно.

Поделиться сообщением


Ссылка на сообщение
Поделиться на других сайтах

задался такой мыслью

допустим знаем только координаты искаженные u,v для n точек, но так же знаем, что расстояния между точками для x,y всегда постоянны.

т.е. для n точек мы можем получить n-1 доп уравнений.

составил уравнения и получил такую формулу

k=4/(n-1)+1

т.е. например надо всего 2 набора по 5 точек или например 5 наборов по 2 точки.

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

Поделиться сообщением


Ссылка на сообщение
Поделиться на других сайтах

Можно и больше уравнений получить, ведь система нелинейна, а значит то, что на линейной системе равномерно, здесь не равномерно.

И расстояние между любой парой точек является линейно-независимым от любой другой пары точек.

Поделиться сообщением


Ссылка на сообщение
Поделиться на других сайтах

кажется я догадался надо взять и переписать эти 30 как A*X=B

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

Поделиться сообщением


Ссылка на сообщение
Поделиться на других сайтах

нет ошибся так не получится потому что система не линейная.

проблема в том, что LeastSquares в mathematica принимает только матрицы ,т.е. я так предполагаю, что линейную систему.

Поделиться сообщением


Ссылка на сообщение
Поделиться на других сайтах

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

http://xenia.media.mit.edu/~cwren/interpolator/

http://www.robots.ox.ac.uk/~vgg/presentations/bmvc97/criminispaper/node1.html#SECTION00010000000000000000

похоже что задача довольно похожа

но я например не понял начиная с момента чего они хотят добиться

eq6.gif

eq11.gif

Поделиться сообщением


Ссылка на сообщение
Поделиться на других сайтах

Просто представили предыдущие расчеты в виде произведения вектора и матрицы.

Написали формальную матрицу для подстановки точек.

Поделиться сообщением


Ссылка на сообщение
Поделиться на других сайтах

ну как матрицу они записали первую я понял, а что значат точки? (типа для произвольного числа уравнений точек?)

и опять же непонятно у них получается неизвестные и слева и справа.

а хотя всё понял, у них как раз неизвестные в матрице гомографии , это у меня неизвестны еще x,y :)

Поделиться сообщением


Ссылка на сообщение
Поделиться на других сайтах

может попробовать решить в общем виде, часть уравнений как линейное уравнение относительно m(i,j) (представить что все остальные данные даны), а потом уже дорешать с остальными уравнениями или так не прокатит?

Поделиться сообщением


Ссылка на сообщение
Поделиться на других сайтах

Ага, сейчас к методу похожему на maximum likelihood estimation придем :)

Там как раз сначала одну половину переменных принимают известной (параметры), по ней ищут приближение второй половины (скрытых переменных), а затем, по найденному приближению скрытых переменных ищут приближение параметров. И по кругу.

Поделиться сообщением


Ссылка на сообщение
Поделиться на других сайтах

попробовал приравнять m32=m31=0 чтобы хотя бы для афинных решить и выдает, что нет решений.

может там где то линейная зависимость в уравнениях, а я её просто не вижу?

может можно как то программно проверить?

Поделиться сообщением


Ссылка на сообщение
Поделиться на других сайтах

Создайте учётную запись или войдите для комментирования

Вы должны быть пользователем, чтобы оставить комментарий

Создать учётную запись

Зарегистрируйтесь для создания учётной записи. Это просто!

Зарегистрировать учётную запись

Войти

Уже зарегистрированы? Войдите здесь.

Войти сейчас


  • Сейчас на странице   0 пользователей

    Нет пользователей, просматривающих эту страницу

×