Euler 0 Жалоба Опубликовано July 24, 2012 я уже не помню как такие уравнения решать, что то типа замены, вообщем какие то трюки. Решаем первое относительно игрека- это квадратное уравнение, значит два корня. Затем решаем второе относительно игрека, там уравнение третьей степени и три корня. Потом приравниваем результаты друг к другу, это будет 6 уравнений сводящихся к уравнениям шестой степени. Они в радикалах не решаются, это было доказано Абелем ещё в 19-ом веке. как это сдетектировать что нельзя? Не понял вопроса. я же говорю, что точек до у меня нету. Тогда ничего не поделаешь. только не точек а уравнений, 1 точка даёт 2 уравнения. Именно точек. Одна точка даст одно уравнение с 8-ю неизвестными(для гомографии). да я что то не пойму как. вот есть у нас изображение эталон A первое изображение будет Dist[A] второе Dist[Rot_90[A]] я могу получить преобразование из Dist[A] в Dist[Rot_90[A]] по 4 точкам, а мне надо получить само преобразование Dist. Чего-то я не понимаю- ты можешь выполнить Dist[A] для произвольного набора точек? Если "нет", то как ты получишь Dist[Rot_90[A]], а если "да", то у тебя есть неограниченный набор точек до и после преобразования и можно составлять систему. Поделиться сообщением Ссылка на сообщение Поделиться на других сайтах
mrgloom 242 Жалоба Опубликовано July 24, 2012 Они в радикалах не решаются, это было доказано Абелем ещё в 19-ом веке. спасибо не знал. это и был ответ на следующий вопрос. только не точек а уравнений, 1 точка даёт 2 уравнения. Именно точек. Одна точка даст одно уравнение с 8-ю неизвестными(для гомографии). u= (m11*x + m12*y + m13)/(m31*x + m32*y + m33); v= (m21*x + m22*y + m23)/(m31*x + m32*y + m33); помоему все таки точек. 1 точка->2 уравнения, 8 параметров(m33==1). Чего-то я не понимаю- ты можешь выполнить Dist[A] для произвольного набора точек? Если "нет", то как ты получишь Dist[Rot_90[A]], а если "да", то у тебя есть неограниченный набор точек до и после преобразования и можно составлять систему. беру объект кладу под камеру снимаю, поворачиваю объект на 90 градусов опять снимаю, какое преобразование не знаю, но могу указать 4 пары точек соответствия на этих 2-х снимках. в итоге получается матричное уравнение Rot*M=M*L Поделиться сообщением Ссылка на сообщение Поделиться на других сайтах
Euler 0 Жалоба Опубликовано July 24, 2012 помоему все таки точек. 1 точка->2 уравнения, 8 параметров(m33==1). Ну вообще логично . Но подставив точку мы узнаем значения x, y, u и v, поэтому прибавив первое ко второму мы получим одно уравнение с 8-ю неизвестными. беру объект кладу под камеру снимаю, поворачиваю объект на 90 градусов опять снимаю, какое преобразование не знаю, но могу указать 4 пары точек соответствия на этих 2-х снимках. Ну дык сфотографируй координатную плоскость, получишь все нужные данные в избытке . За координатную плоскость вполне сойдёт ровный тетрадный лист в клеточку. Поделиться сообщением Ссылка на сообщение Поделиться на других сайтах
mrgloom 242 Жалоба Опубликовано July 24, 2012 если я правильно понимаю, то получается уравнение матрица поворота на 90 град 0 1 0 1 0 0 0 0 1 матрица которую надо найти m11 m12 m13 m21 m22 m23 m31 m32 m33 матрица перехода от одного изображения к другому (известная) b11 b12 b13 b21 b22 b23 b31 b32 b33 Rot*M=M*B получается система уравнений если я правильно перемножил(где её можно решить не в ручную?) m21=m11*b11+m12*b21+m13*b31 m22=m11*b12+m12*b22+m13*b32 m23=m11*b13+m12*b23+m13*b33 m11=m21*b11+m22*b21+m23*b31 m12=m21*b12+m22*b22+m23*b32 m13=m21*b13+m22*b23+m23*b33 m31=m31*b11+m32*b21+m33*b31 m32=m31*b12+m32*b22+m33*b32 m33=m31*b13+m32*b23+m33*b33 Поделиться сообщением Ссылка на сообщение Поделиться на других сайтах
mrgloom 242 Жалоба Опубликовано July 24, 2012 Ну дык сфотографируй координатную плоскость, получишь все нужные данные в избытке . За координатную плоскость вполне сойдёт ровный тетрадный лист в клеточку. у меня нет такой возможности. Поделиться сообщением Ссылка на сообщение Поделиться на других сайтах
Euler 0 Жалоба Опубликовано July 24, 2012 если я правильно понимаю, то получается уравнение матрица поворота на 90 град 0 1 0 1 0 0 0 0 1 матрица которую надо найти m11 m12 m13 m21 m22 m23 m31 m32 m33 матрица перехода от одного изображения к другому (известная) b11 b12 b13 b21 b22 b23 b31 b32 b33 А откуда взялась матрица b? получается система уравнений если я правильно перемножил(где её можно решить не в ручную?) 9 уравнений? Там точно что-то не так . Я использую математику, но вообще с линейной системой любой мат. пакет справится. у меня нет такой возможности. Ну тогда монитор с нужной картинкой. Или на тот прямоугольник линейкой разметку нанести, просто 4-ёх точек мало. Поделиться сообщением Ссылка на сообщение Поделиться на других сайтах
mrgloom 242 Жалоба Опубликовано July 24, 2012 А откуда взялась матрица b? ну я так обозначил ту, что раньше была L. Ну тогда монитор с нужной картинкой. Или на тот прямоугольник линейкой разметку нанести, просто 4-ёх точек мало. я ж говорю, что нет такой возможности по условию. представьте что камера в космосе. может вы загоните в свою математику? а то прям щас нету возможности скачать\поставить. Поделиться сообщением Ссылка на сообщение Поделиться на других сайтах
Smorodov 579 Жалоба Опубликовано July 24, 2012 Opencv - шный solve, решает такие системы, и никакого пакета не нужно. MathCad - решает, он у Вас, судя по картинке установлен. Я использую Maple, т.к. он может большие формулы в оптимизированный C-шный код переводить. А вообще, да, любой мат. пакет с этим справится. Сейчас нет времени, чуть попозже прогоню Вашу систему через Maple. прогнал, хотя можно и без этого видеть, что решение: m_x_x=0. Поделиться сообщением Ссылка на сообщение Поделиться на других сайтах
Euler 0 Жалоба Опубликовано July 24, 2012 ну я так обозначил ту, что раньше была L. А мы разве не её ищем? Короче для вычисления проективного преобразования нужны 4 пары точек, без них ничего не выйдет. может вы загоните в свою математику? а то прям щас нету возможности скачать\поставить. Я и без математики скажу, что тут решением будут нули. При правильном решении будет 8 уравнений и 9 неизвестных(преобразование находится с точностью до множителя). Либо можно сразу приравнять m33 к единице, тогда решение будет однозначным. Поделиться сообщением Ссылка на сообщение Поделиться на других сайтах
mrgloom 242 Жалоба Опубликовано July 25, 2012 ну скачал вчера математику, выглядит приятно,раньше ничем таким не пользовался, но насслышан. попробовал для уравнения бочки. зависимость от 4 параметров, т.е. известны соответсвия 2-х точек(x,y)->(u,v) т.е. каким то образом надо выбрать точки(u,v есть надо выбрать x,y) чтобы откатить бочку, я щас не очень представляю, что будет если их выбрать произвольно, скорее всего ничего хорошего, возможно надо наложить какие либо доп.ограничения? 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))}} а потом я попробовал для 6 неизвестных. Solve[y1 == c2 + (h/2) (1 + k2 ((w/2)^2 + (h/2)^2)) && y2 == c2 + (-h/2) (1 + k2 ((w/2)^2 + (h/2)^2)) && y3 == c2 + (h/2) (1 + k2 (h/2)^2) && x1 == c1 + (w/2) (1 + k1 ((w/2)^2 + (h/2)^2)) && x2 == c1 + (-w/2) (1 + k2 ((w/2)^2 + (h/2)^2)) && x3 == c1 + (w/2) (1 + k2 (w/2)^2), {c1, c2, k1, k2, w, h}] и ответ выдало, но какой то совершенно дикий на несколько листов и непонятно это ответ всё таки точный и как его можно по быстрому перевести в код(maplе может?) или хотя бы в функции в самой математике чтобы было v(x,y) u(x,y) или наоборот x(u,v) y(u,v) или не в параметрическом виде, хотя вроде выше писалось, что это невозможно. Поделиться сообщением Ссылка на сообщение Поделиться на других сайтах
mrgloom 242 Жалоба Опубликовано July 25, 2012 А мы разве не её ищем? Короче для вычисления проективного преобразования нужны 4 пары точек, без них ничего не выйдет. Я и без математики скажу, что тут решением будут нули. При правильном решении будет 8 уравнений и 9 неизвестных(преобразование находится с точностью до множителя). Либо можно сразу приравнять m33 к единице, тогда решение будет однозначным. прогнал, хотя можно и без этого видеть, что решение: m_x_x=0. ну так и о чём это говорит? что невозможно восстановить или я неправильно составил уравнение. кстати почему то даже если поставить m33==1, то математика не выдает ответа, это значит нет решений? Поделиться сообщением Ссылка на сообщение Поделиться на других сайтах
mrgloom 242 Жалоба Опубликовано July 25, 2012 Solve[v == c2 + y (1 + k2 (x^2 + y^2)) && u == c1 + x (1 + k2 (x^2 + y^2)), {x, y}] если попробовать решить относительно, x,y то тоже выдаёт много. Поделиться сообщением Ссылка на сообщение Поделиться на других сайтах
mrgloom 242 Жалоба Опубликовано July 25, 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}] {{m11 -> 0, m12 -> 0, m13 -> 0, m21 -> 0, m22 -> 0, m23 -> 0, m31 -> 0, m32 -> 0, m33 -> 0}} даже в общем виде выдает нули так почему не верно соотношение A*M=M*B1 или же можно еще и так записать A*M*B2=M? попробовал еще Solve[u == (m11*x + m12*y + m13)/(m31*x + m32*y + m33) + c1 + x + k1*(x*x + y*y)*x && v == (m21*x + m22*y + m23)/(m31*x + m32*y + m33) + c2 + y + k2*(x*x + y*y)*y, {x, y}] съело 2 гига оперативки и я так и не дождался конца. Поделиться сообщением Ссылка на сообщение Поделиться на других сайтах
Pavia00 32 Жалоба Опубликовано July 25, 2012 mrgloom может задачу уже наконецто опишешь? А то математик из тебя некудышный. Откуда у тебя взялись точки? Почему дисторсия и перспектива. Откуда вращение? И вообще почему матрицы 3х3? Начни с того что известно, а чего нет. Поделиться сообщением Ссылка на сообщение Поделиться на других сайтах
mrgloom 242 Жалоба Опубликовано July 25, 2012 есть изображение с дисторсией в общем случае перспективное искажение + искажение бочка. на изображении есть объект предположительно прямоугольной формы(когда он без искажений), на реальном изображении он искажён. точки я расставил вручную по периметру контура объекта. надо найти все параметры искажений u = (m11*x + m12*y + m13)/(m31*x + m32*y + m33) + c1 + x + k1*(x*x + y*y)*x; v = (m21*x + m22*y + m23)/(m31*x + m32*y + m33) + c2 + y + k2*(x*x + y*y)*y; кстати я тут походу ошибся,ибо если искажения накладываются, то это значит они применяются последовательно и это не означает, что в формуле для координаты они просто складываются. еще не очень помню как это в терминологии, но неплохо было бы еще определить влияет ли последовательность операций f(g(x))!=g(f(x)) т.е. например для сдвига и поворота без разницы в каком порядке их применять, а вот если 2 перспективных искажения подряд, то похоже уже разница есть. И вообще почему матрицы 3х3? это уже из немного другой задачи, берем изображение снимаем его фото, потом (в реальном мире) крутим его на градус(например 90) и опять снимаем, имеем 2 изображения и допустим имеем матрицу L 3х3 перехода между ними, значит (x,y)*Rot*M=(x,y)*M*L ну и получаем матричное уравнение Rot*M=M*L где Rot-матрица поворота,M-неизвестное искажение которое нам надо найти,L-преобразование между двумя картинками. или тут 2 формы Rot*M=M*L1 Rot*M*L2=M смотря трансформация какого относительно какого изображения L1 или L2. Поделиться сообщением Ссылка на сообщение Поделиться на других сайтах
Smorodov 579 Жалоба Опубликовано July 25, 2012 Лучше когда матрицы 4х4 (однородные координаты), это позволяет избежать делений на 0. Поделиться сообщением Ссылка на сообщение Поделиться на других сайтах
mrgloom 242 Жалоба Опубликовано July 25, 2012 эмм я не знаю как в такой системе координат работать и я не понял даже как это поможет. может прогоните в мапл это уравнение, чтобы там выдало оптимальную форму? 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))}} Поделиться сообщением Ссылка на сообщение Поделиться на других сайтах
Pavia00 32 Жалоба Опубликовано July 25, 2012 это уже из немного другой задачи, берем изображение снимаем его фото, потом (в реальном мире) крутим его на градус(например 90) и опять снимаем, имеем 2 изображения и допустим имеем матрицу L 3х3 перехода между ними, значит (x,y)*Rot*M=(x,y)*M*L ну и получаем матричное уравнение Rot*M=M*L где Rot-матрица поворота,M-неизвестное искажение которое нам надо найти,L-преобразование между двумя картинками. или тут 2 формы Rot*M=M*L1 Rot*M*L2=M (x,y,z)*Rot*M=(x,y,z)*M*L M это преобразование камеры. (x,y,z)*M Это то что мы получили после снимка. распечатали повернули и снова сняли ((x,y,z)*M)*Rot*M И по точкам нашли переход между этими двумя снимками (x,y,z)*M*L=((x,y,z)*M)*Rot*M ((x,y,z)*M)^-1*((x,y,z)*M)*L=((x,y,z)*M)^-1*((x,y,z)*M)*Rot*M L=Rot*M - это в принципе и так из условия было ясно M=Rot^-1*L на изображении есть объект предположительно прямоугольной формы(когда он без искажений), на реальном изображении он искажён. точки я расставил вручную по периметру контура объекта. надо найти все параметры искажений Вначале устрани дисторсию а потом уже перспективу. Берешь точки на кривой строишь премую искаженную дисторсией вычитаешь и минимизируешь ошибку. МНК тут плохо будет рабоать, но пойдёт. А затем уже перспективу. Поделиться сообщением Ссылка на сообщение Поделиться на других сайтах
mrgloom 242 Жалоба Опубликовано July 25, 2012 распечатали повернули и снова сняли так нету возможности распечатать же есть только возможность повернуть на произвольный известный градус в реальном мире предмет относительно неподвижной камеры. т.е. не (x,y,z)*M*L=((x,y,z)*M)*Rot*M а (x,y,z)*M*L=(x,y,z)*Rot*M Поделиться сообщением Ссылка на сообщение Поделиться на других сайтах
mrgloom 242 Жалоба Опубликовано July 25, 2012 Вначале устрани дисторсию а потом уже перспективу. опять же от порядка наложения ведь зависит и порядок устранения? как обычно бывает в оптических системах? Берешь точки на кривой строишь премую искаженную дисторсией вычитаешь и минимизируешь ошибку. МНК тут плохо будет рабоать, но пойдёт. А затем уже перспективу. ну а саму то функцию которую приближаем как задаём? по идее это должно быть u(v) но как решить? u = с1 + x + k1*(x*x + y*y)*x; v = c2 + y + k2*(x*x + y*y)*y; или можно мнк прямо как то в параметрическом виде решать? Поделиться сообщением Ссылка на сообщение Поделиться на других сайтах
Smorodov 579 Жалоба Опубликовано July 25, 2012 а чем параметрический вид не устраивает? Если ищем коэффициенты, какая разница что за уравнения мы ему даем? Заменяем u и v на цифры, заменяем x,y на цифры, получаем обычную переопределенную СЛАУ относительно коэффициентов. Решаем через SVD (получаем решение МНК). Поделиться сообщением Ссылка на сообщение Поделиться на других сайтах
mrgloom 242 Жалоба Опубликовано July 26, 2012 во-первых что имеется ввиду под МНК? такое? http://www.exponenta.ru/educat/class/courses/vvm/theme_7/example.asp во-вторых я ж говорю что (x,y) нету есть только (u,v). допустим я могу взять v(u) одну дугу, но как потом восстановить по ней коэффициенты дисторсии в общем уравнении непонятно. про SVD можете еще пояснить? какое-то уравнение записывается в матричной форме? http://ru.wikipedia.org/wiki/%D0%A1%D0%B8%D0%BD%D0%B3%D1%83%D0%BB%D1%8F%D1%80%D0%BD%D0%BE%D0%B5_%D1%80%D0%B0%D0%B7%D0%BB%D0%BE%D0%B6%D0%B5%D0%BD%D0%B8%D0%B5 через свд получаем псевдообратную матрицу и она как то используется в мнк? что то я не осилил. Поделиться сообщением Ссылка на сообщение Поделиться на других сайтах
Smorodov 579 Жалоба Опубликовано July 26, 2012 Под МНК имеется ввиду метод, минимизирующий квадрат ошибки (X-X_ground_true)^2 . Ну если нет x,y то я тут умыл руки Про SVD все просто. Там действительно решается система линейных уравнений, записанная в матричной форме Ax=b. Обратить классическим методом можно только квадратную матрицу. Если А квадратная, то получим точное решение. Существует метод решать переопределенные системы (как у Вас). При этом, мы не можем найти точное решение. И матрица А у нас прямоугольная. Соответственно используем метод SVD и получаем, псевдообратную матрцу. Решение системы для квадратной и прямоугольной матриц будет записываться одинаково: x=A.inv()*b; Но в случае с прямоугольной матрицей это будет не точное решение (и инверсия будет псевдоинверсией), а приближенное минимизирующее квадрат ошибки, так как точное в таком случае найти нельзя. Вот, для пояснения, реализация псевдоинверсии (для библиотеки Eigen, которая входит в opencv). Opencv и так умеет её искать (cv::invert с опцией DECOMP_SVD). И системы переопределенные решать умеет (см. cv::solve там есть флаги: DECOMP_SVD Singular value decomposition (SVD) method; the system can be over-defined and/or the matrix src1 can be singular и DECOMP_QR QR factorization; the system can be over-defined and/or the matrix src1 can be singular ) //----------------------------------------------------------------------------------------- // Псевдоинверсия матрицы Moore-Penrose //----------------------------------------------------------------------------------------- bool pinv(MatrixXf& a,MatrixXf& a_pinv) { //if(a.rows()<a.cols()){return false;} JacobiSVD<MatrixXf> svd(a, ComputeThinU | ComputeThinV); MatrixXf vSingular; MatrixXf vPseudoInvertedSingular(svd.matrixV().cols(),1); vSingular=svd.singularValues(); for (int iRow =0; iRow<vSingular.rows(); iRow++) { if ( fabs(vSingular(iRow))<=FLT_EPSILON ) { vPseudoInvertedSingular(iRow,0)=0.0f; } else { vPseudoInvertedSingular(iRow,0)=1./vSingular(iRow); } } a_pinv = (svd.matrixV()*vPseudoInvertedSingular.asDiagonal()*svd.matrixU().transpose()); return true; } [/code] Поделиться сообщением Ссылка на сообщение Поделиться на других сайтах
mrgloom 242 Жалоба Опубликовано July 26, 2012 решил попробовать написать всё свое через remap, чтобы проверить. перспективное искажение IplImage* perspective_transfrom(IplImage* img, CvMat* warp_matrix) { double m11= cvmGet(warp_matrix,0,0); // Set M(i,j) double m12= cvmGet(warp_matrix,0,1); // Set M(i,j) double m13= cvmGet(warp_matrix,0,2); // Set M(i,j) double m21= cvmGet(warp_matrix,1,0); // Set M(i,j) double m22= cvmGet(warp_matrix,1,1); // Set M(i,j) double m23= cvmGet(warp_matrix,1,2); // Set M(i,j) double m31= cvmGet(warp_matrix,2,0); // Set M(i,j) double m32= cvmGet(warp_matrix,2,1); // Set M(i,j) double m33= cvmGet(warp_matrix,2,2); // Set M(i,j) //прямое преобразование //simple Solve /*u= (m13 + m11*x + m12*y)/(m33 + m31*x + m32*y); v= (m23 + m21*x + m22*y)/(m33 + m31*x + m32*y);*/ 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++) { for (int x = 0; x < w; x++) { float u= (m13 + m11*x + m12*y)/(m33 + m31*x + m32*y); *pbuf = u; ++pbuf; } } pbuf = (float*)mapy->imageData; for (int y = 0; y < h; y++) { for (int x = 0; x < w; x++) { float v= (m23 + m21*x + m22*y)/(m33 + m31*x + m32*y); *pbuf = v; ++pbuf; } } IplImage* temp = cvCloneImage(img); cvRemap( temp, img, mapx, mapy ); cvReleaseImage(&temp); cvReleaseImage(&mapx); cvReleaseImage(&mapy); return img; } IplImage* un_perspective_transfrom(IplImage* img, CvMat* warp_matrix) { double m11= cvmGet(warp_matrix,0,0); // Set M(i,j) double m12= cvmGet(warp_matrix,0,1); // Set M(i,j) double m13= cvmGet(warp_matrix,0,2); // Set M(i,j) double m21= cvmGet(warp_matrix,1,0); // Set M(i,j) double m22= cvmGet(warp_matrix,1,1); // Set M(i,j) double m23= cvmGet(warp_matrix,1,2); // Set M(i,j) double m31= cvmGet(warp_matrix,2,0); // Set M(i,j) double m32= cvmGet(warp_matrix,2,1); // Set M(i,j) double m33= cvmGet(warp_matrix,2,2); // Set M(i,j) //обратное преобразование //x= -((-(m13*m22) + m12*m23 - m23*m32*u + m22*m33*u + m13*m32*v - m12*m33*v)/ // (m12*m21 - m11*m22 + m22*m31*u - m21*m32*u - m12*m31*v + m11*m32*v)); //y= -((m13*m21 - m11*m23 + m23*m31*u - m21*m33*u - m13*m31*v + m11*m33*v)/ // (m12*m21 - m11*m22 + m22*m31*u - m21*m32*u - m12*m31*v + m11*m32*v)); 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 v = 0; v < h; v++) { for (int u = 0; u < w; u++) { float x= x= -((-(m13*m22) + m12*m23 - m23*m32*u + m22*m33*u + m13*m32*v - m12*m33*v)/ (m12*m21 - m11*m22 + m22*m31*u - m21*m32*u - m12*m31*v + m11*m32*v)); *pbuf = x; ++pbuf; } } pbuf = (float*)mapy->imageData; for (int v = 0;v < h; v++) { for (int u = 0; u < w; u++) { float y= -((m13*m21 - m11*m23 + m23*m31*u - m21*m33*u - m13*m31*v + m11*m33*v)/ (m12*m21 - m11*m22 + m22*m31*u - m21*m32*u - m12*m31*v + m11*m32*v)); *pbuf = y; ++pbuf; } } IplImage* temp = cvCloneImage(img); cvRemap( temp, img, mapx, mapy ); cvReleaseImage(&temp); cvReleaseImage(&mapx); cvReleaseImage(&mapy); return img; } но почему то прямое и обратное на самом деле действуют наоборот.(т.е. как бы поменяны местами) т.е. un_perspective_transfrom соотвествует следующему коду через опенцв IplImage* warp_rect(IplImage* src, CvMat* warp_matrix) { IplImage* dst=cvCloneImage(src); cvWarpPerspective( src, dst, warp_matrix); return dst; } а если сделать cvSaveImage("test_matrix_.png",perspective_transfrom(un_perspective_transfrom(src, warp_matrix),warp_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) в чем ошибка в логике? Поделиться сообщением Ссылка на сообщение Поделиться на других сайтах
mrgloom 242 Жалоба Опубликовано July 26, 2012 опять же если взять дисторсию бочка и предположить, что значение пар точек до и после известны и попробовать по этим точкам восстановить параметры дисторсии. запихиваем в математику (пробовал 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; } в итоге параметры определяются не правильно, возможно математика выдала какие то не те ответы? или из-за того что точки я отмечал вручную,т.е. это приблизительные значения, а не реально известные теоретические ,всё расходится? Поделиться сообщением Ссылка на сообщение Поделиться на других сайтах