Jump to content
Compvision.ru
Sign in to follow this  
mrgloom

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

Recommended Posts

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

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

double Cx=100.0;

double Cy=120.0;

double x= 100;

double y= 120;

то выдаст

x=-1.#IND

Share this post


Link to post
Share on other sites

что то я не понял, как это использовать для развертывания бочки обратно, в начале мы имеем х,у которые как бы хотим определить получаем по ним 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;

}

Share this post


Link to post
Share on other sites

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

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

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

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

Share this post


Link to post
Share on other sites

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

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

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

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;

}

Share this post


Link to post
Share on other sites

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

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

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

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

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

Share this post


Link to post
Share on other sites

Вот держите (работает для положительных 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

Share this post


Link to post
Share on other sites

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

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

Share this post


Link to post
Share on other sites

а если мы знаем где находится центр, обычно в центре изображения, т.е. знаем с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

Share this post


Link to post
Share on other sites

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

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

Share this post


Link to post
Share on other sites
numerical recipes in c

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

Share this post


Link to post
Share on other sites

Библиотек полно. Например: 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) имеет минимум.

Share this post


Link to post
Share on other sites

решил по 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]

Share this post


Link to post
Share on other sites

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

Share this post


Link to post
Share on other sites

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

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

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

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

k=4/(n-1)+1

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

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

Share this post


Link to post
Share on other sites

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

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

Share this post


Link to post
Share on other sites

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

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

Share this post


Link to post
Share on other sites

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

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

Share this post


Link to post
Share on other sites

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

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

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

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

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

eq6.gif

eq11.gif

Share this post


Link to post
Share on other sites

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

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

Share this post


Link to post
Share on other sites

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

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

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

Share this post


Link to post
Share on other sites

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

Share this post


Link to post
Share on other sites

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

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

Share this post


Link to post
Share on other sites

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

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

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

Share this post


Link to post
Share on other sites

Create an account or sign in to comment

You need to be a member in order to leave a comment

Create an account

Sign up for a new account in our community. It's easy!

Register a new account

Sign in

Already have an account? Sign in here.

Sign In Now

Sign in to follow this  

  • Recently Browsing   0 members

    No registered users viewing this page.

×