ITERATIVE ALGORITHM FOR OFFSETS, ROTATION AND SCALE ESTIMATION FOR CONTOUR SUPERPOSITION

Full Text

Введение Совмещение контуров, у которых определено A Ex Zxy Sx B Zxy Ey Sy 0 0 0 0 0 0 соответствие между точками, является достаточ- C Sx Sy N 0 0 0 но известной задачей, которая очень часто реша- D 0 0 0 Ex Zxy Sx ется методом наименьших квадратов. Преобразование, связывающее два контура, E 0 0 0 Zxy Ey Sy содержит шесть неизвестных параметров. И для оценки этих параметров необходимы три и более точки. В общем случае преобразование координат то- F 0 0 0 Zux Zuy Sx Sy N (4) чек описывается выражением (1): Su , u1 A B C xi Zwx Zwy w1 D E F yi , (1) 1 0 0 1 1 где Sw где (xi , yi ), (ui , wi ), i = 1…N - это совмещаемые N N N контуры с количеством точек, равные N; A, B, C, D, E, F - параметры совмещения. 2 Ey yi , i 1 N Sx xi , i 1 N Sy yi , i 1 N Алгоритм оценки параметров совмещения по методу наименьших квадратов определяется вы- Zxy xi yi , i 1 Zux ui xi , i 1 N Zwx wi xi , i 1 N ражениями (2) и (3): N i i i Zuy i 1 ui yi , N Zwy i 1 N wi yi , G i 1 ( Ax By C - u )2 (2) , i (Dx Ey F - w )2 min, Su ui i 1 Sw w . i 1 i i i Однако в ряде случаев на данную задачу на- G 0, A G 0, …, B G 0, F (3) кладываются ограничения, что совмещаемые контуры отличаются: Раскрывая выражения под знаком квадрата в формуле (2) и вычисляя частные производные в формуле (3), определяем систему линейных урав- нений (4), которая может быть решена любым известным способом (методом Крамера, матрич- ным способом и т. д.): смещением, углом поворота, масштабом. В этом случае использование формул (1)-(4) приводит к появлению неоднозначности в интер- претации параметров. «Infokommunikacionnye tehnologii» 2021, Vol. 19, No. 3, pp. 283-292 Сравнивая между собой выражения (1) и (7), можно написать, что A cos( ), C h и т. д. A - sin( ), Однозначная интерпретация возможна, если A E cos( ), D -B sin( ). В этом случае выражение (1) запишется в виде выражения (8): ui A -B C xi wi B A F yi . (8) 1 0 0 1 1 А система линейных уравнений (4) запишется в виде выражения (9): а A ExEy 0 Sx Sy b1 B 0 ExEy mSx Sx b2 , (9) C Sx mSy 1 0 b3 где F Sy Sx 0 1 b4 ExEy N i 1 i i x 2 y 2 , Sx N N i 1 xi , Sy N i 1 yi , mSy -Sy, N b1 ui xi wi yi , i 1 N N b2 -ui yi wi xi , i 1 b3 ui , i 1 b4 wi . i 1 А параметры h, p, , б формуле (10): будут определяться по Рисунок 1. Контуры головок рельсов: исходные данные (а), совмещение (б) h C, p F , A2 B2 , B (10) Ниже представлены формулы, которые покаarcsin . A2 B2 зывают неоднозначность интерпретации. Пусть матрица преобразования описывается выражениями (5) и (6): При рассмотрении контуров, которые имеют одну и ту же форму, но не имеют однозначного соответствия между точками и отличаются сме- 1 0 h 0 0 вать формулы (8)-(10) нельзя. 0 1 p , M 0 2 0 , На рисунке 1 приведены пример конкуров - M1 V M1M 2 M 3 , (5) щением, углом поворота и масштабом, использо- 0 0 1 cos( ) 0 0 1 - sin( ) 0 головок рельсов, которые измеряются триангуляционным оптическим датчиком. На рисунке 1, а показаны исходные данные. На рисунке 1, б содержится результат совме- M1 sin( ) cos( ) 0 . (6) 0 0 1 щения контуров. Из рисунка 1, б видно, что точки не совпадают друг с другом. В таблице 1 и 2 по- В этом случае контуры (xi , yi ), связаны выражением (7): (ui , wi ) будут казаны координаты точек контуров, представлен- ных на рисунке 1, а. Эти данные будут использо- ui cos - sin h xi ваны для проведения численного моделирования для оценки погрешности параметров совмеще- wi sin cos p yi , (7) 1 0 0 1 1 ния. Параметры преобразования для перевода контура (xi , yi ), i 1..N в (u j , wj ), j 1..K слегде h, p - смещения, - угол, - масштаб. дующие: h 10 мм, p 10 мм, 30 , 0, 7. Таблица 1. Координаты первого контура (xi , yi ) № координата, мм № координата, мм 1 (24,860; -41,530) 21 (-15,540; -10,715) 2 (22,685; -37,945) 22 (-17,560; -10,540) 3 (20,195; -34,755) 23 (-19,580; -10,450) 4 (18,465; -30,745) 24 (-21,600; -10,330) 5 (17,025; -26,465) 25 (-23,620; -10,255) 6 (15,600; -22,230) 26 (-25,640; -10,125) 7 (13,600; -18,650) 27 (-27,660; -10,040) 8 (11,025; -15,695) 28 (-29,680; -10,010) 9 (8,700; -14,290) 29 (-31,700; -10,035) 10 (6,680; -13,495) 30 (-33,720; -10,185) 11 (4,660; -12,945) 31 (-35,740; -10,575) 12 (2,640; -12,540) 32 (-37,760; -11,110) 13 (0,620; -12,150) 33 (-39,780; -11,915) 14 (-1,400; -11,870) 34 (-42,650; -13,615) 15 (-3,420; -11,605) 35 (-45,155; -16,330) 16 (-5,440; -11,455) 36 (-47,155; -19,680) 17 (-7,460; -11,265) 37 (-48,235; -24,200) 18 (-9,480; -11,155) 38 (-48,795; -29,445) 19 (-11,500; -11,015) 39 (-49,340; -34,815) 20 (-13,520; -10,765) 40 (-49,560; -40,705) Таблица 2. Координаты первого контура (u j , wj ) № координата, мм № координата, мм 1 (37,808; -3,460) 17 (4,587; 0,249) 2 (32,980; -0,883) 18 (3,221; -0,422) 3 (30,437; 1,231) 19 (2,076; -0,997) 4 (28,373; 3,061) 20 (0,553; -1,801) 5 (26,325; 4,587) 21 (-0,515; -2,348) 6 (24,106; 5,726) 22 (-1,939; -3,076) 7 (21,665; 6,374) 23 (-3,915; -4,116) 8 (18,683; 6,137) 24 (-6,423; -5,651) 9 (15,748; 5,249) 25 (-8,769; -7,705) 10 (13,546; 4,425) 26 (-10,547; -10,141) 11 (11,910; 3,737) 27 (-11,479; -12,722) 12 (10,641; 3,127) 28 (-11,725; -15,339) 13 (9,687; 2,711) 29 (-11,333; -17,866) 14 (8,224; 1,964) 30 (-10,491; -20,305) 15 (7,144; 1,465) 31 (-9,370; -22,748) 16 (5,805; 0,878) 32 (-8,185; -25,488) Существующие работы В работе [1] предлагается метод, состоящий из двух этапов. На первом этапе проводятся грубая оценка смещения путем вычисления центра масс контуров (медианы вектора, содержащего коор- динаты точек контуров) и оценка угла поворота, вычисленная по прямой линии, которая проходит через центр масс и наиболее удаленную от него точку контура. На втором этапе осуществляется точная оценка параметров путем оптимизации (вычисление дополнительных точек на контуре для совмещения, выбор такого преобразования, при котором минимизируется евклидово рассто- яние между контурами). Алгоритм предназначен для совмещения контуров на изображениях, при этом не учитываются масштабные искажения. В работе [2] предлагается алгоритм совмеще- ния на основании выделения характерных точек контура (пересечения прямых линий), сопостав- ления их друг с другом и оценки параметров по двум наборам данных. Ограничение алгоритма заключается в наличии таких характерных точек на контуре. Если их нельзя выделить (на рисун- ке 1 представлен именно такой контур), то и нель- зя вычислить параметры совмещения. В работах [3; 4] предложено совмещать кон- туры путем предварительного разделения точек одного контура на подмножество точек. Этим подмножествам точек в соответствие ставятся линии, которые описываются аналитически. Эти Пусть задана точка x0 , y0 , относительно ко- торой произведем переход из декартовой систе- мы координат точки x, y в логарифмически по- лярную систему согласно выражениям (13)-(15): линии описывают второй контур. Параметры со- x - x 2 y - y 2 , (13) вмещения определяются путем совмещения под- множества точек и этих линий. Этот алгоритм также нельзя применить для рассматриваемой в x 0 0 rx log2 ( x ), x atan 2(x - x0 , y - y0 ), (14) (15) статье задачи. В работе [5] предлагается использовать специгде x , rx - координаты точки x, y в логариф- мически полярной системе координат. альное кодирование контуров на изображении в Примечание: atan 2 это функция, вычислявиде вектора комплексных чисел. Параметры совмещения вычисляются путем нахождения мак- симума скалярного произведения. Однако этот способ имеет ряд ограничений: контуры должны быть одного размера (нет масштабных искажеющая арктангенс отношения y - y0 / x - x0 , выраженного в радианах и с областью значений - ; , в отличие от арктангенса, у которого об- ласть значений - / 2; / 2 . Координаты ний), замкнутыми, количество точек в контурах u0 , w0 определяются согласно должно совпадать. По этой причине этот под- ход не подходит для решения рассматриваемой задачи. В работах [6-8] рассматривается задача соформулам (11) и (12), но вместо x, y в выраже- ниях будет стоять x0 , y0 . Рассмотрим выражение: u - u0 cos( )x - sin( ) y h - вмещения контуров, представляющих собой эл- - ( cos( )x - sin( ) y h) липс. Однако его нельзя модифицировать на слу- чай произвольных контуров. В работах [9-11] предложен способ совмеще- ния, который не касается напрямую контуров, а используется для совмещения изображений. Однако в нем закладывается идея раздельной оценки двух групп параметров: смещений в де- картовой системе, поворота и масштаба в лога- рифмически полярной системе координат. Зная приближенное значение одной группы параме- тров, можно оценить параметры второй группы. Если этот процесс можно повторить несколько раз, то погрешность оценки резко уменьшается с каждой новой итерацией. Этот итерационный подход можно распро- 0 0 (cos( )(x - x0 ) - sin( )( y - y0 )). Тогда: 0 0 (u - u )2 2 (cos2 ( )(x - x )2 0 sin2 ( ))( y - y )2 - - 2 cos( )(x - x0 ) sin( )( y - y0 ). Аналогично: 0 0 (w - w )2 2 (cos2 ( )(x - x )2 0 sin2 ( ))( y - y )2 - - 2 cos( )(x - x0 ) sin( )( y - y0 ). Тогда: 0 0 (u - u )2 (w - w )2 2 2 2 странить на обработку контуров. ((x - x0 ) ( y - y0 ) ). В статье представлено исследование по раз- работке итерационного алгоритма совмещения контуров. По аналогии с выражениями (13) и (14) можно записать выражения (16) и (17): Оценка угла поворота и масштаба u u - u 0 0 2 w - w 2 Угол поворота и масштаб при совмещении 2 x - x0 2 2 y - y0 контуров могут быть оценены в логарифмически полярной системе координат. Рассмотрим следующий пример. Пусть точка x - x0 0 2 y - y 2 u x . x x . (16) x, y первого контура связана с точкой u, w r log ( ) log ( ) u 2 u 2 x первого контура. Связаны они между собой формулами (11) и (12): log2 ( ) log2 ( x ) log2 ( ) rx , u cos( )x - sin( ) y h, w sin( )x - cos( ) y p. (11) (12) ru log2 ( ) rx . Рассмотрим выражение: (17) Рисунок 2. Оценка угла поворота и масштаба в логарифмически полярной системе tg w - w0 В качестве примера рассмотрим контур u (xi , yi ), i 1..N (см. таблицу 1). Преобразуем u - u0 его в контур (ui , wi ), i 1..N согласно формулам sin x - x0 cos y - y0 . и (12) с параметрами h 10 мм, p 10 мм, cos x - x0 - sin y - y0 30 , 0, 7, Умножим делитель и знаменатель дроби на Произведем переход контура (xi , yi ) из декар- товой системы в логарифмически полярную сивыражение 1 . cos x - x0 стему относительно точки x0 , y0 -20; -50 (примеч.: координаты точки были выбраны произвольно). sin x - x0 cos y - y0 Аналогичную операцию проведем для конu tg 0 cos x - x 0 cos x - x тура (ui , wi ) относительно точки u0 , w0 cos x - x0 sin y - y0 - 15,376; -25,311 (примеч.: координаты точки 0 cos x - x 0 cos x - x были рассчитаны по формулам (11) и (12), но вместо (x, y) в выражениях стояли (x0 , y0 )). tg y - y0 На рисунке 2 показаны контуры (xi , yi ) и x - x0 tg tg x (ui , wi ) в логарифмически полярной системе и 1 - tg y - y0 x - x0 1 - tg tg x параметры и . Принцип совмещения контуров tg x , Принцип совмещения контуров (xi , yi ), tg u tg x . i 1..N и (u j , wj ), j 1..K заключается в пред- Таким образом, можем записать выражение (18): варительном выборе точек x0 , y0 и u0 , w0 , относительно которых производится переход из u x . (18) декартовой системы с логарифмически поляр- ную систему (см. рисунок 3, а). где u , ru - координаты точки u, w в логарифмически полярной системе координат. Фактически точки x0 , y0 и u0 , w0 опреде- Сравнивая между собой выражения (14) и (17), а также (15) и (18), можно сделать вывод, ляют смещения h и p. Чем точнее будут рассчитаны эти точки, тем более похожи будут графики контуров в логарифмически полярной системе, что координаты u , ru связаны с координаподобно тем, что представлены на рисунке 2. тами x , rx log2 и . через аддитивные составляющие Используя графики контуров, оцениваются параметры и (см. рисунок 3, б). С исполь- Таким образом, совмещая контуры в логариф- мически полярной системе координат, произво- дим оценку угла поворота и масштаба . зованием параметров и графики совмеща- ются, поле чего передискретизируются на один и тот же шаг дискретизации (см. рисунок 3, в). а б в г Рисунок 3. Этапы совмещения Обратный переход из логарифмически поляр- N N ной системы в декартовую систему образует два контура, у которых точки имеют однозначное со- xi x i 1 , y yi i 1 , ответствие: точка 1 первого контура соответству- ет точке 1’ второго контура, точка 2 - точке 2’ 0 N 0 N K K u w (19) и т. д. (см. рисунок 3, г). Оценить параметры по двум последователь- j u j 1 , 0 K j w j 1 . 0 K ностям точек с однозначным соответствием по- зволяют формулы (9) и (10). Описанную процедуру оценки параметров можно повторить несколько раз. При этом коор- При этом полагается, что выражения (11) и будут выполняться приближенно для полу- ченных точек. динаты точки x0 , y0 , выбранные для первого Алгоритм совмещения приближения, остаются неизменными, а u0 , w0 Исходные данные: могут быть приравнены к x0 , y0 на второй и последующих итерациях. (xi , yi ), i 1..N первый контур, С каждой новой итерацией погрешность оцен- (u j , wj ), j 1..K второй контур. ки будет становиться меньше. Выходные данные: Предварительный выбор точек hˆ, pˆ ˆ , ˆ смещения вдоль координатных осей; угол поворота и масштаб. Координаты точек x0 , y0 и u0 , w0 для пер- Алгоритм состоит из следующих шагов. вой итерации совмещении будут определяться Загрузка контуров (xi , yi ) и (u j , wj ). как средние значения от координат всех точек контура согласно выражению (19): Вычисление координат x0 , y0 по формуле (19). и u0 , w0 3. Выбор количества итераций L (примеч.: в работе L = 9). Процедура оценки параметров п. 5.3 описана ниже. Инициализация матрицы преобразования в Входные данные: (xi , yi ), i 1..N - персоответствии с количеством итераций: вый контур, (u j , wj ), j 1..K - второй контур, 1 0 0 1 0 0 x0 , y0 , u0 , w0 - точки для преобразования из 1 M 0 1 0 , 0 0 1 2 M 0 1 0 , 0 0 1 декартовой в логарифмически полярную систему. Выходное значение: M - матрица преобразования размером 3x3. 1 0 0 Преобразование контура (xi , yi ) в логариф- L …, M 0 1 0 . 0 0 1 мически полярную систему относительно точки x0 , y0 : Цикл по количеству итераций, k 1. Определение промежуточной матрицы rxi log2 i 0 x - x 2 i 0 y - y 2 , (25) преобразования Mres: если k 1, то Mres M1 , xi atan 2(xi - x0 , yi - y0 ). Преобразование контура (u j , wj ) (26) в логаесли k 2, … то Mres M 2 M1 , рифмически полярную систему относительно точки u0 , w0 : если k L, то Mres ML ..M 2 M1. Вычисление положение контура (u j , wj ) ru j log 2 u j 0 w 2 j - u - w0 2 , (27) в соответствии с матрицей преобразования Mres: uj atan 2(u j - u0 , wj - w0 ). (28) u j u j Mres 1,1 (20) Дискретизация сигнала xi , rxi на равно- wj Mres 1, 2 Mres 1,3 , мерную сетку аргумента с шагом step: xn x1 , w j u j Mres 2,1 x1 step, .. xN . Формирование вектора rxn. wj Mres 2, 2 Mres 2,3 . (21) Примеч.: для обозначения векторов использовался полужирный шрифт, чтобы не вводить дополнительных индексов и не затруднять вос- Вычисление матрицы Mk по данным приятие, ранее при обозначении векторов исполь- (xi , yi ), x0 , y0 , (u j , wj ), u0 , w0 согласно прозовались индексы, как, например, «rx », « x », i i цедуре оценки параметров, описанной ниже. чтобы показать связь с исходными данными « xi », Приравнивание точки u0 , w0 к точке « yi ». x0 , y0 (примеч.: если первоначальное прибли- Дискретизация сигнала uj , ru j на равножение по п. 5.3 было выбрано верно, то в дальнейшем нет необходимости делать оценки этой мерную сетку аргумента с шагом step: un u1 , координаты, так как после первой итерации будут u1 step, .. uK . Формирование вектора run. вычислены достаточно точные значения параме- тров преобразования и на последующих итераци- Вычисление первоначального смещения вдоль оси «угол полярной системы координат»: d 0 k - u . (29) ях можно считать, что точка u0 , w0 x0 , y0 ). 1 1 Проверка на выход из цикла: k k 1, Совмещение векторов rxn и run с испольесли k L, то выход из цикла, в противном слузованием корреляционно-экстремального алгочае переход к п. 5.1. Оценка параметров Mres: hˆ, pˆ , ˆ , ˆ по матрице ритма (по максимуму коэффициента корреляции) [12]. Результатом будет являться значение смеще- ния d (в отсчетах) между векторами rxn и run. ˆ hˆ Mres 1,3 , pˆ Mres 2,3 , Mres2 1,1 Mres2 1, 2 , (22) (23) Компенсация смещения: rum = run, um = un d 0 d step, Дискретизация сигнала ( um, rum) (30) (31) на ˆ arcsin Mres 2,1 180 . (24) сетку аргумента xi. Формирование сигнала Mres2 1,1 Mres2 1, 2 uki , ruki , uki xi . Выход из программы. Преобразование в декартову систему коор- динат: а б в Рисунок 4. Погрешность оценки параметров в зависимости от СКО помехи i i u 2ruki cos uk - d 0 d step u0 , (32) где a - действительное значение параметра; ai - измеренное значение параметра в i-м опыте; u 2 ruki sin uk - d 0 d step w0 . (33) n - количество опытов. i i Оценка матрицы преобразования M соглас- но формулам (7)-(9) по точкам контуров xi , yi , (ui , wi ). Численное моделирование Для определения погрешности оценки пара- метров совмещения было проведено численное моделирование, которое заключалось в добавле- нии помехи с гауссовым распределением к аб- сциссе и ординате каждой точки контуров. При выбранном значении мощности (диспер- По результатам численного моделирования был построен график зависимости RMS от средне- квадратичного отклонения (СКО) помехи ( n ). Результаты численного моделирования приведены на рисунке 4. Разработанный алгоритм сравнивался с алго- ритмом полного перебора (проверка гипотез от- носительно измеряемых параметров с выбором той гипотезы, которая определяет минимальное евклидово расстояние между контурами). По- грешности обоих алгоритмов совпадают. Но у разработанного алгоритма скорость вычисления сии помехи, n 2 ) проводилась серия опытов, по приблизительно в 140 раз выше. которым оценивалась погрешность оценки пара- метра. Как видно из графиков (см. рисунок 4), значение погрешности практически не меняется на ин- В качестве погрешности выступало значение тервале n 10-2 ;10-1 (мм). Это связано с тем, что среднеквадратичной величины (RMS): n линейная интерполяция при дискретизации контура приводит к появлению погрешности коор- RMS i 1 ai - a 2 / n, (4) динат точек контура, которая по мощности выше, чем мощность аддитивной помехи. И только при превышении СКО помехи более 10-1 мм стано- вится заметно ее влияние на погрешность оценки параметров совмещения. Таким образом, разработанный алгоритм ха- рактеризуется низкой погрешностью оценки параметров, совпадающей с погрешностями по алгоритму полного перебора, и высоким быстро- действием. Заключение Разработанный алгоритм совмещения апроби- рован на данных, полученных в измерительных системах вагонов-путеизмерителей. Кроме железнодорожной сферы разработан- ный алгоритм может быть применен для пре- цизионного совмещения медицинских данных (гистологические измерения, трехмерная рекон- струкция), для данных, полученных при аэро- съемке (совмещение подстилающей поверхно- сти), а также в других промышленных сферах. Погрешность оценки параметров у разрабо- танного алгоритма сравнима с погрешностью по алгоритму полного перебора, при этом скорость обработки выше более чем на два порядка. Это преимущество позволяет применить раз- работанный алгоритм в измерительных системах, ориентированных на работу в режиме реального времени.

About the authors

R. R Diyazitdinov

Povolzhskiy State University of Telecommunications and Informatics

Email: rinat.diyazitdinov@gmail.com
Samara, Russian Federation

References

Supplementary files