Smorodov 579 Жалоба Опубликовано July 26, 2012 Ух,ниасилил я обычно просто обратную подстановку использую. Подставляю именно те значения, по которым получал зависимость, и смотрю верно ли равентсво в уравнениях. то есть подставьте все в u = с1 + x + k1*(x*x + y*y)*x; v = c2 + y + k2*(x*x + y*y)*y; и посмотрите равны ли правые части левым. Поделиться сообщением Ссылка на сообщение Поделиться на других сайтах
mrgloom 242 Жалоба Опубликовано July 26, 2012 кажется понял в чем ошибка, забыл делать смещение до центра (с1,с2) но в такой форме если решать, то математика почему то опять зависает и не развисает. ну тогда я попробовал решить это дело относительно k1,k2 , потому что центр с1,с2 мы можем и так найти вручную путем нахождения максимально удаленных точек на противоположных дугах и проведения по этим точкам прямым, а на пересечении и будет искомый центр. проверил на ассиметричной бочке.в итоге ошибся на 1 пиксель по х. 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, возможно не точно взял точку на рисунке или просто нельзя добиться большей точности из-за вычислений? про подстановку не разобрался пока как это сделать в математике. Поделиться сообщением Ссылка на сообщение Поделиться на других сайтах
Smorodov 579 Жалоба Опубликовано July 26, 2012 подстановку можно и на бумажке сделать. допустим алгоритм для 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; получаем что левая часть не равна правой, значит ошиблись. Поделиться сообщением Ссылка на сообщение Поделиться на других сайтах
mrgloom 242 Жалоба Опубликовано July 26, 2012 написал для тестирования точки 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}] Поделиться сообщением Ссылка на сообщение Поделиться на других сайтах
Smorodov 579 Жалоба Опубликовано July 26, 2012 Самый простой способ - соорудить итеративную схему, как в школе учили // считаем 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] Поделиться сообщением Ссылка на сообщение Поделиться на других сайтах
mrgloom 242 Жалоба Опубликовано July 27, 2012 попробую для 1 точки эту итеративную схему, но что то мне кажется что для разных точек могут определится разные параметры. еще по поводу уравнения A*X=X*B нашел в книжке про решение такого уравнения там есть такой комментарий только всё же не понятно, почему выдало всё 0, если уравнения были записаны в общем виде с переменными a,b. произвольный объект и произвольное перспективное искажение для иллюстрации задачи. я всё таки так и не пойму реально ли восстановить искажение или нет. произвольный объект X произвольный объект повернутый на 90 градусов X*Rot произвольный объект с искажением X*M произвольный объект повернутый на 90 с искажением X*Rot*M Поделиться сообщением Ссылка на сообщение Поделиться на других сайтах
mrgloom 242 Жалоба Опубликовано July 27, 2012 попробовал еще в гимпе откатить свою бочку, но там какие то странные параметры. Filters->Distorts->Lens Distortion есть ползунок Main, а есть Edge но зато можно попробовать найти сорцы вроде как еще можно откатить бочку через cvInitUndistortMap Поделиться сообщением Ссылка на сообщение Поделиться на других сайтах
Smorodov 579 Жалоба Опубликовано July 27, 2012 У этой системы 6 решений, итеративная схема дает одно, оно зависит от начального приближения, но в большинстве случаев это как раз то решение, которое мы ищем. Я пробовал (в куске кода выше) посчитать координаты точек, все сходится к правильному результату. Поделиться сообщением Ссылка на сообщение Поделиться на других сайтах
mrgloom 242 Жалоба Опубликовано July 27, 2012 а как посчитать матрицу гомографии для поворота на произвольный угол? тут пишут 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 но не написано как найти Поделиться сообщением Ссылка на сообщение Поделиться на других сайтах
Smorodov 579 Жалоба Опубликовано July 27, 2012 Аффинные преобразования: 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 Поделиться сообщением Ссылка на сообщение Поделиться на других сайтах
mrgloom 242 Жалоба Опубликовано July 27, 2012 пробую все таки решить это матричное уравнение если записать полностью, выдает все нули 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); } Поделиться сообщением Ссылка на сообщение Поделиться на других сайтах
mrgloom 242 Жалоба Опубликовано July 27, 2012 Аффинные преобразования: 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 а что под родригисом?) там что то про калибровку вроде. Поделиться сообщением Ссылка на сообщение Поделиться на других сайтах
mrgloom 242 Жалоба Опубликовано July 30, 2012 посмотрел сорцы и там такая интересная функция была 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 по часовой, а еще и отражает. Поделиться сообщением Ссылка на сообщение Поделиться на других сайтах
mrgloom 242 Жалоба Опубликовано July 30, 2012 еще про уравнение A*X=X*B учитывая ограничения налагаемые на матрицы А и М А- матрица поворота,значит у неё последняя строка = 0 0 1 M- матрица перспективного искажения у нее m33==1 ( и b33==1 кстати тоже) получаем 3 последних уравнения из системы т.е. получается, что решений нет или это просто даёт доп. ограничения на b(i,j)? Поделиться сообщением Ссылка на сообщение Поделиться на других сайтах
mrgloom 242 Жалоба Опубликовано July 30, 2012 попробовал выразить теоретически 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)); Поделиться сообщением Ссылка на сообщение Поделиться на других сайтах
mrgloom 242 Жалоба Опубликовано July 30, 2012 наблюдение о том т.е. 2 последовательные операции невозможно проделать за 1 операцию, заменив произведением матриц и ваш вопрос навёл меня на мысль о том, что матричное уравнение составлено не верно. т.е. тут не подходит перемножение матриц (оператор как бы не линейный (не знаю как это звучит строго математически)) т.е. в общем виде уравнение надо записывать так а преобразование записываются в такой форме прямое и обратное можно составить уравнение для 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}] Поделиться сообщением Ссылка на сообщение Поделиться на других сайтах
Smorodov 579 Жалоба Опубликовано July 30, 2012 Свойства матричного умножения: 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 ) ; Поделиться сообщением Ссылка на сообщение Поделиться на других сайтах
mrgloom 242 Жалоба Опубликовано July 30, 2012 хорошо тогда непонятно как наша координата (x,y) умножается на матрицу 3х3, так что получается формула вроде как то так но непонятно почему так. 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 делится. Поделиться сообщением Ссылка на сообщение Поделиться на других сайтах
Smorodov 579 Жалоба Опубликовано July 30, 2012 А эта формула и не является результатом умножения матриц, как Вы правильно заметили, она нелинейна. Поделиться сообщением Ссылка на сообщение Поделиться на других сайтах
mrgloom 242 Жалоба Опубликовано July 31, 2012 ну получается перспективные искажения (и вроде даже аффинные) нельзя записать как произведения матриц и получается матричное уравнение не правильно. т.е. матрица гомографии к матричным операциям никак не относится?) но тогда встает вопрос как решать не матричное уравнение, опять итеративным методом? там получается 8 неизвестных, надо оптимизировать по каждой переменной отдельно, остальные фиксируя или как? как вообще эти методы называются? Поделиться сообщением Ссылка на сообщение Поделиться на других сайтах
Smorodov 579 Жалоба Опубликовано July 31, 2012 В книжке 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 Поделиться сообщением Ссылка на сообщение Поделиться на других сайтах
mrgloom 242 Жалоба Опубликовано July 31, 2012 я видел первый документ, точнее не точно этот, но с теми же картинками, они там использую вероятности нахождения линий как то так, я щас рассматриваю более простой случай, если известны 4 точки углов и объект предполагается прямоугольным (или например что вроде эквивалентно - известны по 2 точки на 2-х прямых одна из которых без дисторсии горизонтальна, а другая вертикальна).или же я могу разметить стороны прямоугольника точками или некоторые вертикальные\горизонтальные предполагаемые линии. если нам известны 2 пары точек до и после = 4 уравнения, то мы можем восстановить с1,c2,k1,k2. это я правда делал в математике, пока как сделать на с++ не очень понятно. Но опять же у нас нету x,y до преобразования, есть только информация о том что точки лежат на какой то вертикальной или горизонтальной прямой, т.е. можно взять какой то примерный прямоугольник и его использовать(или же наложить доп. ограничения уравнения на точки например х1=х2=...х6), но непонятно правильно ли восстановимо будет искажение для всех точекю т.е. вопрос сводится к решению системы уравнений для с1,c2,k1,k2 , причем они выражены не явно. второй вопрос если считать, что х,у неизвестны как найти их примерные значения. кстати с матрицами я немного ошибся все у которых последняя строка = 0 0 1 являют собой линейную операцию. Поделиться сообщением Ссылка на сообщение Поделиться на других сайтах
mrgloom 242 Жалоба Опубликовано July 31, 2012 что то я всё усложнил, если точки известны то параметры дисторсии выводятся 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))); Поделиться сообщением Ссылка на сообщение Поделиться на других сайтах
mrgloom 242 Жалоба Опубликовано August 1, 2012 еще насчет итеративного метода, если брать формулу со смещением в центр, то что то ничего не сходится. и опять же я не понял как выводится формула. допустим в обычном виде, без смещения 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; Поделиться сообщением Ссылка на сообщение Поделиться на других сайтах
Smorodov 579 Жалоба Опубликовано August 1, 2012 До этого был метод простых итераций Здесь приведен Метод Ньютона: 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] 2 Поделиться сообщением Ссылка на сообщение Поделиться на других сайтах