Точные решения вариационной задачи алгоритма ICP в классе аффинных преобразований Текст научной статьи по специальности «Математика»

Челябинский физико-математический журнал. 2017. Т. 2, вып. 3. С. 282-294.

УДК 517.972.9

ТОЧНЫЕ РЕШЕНИЯ ВАРИАЦИОННОЙ ЗАДАЧИ АЛГОРИТМА ICP В КЛАССЕ АФФИННЫХ ПРЕОБРАЗОВАНИЙ

А. Ю. Маковецкий", С. М. Воронин6, Д. В. Тихонькихс, М. Н. Алексеев^

Челябинский государственный университет, Челябинск, Россия "artemmac@csu.ru, bvoron@csu.ru, csparsemind@gmail.com, dalexeev@csu.ru

Наиболее используемым методом регистрации облаков точек в трёхмерном пространстве является «итерационный алгоритм ближайших точек» (iterative closest points, ICP). Цель работы алгоритма — вычисление оптимального относительно заданной метрики геометрического преобразования, совмещающего два данных облака. Важным этапом алгоритма ICP является решение задачи минимизации функционала, соответствующего данной метрике для данного класса геометрических преобразований. В работе представлен метод решения вариационной задачи алгоритма ICP для метрики «соответствие типа точка-точка» (point-to-point) в классе аффинных преобразований. С помощью компьютерного моделирования демонстрируется корректность работы алгоритма.

Ключевые слова: трёхмерная реконструкция, регистрация облаков точек, локализация, аффинные преобразования.

Введение

Алгоритм ICP (Iterative Closest Points) является основным среди методов «выравнивания» трёхмерных моделей, основанных на использовании исключительно геометрических характеристик. Под выравниванием понимается определение геометрического преобразования, наилучшим относительно нормы L2 образом связывающего два данных набора (облака) точек в R3. Алгоритм широко используется для регистрации данных, получаемых с помощью SD-сканеров. Первоначально описанный в работах Чена и Медиони [1] и Бесла и Маккея [2] алгоритм ICP состоит из следующих итеративно применяемых шагов:

1. Выбор некоторого подмножества точек в обоих облаках.

2. Определение соответствия между точками выбранных подмножеств.

3. Сопоставление весовых коэффициентов полученным парам.

4. Отбрасывание некоторых пар на основе различных критериев.

5. Выбор метрики ошибки для пар точек.

6. Минимизация метрики ошибки (вариационная подзадача алгоритма ICP).

Существуют два основных подхода к выбору метрики ошибки для пар точек. В рамках первого подхода (point-to-point) [2] используется расстояние между элементами пары в R3, в рамках второго подхода (point-to-plane) [1] учитывается расстояние между точкой первого облака и касательной плоскостью к соответствующей точке второго облака.

Ключевым элементом [3] алгоритма ICP является поиск ортогонального или аффинного преобразования, наилучшим образом в смысле квадратичной метрики

совмещающего два облака точек с заданным соответствием между точками (вариационная подзадача алгоритма ICP). Для ортогональных преобразований решение этой задачи в замкнутой форме получено Б. Хорном в работах [4] и [5]. В первой из них решение основано на использовании кватернионов, во второй — на использовании ортогональных матриц. Решения являются линейными по времени по отношению к числу пар точек.

Оригинальный алгоритм ICP широко применяется для регистрации жёстких объектов, но плохо работает в случае регистрации нежёстких объектов. В работе [6] предложено расширение алгоритма ICP при использовании масштабирования помимо поворота и параллельного переноса. Обобщение алгоритма на случай произвольного аффинного преобразования сделано в работах [7; 8].

Описанные выше подходы к решению вариационной подзадачи алгоритма ICP относятся к point-to-point методам. Для метода point-to-plane аналитическое решение неизвестно. Для решения задачи используется либо итерационный метод Левенберга — Марквардта, либо линеаризация проблемы [9].

В данной работе описываются точные решения вариационной подзадачи алгоритма ICP для случая аффинных преобразований с использованием метода point-to-point. В отличие от работы [7] представляемый подход позволяет получать точные решения для всех возможных вырожденных случаев. Частично метод был описан в работе [10]. Компьютерное моделирование содержит результаты вычислительных экспериментов, основанных на представляемом методе решения задачи минимизации.

1. Формулировка задачи минимизации

Обозначим через X = ..., жп} первое облако точек, через У = ..., уп} — второе облако точек в К3. Предположим, что известно такое соответствие между облаками X и У, что каждой точке ж первого облака сопоставлена точка у второго облака. Во многих работах, посвящённых алгоритму 1СР, рассматривается ортогональное геометрическое преобразование, связывающее X и У: Кж + где К — матрица поворота, £ — вектор параллельного переноса, г = 1, 2,..., п. В рамках алгоритма Б-ГСР рассматривается геометрическое преобразование следующего вида: КБж + где Б — матрица масштабирования.

Аффинное преобразование общего вида в К3 задаётся функцией от 12 переменных. Рассмотрим вариационную подзадачу алгоритма 1СР для произвольного аффинного преобразования. Обозначим через (жн, ж2г, и (ун, у2г, у3г) координаты точек ж и у в данной системе координат. Обозначим через 3(А,£) следующий функционал:

J(A,t) = ^ \\Axi + t - y

i=0

A

aii ai2 ai3 \ ti

a2i a22 a23 1 , t = | t2

a3i a32 a33 / t3

xi

i = 1,2, ...,n. Вариационная задача алгоритма ICP заключается в поиске arg min J(A,t), где функционал J(A,t) имеет вид

J(A,t) = E(«iixii + ai2X2i + ais^si + ti - yii)2+

i=1

+ (a2iXii + a22X2i + a23X3i + t2 - y2i)2 + (a3iXii + аз2^ + азз^ + t3 - y3i)2.

2

y

2. Решение задачи минимизации

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

1 П 1 П 1 п

Хц = Хц ^ ^ х1.7 , х2г = Х2г ^ ^ х2., х3г = х3г ^ ^ х3.

7=1 7=1 7=1

а к точкам второго облака — параллельный перенос

1 п 1 п 1 п

у!г = Ун - П ^ Уи, у2г = У2г - ^ ^ , у3г = У3г - ^ ^

7=1 7=1 .7 = 1

В новых координатах получаем

п

J(А) = Е^цх'^ + Я12х2г + а13х3г - у!г)2 + ^х'^ + а22х2г + а23х3 - у2г)2+

г=1

+(а31х1г + а32Х;2г + а33х3 - у3г)2.

Таким образом, в новых координатах {х1, х'2, х'3| и {у^,у2,у3} функционал 3 не зависит от компонент параллельного переноса. Далее штрихованные координаты записываются без штрихов.

Рассмотрим различные вырожденные случаи взаимного расположения точек первого облака X.

2.1. Пусть все точки лежат на прямой Ох, т. е. х2г = х3г = 0 при всех г = 1, 2,... , п. Функционал 3(А) принимает вид

3(А) = ^(аИХ1г - у1г)2 + (а21Х1г - у2г)2 + (а31ХН - у3г)2.

г=1

Вычислим компоненты градиента V/(А): д3 (А)

Е у1гХ1г

2 У^(ацХ1г - у1г)Х1г = 0, ап = г=1-, У^х1г = 0

да-м — 2 —

11 г=1 ^Х?г г=1

г=1

Аналогично

Е у2гХ1г Е у3гХ1г

_ г=1 _ г=1 а21 = -П-, а31 = -П-.

Е Х1г Е Х1г

г=1 г=1

Остальные элементы матрицы а., г = 1, 2, 3, ] = 2, 3, принимают любые значения. Так как рассматриваемые координаты являются штрихованными, заменим их на исходные:

п

Х1к = х1к + ПЕ Х1., к = 1, 2, 3, п . =1

п

уи = У'И + ПЕ у., I = 1, 2, 3. . =1

Вычислим компоненты вектора используя найденную матрицу А:

1п

¿й = - У^(уйг - (а&1Х1г + а^Х2г + а^Х3г)) = 0, к = 1, 2, 3. п

г=1

Замечание 1. Значения элементов матрицы а&2, а&3, к = 1, 2, 3, не влияют на результат суммирования, так как ж2^ = ж3^ = 0, г = 1, 2,... , п.

2.2. Пусть все точки первого облака X лежат на прямой

Ж1 - Рх _ Ж2 - Ру _ Ж3 - Р*

VI

^2

^3

где V = (^1,^2,^3) — направляющий единичный вектор прямой. Рассмотрим орто-нормальную систему координат {ж1, ж'2, ж'3}:

ж1 \ 1 1 ( щ ^2 ^3

ж2 1 = 1

ж3 У I 1 V Ш 0 7)2

1

Х1 - Рх Ж2 - Ру Ж3 - Р*

Направляющий вектор оси рж^ этой системы является вектором V, центр системы — точка р = (рх,ру,рг). В системе координат {ж1,ж2, ж'3} задача сводится к п. 2.1. Запишем последнее равенство в виде

ж' = М—1 (ж - р). (1)

Используя вычисленные матрицу К и вектор получаем итоговое преобразование

КМ-1ж - Яр + (2)

2.3. Пусть все точки первого облака X лежат в плоскости Ожу. После параллельного переноса облака точек функционал 3(А) принимает вид

п

3(А) = ^(ацжн + Й12ж2г - У1г)2 + («21 жи + Й22ж2г - У2г)2 + («31 жи + Й32ж2г - У3г)2.

г=1

Имеем

(А) 5а"

п п п

а11 ^^ ж1г = ^^ УижИ - а12 ^^ ж2гж1г, а11 г=1 г=1 г=1

2 ^(аижн + Й12ж2г - Ун)жн = 0,

пп

Е У^жн - а12 Е ж2гж1г

г=1

г=1

г=1

ж2

Е жи = 0.

1г

г=1

г=1

Функционал 3(А) принимает вид

3 (А) = Е

г=1

п

ж^ Е ж2^- ж^

«12

V

¿=1

ж2г —

V

ж2

¿=1

У'

/

Е У1^- жу-

Уи - Ч-ж1»

V

ж2

¿=1

У

+...

Обозначим через ^и и 72» следующие выражения:

71г = ж2г

Е ж2Уж1У ¿=1

ж2

-ж1г, 72г = У1г

Е У1У ж1У ¿=1

-ж1г,

¿=1

1У

ж2

¿=1

У

2

где ЕП=1 Ху = 0. Тогда 3(А) = Еп=1(а1271г - 72г)2 + ...,

п

д 3 (А) п ? ^2г^1г

- 2 ^(а12Т1г - 72г)71г = 0, а12 = г=1;

да12

г=1 Ет?г

г=1

2.3.1. Пусть

п Е Х2.х1. п п

71г = 0, х2г Л х1г = 0, х2г ^ ] х2. = х1г ^ ] х2.х1..

г=1 Е Х2. 7= 7=

7=1

Обозначим А = Еп=1 х2., В = Еп=1 х2.Х17, тогда х2г = х1гВ/А. Таким образом, для точек Рг = (х1г,х2г, 0) соответствующие им радиус-векторы ОРг = (х1г,х2г, 0) = (х1г,х1гВ/А, 0) = х1г(1,В/А, 0) коллинеарны, г = 1,2,..., п. Следовательно, все точки Рг лежат на одной прямой с направляющим вектором V = (1, В/А, 0), проходящей через центр О. Задача сводится к п. 2.1.

2.3.2. Пусть Еп=1 Тн = 0. Рассмотрим первую строку матрицы А:

п п п

Е 72г71г Е Х27х17 Е у.Х17 а12 = ЦП-, 71 г = Х2г - 7—Л-Х1г, 72г = у1г - 7=П-Х1г,

Е 1и Е х2. Е х1.

г=1 7=1 7=1

пп

Е у1гХ1г - ац Е Х2гХ1г

г=1 г=1 п ац =-п-, а13 = 0.

п

г1г

Ех?

г=1

Для второй и третьей строк получим

п п п

Е 72г71г Е Х27Х17 Е у2.х17

а22 = ЦП-, 71 г = Х2г - 7—Л-Х1г, 72г = у2г - 7=П-Хц

Е 7н Е х27 Е х17

г=1 7=1 7=1

пп

Е у2гХ1г - а22 Е Х2гХ1г

г=1 г=1 п

ац =-п-, а23 = 0,

Х21 г

г=1

п п п

Е 72г71г Е Х27Х17 Е у37Х17 а32 = ЦП-, 71 г = Х2г - 7—Л-Х1г, 72г = у3г - 7=П-ХН

Е 7н Е х7 Е х17

г=1 7=1 7=1

пп

у3гХ1г - а32 Х2гХ1г

г=1 г=1 п

а31 =-п-, а33 = 0.

п

1г

Х21

г=1

Вернёмся к исходной системе координат и вычислим компоненты вектора используя найденную матрицу А по формулам (1) и (2).

2.4. Пусть все точки первого облака X лежат в плоскости Ожг или Оуг. Переобозначением осей системы координат сведём задачу к п. 2.3. Используя вычисленные матрицу К и вектор получаем итоговое преобразование КМ-1ж + где

1 0 0 0 0 1

М = I 001 I или М = I 100

0 1 0 0 1 0

соответственно.

2.5. Пусть £П=1 ж?, = 0, ЕП=1 ж2г = 0 и ЕП=1 ж3г = 0, тогда

д/(А) = дат к ^

2(аижн + Й12ж2г + а13ж3, - Ун)жь = 0, к =1, 2, 3,

,=1

а^Е жк, жт, + а^жы - У1г)жкг = 0, т, п =1, 2, 3, т, п = к,

,=1 ,=1

п п Е(у1г а1тжт, а1пжпг)жкг

а1к ^ жк ^ ^(у1, а1тжт, а1nжnг)жfcг, а1к "

,= 1 ,= 1

ж2

к,

,= 1

3 (А) = Е

,=1

(у1? а1та1пжп^'

¿=1

а1тжт, + п + а1пжп, у1,

V

ж2

¿=1 2

к?

+

/

+ (а21ж1, + а22ж2г + а23Ж3, - У2,)2 + (а31жн + а32ж2, + а33ж3г - У3,)2. Рассмотрим отдельно

(у17 a1mжmj a1nжnj)жfcj

¿=1

а1тжт, + п + а1пжп, у1,

ж2

¿=1

fcJ

¿=1

^ ^ ylJ жfcJ' а1т жmJ' жfcJ' а1п ^ жnJ'I + а1тжт, + а1пжп, у1 ¿=1

¿=1

а1т

п _

Е Жk7 ¿=1 ¿=1 "

Е

ЖmJ' ЖfcJ'

+ а 1п

жп,

п _

Е Ж/Ь ¿=1 ¿=1 "

Е

ЖnJ' ЖfcJ'

У1,

п _

Е Ж I' ¿=1

Е

ylJ жfc7

2

Подставим полученное выражение в 3(А):

(

3 (А) = Е

г=1

(

а1т

хтг ~Л / J Хт7 ХЙ7

Е хк7 7=1

7=1 (

у1г

+ а 1п

Хпг

^пг Л / у Хп7 Хк7

Е хк7 7=1

7=1

V

хкг п _

Е Хк7 7=1

Е

у17 хк7

7=1

/

+

/

+ (а21Х1г + а22Х2г + а23Х3г - у2г)2 + (а31Х1г + а32Х2г + а33Х3г - у3г)2. Обозначим

7тг х

хкг

тг п

Еп

7=1

Е Х/г7 7=1

^ ^ хт7 ХЙ7,

7гг = х.

хкг

пг - п

Еп

7=1

Е Х/г7 7=1

^ ] хп7 х&7,

7кг = у1г

хкг

Е х17 7=1

^^ у 17хк7.

7=1

Тогда

3(А) = У^(а1тХтг + а^г - 7&г)2 + (а21Х1г + а22Х2г + а23Х3г - у2г)2+

г=1

+(а31Х1г + а32Х2г + а33Х3г - у3г)2.

2.5.1. Если все 7тг = 0, г =1, 2,..., п, то в качестве переменной дифференцирования вместо т берётся 1. Предположим, что все коэффициенты 7^г, г = 1, 2,... , п, равны нулю. Тогда если все 7тг = 0 и 7пг = 0, то

хтг ХЙ7 Хкг ^ ^ хт7ХЙ7, Хпг ^ ХЙ7 Хкг ^ ^ Х7Хк7, г 1, 2, . . .

7=1 7=1 7=1 7=1

п.

Пусть ЕЛ=1 Х/г7 = А, ЕЛ=1 Хт7Хй7 = В, ЕЛ=1 Х7х&7 = С, тогда получим при А = 0

_ ХкгВ _ ХкгС • _ , ^

Хтг ] , Х1г ' , г 1, 2, . . . , п.

А

А

Таким образом, Р = (хтг,хн,х^), ОРг = (хтг,Хн,х^) = Хйг(В/А,С/А, 1). Следовательно, все векторы ОРг коллинеарные, г = 1, 2,... , п, поэтому все точки Рг лежат на одной прямой с направляющим вектором V = (В/А, (С)/А, 1), проходящей через центр О. Далее применим рассуждения из п. 2.2.

2.5.2. Предположим, что не все коэффициенты 7^г, г = 1, 2,... , п, равны нулю, в таком случае

п

3(А) = У^(а1т7тг + а1п7гг - 7йг)2 + (ацх^ + а22Х2г + а23Х3г - у2г)2+

г=1

2

+ (а31Ж1, + а32Ж2г + а33Ж3г - У3,)2

д3 (А)

да

1т

2У^(а1тТтг + «^71» - 7&г)7тг = 0,

г=1

Е «1т7тг + а1п7г,7т, - 7тг = «1т ^ ^тг + «11 ^ 7н7тг - ^ 7к»7тг = 0,

г=1

г=1

г=1

г=1

п п п

2

«1т 7тг ^ ^ 7кг7тг «11 ^ у 7ií7mг, «1т

г=1 г=1 г=1

Е Тк^ аmJ• - Е Тг^ 7mJ• ¿=1 ¿=1

Е72

¿=1

3 (А) = Е

пп

Тк^ 7mJ' «11 ^ ! 7гJ 7mJ' ¿=1 ¿=1

г=1

\

V

«2

7тг + «1г7гг - 7кг

¿=1

+

/

+ (a2lЖlг + «22Ж2г + «23Ж3г - У2г)2 + ^Жн + «32Ж2г + «33Ж3г - У3г)2,

Тк^ 7mJ' «11 ^ ! 7гJ 7mJ' ¿=1 ¿=1

«2

7тг + «1г7гг - 7кг

¿=1

7тг 7mJ' «1п7тг 53 7гJ 7mJ'

¿=1 ¿=1

п

+ «1п7гг - 7кг

Е72

¿=1

п

«1п

7г»

7тг Е Тг^ Тт^ ¿=1

п

Е72

¿=1

п

7тг 53 Тк^ 7mJ'

+ ¿=1__

+ п 7кг

Е72

¿=1

«1п

7г»

V

7тг 53 Тг^ 7mJ' ¿=1 п

Е72 ■

/mJ

¿=1 " / \ ¿=1 ( ( „ ^ Л /

7кг

7тг 53 Тк^ 7mJ' ¿=1

V

Е72 ■

3 (А) = Е

г=1

«1г

V

71г

7тг Е 7т^ ¿=1

V

Е 7

2

тт/

¿=1 " / \ ¿=1 22

7кг

7тг 53 7fcJ 7mJ' ¿=1

V

53 7т /

/

+

/

+ (a2lЖlг + «22Ж2г + «23Ж3» - У2г) + ^Жн + «32Ж2г + «33Ж3» - У3»)

Обозначим через и следующие выражения:

7тг 53 7гJ 7mJ' 7тг 53 7fcJ 7mJ' ¿=1 , ¿=1 ^г = 71г--п-, г» = 7кг--п-•

Е72

¿=1

Е72

¿=1

2

2

Тогда получим

п

3 (А) = ^(а^г - "г)2 + (а^х^ + а22Х2г + а23Х3г - у2г)2 + (а31Х1г + а32Х2г + а33Х3г - у3г )2

д3 (А)

г=1

- "г)^г = 0,

дан ^

г=1

п

п п п Е №"

Е(а1п^г - "г)^г = а^ ^ ^ - Е = 0, а^ = к=П-.

г=1 г=1 г=1 Е

к=1

п2

пРи ЕЛ=1 = 0 пусть

А = ^^ Х/г7 = 0 В = ^^ Хт7хк7 , С = ^^

7=1 7=1 7=1

в таком случае

хкг ^^ __В __хкг ^^ _ С

' ~Л Хт7хк7 хт7 AХfcг, 7пг хпг л / у хп7хк7 хп7 А"

Е хк7 7=1 Е хк7 7=1

7=1 7=1

Для Е = ЕЛ=1 7^7, В = ЕЛ=1 7л7 7т7 имеем

п

7таг 7л7 7т7

7=1

В

^г 7пг Л 7пг ^ 7тг

/т7

Е7

7=1

2 Е

С \ В / В В В - С _

хпг АХк/ ^Е V Хтг АХк/ хпг Е^тг + А хкг ° г -1,2,... , n,

А А В В В - С

-у, . —— _-у* . ___-у . -у* . —— _-у* . _ _-у, .

г С В пг С В Е т^ пг Е А кг?

Р = ( ) = ^ В _ В - С \

Р г (хтг, хпг, хкг) I хтг, е хтг а хкг, хкг I ,

пд = ^ В _ В - С \ = / В Л ^ 0 С - В '

ОРг ( xmг, ЕХтг а хкг,хк:г I I xmг, EХmг, 0 I + I ° а хкг,хкг

. В \ / С - В . _ __

ХтН 1, Е, 0) + ХкИ 0, —А—, 1 I = Хтг^1 + Хкг^2, г>1

Таким образом, все ОРг лежат в плоскости (О,^1 ,^2). Далее применим рассуждения из п. 2.3.

Если ЕЛ=1 ^ = 0, то

п п п п

Е Тк77т7 - аи Е 7^7Тт7 Е №" Е(у1г - а1тХтг а11 хпг)хкг _ 7=1 7=1 _ й=1 _ г=1 а1т = -Л-, а11 = -Л-, а1к = -Л-,

Е 7т7 Х>к Е хкг

7=1 к=1 г=1

где

п N п

жкг \ Л жкг \ Л жкг

__жкг \ Л __жкг \ Л __жкг \ Л

7тг жтг п / у Жт/жк/, Т1г ж1г п / у Жг/Жк/, 7кг у1г п / у у17 Жк/,

Е Жк/ ¿=1 Е Ж/г7 ¿=1 Е Жк/ ¿=1

¿=1 ¿=1 ¿=1

пп

7тг 7г/ 7т/ 7тг 53 7к/ 7т/ ¿=1 , ¿=1 = 7г»--п-, г» = 7кг--п-.

53 7т/ 53 7т/

¿=1 ¿=1

Для второй строки:

п п п

53 7к/7т/ - «2п 53 7т/ 53 "к

_ 7=1 7=1 _ к=1

«2т = п , «2г = п ,

Е тт? Е ^к

7=1 к=1

п

Е(У2г - «2тЖтг - «2гЖгг)Жкг

_ г=1

«2к = -п-,

Е Жкг

г=1

где

жкг \ Л жкг

7тг жтг п / у Жт/жк/, 7гг ж1» п N

Е жк? ¿=1 53 ж)у 53 жг?жк/

7=1 7=1 7 7=1

п

жкг

жкг

7кг = У2г - "п- У27Жк7

Е Ж|/ 7=1 ¿=1

пп

7тг 53 7г/ 7т/ 7тг 53 7к/ 7т/

7=1 , 7=1

= 7гг--п-, г» = 7кг--п-.

53 7т/ 53 7т/

7=1 7=1

Для третьей строки:

п п п п

Етк/7т7 - <03г Е 7г/7т7 Е №"к Е(У3г - «3тЖтг «3г жпг)жкг

_ 7 = 1 7=1 _ к=1 _ г=1

«3т = -п-, «3г = -п-, = -п-;

Е 7т7 Е^к Е жк

где

/т/ Г к ¿-^ "^кг

7=1 к=1 г=1

пп жкг жкг

7тг жтг ""п Л, Жт/ жк/, 7гг ж1» ""п / у жг/ жк/,

Е Жк/ 7=1 Е Жк/ 7=1

¿=1 ¿=1

п

жкг

7кг = У3г - "п- У37 Жк7 ,

Е жк? /=1

7=1

п п

У] Т. Тт. У] Тк. Тт. 7=1 , 7=1 = 7Н--П-, ^ = 7кг--П-.

Тт. Тт.

7=1 7=1

Вернёмся к исходной системе координат и вычислим компоненты вектора используя найденную матрицу А, по формулам (1) и (2).

3. Компьютерное моделирование

Пусть облако X состоит из 10 000 точек. Координаты х и у точек случайно выбраны (с помощью равномерного распределения). Координаты относительно оси Ог вычисляются по формуле г = х2-у2. Облако У получено из облака X с помощью геометрического преобразования У = Я * X + где Я и £ описаны ниже:

/ 0.5 0 0.866025 \ / 10 \

Я = I 0.433013 0.866025 -0.25 I , £ = I 20 I . \ -0.75 0.5 0.433013 у \ 7 /

На рис. 1 показано исходное взаимное расположение облаков X и У. Здесь и далее облако X светлое, а облако У — тёмное. Результатом работы алгоритма является

Рис. 1. Исходное расположение облаков Рис. 2. Результат работы алгоритма

следующие матрица поворота и вектор параллельного переноса:

/ 0.5 1.74302е - 08 0.866025 \ / 10 \

Я= I 0.433012 0.866025 -0.25 I , ¿= I 20 I . -0.75 0.5 0.433012 7

Рис. 2 показывает взаимное расположение облаков X и У, полученных в результате применения алгоритма.

Рис. 3. Исходное расположение облаков Рис. 4. Результат работы алгоритма

На рис. 3 изображено исходное взаимное расположение ещё одной пары облаков X и У .С помощью алгоритма получаются следующие матрица поворота и вектор

параллельного переноса:

R

0.5 1.223e - 08 0.866025 0.433013 0.866025 -0.25 -0.75 0.5 0.433012

10 20 7

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

Список литературы

1. Chen, Y. Object modeling by registration of multiple range images / Y. Chen, G. Medioni // Image and Vision Computing. — 1992. — Vol. 10, no. 3. — P. 145-155.

2. Besl, P. A method for registration of 3-D shapes / P. Besl, N. McKay // IEEE Transactions of Pattern Analysis and Machine Intelligence. — 1992. — Vol. 14, no. 2. — P. 239-256.

3. Turk, G. Zippered polygon meshes from range images / G. Turk, M. Levoy // Computer Graphics Proceedings. Annual Conference Series, ACM SIGGRAPH. — 1994. — P. 311-318.

4. Horn, B. Closed-form solution of absolute orientation using unit quaternions / B. Horn // Journal of the Optical Society of America. Series A. — 1987. — Vol. 4, no. 4. — P. 629-642.

5. Horn, B. Closed-form solution of absolute orientation using orthonormal matrices / B. Horn, H. Hilden, S. Negahdaripour // Journal of the Optical Society of America. Series A. — 1988. — Vol. 5, no. 7. — P. 1127-1135.

6. Du, S. An extension of the ICP algorithm considering scale factor / S. Du, N. Zheng, S. Ying, Q. You, Y. Wu // Proceedings of the 14th IEEE International Conference on Image Processing. — 2007. — P. 193-196.

7. Du, S. Affine registration of point sets using ICP and ICA / S. Du, N. Zheng, G. Meng, Z. Yuan // IEEE Signal Processing Letters. — 2008. — Vol. 15. — P. 689-692.

8. Du, S. Affine iterative closest point algorithm for point set registration / S. Du, N. Zheng, S. Ying, J. Liu // Pattern Recognition Letters. — 2010. — Vol. 31. — P. 791-799.

9. Rusinkiewicz, S. Efficient variants of the ICP algorithm / S. Rusinkiewicz, M. Levoy // Proceedings of the International Conference on 3-D Digital Imaging and Modeling. — 2001. — P. 145-152.

10. Tihonkih, D. A modified iterative closest point algorithm for shape registration / D. Tihonkih, A. Makovetskii, V. Kuznetsov // Proceedings SPIE Applications of Digital Image Processing XXXIX. — 2016. — Vol. 9971. — 99712D.

Поступила в 'редакцию 28.09.2017

После переработки 16.10.2017

Маковецкий Артём Юрьевич, кандидат физико-математических наук, доцент кафедры вычислительной механики и информационных технологий, Челябинский государственный университет, Челябинск, Россия; e-mail: artemmac@csu.ru.

Воронин Сергей Михайлович, доктор физико-математических наук, доцент, профессор кафедры математического анализа, Челябинский государственный университет, Челябинск, Россия; e-mail: voron@csu.ru.

Тихоньких Дмитрий Вадимович, ассистент кафедры вычислительной механики и информационных технологий, Челябинский государственный университет, Челябинск, Россия; e-mail: sparsemind@gmail.com.

Алексеев Михаил Николаевич, кандидат педагогических наук, доцент кафедры вычислительной механики и информационных технологий, Челябинский государственный университет, Челябинск, Россия; e-mail: alexeev@csu.ru.

Сведения об авторах

Chelyabinsk Physical and Mathematical Journal. 2017. Vol. 2, iss. 3. P. 282-294.

CLOSED FORM SOLUTION OF ICP ERROR MINIMIZATION PROBLEM FOR AFFINE TRANSFORMATIONS

A.Yu. Makovetskiia, S.M. Voroninb, D.V. Tihonkihc, M.N. Alekseevd

Chelyabinsk State University, Chelyabinsk, Russia

aartemmac@csu.ru, bvoron@csu.ru, csparsemind@gmail.com, dalexeev@csu.ru

The iterative closest point (ICP) algorithm is one of the most popular approaches to shape registration. The aim of the algorithm is to calculate the optimal geometric transformation relative to the given metric, combining the two given clouds. An important step in the ICP algorithm is the solution of the problem of minimizing the functional corresponding to a given metric for a given class of geometric transformations. In this paper, a method is presented for solving the variational problem of the ICP algorithm for the point-to-point metric in the class of affine transformations. With the help of computer simulation, the correctness of the proposed method is demonstrated.

Keywords: 3D reconstruction, registration of point clouds, localization, affine transformation

References

1. Chen Y., Medioni G. Object modeling by registration of multiple range images. Image and Vision Computing, 1992, vol. 10, no. 3, pp. 145-155.

2. Besl P., McKay N. A method for registration of 3-D shapes. IEEE Transactions of Pattern Analysis and Machine Intelligence, 1992, vol. 14, no. 2, pp. 239-256.

3. Turk G., Levoy M. Zippered polygon meshes from range images. Computer Graphics Proceedings, Annual Conference Series, ACM SIGGRAPH, 1994, pp. 311-318.

4. Horn B. Closed-form solution of absolute orientation using unit quaternions. Journal of the Optical Society of America. Series A, 1987, vol. 4, no. 4, pp. 629-642.

5. Horn B., Hilden H., Negahdaripour S. Closed-form solution of absolute orientation using orthonormal matrices. Journal of the Optical Society of America. Series A, 1988, vol. 5, no. 7, pp. 1127-1135.

6. Du S., Zheng N., Ying S., You Q., Wu Y. An extension of the ICP algorithm considering scale factor. Proceedings of 14th IEEE International Conference on Image Processing, 2007, pp. 193-196.

7. Du S., Zheng N., Meng G., Yuan Z. Affine registration of point sets using ICP and ICA. IEEE Signal Processing Letters, 2008, vol. 15, pp. 689-692.

8. Du S., Zheng N., Ying S., Liu J. Affine iterative closest point algorithm for point set registration. Pattern Recognition Letters, 2010, vol. 31, pp. 791-799.

9. Rusinkiewicz S., Levoy M. Efficient variants of the ICP algorithm. Proceedings of the International Conference on 3-D Digital Imaging and Modeling, 2001, pp. 145-152.

10. Tihonkih D., Makovetskii A., Kuznetsov V. A modified iterative closest point algorithm for shape registration. Proceedings SPIE Applications of Digital Image Processing XXXIX, 2016, vol. 9971, 99712D.

Accepted article received 28.09.2017 Corrections received 16.10.2017