CC BY
6
Итерационный алгоритм совмещения контуров с неравномерным шагом дискретизации
Р. Р. Диязитдинов1
1ПГУТИ - Поволжский государственный университет телекоммуникаций и информатики, 443010, Россия, г. Самара, ул. Льва Толстого, д. 23
Аннотация
В статье представлен итерационный алгоритм для совмещения контуров. В рассматриваемой задаче совмещаемые контуры имеют одинаковую форму, но шаг дискретизации является неравномерным и между точками контуров нет соответствия. По этой причине применение методов, связывающих между собой определенные точки контуров через уравнения, невозможно. В работе представлен алгоритм, который проводит раздельную оценку параметров: смещений вдоль осей координат и угла поворота для таких контуров. Идея алгоритма заключается в итерационном уточнении параметров. Оценка параметров смещения используется для расчета угла поворота, и оценка угла поворота используется для расчета смещения. Алгоритм характеризуется более высокой скоростью обработки, чем алгоритм полного перебора, и меньшей погрешностью совмещения по сравнению с алгоритмами, основанными на вычислении макропараметров контура.
Ключевые слова: совмещение, итерационный, пространственно-временной, контур.
Цитирование: Диязитдинов, Р.Р. Итерационный алгоритм совмещения контуров с неравномерным шагом дискретизации / Р.Р. Диязитдинов // Компьютерная оптика. - 2023. - Т. 47, № 1. -С. 102-111. - DOI: 10.18287/2412-6179-CO-1123.
Citation: Diyazitdinov RR. Iterative algorithm for accurate superposition of contours with nonuniform sampling step. Computer Optics 2023; 47(1): 102-111. DOI: 10.18287/2412-6179-CO-1123.
Введение
Совмещение контуров представляет собой актуальную тему исследования при решении прикладных задач [1]. Совмещение изображений по контурам применяется в медицине (совмещение гистологических изображений, трехмерная реконструкция органов и т.д.) [2], совмещение разнородных снимков подстилающей поверхности при аэросъемке [3, 4], для измерения специфических параметров и восстановления профиля рельса на железнодорожных диагностических комплексах [5, 6] и в других наукоемких прикладных сферах.
Одной из прикладных задач при оценке параметров железнодорожной колеи, в которой необходимо совмещать контуры, является измерение износа рельса. Этот параметр играет важную роль с точки зрения обеспечения безопасности движения. По результатам измерения износа принимается решение о необходимости ремонта пути, а иногда и об экстренном закрытии движения. Параметр износа определяется в результате совмещения контуров рельсов.
Для оценки износа необходимо совмещать эталонный профиль с измеренным. Несоответствие профилей в контрольных точках будет определять боковой износ [5]. Однако погрешность оценки бокового износа на существующих диагностических комплексах превышает значения, приведенные в стандарте [7] (примеч.: согласно стандарту погрешность оценки бокового износа не должна превышать 0,25 мм).
Измерения на диагностических комплексах проводятся в режиме реального времени для обеспече-
ния возможности принятия оперативных решений. В первую очередь о необходимости закрытия движения.
Для измерения износа в режиме реального времени алгоритм совмещения контуров должен обладать соответствующим быстродействием.
Таким образом, алгоритм совмещения контуров, предназначенный для программного обеспечения диагностического комплекса, должен удовлетворять требованиям по погрешности измерений и по быстродействию.
Однако достижение точности, предусмотренной стандартом, требует вычислительных затрат, не удовлетворяющих требованиям реального времени. Одним из возможных способов разрешения этого противоречия является оценка в реальном времени идентичности профилей, измеренных на небольшом участке пути (примеч.: так как форма профилей на небольшом участке пути меняется незначительно), с последующей оценкой износа одного из них по эталону. Это направление в части оценки идентичности контуров рельсов исследовано в представленной статье.
Существуют различные способы совмещения контуров. Способ совмещения во многом определяется особенностями контуров. Но какими бы особенностями не обладали контуры, общей идеей всех методов и алгоритмов является:
- сопоставление точек;
- вычисление параметров совмещения по этим
точкам.
На рис. 1 представлен наиболее простой пример совмещения, когда контуры представляют собой по-
следовательность точек, между которыми есть однозначное соответствие.
Точка первого контура, обозначенная «1», соответствует точке «1'» второго контура, точки «2» и «3» соответствуют точкам «2'» и «3'».
\у ¡
о2 X О 5'
"Г о 2' о
3 о 1' о
Рис. 1. Пример задачи совмещения контуров
Если обозначить первый контур как последовательность точек (х,, у,), а второй контур как (и,, w ), , = 1 .. N где N - количество точек, то можно записать уравнения, связывающие их друг с другом:
u¡ = x¡ • cos(ф)- y¡ • sin(ф) + h, w¡ = x¡ • sin (ф) + y ¡ • cos(ф) + p,
(1) (2)
где ф - это угол поворота,
И, р - смещения вдоль осей координат.
Оценку параметров (ф, И, р) можно найти по методу наименьших квадратов:
(ф, h, p) = argmin I V(x, • cos (ф)- y • sin (ф)+ h - и, )2
N 2 ]
• sin(ф) + y • cos(ф) + p-w) |.
(3)
В работе [8] описан алгоритм вычисления параметров.
Вводится замена тригонометрических функций на линеаризованный эквивалент: sin (ф) ~ ф, cos (ф) ~ 1- 0,5ф2.
Выражение (3) в этом случае перепишется как линеаризованная модель:
h, p ) =
I N 2
= argminI E(x •(l-0,5ф2)-yt 4 + h-ut) + (4)
Ф.А .p V ¡=i
+E(x,. 4 + y •(l-0,5ф2) + p - w)2 |.
Упрощая выражения и вычисляя частные производные по искомым параметрам, будет определена система уравнений:
-2p • Sx + 2ф^ Zu x + 2ф-Zw y + 2ф^ p • Sy + 2ф-h • Sx + +2• Zu y -2• Zw x + 2• h • Sy +ф3Еу +ф3Ex = 0; -2ф^ Sx +ф2Sy + 2 • Sw - 2 • Sy + 2p • N = 0; ф2 • Sx + 2ф-Sy + 2h • N+2Su - 2Sx = 0,
(5)
где
Ex=Ex2, Ey=Ey¡2, Sx=Ex¡, Sy=Ey<
¡=1 i=1 i=1 ¡=1 N N N
Su = E u¡, Sw = E w, Zu_x = E u¡ • x¡,
Zw _ y = E w • y ¡, Zu _ y = E и ¡ • У ¡ , Zw _ x = E
w • x¡.
В системе (5) можно выделить параметры И и р как функции параметра ф. Подставляя выражения для И и р в первое уравнение системы, определяется уравнение третьей степени от ф:
А-ф3 + В-ф2 + С-ф + Б = 0, (6)
где
А = N - Ех + N - Еу - Бх2 -Бу2; В = 0;
С = 2 - 2„_х - N - 2 - Б - Би - 2 - Бу - + 2 - _у - N;
Б = 2 - 2„_у - N - 2 - Бу - Би + 2 - Бх - - 2 - _х - N.
Решение уравнения третьей степени можно найти аналитически. Два корня в этом уравнении будут являться комплексными числами, и они не рассматриваются как решения. А третье число является оценкой угла ф. Параметры И и р вычисляются по второму и третьему уравнению системы (5).
Матрица преобразования, связывающая между собой точки (х,,у,) и (и,, V,), , = 1 .. N, имеет вид:
(
M =
cos (ф) - sin (ф) h sin (ф) cos (ф) p 0 0 1
Л
(7)
А формулы (1) и (2) можно переписать в виде:
I ut > cos (ф) - sin (ф) hл ' x¡ Л
w = sin (ф) cos (ф) p y¡
V 1 V 0 0 V 1 V V 1 J
(8)
В том случае, если угол поворота не является малой величиной (при абсолютном значении угла ф > 15° «sin (ф)» значительно отличается от «ф», а «cos (ф)» от «1 - 0,5ф2»), то замена на линеаризованную модель не является правомерной. В этом случае поворот и смещение будут неверными.
Решением данной проблемы является итерационное оценивание. На каждой итерации по параметрам (ф, h, p) рассчитывается (x¡,y¡) (примеч.: (x¡,y¡) определяются по формуле (8), в которой вместо (и, w¡) стоят (x¡, y¡)). После этого процедура оценки повторяется по данным (x', y'¡) и (и¡, w,).
Матрица преобразования будет определяться выражением:
M=MlML -1 •.. •MÍ •.. • Mb (9)
где M - матрица преобразования, вычисленная на -й итерации,
¡=1
Ь - количество итераций.
Оценка параметров по формулам (4 - 9) является частью алгоритма рассматриваемой задачи.
Возвращаясь к вопросу об особенностях контуров, на практике встречаются такие контуры, которые имеют одну и ту же форму, но при этом характеризуются:
- неравномерным шагом дискретизации;
- отсутствием соответствия между точками;
- отсутствием характерных точек или фрагментов, по которым можно совместить контуры.
На рис. 2 представлен пример контуров головок рельсов, измеренных диагностическим комплексом с помощью триангуляционного датчика. Такие ситуации возникают из-за изменения положения триангуляционного датчика относительно рельса при движении диагностического комплекса на различных скоростях.
На рис. 2а показаны исходные данные. Из них видно, что контуры сильно отличаются друг от друга как по смещению, так и по углу наклона. На рис. 2б показан результат совмещения, на котором видно, что точки принадлежат идентичным контурам, но имеют на них разное положение и в этом смысле «между ними нет соответствия».
Совместить подобные контуры в режиме реального времени с высокой точностью является сложной задачей.
ю о
-ю -20 -зо
-40 -50 -60
.у, мм
о первый .контур. • второй контур х
а) ~~ "40 -30 -20 -10 0 10 20 30 40 50 ^ <у, мм
-40 -30 -20 -10 0 10 20 30 40 50 Рис. 2. Контуры головок рельсов до (а) и после (б) совмещения
Решению задачи совмещения контуров, подобных тем, что представлены на рис. 2, посвящена данная статья.
Для совмещения этих контуров необходимо оценивать три параметра: два смещения вдоль осей ко-
ординат и угол поворота. Более общая задача, включающая масштабирование, в данной статье не рассматривается.
В табл. 1 и 2 приведены координаты точек контуров, представленных на рис. 2а. Эти данные используются для проведения численного моделирования и для описания примеров.
Табл. 1. Координаты первого контура (х, у)
№ Координата, мм № Координата, мм
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. Координаты второго контура (щ, wj)
№ Координата, мм № Координата, мм
1 (49,726; - 10,087) 13 (- 2,328; - 7,081)
2 (42,829; - 6,404) 14 (- 7,056; - 9,538)
3 (39,195; - 3,383) 15 (- 13,462; - 13,216)
4 (36,248; - 0,769) 16 (- 16,814; - 16,150)
5 (33,321; 1,411) 17 (- 19,353; - 19,631)
6 (26,664; 3,963) 18 (- 20,685; - 23,317)
7 (22,405; 3,624) 19 (- 21,036; - 27,056)
8 (18,211; 2,356) 20 (- 20,476; - 30,665)
9 (12,728; 0,196) 21 (- 19,273; - 34,151)
10 (8,428; - 1,792) 22 (- 17,672; - 37,641)
11 (4,612; - 3,599) 23 (- 15,979; - 41,554)
12 (0,994; - 5,386) 24 (- 13,918; - 45,722)
Существующие работы
В области пространственно-временной обработки для совмещения контуров исследователями рассматривались различные задачи и предлагались алгоритмы оценки параметров совмещения.
В работах [1, 2, 9, 10] в качестве предварительного «грубого» совмещения предлагалось использовать «центр масс» (точку, рассчитанную как медиана всех точек контура) для оценки смещений по координатным осям. А линия, соединяющая центр масс и наиболее удаленную от нее точку контура, используется для оценки угла поворота. Более точное совмещение определяется дополнительными оптимизационными процедурами. Описанные алгоритмы характеризуются высокой скоростью обработки данных, но погрешность оценки параметров совмещения, даже
при использовании оптимизации, не удовлетворяет требованиям, которые приведены выше.
В работах [3, 4] описана задача, в которой контуры содержат характерные точки, по которым их можно совместить. Характерные точки на рассматриваемых контурах измеряются с большой погрешностью, так как контур меняется плавно и не содержит ярко выраженных перегибов, которые определяют характерные точки. Погрешность параметров совмещения, вычисленных по этим точкам, не удовлетворяет требованиям, которые приведены во введении.
В работах [5, 6, 8] предлагается способ совмещения, где один из контуров описывается набором уравнений (прямых и окружностей), а второй контур делится на соответствующие фрагменты, которые сопоставляются этим уравнениям. Исходные данные отличаются от тех, что рассматриваются в данной статье (см. рис. 2).
В работах [11, 12] рассматриваются контуры изображений и для совмещения используется их представление в комплексном виде. Параметры совмещения определяются по максимуму скалярного произведения. Алгоритм можно использовать для контуров изображений при условии, что количество точек в совмещаемых контурах будет одинаковым, а сами контуры являются замкнутыми. Данные требования трудно выполнить для рассматриваемого случая (см. рис. 2).
В работах [13 - 16] рассматриваются контуры на плоскости, представляющие собой эллипсы. Совмещение определяется как параметры эллипсов. Естественным ограничением является узкая область применения подобных алгоритмов.
Метод полного перебора [17, 18] не рассматривает именно данную задачу, но закладывает общую идею -перебор всех возможных параметров совмещения и выбор такого набора, при котором гарантируется минимизация метрики (расстояние по Евклиду между двумя контурами). Одним из возможных вариантов реализации метода полного перебора для совмещения контуров является проверка возможных углов поворота с последующей равномерной дискретизацией контуров и оценкой смещений корреляционным способом. Однако такой подход связан с большим объемом вычислений, что приводит к снижению скорости обработки. Время обработки, даже при использовании современных процессоров, не удовлетворяет требованиям проведения измерений в режиме реального времени.
Следует отметить работы [19 - 22]. В этих работах заложена идея раздельной оценки параметров совмещения с использованием декартовой и полярной (или логарифмически-полярной) систем координат. Этот способ позволяет совместить изображения. Однако в том виде, что он описан в работах, он не применим для контуров.
Разработанный алгоритм совмещения контуров использует прием оценки параметров [19 - 22] для
выполнения, с одной стороны, требований погрешности, не превышающей заданную, а с другой стороны, требований обработки в режиме реального времени.
Описание решаемой задачи
Пусть заданы два контура (x¡, y,), i = 1 .. N, (uj, wj), j = 1 .. K, N и K - количество точек в первом и втором контурах соответственно. Оба контура описывают кривую на плоскости одной и той же формы. Обозначим кривые как f(x, y) = 0 и g(u, w) = 0. Между точками (x,, y¡) и (uj, wj) нет соответствия, как в примере, показанном на рис. 1, но аргументы функций f(x, y), g(u, w) связаны между собой:
u = x • cos (ф)-y • sin (ф) + h, (10)
w = x • sin (ф) + y • cos (ф) + p. (11)
Выражения (10) и (11) совпадают с выражениями (1) и (2) с точностью до индексов. Отсутствие соответствия между отдельными точками контуров является основной трудностью при решении задачи совмещения.
Оценка угла поворота в полярной системе координат
Угол поворота при совмещении контуров может быть оценен в полярной системе координат.
В качестве примера рассмотрим контур (x,, y¡), i = 1, .. N (см. табл. 1). Преобразуем его в контур (uí, w,), i = 1, .. N по формулам (10) и (11). Для примера зададим h = 5 мм, p = - 8 мм, ф = 15°. И рассмотрим два варианта точки (x0, y0) (примеч.: значения выбраны произвольно):
- (x0, y0) = (-20, - 50);
- (x0, y0) = (0, - 60).
Осуществим переход контура (x¡, y,) из декартовой системы в полярную систему относительно точки (x0, y0). Аналогичную операцию проведем для контура (u,, wi) относительно точки (u0, w0) (примеч.: (u0, w0) рассчитываются с использованием (x0, y0) по формулам (10) и (11), точке (x0, y0) = (- 20, -50) соответствует (u0, w0) = (-1,378; - 61,473), а точке (x0,y0) = (0, - 60) соответствует (u0, w0) = (20,529; - 65,956)).
На рис. 3 показано представление контуров (x , y ) и (u,, w,) в полярной системе координат.
Как можно видеть из рис. 3, выбор положения точки влияет на график сигнала в полярной системе координат.
Контуры, имеющие одну и ту же форму, могут иметь одинаковый график в полярной системе координат и отличаться только на смещение вдоль оси «угол полярной системы координат», если особым образом выбрать точку, относительно которой производится переход из декартовой системы в полярную.
Смещение между графиками в полярной системе координат будет численно равно углу поворота ф.
ь I радиус полярной СК, мм
а)
к радиус полярной СК, мм
угол полярной СК, граду ;
б)
Рис. 3. Представление контуров (х, у) и "м) в полярной системе координат (СК) относительно точки (х0,у0) = (-20, -50), (М, -н>0) = (-1,378;-61,473) (а) и (х0, у0) = (- 30, - 40), (М, -н>0) = (20,529; - 65,956) (б)
Оценка смещений в декартовой системе координат
Если угол поворота известен, то для оценки смещений между контурами достаточно знать соответствие между точками, и в этом случае смещение будет численно равно средней разности между координатами соответствующих точек (см. рис. 4):
h = -
Е (х< - ч) Е (у< - w')
N
Р ="
N
(12)
где (ui', w/) - координаты второго контура после компенсации угла поворота.
пу
(«Л:
: *
у
¿4 Д дал ).т
ч \
■ ч\ : \\
1 \(*зЬз)
Рис. 4. Оценка смещений вдоль осей координат
Однако, как утверждалось в самом начале статьи, точки контуров не имеют однозначного соответствия. Для исправления этого недостатка будет использован прием дискретизации сигнала (примеч.: в работе используется дискретизация сигнала с линейной интерполяцией).
На рис. 5 условно показаны этапы обработки, которые позволяют провести оценку смещений вдоль осей координат:
- исходные контуры (рис. 5а);
- представление контуров в полярной системе,
оценка угла поворота (рис. 5б);
- совмещение и дискретизация второго контура на
шаг первого контура (рис. 5в);
- представление контуров в декартовой системе с однозначным соответствием между точками (1 - 1', 2 - 2', .. , N- Щ (рис. 5г).
а)
о первый: контур второй контур
б)
радиус полярной СК, мм
полярной СК, градус
в)
радиус полярной СК, мм
угол полярной СК, градус
А У
первый'контур у второй вонтур,: дискретизированный по отсчетам первого ■ контура : :
г)
Рис. 5. Оценка смещений
Алгоритм оценки параметров совмещения по методу полного перебора
Вывод, который можно сделать из параграфа об оценке угла поворота, заключается в следующем: если точки, относительно которых происходит преобразование в полярную систему, связаны между собой формулами (10) и (11), то по графикам в полярной системе координат можно оценить угол поворота ф.
Для алгоритма полного перебора точка (х0, у0) для контура (хиу), i = 1, .. N выбирается произвольным образом (примеч.: в работе использовалось среднее значение координат от всех точек контура).
А точка ^0, w0) для контура Wj), ] = 1, .. К, определяется переборным способом. В работе ис-
пользовались значения из области [ит - иТ, ит + иТ], \у>т - ШТ, »т + Шт] с шагом (ит, »т) - это среднее значение координат от всех точек контура, ит = 3 мм, Шт = 3 мм, = 0,1 мм. Примеч.: полагается, что в указанной области есть такая точка (и0, »0), которая будет связана с (х0,у0) формулами (10) и (11); обоснование значений параметров Пт, Шт представлено в параграфе «Обоснование выбора значений параметров».
На рис. 6 показано схематично точки и области. АУ
о первый: контур второй контур
Рис. 6. Выбор точки (х0, у0) и (и0, »0) для алгоритма полного перебора
Для совмещения графиков в полярной системе используется корреляционный способ [23]. Для каждого положения (и0, »0) будет определен угол ф.
В ходе совмещения будут сформированы два контура, у которых определено однозначное соответствие между точками (см. рис. 5). Зная однозначное соответствие, проводится оценка параметров совмещения < И, р, ф > по (4 - 9).
Для того чтобы выбрать наилучшую оценку, используется метрика, которая представляет собой евклидово расстояние. Минимальное значение метрики будет соответствовать наилучшему совпадению контуров и, соответственно, минимальной погрешности оценки параметров.
Описанный алгоритм достаточно прост в реализации, однако обладает существенным недостатком, связанным с большим временем обработки. Например, при выбранных значениях иТ, Жт, количество проверок будет равно 612 = 3721 (примеч.: параметры ит = 3,0 мм и = 0,1 мм, определяют координаты: ит -3,0; ит - 2,9 .. ит + 3,0, то есть количество точек вдоль оси абсцисс будет 61, аналогичным образом рассчитывается количество точек вдоль оси ординат).
Чтобы уменьшить время обработки, был разработан алгоритм совмещения, в котором не требуется такого огромного количества проверок.
Разработанный алгоритм оценки параметров совмещения
Идея разработанного алгоритма заключается в итерационном совмещении контуров. Если для точек, относительно которых происходит преобразование в полярную систему координат, на первой итерации выражения (10) и (11) выполняются приближённо (примеч.: в ходе исследований было определено, что
отклонение может достигать 1,0 мм), то на каждой последующей итерации погрешность оценки параметров совмещения будет становиться меньше.
В разработанном алгоритме точка (x0,y0) рассчитывается, как и в алгоритме полного перебора, и точка (u0, w0) также лежит в области [um- UT, um + UT], [wm - WT, wm + WT]. Однако шаг st можно выбрать много больше. Даже при относительно большом шаге st наблюдается уменьшение погрешности с увеличением числа итераций.
При проведении сравнительного анализа значение st = 1,0 мм обеспечивало совпадение результатов оценки параметров разработанного алгоритма с результатами алгоритма полного перебора.
При выбранных значениях UT, WT, st количество проверок будет равно 72 = 49 (примеч.: параметры UT = 3,0 мм и st = 1,0 мм определяют координаты: um -3,0; um - 2,0 .. um + 3,0, то есть количество точек вдоль оси абсцисс будет 7, аналогичным образом рассчитывается количество точек вдоль оси ординат, обоснование значений параметров Ut, Wt представлено в параграфе «Обоснование выбора значений параметров»).
Таким образом, число проверок уменьшилось с 612 до 72 (с 3721 до 49).
Итерационная процедура совмещения представлена ниже.
Исходными данными являются:
1. (xi,y), i = 1 .. N - первый контур;
2. (u,, Wj), j = 1 .. K - второй контур;
3. (x0,y0), (u0, w0) - точки, относительно которых проводится преобразование из декартовой системы в полярную систему.
Выходными данными являются:
1. Mres - матрица [3*3] с параметрами преобразования;
2. (h,p, ф) - параметры совмещения;
3. metric - метрика (евклидово расстояние между контурами).
Алгоритм оценки параметров совмещения имеет следующий вид:
1. Загрузка данных: (xi,yi), (u,, w,), (x0,У0), (u0, W0).
2. Выбор количества итераций L.
3. Инициализация матрицы преобразования в соответствии с количеством итераций:
M1 =
(100 ^ 0 1 0 0 0 1
M2 = Mi, ...ML = Mi.
4. Цикл по количеству итераций, к = 1.
4.1. Определение промежуточной матрицы преобразования Мге5: если к = 1, то Мге5 = М1, если к = 2, то Мтея=М2М1,
если к=Ь, то Мтея=Мь .. М2М1.
4.2. Вычисление положения контура (и/, м/') в соответствии с матрицей преобразования Мге5:
и = и -Мтез(1,1) + » -Мге,5(1,2)+ Мтез(1,3), (13)
м>) = и/ -Ме(2,1) + -Ме(2,2)+ Ме(2,3). (14)
4.3. Вычисление матрицы Мк, параметров (Ик,рк, фк), метрики те^вк по данным (х,,у,), (х0,у0), (и/', »/'), (и0, »0) согласно процедуре оценки параметров, описанной ниже.
4.4. Приравнивание точки (и0, »0) к точке (х0, у0) (примеч.: если первоначальное приближение по п. 4.3 было выбрано верно, то в дальнейшем нет необходимости делать оценку этой координаты, так как после первой итерации будут вычислены достаточно точные значения параметров преобразования, и на последующих итерациях можно считать, что точка (и0, »0) = (х0, у0)).
4.5. Проверка на выход из цикла: если (|Ик - Ик-1| < Мг_И и рк -рк-11 < Лг_р и |фк - фк-11 < Лг_ф) или к+1 > Ь, то выход из цикла, в противном случае к = к +1 и переход к п. 4.1 (примеч.: №г_И = 0,1 мм, 1Ит_р = 0,1 мм, IИг_ф = 0,1°, И0 = 0,0 мм, р0 = 0,0 мм, ф0 = 0,0°).
5. В результате работы алгоритма будет определена матрица Мге^\ из которой можно оценить пара-
метры (ф, h, p):
ф = arcsin (Mres (2,1)); (15)
h = Mres (1,3); (16)
p = Mres (2,3). (17)
metric = metrick (18)
6. Выход из программы.
Процедура оценки параметров п. 4.3 описана ниже. Входные данные:
1. (xi,yi), i = 1.. N - первый контур,
2. (uj, Wj'), j = 1.. K - второй контур,
3. (x0, У0), (u0, W0) - точки для преобразования из декартовой в полярную систему.
Выходные данные:
1. M- матрица [3 *3] с параметрами преобразования;
2. (h, p, ф) - параметры совмещения;
3. metric - метрика (евклидово расстояние между контурами).
1. Преобразование контура (x, yi) в полярную систему относительно точки (x0, y0):
px, =^/(x, - x„ )2 + (y - y0 )2, (19)
фxI■ = atan2 (x,- - x0, yi -y>). (20)
Ко всем отрицательным значениям фxi добавляется значение 2п, чтобы параметр принимал только положительные значения.
2. Преобразование контура (u/, Wj') в полярную систему относительно точки (u0, w0):
puj =^(u'j - u0 )2 +(w'j - W0 )2, (21)
ф^- = atan2(u/ - u0, wj' - w0). (22)
Ко всем отрицательным значениям ф^ добавляется значение 2п, чтобы параметр принимал только положительные значения.
3. Дискретизация сигнала ^x¡, px,) на равномерную сетку аргумента с шагом step: фxnii = фx1, фГ1 + step, .., фXNN, ii = 1, 2 .. NN, NN = (9xnn -фx1) /step (примеч.: step = 0,3°, обоснование значения параметра step представлено в параграфе «Обоснование выбора значений параметров»). Формирование вектора pxn,¡.
4. Дискретизация сигнала (щ-, puj на равномерную сетку аргумента с шагом step: фunj- = фu1, фщ + step, .., фuк, jj = 1, 2 .. KK, KK = (фuк -ф^) /step. Формирование вектора pu%.
5. Вычисление первоначального смещения вдоль оси «угол полярной системы координат»:
d0 = фXl - фul. (23)
6. Совмещение векторов pxn,, и pu% с использованием корреляционно-экстремального алгоритма (по максимуму коэффициента корреляции) [23]. Результатом будет являться значение смещения d (в отсчетах) между векторами pxn,, и pu%.
7. Компенсация смещения:
pum=puj (24)
фumjj = фunjj + d0 + d^step. (25)
8. Дискретизация сигнала ^umj-, pumjj) на сетку аргумента фXi. Формирование сигнала (uk, puki), фuki = фxi.
9. Преобразование в декартовую систему координат:
u¡ = puk¡ • cos (uk, - (d0 + d • step)) + u0, (26)
wj = puk¡ • sin (uk, - (d0 + d • step)) + w0. (27)
10. Оценка матрицы преобразования M согласно (4 - 9) по точкам контуров (x¡, y,), (u ,'', w ,■").
ф = arcsin (M(2,1)); h = M(1,3); p = M(2,3);
p-2-2
metric = JE(x¡ - u-")2 + ( - w")2.
Как и в случае с алгоритмом полного перебора, минимальное значение метрики используется для вынесения решения о параметрах (ф, h, p).
Обоснование выбора значений параметров
В итерационной процедуре оценки введен параметр L (количество итераций). Чтобы выбрать параметр L, были построены зависимости среднего значения метрики от числа итераций при различном среднеквадратичном отклонении аддитивной помехи (см. рис. 7). Описание моделирования подробно приведено в следующем параграфе.
Как можно видеть из рис. 7, при L > 5 значение метрики практически не меняется. По этой причине
можно ограничиться количеством итераций Ь = 5, так как с дальнейшим ростом количества итераций будет увеличиваться время обработки, а погрешность оценки параметров уменьшаться не будет.
' метрика, мм
1(Г 1<Г
Рис. 7. Выбор числа итераций L
Аналогичным образом выбирались значение шага дискретизации угла Ut. Wt (см. рис. 9) и диапазоны поиска step (см. рис. 8).
10J 102 10' 10°
Рис. 8. Выбор параметров диапазона поиска
метрика, мм
103 10"2 10"' 10° Рис. 9. Выбор параметра шага дискретизации угла
Численное моделирование
Для определения погрешности оценки параметров совмещения было проведено численное моделирование, которое заключалось в добавлении помехи с Гауссовым распределением к абсциссе и ординате каждой точке контуров.
При выбранном значении мощности (дисперсии помехи, си2) проводилась серия моделирований (опытов), по которым оценивалась погрешность оценки параметра.
В качестве погрешности выступало значение среднеквадратичной величины (RMS):
RMS =
" a )2
(28)
где a - действительное значение параметра, at - измеренное значение параметра в i-м опыте; n - количество опытов.
По результатам численного моделирования был построен график зависимости RMS от среднеквадратичного отклонения (СКО) помехи (cn).
На рис. 10 приведен результат совмещения профилей с помехой cn = 1,0 мм. Несмотря на большие искажения, результат совмещения выглядит корректно.
10
-10 -20 -30 -40 -50
у, мм
перцыи контур второй кшггур..
: i X,
ММ
-30 -20 -10 0 10 20 30 40 50 60
Рис. 10. Пример совмещения контуров при Оп = 1,0 мм
На рис. 11 показаны результаты оценки параметров от количества итераций Ь.
Как можно видеть из рис. 12, при Ь > 5 значение погрешности практически не меняется, что еще раз свидетельствует в пользу выбора количества итераций, равного пяти.
Результаты погрешностей совмещения (рис. 10) сравнивались с методами:
- полного перебора [17];
- макропараметров [9];
- макропараметров с оптимизацией [2].
Как можно видеть из рис. 12, погрешность предложенного алгоритма и алгоритма по методу полного перебора совпадают, а погрешность алгоритмов с использованием макропараметров значительно выше.
1=1
n
к VRMS <p, рад
б)
Рис. 11. Значение погрешности оценки параметров от количества итераций Ь
рад
..... полный перебор [17]
■- разработанный алгоритм
--алгоритм по макропараметрам [9]
--алгоритм по макропараметрам
с оптимизацией [2]
а)
12 10 8 6 4 2 0
t iRMS h(p), мм
полный перебор [17] разработанный алгоритм^ - - : алгоритм по макропараметрам [9] алгоритм по макропараметрам с оптимизацией [2] :
б)
Рис. 12. Результаты численного моделирования
Сравнительный анализ показал, что скорость обработки данных в алгоритмах с макропараметрами [2,
9] приблизительно в 100 - 200 раз быстрее, чем в предложенном алгоритме. А алгоритм полного перебора [17], напротив, приблизительно в 50 раз медленнее предложенного.
Несмотря на то, что скорость обработки предложенного алгоритма много меньше, чем у алгоритмов с макропараметрами [2, 9], на современных процессорах алгоритм выполняет требование работы в режиме реального времени (примеч: при работе триангуляционного датчика на частоте 1000 Гц (1 мс) время обработки данных составляет ~0,5 мс на процессоре Intel Core i7-3770). Алгоритм полного перебора [17] работает значительно медленнее, чем разработанный алгоритм, и даже на современных процессорах обработка не может быть выполнена в реальном масштабе времени.
С точки зрения погрешности совмещения можно ориентироваться на измерения при СКО помехи си = 0,1 мм (примеч. : современные триангуляционные датчики характеризуются помехами менее 0,1 мм) [24]. Алгоритмы с макропараметрами [2, 9] при си=0,1 мм имеют погрешности равными 10,0 мм и 1,2 мм по смещению и 0,048 рад и 0,012 рад по углу соответственно. А алгоритмы полного перебора [17] и разработанный алгоритм - 0,1 мм по смещению и 0,003 рад по углу.
Предложенный алгоритм был реализован в виде программного обеспечения для измерения бокового износа на диагностическом комплексе. Внедрение предложенного алгоритма повысило точность измерения этого параметра. До внедрения погрешность оценки бокового износа составляла owear = 0,38 мм, после внедрения cwear = 0,22 мм.
Заключение
В работе представлен алгоритм совмещения контуров, которые характеризуются неравномерностью шага дискретизации и отсутствием соответствия между точками контуров.
Разработанный алгоритм представляет собой решение, обеспечивающее скорость обработки данных в режиме реального времени и позволяющее уменьшить погрешность оценки бокового износа при его использовании в составе диагностического комплекса на железной дороге. Также алгоритм может быть использован при разработке систем анализа и для измерительных систем авиационных предприятий, а также при разработке медицинского оборудования.
References
[1] Soifer VA, ed. Methods for computer image processing [In Russian]. Moskow: "Fizmatlit" Publisher; 2003.
[2] Sungatullina DI, Krilov AV. Fast algorithm for image contour superposition linking isotropic affine transformation [In Russian]. GRAFIKON 2014; 92-95.
[3] Efimov AI. Developing and researching image superposition algorithm for video sensors with a virtual terrain model [In Russian]. Ryazan: "Ryazan State Radio Engineering University" Publisher; 2016.
[4] Efimov AI, Novikov AI. An algorithm for multistage projective transformation adjustment for image superimposition [In Russian]. Computer Optics 2016; 40(2): 258265. DOI: 10.18287/2412-6179-2016-40-2-258-265.
[5] Vasin NN, Diyazitdinov RR. Processing of triangulation scanner data for measurements of rail profiles [In Russian]. Computer Optics 2018; 42(6): 1054-1061. DOI: 10.18287/2412-6179-2018-42-6-1054-1061.
[6] Diyazitdinov RR. Video signal recovery in measuring systems with optical triangulation sensors [In Russian]. Infocommunication Technologies 2019; 17(3): 324-331. DOI: 10.18469/ikt.2019.17.3.09.
[7] DIN EN 13674-1-2011. Railway applications - Track -Rail - Part 1 : Vignole railway rails 46 kg/m and above; German version EN 13674-1:2011.
[8] Diyazitdinov R. Iterative algorithm of optical triangulation sensors signals superposition for measuring solid deformation. CEUR Workshop Proceedings 2020; 2665: 93-99.
[9] Wang W, Jiang Y, Xiong B, Zhao L. Contour matching using the affine-invariant support point set. IET Computer Vision 2014; 8: 35-44. DOI: 10.1049/iet-cvi.2013.0031.
[10] Efimov AI, Novikov AI. Software and algorithmic system for image superposition in aircraft vision systems [In Russian]. The III Int Conf on Information Technology and Nanotechnology (ITNT-2017) 2017: 400-409.
[11] Furman YaA, Kreversky AV, Peredreyev AK, Rozentsov AA, Hafizov RG, Yegoshina IL, Leukhin AN. Contour analysis and its image and signal processing application [In Russian]. Moscow: "Fizmatlit" Publisher; 2003. ISBN: 5-9221-0374-1.
[12] Makarov MA. Contour analisys in the problems of description and classification of objects [In Russian]. Modern Problems of Science and Education. Surgery 2014; 3: 38-38.
[13] Ellis T, Abbood A, Brillault B. Ellipse detection and matching with uncertainty. Image Vis Comput 1992; 10(5): 271-276. DOI: 10.1016/0262-8856(92)90041-Z.
[14] Fitzgibbon A, Fisher R. A bayer's guide to conic fitting. Proc 6th British Conf on Machine Vision 1995; 2: 513522.
[15] Gander W, Golub GH, Strebel R. Least-square fitting of circles and ellipses. BIT Numer Math 1994; 34(4): 558578. DOI: 10.1007/BF01934268.
[16] Bookstein FL. Fitting conic sections to scattered data. Comput Graph Image Process 1979; 9(1): 56-71. DOI: 10.1016/0146-664X(79)90082-0.
[17] Linnik YuV. Least Squares method and base of mathematical and statistical theory of data processing [In Russian]. Moscow: "Fizmatgiz" Publisher; 1958.
[18] Baklitckiy VK. The method of signal filtering in correlation-extreme navigation systems [In Russian]. Tver: "Knigniy klub" Publisher; 2009.
[19] Myasnikov EV. Determination of parameters of geometric transformation to combine portrait images. Computer Optics 2007; 31(3): 77-82.
[20] Reddy B, Chatterji B. An FFT-based technique for translation, rotation, and scale-invariant image registration. IEEE Trans Image Process 1996; 5(8): 1266-1271. DOI: 10.1109/83.506761.
[21] Alba A, Aguilar-Ponce R, Vigueras-Gomez J, Arce-Santana E. Phase correlation based image alignment with subpixel accuracy. In Book: Batyrshin I, Mendoza MG, eds. Advances in artificial intelligence. 11th Mexican Int Conf on Artificial Intelligence (MICAI 2012), Part 1 2012: 171-182. DOI: 10.1007/978-3-642-37807-2_15.
[22] Evangelidis G, Psarakis E. Parametric image alignment using enhanced correlation coefficient maximization. IEEE Trans Pattern Anal Mach Intell 2008; 30(10): 18581865. DOI: 10.1109/TPAMI.2008.113.
[23] Kuzmin SV. Scale-invariant delay estimation between two one-dimensional digital signals [In Russian]. Info-communication Technologies 2011; 9(2): 7-10.
[24] Catalog optoNCDT laser sensors (Laser displacement sensors - triangulation). Source: (https://www.micro-epsilon.ru/download/products/cat--optoNCDT--en.pdf).
Сведения об авторе
Диязитдинов Ринат Радмирович, 1987 года рождения, в 2009 году окончил Поволжский государственный университет телекоммуникаций и информатики (ПГУТИ) по специальности 210709 «Многоканальные телекоммуникационные системы», работает доцентом в ПГУТИ. Область научных интересов: компьютерная оптика, обработка изображений, триангуляционные системы измерения. E-mail: rinat.diyazitdinov@smail.com .
ГРНТИ: 28.23.15
Поступила в редакцию 15 марта 2022 г. Окончательный вариант - 26 июля 2022 г.
Iterative algorithm for accurate superposition of contours with non-uniform sampling step
R.R. Diyazitdinov1
1PSUTY - Povolzhskiy State University of Telecommunications and Informatics, 443010, Samara, Russia, Tolstoy str., 23
Abstract
In this article, we describe an iterative algorithm for accurate superposition of contours with non-uniform sampling step. The processing contours are characterized by the same shape, but the sampling step is non-uniform, with no matching between points of the superposed contours. This makes impossible the use of methods for estimating superposition parameters by matching points. The algorithm proposed herein allows estimating the offsets and rotation angle separately. The idea of the algorithm is to perform the iterative correction of parameters. An estimate of the offsets is used to estimate the rotation angle and, vice versa, an estimate of the rotation angle is used to estimate the offsets. The proposed algorithm is characterized by a higher speed of processing than a brute force algorithm and a lower estimation error than algorithms that analyze contour macro-parameters.
Keywords: superposition, iterative, space-time, contour, accuracy.
Citation: Diyazitdinov RR. Iterative algorithm for accurate superposition of contours with nonuniform sampling step. Computer Optics 2023; 47(1): 102-111. DOI: 10.18287/2412-6179-CO-1123.
Author's information
Rinat Radmirovich Diyazitdinov (b. 1987) graduated from Povolzhskiy State University of Telecommunications and Informatics in 2009 (PSUTI), Multychannel Telecommunication System sub-department. He works as the assistant professor. His research interests are currently focused on computer optics, image processing, triangulation measurements system. E-mail: rinat.diyazitdinov@smail.com .
Received March 15, 2022. The final version - July 26, 2022.
