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

[mrgloom 242]

есть набор точек типа такого.

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

u= (m11*x + m12*y + m13)/(m31*x + m32*y + m33)+С1+x+К1*(x*x+y*y)*x;

v= (m21*x + m22*y + m23)/(m31*x + m32*y + m33)+C2+y+К2*(x*x+y*y)*y;

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

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

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

[Smorodov 578]

Если там прямоугольник, то я бы, при отсутствии x,y, делал так:

1) нашел аппроксимирующие прямые для каждой грани, методом наименьших квадратов.

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

3) при заданных x,y и u,v, получаем обычную СЛАУ, которую и решаем через SVD, т.к. система получится переопределенной (кол-во уравнений > кол-ва неизвестных).

[mrgloom 242]

там прямоугольник

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

дуги из которых состоит фигура на картинке я смог аппроксимировать полиномами 3-го порядка, но вот как это потом объединить в общую формулу?

т.е. должна быть двумерная какая то функция, заданная из центральной точки.

кстати центральную точку(т.е. коэффициенты С1,С2), можно найти

правда ,если упростить формулу до

u= С1+x+К1*(x*x+y*y)*x;

v= C2+y+К2*(x*x+y*y)*y;

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

нашел аппроксимирующие прямые для каждой грани, методом наименьших квадратов.

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

это неверно

посмотрите на рисунок.

в общем случае

при заданных x,y и u,v, получаем обычную СЛАУ, которую и решаем через SVD, т.к. система получится переопределенной (кол-во уравнений > кол-ва неизвестных).

ну по x,y данных нету.

и формула не линейная же

u= (m11*x + m12*y + m13)/(m31*x + m32*y + m33)+С1+x+К1*(x*x+y*y)*x;

v= (m21*x + m22*y + m23)/(m31*x + m32*y + m33)+C2+y+К2*(x*x+y*y)*y;

[mrgloom 242]

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

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

[Smorodov 578]

Если мы ничего не знаем, то что же мы хотим от компа?

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

Формула линейная, относительно коэффициентов.

Мы же ищем коэффициенты?

При калибровке обычно известны искаженные и неискаженные координаты точек.

ЗЫ: хранилище на 4shared соорудил, посмотрим нужно - ли оно.

[mrgloom 242]

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

т.е.

u= (m11*x + m12*y + m13)/(m31*x + m32*y + m33);

v= (m21*x + m22*y + m23)/(m31*x + m32*y + m33);

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

перспективные искажения восстанавливаются по 4-м точкам.

и

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

u= С1+x+К1*(x*x+y*y)*x;

v= C2+y+К2*(x*x+y*y)*y;

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

про аппроксимацию дуги

аппроксимирую дугу и получаю например

v= a0+a1*u+a2*u^2+a3*u^3

и имею

u= С1+x+К1*(x*x+y*y)*x;

v= C2+y+К2*(x*x+y*y)*y;

может это как то можно решить?

[Smorodov 578]

Ошибся я, преобразование относително коэффициентов квадратичное (из-за преобразования гомографии).

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

Относительно второй части: если мы приняли допущение, что то что мы видим это прямоугольник изображения с камеры (с размерами равными разрешению камеры по вертикали и горизонтали), то мы знаем координаты x,y, и можем найти коэффициенты.

Без такого рода допущений не получится, так как не хватает информации (как не преобразовывай. Это что-то типа сохранения энергии) если линейно независимых уравнений меньше количества неизвестных, то ничего не получится.

[Smorodov 578]

Есть еще вариант:

Берем в качестве начальных приближений (x,y)=(u,v), добавляем уравнение, которое обеспечивает равномерность сетки (равенство расстояний между узлами). Уравнение, которое задает прямоугольность решетки.

Впихнуть это в итеративный цикл (выразить ур-ния в форме: x=f(x). Есть масса готовых процедур, решающих уравнения в такой форме.), теоретически должно сходиться к ближайшему локальному минимуму. Вопрос только это ли нам нужно?

[mrgloom 242]

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

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

есть еще идея фотографировать прямоугольный объект, а потом поворачивать его на 90 градусов и опять фотографировать.

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

а насчет набора точек

а если мы получим по ним, что то типа "функции поля"

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

[Smorodov 578]

По поводу поворота объекта и повторного фотографирования, правильно, это плюс еще одно уравнение, можно будет найти еще одну переменную.

Если получится потенциальная поверхность ("функция поля"), то можно произвести оптимизацию (найти минимум или максимум параметра (какого?) ), при помоци градиентного спуска (и не только). Решить задачу в такой форме просто, но вот как ее перевести в эту форму, большой вопрос.

[mrgloom 242]

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

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

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

черный это реальный

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

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

синий подобен черному.

[Smorodov 578]

Полином имеется ввиду что то вроде полинома Эрмита (сплайновая поверхность)? В принципе возможно задать сплайн на решетке, а высотой считать отклонение положения узла на искаженном изображении от его положения на неискаженном, но мы уткнулись опять в ту-же задачу, сформулированную другими словами.

По поводу пропорций, неудачный рисунок.

Пропорции поворотом воруг осей можно получить любые.

Поверните хотя бы только вокруг оси x например на 45 градусов, или на 90.

[mrgloom 242]

полинома Эрмита

не знаю. я не очень понимаю как вообще задать ,т.е. если взять типа y(x)=a0+a1*x+a2*x^2+a3*x^3 ,то так можно описать только 1 дугу и это не подходит.

т.е. не очень понятно как это задавать, возможно можно в параметрическом виде u(x,y) v(x,y).

Пропорции поворотом воруг осей можно получить любые.

Поверните хотя бы только вокруг оси x например на 45 градусов, или на 90.

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

[Smorodov 578]

красный пропорции синего (повернутого черного) вокруг оси x. Подобия нет.

ЗЫ: Так может не выдумывать велосипед, и использовать эти уравнения:

u= С1+x+К1*(x*x+y*y)*x;

v= C2+y+К2*(x*x+y*y)*y;

это тоже полиномы 3 степени.

[mrgloom 242]

красный пропорции синего (повернутого черного) вокруг оси x.

не осилил такую конструкцию)

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

т.е.

w1=k1*w2

h1=k2*h2

или имеется ввиду подобие это когда k1=k2?

[Smorodov 578]

Ух, прочитал, что написал ))

имелось ввиду

красный - пропорции синего (черного, повернутого вокруг оси x).

Да, подобие сохраняет пропорции объектов.

http://ru.wikipedia.org/wiki/%D0%9F%D0%BE%D0%B4%D0%BE%D0%B1%D0%B8%D0%B5

То есть k1=k2.

[mrgloom 242]

попробовал выправить по 4-м точкам

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

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

т.е. видимо одному и тому же передвижению точки (x,y)->(u,v) могут соответствовать разные комбинации.

[Smorodov 578]

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

[mrgloom 242]

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

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

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

[Euler 0]

есть набор точек типа такого.

А есть набор точек до и после преобразования? Если есть, то хотелось бы на него взглянуть.

[mrgloom 242]

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

u= С1+x+К1*(x*x+y*y)*x;

v= C2+y+К2*(x*x+y*y)*y;

т.е. из этого получить v(u)

данные есть только после.

Номер Коорд. X Коорд. Y

1 287 46,1

2 286 73,8

3 284 106

4 279 171

5 274 231

6 268 337

7 265 395

8 262 453

9 258 513

10 254 571

11 255 619

12 253 674

13 252 732

14 251 786

15 250 845

16 250 902

17 248 957

18 249 1010

19 249 1070

20 251 1130

21 252 1180

22 254 1240

23 256 1300

24 258 1350

25 261 1410

26 264 1470

27 268 1530

28 269 1550

29 273 1600

30 282 1700

31 285 1720

32 271 279

33 277 1650

34 318 42,6

35 350 39,8

36 392 35,9

37 422 33,2

38 456 30,7

39 488 28,2

40 526 25,7

41 556 24

42 585 22,2

43 619 20,7

44 651 19,5

45 680 18,5

46 710 17,5

47 738 16,6

48 769 16,3

49 799 15,5

50 832 15,4

51 861 15,1

52 886 14,9

53 915 14,9

54 943 14,9

55 966 15,1

56 985 15,2

57 1010 15,1

58 1030 15,3

59 1050 16

60 1070 16,1

61 1100 16,9

62 1120 17,4

63 1140 17,7

64 1160 18,3

65 1180 19,1

66 1200 19,4

67 1220 20,5

68 1240 21

69 1260 21,9

70 1280 23

71 1300 23,5

72 1320 24,5

73 1340 25,6

74 1360 27

75 1380 28,3

76 1400 28,8

77 1420 30

78 1450 31,3

79 1470 32,9

80 1490 34,4

81 1520 36,4

82 1540 38,4

83 1560 40,2

84 1580 41,3

85 1600 43,1

86 1630 45,6

87 1650 47,5

88 1680 50,5

89 1710 52,8

90 1730 54,9

91 1750 56,8

92 1760 57,7

93 1770 78,3

94 1770 103

95 1770 129

96 1770 153

97 1780 178

98 1780 208

99 1780 228

100 1790 252

101 1790 275

102 1790 299

103 1790 324

104 1790 345

105 1790 367

106 1800 388

107 1800 410

108 1800 431

109 1800 451

110 1800 468

111 1800 489

112 1800 513

113 1800 531

114 1810 553

115 1810 571

116 1810 591

117 1810 610

118 1810 630

119 1810 647

120 1810 665

121 1810 687

122 1810 704

123 1810 719

124 1810 744

125 1810 764

126 1810 794

127 1810 816

128 1810 835

129 1810 854

130 1810 873

131 1810 894

132 1810 916

133 1810 936

134 1810 956

135 1810 982

136 1810 1010

137 1810 1030

138 1810 1060

139 1810 1080

140 1810 1110

141 1810 1130

142 1810 1160

143 1810 1180

144 1810 1210

145 1810 1230

146 1810 1240

147 1810 1260

148 1810 1280

149 1800 1300

150 1800 1330

151 1800 1350

152 1800 1380

153 1800 1400

154 1800 1420

155 1800 1440

156 1800 1460

157 1790 1480

158 1790 1500

159 1790 1510

160 1790 1540

161 1790 1560

162 1790 1590

163 1790 1610

164 1780 1630

165 1780 1650

166 1780 1680

167 1780 1700

168 1750 1710

169 1720 1710

170 1690 1710

171 1650 1720

172 1630 1720

173 1600 1720

174 1580 1730

175 1550 1730

176 1520 1730

177 1490 1730

178 1460 1730

179 1440 1740

180 1420 1740

181 1400 1740

182 1380 1740

183 1350 1740

184 1330 1740

185 1310 1740

186 1280 1740

187 1260 1750

188 1240 1750

189 1220 1750

190 1200 1750

191 1180 1750

192 1150 1750

193 1130 1750

194 1110 1750

195 1080 1750

196 1060 1750

197 1040 1750

198 1010 1750

199 990 1750

200 957 1750

201 922 1750

202 879 1750

203 901 1750

204 850 1750

205 819 1750

206 783 1750

207 758 1750

208 729 1750

209 699 1750

210 674 1750

211 649 1750

212 627 1750

213 603 1750

214 579 1740

215 549 1740

216 523 1740

217 494 1740

218 467 1740

219 434 1740

220 404 1730

221 379 1730

222 350 1730

223 322 1730

224 299 1730

но можно сгенерировать впринципе синтетические данные

IplImage* warp_rect(IplImage* src, std::vector<CvPoint2D32f> vec)

{

CvPoint2D32f srcQuad[4], dstQuad[4];

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

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

srcQuad[0].y = vec[0].y;

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

srcQuad[1].y = vec[1].y;

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

srcQuad[2].y = vec[2].y;

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

srcQuad[3].y = vec[3].y;

//находим описанный прямоугольник

dstQuad[0].x = MIN(vec[0].x,vec[3].x); //Top left

dstQuad[0].y = MIN(vec[0].y,vec[1].y);

dstQuad[1].x = MAX(vec[1].x,vec[2].x); //Top right

dstQuad[1].y = MIN(vec[0].y,vec[1].y);

dstQuad[2].x = MAX(vec[1].x,vec[2].x); //Bottom left

dstQuad[2].y = MAX(vec[2].y,vec[3].y);

dstQuad[3].x = MIN(vec[0].x,vec[3].x); //Bottom right

dstQuad[3].y = MAX(vec[2].y,vec[3].y);

if (vec.size()!=4) //проверка вектора на содержание 4 углов

return NULL;

cvGetPerspectiveTransform(srcQuad,dstQuad,warp_matrix);

IplImage* dst=cvCloneImage(src);

cvWarpPerspective( src, dst, warp_matrix );

cvRectangle(dst, cvPoint(dstQuad[0].x,dstQuad[0].y), cvPoint(dstQuad[2].x,dstQuad[2].y), cvScalar(0,0,255));

cvReleaseMat(&warp_matrix);

return dst;

}

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;

}

[mrgloom 242]

что если задачу поставить так.

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

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

можем мы восстановить реальную картину?

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

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

[Euler 0]

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

Нужно решить систему относительно x и y. Указанный пример сведётся к уравнениям шестой степени, т.е. решение не получится выразить в элементарных функциях.

данные есть только после.

...

но можно сгенерировать впринципе синтетические данные

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

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

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

можем мы восстановить реальную картину?

Да, можем любое афинное преобразование вычислить. Для гомографии нужно уже 8 точек.

[mrgloom 242]

Нужно решить систему относительно x и y. Указанный пример сведётся к уравнениям шестой степени, т.е. решение не получится выразить в элементарных функциях.

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

т.е. решение не получится выразить в элементарных функциях

как это сдетектировать что нельзя?

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

я же говорю, что точек до у меня нету.

Да, можем любое афинное преобразование вычислить. Для гомографии нужно уже 8 точек.

только не точек а уравнений, 1 точка даёт 2 уравнения.

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

вот есть у нас изображение эталон A

первое изображение будет Dist[A]

второе Dist[Rot_90[A]]

я могу получить преобразование из Dist[A] в Dist[Rot_90[A]] по 4 точкам, а мне надо получить само преобразование Dist.

[mrgloom 242]

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

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

L- это матрица перехода между нашими изображениями(мы ее можем легко найти)

M- матрица которую требуется найти

Rot - матрица поворота (её мы тоже знаем)