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

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

Recommended Posts

Ух,ниасилил :) я обычно просто обратную подстановку использую.

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

то есть подставьте все в

u = с1 + x + k1*(x*x + y*y)*x;

v = c2 + y + k2*(x*x + y*y)*y;

и посмотрите равны ли правые части левым.

Share this post


Link to post
Share on other sites

кажется понял в чем ошибка, забыл делать смещение до центра (с1,с2)

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

ну тогда я попробовал решить это дело относительно k1,k2 , потому что центр с1,с2 мы можем и так найти вручную путем нахождения максимально удаленных точек на противоположных дугах и проведения по этим точкам прямым, а на пересечении и будет искомый центр.

проверил на ассиметричной бочке.в итоге ошибся на 1 пиксель по х.

assimetric_barrel.png

CForm[k1 /.

Solve[v1 == c2 + (y1 - c2) (1 + k2 ((x1 - c1)^2 + (y1 - c2)^2)) &&

u1 == c1 + (x1 - c1) (1 + k1 ((x1 - c1)^2 + (y1 - c2)^2)), {k1,

k2}]]

CForm[k2 /.

Solve[v1 == c2 + (y1 - c2) (1 + k2 ((x1 - c1)^2 + (y1 - c2)^2)) &&

u1 == c1 + (x1 - c1) (1 + k1 ((x1 - c1)^2 + (y1 - c2)^2)), {k1,

k2}]]

void barrel_parameters(double x1,double y1,double u1,double v1,double c1,double c2)

{

	double k1= (-u1 + x1)/((c1 - x1)*(pow(c1,2) + pow(c2,2) - 2*c1*x1 + pow(x1,2) - 2*c2*y1 + pow(y1,2)));

	double k2= (-v1 + y1)/((c2 - y1)*(pow(c1,2) + pow(c2,2) - 2*c1*x1 + pow(x1,2) - 2*c2*y1 + pow(y1,2)));

}

пробую посчитать недостающие параметры и получаю 9.4*10^-6 хотя на самом деле 2*10^-5, возможно не точно взял точку на рисунке или просто нельзя добиться большей точности из-за вычислений?

про подстановку не разобрался пока как это сделать в математике.

Share this post


Link to post
Share on other sites

подстановку можно и на бумажке сделать.

допустим алгоритм для x=1,y=2,u=3,v=4 выдает коэффициенты k1=5 k2=6, с1=7 с2=8.

Берем формулу и тупо вбиваем цифры:

3 = 8 + 1 + 5*(1*1 + 2*2)*1;

4 = 7 + 2 + 6*(1*1 + 2*2)*2;

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

Share this post


Link to post
Share on other sites

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

void test_barrel(double x,double y,double c1,double c2,double k1,double k2)

{

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

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

	double u= c1+(x-c1)*(1+tx);

	double v= c2+(y-c2)*(1+ty);

}

для точки (0,0) выдает отрицательные координаты, что странно.

хотя может это получаются уже относительно центра?

потом еще непонятен вопрос как потом сделать обратное преобразование, т.к.

barrel_pincusion_dist(barrel_pincusion_dist(src, с1,с2,к1,к2)с1,с2,-к1,-к2)

не даст обратного преобразования.

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

Solve[v == c2 + (y - c2) (1 + k2 ((x - c1)^2 + (y - c2)^2)) &&

u == c1 + (x - c1) (1 + k1 ((x - c)^2 + (y - c2)^2)), {x,

y}]

Share this post


Link to post
Share on other sites

Самый простой способ - соорудить итеративную схему, как в школе учили :)

	// считаем u,v по x,y
// здесь все очевидно
double k1=2e-4;
double k2=3e-4;
double c1=100.0;
double c2=100.0;
double x=10;
double y=20;
double u = c1 + x + k1*(x*x + y*y)*x;
double v = c2 + y + k2*(x*x + y*y)*y;
cout << u << endl;
cout << v << endl;


// Вычисление x,y по u,v
// выразили x и y через известные и самих себя
// и крутим цикл пока не сойдется.

double x_prev=0,y_prev=0;
x=0,y=0;
double delta=DBL_MAX;
for(;delta>DBL_EPSILON;)
{
x=u-c1-k1*(x_prev*x_prev + y_prev*y_prev)*x_prev;
y=v-c2-k2*(x_prev*x_prev + y_prev*y_prev)*y_prev;
delta=sqrt(pow((x-x_prev),2)+pow((y-y_prev),2));
x_prev=x;
y_prev=y;
}

cout << x << endl;
cout << y << endl;[/code]

Share this post


Link to post
Share on other sites

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

еще по поводу уравнения A*X=X*B

matrix_eq.gif

нашел в книжке про решение такого уравнения там есть такой комментарий

Image_27-07-2012%20%5B10-34-59%5D.png

только всё же не понятно, почему выдало всё 0, если уравнения были записаны в общем виде с переменными a,b.

произвольный объект и произвольное перспективное искажение для иллюстрации задачи.

я всё таки так и не пойму реально ли восстановить искажение или нет.

произвольный объект X

%D0%BF%D1%80%D0%BE%D0%B8%D0%B7%D0%B2%D0%BE%D0%BB%D1%8C%D0%BD%D1%8B%D0%B9%20%D0%BE%D0%B1%D1%8A%D0%B5%D0%BA%D1%82.PNG

произвольный объект повернутый на 90 градусов X*Rot

%D0%BF%D1%80%D0%BE%D0%B8%D0%B7%D0%B2%D0%BE%D0%BB%D1%8C%D0%BD%D1%8B%D0%B9%20%D0%BE%D0%B1%D1%8A%D0%B5%D0%BA%D1%82_90_%D0%93%D1%80%D0%B0%D0%B4.png

произвольный объект с искажением X*M

%D0%BF%D1%80%D0%BE%D0%B8%D0%B7%D0%B2%D0%BE%D0%BB%D1%8C%D0%BD%D1%8B%D0%B9%20%D0%BE%D0%B1%D1%8A%D0%B5%D0%BA%D1%82_%D0%9C.PNG.png

произвольный объект повернутый на 90 с искажением X*Rot*M

%D0%BF%D1%80%D0%BE%D0%B8%D0%B7%D0%B2%D0%BE%D0%BB%D1%8C%D0%BD%D1%8B%D0%B9%20%D0%BE%D0%B1%D1%8A%D0%B5%D0%BA%D1%82_90_%D0%93%D1%80%D0%B0%D0%B4_%D0%9C.png

Share this post


Link to post
Share on other sites

попробовал еще в гимпе откатить свою бочку, но там какие то странные параметры.

Filters->Distorts->Lens Distortion

есть ползунок Main, а есть Edge

но зато можно попробовать найти сорцы

вроде как еще можно откатить бочку через cvInitUndistortMap

Share this post


Link to post
Share on other sites

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

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

Share this post


Link to post
Share on other sites

а как посчитать матрицу гомографии для поворота на произвольный угол?

тут пишут

http://tech.groups.yahoo.com/group/OpenCV/message/84530

H =

[scos(\theta) -ssin(\theta) t_x;

ssin(\theta) scos(\theta) t_y;

0 0 1]

но я проверил это не работает.

например по часовой на 90 не соответствует формуле

cvmSet(warp_matrix,0,0,0); // Set M(i,j)

cvmSet(warp_matrix,0,1,1); // Set M(i,j)

cvmSet(warp_matrix,0,2,0); // Set M(i,j)

cvmSet(warp_matrix,1,0,1); // Set M(i,j)

cvmSet(warp_matrix,1,1,0); // Set M(i,j)

cvmSet(warp_matrix,1,2,0); // Set M(i,j)

cvmSet(warp_matrix,2,0,0); // Set M(i,j)

cvmSet(warp_matrix,2,1,0); // Set M(i,j)

cvmSet(warp_matrix,2,2,1); // Set M(i,j)

тут уже другое пишут http://stackoverflow.com/questions/10969170/rotation-angle-of-the-perspective-matrix

но не написано как найти

Share this post


Link to post
Share on other sites

Аффинные преобразования:

http://en.wikipedia.org/wiki/Transformation_matrix#Affine_transformations

И, судя по Вашим ссылкам, полезно будет копнуть под Родригиса:

http://docs.opencv.org/modules/calib3d/doc/camera_calibration_and_3d_reconstruction.html?highlight=rodrigues#cv.Rodrigues2

Share this post


Link to post
Share on other sites

пробую все таки решить это матричное уравнение

если записать полностью, выдает все нули

Solve[a11*m11 + a12*m21 + a13*m31 == m11*b11 + m12*b21 + m13*b31 &&

a11*m12 + a12*m22 + a13*m32 == m11*b12 + m12*b22 + m13*b32 &&

a11*m13 + a12*m23 + a13*m33 == m11*b13 + m12*b23 + m13*b33 &&

a21*m11 + a22*m21 + a23*m31 == m21*b11 + m22*b21 + m23*b31 &&

a21*m12 + a22*m22 + a23*m32 == m21*b12 + m22*b22 + m23*b32 &&

a21*m13 + a22*m23 + a23*m33 == m21*b13 + m22*b23 + m23*b33 &&

a31*m11 + a32*m21 + a33*m31 == m31*b11 + m32*b21 + m33*b31 &&

a31*m12 + a32*m22 + a33*m32 == m31*b12 + m32*b22 + m33*b32 &&

a31*m13 + a32*m23 + a33*m33 == m31*b13 + m32*b23 + m33*b33, {m11,

m12, m13, m21, m22, m23, m31, m32,m33}]

если положить m33==1 и 8 неизвестных пишет нет решений (или не знаю как интерпретировать пустой вывод)

Solve[a11*m11 + a12*m21 + a13*m31 == m11*b11 + m12*b21 + m13*b31 &&

a11*m12 + a12*m22 + a13*m32 == m11*b12 + m12*b22 + m13*b32 &&

a11*m13 + a12*m23 + a13 == m11*b13 + m12*b23 + m13*b33 &&

a21*m11 + a22*m21 + a23*m31 == m21*b11 + m22*b21 + m23*b31 &&

a21*m12 + a22*m22 + a23*m32 == m21*b12 + m22*b22 + m23*b32 &&

a21*m13 + a22*m23 + a23 == m21*b13 + m22*b23 + m23*b33 &&

a31*m11 + a32*m21 + a33*m31 == m31*b11 + m32*b21 + b31 &&

a31*m12 + a32*m22 + a33*m32 == m31*b12 + m32*b22 + b32 &&

a31*m13 + a32*m23 + a33 == m31*b13 + m32*b23 + b33, {m11, m12, m13,

m21, m22, m23, m31, m32}]

если убрать одно уравнение, например первое, то выдает решение примерно на страницу.

попробовал применить его на практике, написал код для тестирования.

но выдает что то не похожее на правду, возможно потому что матрицу поворота и матрицу перехода я получал по точках ,а значит не совсем точную.

void get_transform_matrix(CvMat* A,CvMat* B,CvMat* err)

{

	double a11= cvmGet(A,0,0); // Set M(i,j)

	double a12= cvmGet(A,0,1); // Set M(i,j)

	double a13= cvmGet(A,0,2); // Set M(i,j)

	double a21= cvmGet(A,1,0); // Set M(i,j)

	double a22= cvmGet(A,1,1); // Set M(i,j)

	double a23= cvmGet(A,1,2); // Set M(i,j)

	double a31= cvmGet(A,2,0); // Set M(i,j)

	double a32= cvmGet(A,2,1); // Set M(i,j)

	double a33= cvmGet(A,2,2); // Set M(i,j)


	double b11= cvmGet(B,0,0); // Set M(i,j)

	double b12= cvmGet(B,0,1); // Set M(i,j)

	double b13= cvmGet(B,0,2); // Set M(i,j)

	double b21= cvmGet(B,1,0); // Set M(i,j)

	double b22= cvmGet(B,1,1); // Set M(i,j)

	double b23= cvmGet(B,1,2); // Set M(i,j)

	double b31= cvmGet(B,2,0); // Set M(i,j)

	double b32= cvmGet(B,2,1); // Set M(i,j)

	double b33= cvmGet(B,2,2); // Set M(i,j)


	//matrix A - rotation matrix

	//matrix B - transfer matrix

	//matrix M - matrix of distortion

	//near solve of A*M=M*B matrix equation

	//m33==1


	double m33=1;


	//2 лиска до по вычислению значений матрицы М



	double err11= cvmGet(err,0,0); // Set M(i,j)

	double err12= cvmGet(err,0,1); // Set M(i,j)

	double err13= cvmGet(err,0,2); // Set M(i,j)

	double err21= cvmGet(err,1,0); // Set M(i,j)

	double err22= cvmGet(err,1,1); // Set M(i,j)

	double err23= cvmGet(err,1,2); // Set M(i,j)

	double err31= cvmGet(err,2,0); // Set M(i,j)

	double err32= cvmGet(err,2,1); // Set M(i,j)

	double err33= cvmGet(err,2,2); // Set M(i,j)


//смотрим отклонение от реально примененного (гипотетически неизвестного) искажения

	double t11= abs(m11-err11);

	double t12= abs(m12-err12);

	double t13= abs(m13-err13);

	double t21= abs(m21-err21);

	double t22= abs(m22-err22);

	double t23= abs(m23-err23);

	double t31= abs(m31-err31);

	double t32= abs(m32-err32);

	double t33= abs(m33-err33);

}

void test_hipothesis()

{

	//для начала определяем произвольную матрицу поворота

	//по соответствиям точек на 2-х изображениях

	CvMat* Rot_matrix= cvCreateMat(3,3,CV_32FC1);

	CvPoint2D32f srcQuad[4], dstQuad[4];

	srcQuad[0].x = 32.5;  //src Top left

	srcQuad[0].y = 79;

	srcQuad[1].x = 236;  //src Top right

	srcQuad[1].y = 19.3;

	srcQuad[2].x = 243;  //src Bottom left

	srcQuad[2].y = 150;

	srcQuad[3].x = 80.5;  //src Bot right

	srcQuad[3].y = 247;

	dstQuad[0].x = 79.4;  //Top left

	dstQuad[0].y = 267;  

	dstQuad[1].x = 19.3;  //Top right

	dstQuad[1].y = 63.4;

	dstQuad[2].x = 150;  //Bottom left

	dstQuad[2].y = 57.3;

	dstQuad[3].x = 247;  //Bottom right

	dstQuad[3].y = 220;

	cvGetPerspectiveTransform(srcQuad,dstQuad,Rot_matrix);

	double rot11= cvmGet(Rot_matrix,0,0); // Set M(i,j)

	double rot12= cvmGet(Rot_matrix,0,1); // Set M(i,j)

	double rot13= cvmGet(Rot_matrix,0,2); // Set M(i,j)

	double rot21= cvmGet(Rot_matrix,1,0); // Set M(i,j)

	double rot22= cvmGet(Rot_matrix,1,1); // Set M(i,j)

	double rot23= cvmGet(Rot_matrix,1,2); // Set M(i,j)

	double rot31= cvmGet(Rot_matrix,2,0); // Set M(i,j)

	double rot32= cvmGet(Rot_matrix,2,1); // Set M(i,j)

	double rot33= cvmGet(Rot_matrix,2,2); // Set M(i,j)


	////наш изначальный произвольный объект

	//IplImage* src=cvLoadImage("произвольный объект.PNG");

	////проверяем результат повернутый объект

	//cvSaveImage("test_matrix.png",warp_rect(src, rot_matrix));


	//далее искажаем проивольным искажением которое как бы неизвестно и его надо определить

	CvMat* warp_matrix= cvCreateMat(3,3,CV_32FC1);

	cvmSet(warp_matrix,0,0,0.59); // Set M(i,j)

	cvmSet(warp_matrix,0,1,-0.22); // Set M(i,j)

	cvmSet(warp_matrix,0,2,66); // Set M(i,j)

	cvmSet(warp_matrix,1,0,-0.19); // Set M(i,j)

	cvmSet(warp_matrix,1,1,0.59); // Set M(i,j)

	cvmSet(warp_matrix,1,2,57); // Set M(i,j)

	cvmSet(warp_matrix,2,0,-0.00063); // Set M(i,j)

	cvmSet(warp_matrix,2,1,-0.00073); // Set M(i,j)

	cvmSet(warp_matrix,2,2,1); // Set M(i,j)


	////получаем ищображения

	//IplImage* src=cvLoadImage("произвольный объект.PNG");

	//cvSaveImage("test_matrix.png",warp_rect(warp_rect(src, rot_matrix),warp_matrix));


	//получаем матрицу перехода между двумя искаженными изображениями


	CvMat* L_matrix= cvCreateMat(3,3,CV_32FC1);

	srcQuad[0].x = 74;  //src Top left

	srcQuad[0].y = 106;

	srcQuad[1].x = 240;  //src Top right

	srcQuad[1].y = 28.3;

	srcQuad[2].x = 239;  //src Bottom left

	srcQuad[2].y = 135;

	srcQuad[3].x = 77.3;  //src Bot right

	srcQuad[3].y = 244;

	dstQuad[0].x = 71.3;  //Top left

	dstQuad[0].y = 264;  

	dstQuad[1].x = 68;  //Top right

	dstQuad[1].y = 96.3;

	dstQuad[2].x = 164;  //Bottom left

	dstQuad[2].y = 72.3;

	dstQuad[3].x = 239;  //Bottom right

	dstQuad[3].y = 204;

	cvGetPerspectiveTransform(srcQuad,dstQuad,L_matrix);

	double L_11= cvmGet(L_matrix,0,0); // Set M(i,j)

	double L_12= cvmGet(L_matrix,0,1); // Set M(i,j)

	double L_13= cvmGet(L_matrix,0,2); // Set M(i,j)

	double L_21= cvmGet(L_matrix,1,0); // Set M(i,j)

	double L_22= cvmGet(L_matrix,1,1); // Set M(i,j)

	double L_23= cvmGet(L_matrix,1,2); // Set M(i,j)

	double L_31= cvmGet(L_matrix,2,0); // Set M(i,j)

	double L_32= cvmGet(L_matrix,2,1); // Set M(i,j)

	double L_33= cvmGet(L_matrix,2,2); // Set M(i,j)


	//затем подсчитываем искомое искажение

	get_transform_matrix(Rot_matrix,L_matrix,warp_matrix);

}

Share this post


Link to post
Share on other sites

Аффинные преобразования:

http://en.wikipedia.org/wiki/Transformation_matrix#Affine_transformations

И, судя по Вашим ссылкам, полезно будет копнуть под Родригиса:

http://docs.opencv.org/modules/calib3d/doc/camera_calibration_and_3d_reconstruction.html?highlight=rodrigues#cv.Rodrigues2

нет (по крайней мере в опенцв) матрица задается как то не так

например матрица

0 1 0

1 0 0

0 0 1

повернёт изображение на 90 градусов по часовой стрелке, да и значения в первом квадрате 2х2 могут быть больше 1, так что там не син\кос и даже похоже не син\кос умноженные на константу, как я писал выше, либо непонятно от чего зависит эта константа.

вообщем надо функцию типа повернуть на phi угол

а она выдает матрицу 3х3

а что под родригисом?)

там что то про калибровку вроде.

Share this post


Link to post
Share on other sites

посмотрел сорцы и там такая интересная функция была

cv::Mat cv::getRotationMatrix2D( Point2f center, double angle, double scale )

{

angle *= CV_PI/180;

double alpha = cos(angle)*scale;

double beta = sin(angle)*scale;

Mat M(2, 3, CV_64F);

double* m = (double*)M.data;

m[0] = alpha;

m[1] = beta;

m[2] = (1-alpha)*center.x - beta*center.y;

m[3] = -beta;

m[4] = alpha;

m[5] = beta*center.x + (1-alpha)*center.y;

return M;

}

понял, что матрица

0 1 0

1 0 0

0 0 1

не просто поворачивает на 90 по часовой, а еще и отражает.

Share this post


Link to post
Share on other sites

еще про уравнение A*X=X*B

учитывая ограничения налагаемые на матрицы А и М

А- матрица поворота,значит у неё последняя строка = 0 0 1

M- матрица перспективного искажения у нее m33==1 ( и b33==1 кстати тоже)

получаем 3 последних уравнения из системы

de8ac9d80e48b8a3c848021f94dcbeae82.png

т.е. получается, что решений нет или это просто даёт доп. ограничения на b(i,j)?

Share this post


Link to post
Share on other sites

попробовал выразить теоретически L

M*L=Rot*M => L=M^-1*Rot*M

но похоже проблема в том, что если последовательно применить 2 матрицы искажений, то результат будет один, а если матрицы сначала перемножить, а потом применить результат будет другой.

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

или всё таки у меня ошибка в рассуждениях?

IplImage* warp_rect(IplImage* src, CvMat* warp_matrix) 

{

	IplImage* dst=cvCloneImage(src);

	cvWarpPerspective( src, dst, warp_matrix);

	return dst;

}

	Mat M(3,3,CV_64FC1);

	Mat E(3,3,CV_64FC1);

	Mat L1(3,3,CV_64FC1);

	Mat L2(3,3,CV_64FC1);

	Mat Rot(3,3,CV_64FC1);

	Mat S(3,3,CV_64FC1);

	M.at<double>(0,0)=0.59;

	M.at<double>(0,1)=-0.22;

	M.at<double>(0,2)=66;

	M.at<double>(1,0)=-0.19;

	M.at<double>(1,1)=0.59;

	M.at<double>(1,2)=57;

	M.at<double>(2,0)=-0.00063;

	M.at<double>(2,1)=-0.00073;

	M.at<double>(2,2)=1;

	Rot.at<double>(0,0)=0;

	Rot.at<double>(0,1)=-1;

	Rot.at<double>(0,2)=300;

	Rot.at<double>(1,0)=1;

	Rot.at<double>(1,1)=0;

	Rot.at<double>(1,2)=0;

	Rot.at<double>(2,0)=0;

	Rot.at<double>(2,1)=0;

	Rot.at<double>(2,2)=1;


	L1= M.inv()*Rot*M;

	cout<<L1;

	L2= M.inv()*Rot.inv()*M;

	cout<<L2;


	S=Rot*M;


CvMat* S_matrix= cvCreateMat(3,3,CV_32FC1);

	cvmSet(S_matrix,0,0,S.at<double>(0,0)); 

	cvmSet(S_matrix,0,1,S.at<double>(0,1)); 

	cvmSet(S_matrix,0,2,S.at<double>(0,2));

	cvmSet(S_matrix,1,0,S.at<double>(1,0));

	cvmSet(S_matrix,1,1,S.at<double>(1,1));

	cvmSet(S_matrix,1,2,S.at<double>(1,2));

	cvmSet(S_matrix,2,0,S.at<double>(2,0));

	cvmSet(S_matrix,2,1,S.at<double>(2,1));

	cvmSet(S_matrix,2,2,S.at<double>(2,2));

	CvMat* M_matrix= cvCreateMat(3,3,CV_32FC1);

	cvmSet(M_matrix,0,0,M.at<double>(0,0)); 

	cvmSet(M_matrix,0,1,M.at<double>(0,1)); 

	cvmSet(M_matrix,0,2,M.at<double>(0,2));

	cvmSet(M_matrix,1,0,M.at<double>(1,0));

	cvmSet(M_matrix,1,1,M.at<double>(1,1));

	cvmSet(M_matrix,1,2,M.at<double>(1,2));

	cvmSet(M_matrix,2,0,M.at<double>(2,0));

	cvmSet(M_matrix,2,1,M.at<double>(2,1));

	cvmSet(M_matrix,2,2,M.at<double>(2,2));

	CvMat* L_matrix= cvCreateMat(3,3,CV_32FC1);

	cvmSet(L_matrix,0,0,L1.at<double>(0,0));

	cvmSet(L_matrix,0,1,L1.at<double>(0,1));

	cvmSet(L_matrix,0,2,L1.at<double>(0,2));

	cvmSet(L_matrix,1,0,L1.at<double>(1,0));

	cvmSet(L_matrix,1,1,L1.at<double>(1,1));

	cvmSet(L_matrix,1,2,L1.at<double>(1,2));

	cvmSet(L_matrix,2,0,L1.at<double>(2,0));

	cvmSet(L_matrix,2,1,L1.at<double>(2,1));

	cvmSet(L_matrix,2,2,L1.at<double>(2,2));

	CvMat* Rot_matrix= cvCreateMat(3,3,CV_32FC1);

	cvmSet(Rot_matrix,0,0,Rot.at<double>(0,0));

	cvmSet(Rot_matrix,0,1,Rot.at<double>(0,1));

	cvmSet(Rot_matrix,0,2,Rot.at<double>(0,2));

	cvmSet(Rot_matrix,1,0,Rot.at<double>(1,0));

	cvmSet(Rot_matrix,1,1,Rot.at<double>(1,1));

	cvmSet(Rot_matrix,1,2,Rot.at<double>(1,2));

	cvmSet(Rot_matrix,2,0,Rot.at<double>(2,0));

	cvmSet(Rot_matrix,2,1,Rot.at<double>(2,1));

	cvmSet(Rot_matrix,2,2,Rot.at<double>(2,2));

IplImage* src=cvLoadImage("произвольный объект.PNG");

cvSaveImage("Rot_M_matrix1.png",warp_rect(warp_rect(src,Rot_matrix),M_matrix));

cvSaveImage("Rot_M_matrix2.png",warp_rect(src,S_matrix));

Share this post


Link to post
Share on other sites

наблюдение о том 5057e0c370be3a1af8a4cd7c132e5bda82.png

т.е. 2 последовательные операции невозможно проделать за 1 операцию, заменив произведением матриц

и ваш вопрос навёл меня на мысль о том, что матричное уравнение составлено не верно.

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

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

b72962aafc759ce539415f71981f9cd082.png

а преобразование записываются в такой форме

прямое

5e8e8984a34d8aed41e5b5020de5adf082.png

и обратное

d94f1c566dcff2552ef729bc3e4638ff82.png

можно составить уравнение для 4-х точек, но оно опять же как то не очень хорошо решается.

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

Solve[(b13+(b11 (m13+m11 x1+m12 y1))/(m33+m31 x1+m32 y1)+(b12 (m23+m21 x1+m22 y1))/(m33+m31 x1+m32 y1))/(b33+(b31 (m13+m11 x1+m12 y1))/(m33+m31 x1+m32 y1)+(b32 (m23+m21 x1+m22 y1))/(m33+m31 x1+m32 y1))==

(m13+(m11 (a13+a11 x1+a12 y1))/(a33+a31 x1+a32 y1)+(m12 (a23+a21 x1+a22 y1))/(a33+a31 x1+a32 y1))/(m33+(m31 (a13+a11 x1+a12 y1))/(a33+a31 x1+a32 y1)+(m32 (a23+a21 x1+a22 y1))/(a33+a31 x1+a32 y1))&&

(b23+(b21 (m13+m11 x1+m12 y1))/(m33+m31 x1+m32 y1)+(b22 (m23+m21 x1+m22 y1))/(m33+m31 x1+m32 y1))/(b33+(b31 (m13+m11 x1+m12 y1))/(m33+m31 x1+m32 y1)+(b32 (m23+m21 x1+m22 y1))/(m33+m31 x1+m32 y1))==

(m23+(m21 (a13+a11 x1+a12 y1))/(a33+a31 x1+a32 y1)+(m22 (a23+a21 x1+a22 y1))/(a33+a31 x1+a32 y1))/(m33+(m31 (a13+a11 x1+a12 y1))/(a33+a31 x1+a32 y1)+(m32 (a23+a21 x1+a22 y1))/(a33+a31 x1+a32 y1))&&

(b13+(b11 (m13+m11 x2+m12 y2))/(m33+m31 x2+m32 y2)+(b12 (m23+m21 x2+m22 y2))/(m33+m31 x2+m32 y2))/(b33+(b31 (m13+m11 x2+m12 y2))/(m33+m31 x2+m32 y2)+(b32 (m23+m21 x2+m22 y2))/(m33+m31 x2+m32 y2))==

(m13+(m11 (a13+a11 x2+a12 y2))/(a33+a31 x2+a32 y2)+(m12 (a23+a21 x2+a22 y2))/(a33+a31 x2+a32 y2))/(m33+(m31 (a13+a11 x2+a12 y2))/(a33+a31 x2+a32 y2)+(m32 (a23+a21 x2+a22 y2))/(a33+a31 x2+a32 y2))&&

(b23+(b21 (m13+m11 x2+m12 y2))/(m33+m31 x2+m32 y2)+(b22 (m23+m21 x2+m22 y2))/(m33+m31 x2+m32 y2))/(b33+(b31 (m13+m11 x2+m12 y2))/(m33+m31 x2+m32 y2)+(b32 (m23+m21 x2+m22 y2))/(m33+m31 x2+m32 y2))==

(m23+(m21 (a13+a11 x2+a12 y2))/(a33+a31 x2+a32 y2)+(m22 (a23+a21 x2+a22 y2))/(a33+a31 x2+a32 y2))/(m33+(m31 (a13+a11 x2+a12 y2))/(a33+a31 x2+a32 y2)+(m32 (a23+a21 x2+a22 y2))/(a33+a31 x2+a32 y2))&&

(b13+(b11 (m13+m11 x3+m12 y3))/(m33+m31 x3+m32 y3)+(b12 (m23+m21 x3+m22 y3))/(m33+m31 x3+m32 y3))/(b33+(b31 (m13+m11 x3+m12 y3))/(m33+m31 x3+m32 y3)+(b32 (m23+m21 x3+m22 y3))/(m33+m31 x3+m32 y3))==

(m13+(m11 (a13+a11 x3+a12 y3))/(a33+a31 x3+a32 y3)+(m12 (a23+a21 x3+a22 y3))/(a33+a31 x3+a32 y3))/(m33+(m31 (a13+a11 x3+a12 y3))/(a33+a31 x3+a32 y3)+(m32 (a23+a21 x3+a22 y3))/(a33+a31 x3+a32 y3))&&

(b23+(b21 (m13+m11 x3+m12 y3))/(m33+m31 x3+m32 y3)+(b22 (m23+m21 x3+m22 y3))/(m33+m31 x3+m32 y3))/(b33+(b31 (m13+m11 x3+m12 y3))/(m33+m31 x3+m32 y3)+(b32 (m23+m21 x3+m22 y3))/(m33+m31 x3+m32 y3))==

(m23+(m21 (a13+a11 x3+a12 y3))/(a33+a31 x3+a32 y3)+(m22 (a23+a21 x3+a22 y3))/(a33+a31 x3+a32 y3))/(m33+(m31 (a13+a11 x3+a12 y3))/(a33+a31 x3+a32 y3)+(m32 (a23+a21 x3+a22 y3))/(a33+a31 x3+a32 y3))&&

(b13+(b11 (m13+m11 x4+m12 y4))/(m33+m31 x4+m32 y4)+(b12 (m23+m21 x4+m22 y4))/(m33+m31 x4+m32 y4))/(b33+(b31 (m13+m11 x4+m12 y4))/(m33+m31 x4+m32 y4)+(b32 (m23+m21 x4+m22 y4))/(m33+m31 x4+m32 y4))==

(m13+(m11 (a13+a11 x4+a12 y4))/(a33+a31 x4+a32 y4)+(m12 (a23+a21 x4+a22 y4))/(a33+a31 x4+a32 y4))/(m33+(m31 (a13+a11 x4+a12 y4))/(a33+a31 x4+a32 y4)+(m32 (a23+a21 x4+a22 y4))/(a33+a31 x4+a32 y4))&&

(b23+(b21 (m13+m11 x4+m12 y4))/(m33+m31 x4+m32 y4)+(b22 (m23+m21 x4+m22 y4))/(m33+m31 x4+m32 y4))/(b33+(b31 (m13+m11 x4+m12 y4))/(m33+m31 x4+m32 y4)+(b32 (m23+m21 x4+m22 y4))/(m33+m31 x4+m32 y4))==

(m23+(m21 (a13+a11 x4+a12 y4))/(a33+a31 x4+a32 y4)+(m22 (a23+a21 x4+a22 y4))/(a33+a31 x4+a32 y4))/(m33+(m31 (a13+a11 x4+a12 y4))/(a33+a31 x4+a32 y4)+(m32 (a23+a21 x4+a22 y4))/(a33+a31 x4+a32 y4))&&m33==1,{m11,m12,m13,m21,m22,m23,m31,m32}]

Share this post


Link to post
Share on other sites

Свойства матричного умножения:

1.ассоциативность (AB)C = A(BC);

2.некоммутативность (в общем случае): AB!=BA;

3.произведение коммутативно в случае умножения с единичной матрицей: AI = IA;

4.дистрибутивность: (A+B)C = AC + BC, A(B+C) = AB + AC;

5.ассоциативность и коммутативность относительно умножения на число: (λA)B = λ(AB) = A(λB ) ;

Share this post


Link to post
Share on other sites

хорошо тогда непонятно как наша координата (x,y) умножается на матрицу 3х3, так что получается формула

5e8e8984a34d8aed41e5b5020de5adf082.png

вроде как то так но непонятно почему так.

M*V

V=(X,Y,1)^t

New_X:=m11*X+m12*Y+m13*1

New_Y:=m21*X+m22*Y+m23*1

New_Z:=m31*X+m32*Y+m33*1

u=New_X/New_Z

v=New_Y/New_Z

вообщем в итоге сама то операция получается нелинейная потому что на New_Z делится.

Share this post


Link to post
Share on other sites

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

Share this post


Link to post
Share on other sites

ну получается перспективные искажения (и вроде даже аффинные) нельзя записать как произведения матриц и получается матричное уравнение не правильно.

т.е. матрица гомографии к матричным операциям никак не относится?)

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

Share this post


Link to post
Share on other sites

В книжке Multiview Geometry в качестве метода коррекции радиальной дисторсии ссылаются на алгоритм plumb line correction.

Вот например: http://arxiv.org/pdf/0810.4426v2.pdf

Еще один, более простой: http://www.isprs.org/proceedings/2006/euroCOW06/euroCOW06_files/papers/Tommaselli_Telles_line_calib_eurocow06.pdf

Share this post


Link to post
Share on other sites

я видел первый документ, точнее не точно этот, но с теми же картинками, они там использую вероятности нахождения линий как то так, я щас рассматриваю более простой случай, если известны 4 точки углов и объект предполагается прямоугольным (или например что вроде эквивалентно - известны по 2 точки на 2-х прямых одна из которых без дисторсии горизонтальна, а другая вертикальна).или же я могу разметить стороны прямоугольника точками или некоторые вертикальные\горизонтальные предполагаемые линии.

если нам известны 2 пары точек до и после = 4 уравнения, то мы можем восстановить с1,c2,k1,k2. это я правда делал в математике, пока как сделать на с++ не очень понятно. Но опять же у нас нету x,y до преобразования, есть только информация о том что точки лежат на какой то вертикальной или горизонтальной прямой, т.е. можно взять какой то примерный прямоугольник и его использовать(или же наложить доп. ограничения уравнения на точки например х1=х2=...х6), но непонятно правильно ли восстановимо будет искажение для всех точекю

т.е. вопрос сводится к решению системы уравнений для с1,c2,k1,k2 , причем они выражены не явно.

второй вопрос если считать, что х,у неизвестны как найти их примерные значения.

кстати с матрицами я немного ошибся все у которых последняя строка = 0 0 1 являют собой линейную операцию.

Share this post


Link to post
Share on other sites

что то я всё усложнил, если точки известны то параметры дисторсии выводятся


c1= -((x1*(u2-x2)*(pow(x1,2)+pow(y1,2))-(u1-x1)*x2*(pow(x2,2)+pow(y2,2)))/(-(x1*(pow(x1,2)+pow(y1,2)))+x2*(pow(x2,2)+pow(y2,2))));

c2= v1-y1+((pow(x1,2)*y1+pow(y1,3))*(-v1+v2+y1-y2))/(pow(x1,2)*y1+pow(y1,3)-pow(x2,2)*y2-pow(y2,3));

k1= -((u1-u2-x1+x2)/(-pow(x1,3)+pow(x2,3)-x1*pow(y1,2)+x2*pow(y2,2)));

k2= -((-v1+v2+y1-y2)/(pow(x1,2)*y1+pow(y1,3)-pow(x2,2)*y2-pow(y2,3)));

Share this post


Link to post
Share on other sites

еще насчет итеративного метода, если брать формулу со смещением в центр, то что то ничего не сходится.

и опять же я не понял как выводится формула.

допустим в обычном виде, без смещения

double u= c1+x+k1*((x)*(x)+(y)*(y))*(x);
отсюда просто выражается х? т.е. получается
double x= u-c1-k1*((x)*(x)+(y)*(y))*(x);
потом в правой части заменяем х,y
double x= u-c1-k1*((x_)*(x_)+(y_)*(y_))*(x_);
а если х просто так в явном виде не выводится, то тогда что? из каких соображений выбирается первое приближение? для формулы со смещением координат почему то не работает.
 // считаем u,v по x,y

    // здесь все очевидно

    double k1=2e-4;

    double k2=3e-4;

    double c1=100.0;

    double c2=120.0;


	double x= 10;

	double y= 20;

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

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

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


    cout << u << endl;

    cout << v << endl;



    // Вычисление x,y по u,v

    // выразили x и y через известные и самих себя

    // и крутим цикл пока не сойдется.


    double x_prev=0,y_prev=0;

    double x_curr=0,y_curr=0;

    double delta=DBL_MAX;

	int iter=0;

    for(;delta>DBL_EPSILON;)

    {

		/*x_curr=u-c1-k1*(x_prev*x_prev+y_prev*y_prev)*x_prev;

		y_curr=v-c2-k2*(x_prev*x_prev+y_prev*y_prev)*y_prev;*/

		x_curr=u-k1*((x_prev-c1)*(x_prev-c1)+(y_prev-c2)*(y_prev-c2))*(x_prev-c1);

		y_curr=v-k2*((x_prev-c1)*(x_prev-c1)+(y_prev-c2)*(y_prev-c2))*(y_prev-c2);

		delta=sqrt(pow((x_curr-x_prev),2)+pow((y_curr-y_prev),2));

		x_prev=x_curr;

		y_prev=y_curr;

		iter++;

    }


    cout << x_curr << endl;

    cout << y_curr << endl;

Share this post


Link to post
Share on other sites

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

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

void main(void)
{
// считаем u,v по x,y
// здесь все очевидно
double k1=2e-4;
double k2=3e-4;
double c1=100.0;
double c2=120.0;

double x= 10;
double y= 20;
/*double u= c1+x+k1*((x)*(x)+(y)*(y))*(x);
double v= c2+y+k2*((x)*(x)+(y)*(y))*(y);*/
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);

cout << u << endl;
cout << v << endl;


// Вычисление x,y по u,v
// выразили x и y через известные и самих себя
// и крутим цикл пока не сойдется.

double x_prev=0,y_prev=0;

double delta=DBL_MAX;
int iter=0;
for(;delta>DBL_EPSILON;)
{

double t1 = x - c1;
double t2 = t1 * t1;
double t3 = y - c2;
double t4 = t3 * t3;
double deltax = (k1 * (t2 + t4) * t1 + x - u) / k1 / t3 / t1 / 0.2e1;


t2 = pow(x - c1, 0.2e1);
t3 = y - c2;
t4 = t3 * t3;
double t6 = k2 * (t2 + t4);
double deltay = (t6 * t3 + y - v) / (0.2e1 * k2 * t3 * t3 + t6 + 0.1e1);

x=x-deltax;
y=y-deltay;

//cout << deltax << endl;

delta=sqrt(pow((deltax),2)+pow((deltay),2));
iter++;
//cout << delta << endl;
}

cout << x << endl;
cout << y << endl;

getchar();
}[/code]

  • Like 2

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.

×