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

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

Recommended Posts

я уже не помню как такие уравнения решать, что то типа замены, вообщем какие то трюки.

Решаем первое относительно игрека- это квадратное уравнение, значит два корня. Затем решаем второе относительно игрека, там уравнение третьей степени и три корня. Потом приравниваем результаты друг к другу, это будет 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]], а если "да", то у тебя есть неограниченный набор точек до и после преобразования и можно составлять систему.

Share this post


Link to post
Share on other sites
Они в радикалах не решаются, это было доказано Абелем ещё в 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

Share this post


Link to post
Share on other sites

помоему все таки точек. 1 точка->2 уравнения, 8 параметров(m33==1).

Ну вообще логично :). Но подставив точку мы узнаем значения x, y, u и v, поэтому прибавив первое ко второму мы получим одно уравнение с 8-ю неизвестными.

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

Ну дык сфотографируй координатную плоскость, получишь все нужные данные в избытке :). За координатную плоскость вполне сойдёт ровный тетрадный лист в клеточку.

Share this post


Link to post
Share on other sites

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

матрица поворота на 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

Share this post


Link to post
Share on other sites
Ну дык сфотографируй координатную плоскость, получишь все нужные данные в избытке . За координатную плоскость вполне сойдёт ровный тетрадный лист в клеточку.

у меня нет такой возможности.

Share this post


Link to post
Share on other sites

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

матрица поворота на 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-ёх точек мало.

Share this post


Link to post
Share on other sites
А откуда взялась матрица b?

ну я так обозначил ту, что раньше была L.

Ну тогда монитор с нужной картинкой. Или на тот прямоугольник линейкой разметку нанести, просто 4-ёх точек мало.

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

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

Share this post


Link to post
Share on other sites

Opencv - шный solve, решает такие системы, и никакого пакета не нужно.

MathCad - решает, он у Вас, судя по картинке установлен.

Я использую Maple, т.к. он может большие формулы в оптимизированный C-шный код переводить.

А вообще, да, любой мат. пакет с этим справится.

Сейчас нет времени, чуть попозже прогоню Вашу систему через Maple.

прогнал, хотя можно и без этого видеть, что решение: m_x_x=0.

Share this post


Link to post
Share on other sites

ну я так обозначил ту, что раньше была L.

А мы разве не её ищем? Короче для вычисления проективного преобразования нужны 4 пары точек, без них ничего не выйдет.

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

Я и без математики скажу, что тут решением будут нули. При правильном решении будет 8 уравнений и 9 неизвестных(преобразование находится с точностью до множителя). Либо можно сразу приравнять m33 к единице, тогда решение будет однозначным.

Share this post


Link to post
Share on other sites

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

попробовал для уравнения бочки.

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

Share this post


Link to post
Share on other sites

А мы разве не её ищем? Короче для вычисления проективного преобразования нужны 4 пары точек, без них ничего не выйдет.

Я и без математики скажу, что тут решением будут нули. При правильном решении будет 8 уравнений и 9 неизвестных(преобразование находится с точностью до множителя). Либо можно сразу приравнять m33 к единице, тогда решение будет однозначным.

прогнал, хотя можно и без этого видеть, что решение: m_x_x=0.

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

кстати почему то даже если поставить m33==1, то математика не выдает ответа, это значит нет решений?

Share this post


Link to post
Share on other sites
Solve[v == c2 + y (1 + k2 (x^2 + y^2)) &&

u == c1 + x (1 + k2 (x^2 + y^2)), {x, y}]

если попробовать решить относительно, x,y то тоже выдаёт много.

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}]

{{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 гига оперативки и я так и не дождался конца.

Share this post


Link to post
Share on other sites

mrgloom может задачу уже наконецто опишешь? А то математик из тебя некудышный.

Откуда у тебя взялись точки? Почему дисторсия и перспектива. Откуда вращение? И вообще почему матрицы 3х3?

Начни с того что известно, а чего нет.

Share this post


Link to post
Share on other sites

есть изображение с дисторсией в общем случае перспективное искажение + искажение бочка.

на изображении есть объект предположительно прямоугольной формы(когда он без искажений), на реальном изображении он искажён.

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

надо найти все параметры искажений

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.

Share this post


Link to post
Share on other sites

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

Share this post


Link to post
Share on other sites

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

может прогоните в мапл это уравнение, чтобы там выдало оптимальную форму?

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

Share this post


Link to post
Share on other sites
это уже из немного другой задачи, берем изображение снимаем его фото, потом (в реальном мире) крутим его на градус(например 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

на изображении есть объект предположительно прямоугольной формы(когда он без искажений), на реальном изображении он искажён.

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

надо найти все параметры искажений

Вначале устрани дисторсию а потом уже перспективу.

Берешь точки на кривой строишь премую искаженную дисторсией вычитаешь и минимизируешь ошибку. МНК тут плохо будет рабоать, но пойдёт. А затем уже перспективу.

Share this post


Link to post
Share on other sites
распечатали повернули и снова сняли

так нету возможности распечатать же

есть только возможность повернуть на произвольный известный градус в реальном мире предмет относительно неподвижной камеры.

т.е. не

(x,y,z)*M*L=((x,y,z)*M)*Rot*M

а

(x,y,z)*M*L=(x,y,z)*Rot*M

Share this post


Link to post
Share on other sites
Вначале устрани дисторсию а потом уже перспективу.

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

как обычно бывает в оптических системах?

Берешь точки на кривой строишь премую искаженную дисторсией вычитаешь и минимизируешь ошибку. МНК тут плохо будет рабоать, но пойдёт. А затем уже перспективу.

ну а саму то функцию которую приближаем как задаём?

по идее это должно быть u(v)

но как решить?

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

а чем параметрический вид не устраивает?

Если ищем коэффициенты, какая разница что за уравнения мы ему даем?

Заменяем u и v на цифры, заменяем x,y на цифры, получаем обычную переопределенную СЛАУ относительно коэффициентов.

Решаем через SVD (получаем решение МНК).

Share this post


Link to post
Share on other sites

во-первых что имеется ввиду под МНК?

такое?

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

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

Share this post


Link to post
Share on other sites

Под МНК имеется ввиду метод, минимизирующий квадрат ошибки (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]

Share this post


Link to post
Share on other sites

решил попробовать написать всё свое через 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)

в чем ошибка в логике?

Share this post


Link to post
Share on other sites

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

запихиваем в математику (пробовал 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));

}
для теста использовал картинки %D1%81%D0%B5%D1%82%D0%BA%D0%B0_%D0%B1%D0%BE%D1%87%D0%BA%D0%B0_%D0%BC%D0%B5%D1%82%D0%BA%D0%B0.png%D1%81%D0%B5%D1%82%D0%BA%D0%B0_%D0%BC%D0%B5%D1%82%D0%BA%D0%B0.png вызываю
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;

}

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

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

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.

×