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

ВИЗУАЛЬНАЯ НАВИГАЦИЯ АВТОНОМНОГО НЕОБИТАЕМОГО ПОДВОДНОГО АППАРАТА С ГЛОБАЛЬНЫМ ПОИСКОМ СВЯЗЕЙ МЕЖДУ ИЗОБРАЖЕНИЯМИ

А.Н. Камаев1, Д.А. Карманов1 1 Вычислительный центр ДВО РАН, Хабаровск, Россия

Аннотация

В статье рассматривается задача навигации автономного необитаемого подводного аппарата. Изображения, получаемые бортовой стереокамерой, используются для построения облаков точек, привязанных к конкретным положениям аппарата. Квантованные дескрипторы точек хранятся в метрическом дереве для организации процедуры быстрого поиска с использованием алгоритма best bin first. Связи для нового облака точек ищутся в небольшой группе облаков, которые имеют наибольшее число похожих дескрипторов, сохранённых в дереве. Таким образом, новое облако может быть позиционировано относительно других облаков без какой-либо априорной информации о положении АНПА и точности, с которой это положение известно. Это приводит к повышению надежности системы навигации, делает систему нечувствительной к пропускам данных, однотонным участкам дна и длинным проходам без пересечений собственной траектории. В статье рассматриваются алгоритмы построения облаков точек, сопоставления этих облаков, построения групп потенциально связанных облаков для ускорения процедуры глобального поиска связей. Также предлагается основной алгоритм навигации, состоящий из трёх выполняемых параллельно подпрограмм: добавление изображений в систему, обновление дерева поиска, оптимизация. Предложенный алгоритм навигации тестируется как на реальных, так и на синтетических данных. Тесты на реальных данных показывают, что траектория может быть построена даже для последовательностей с 60 % пропущенных изображений и малым или нулевым перекрытием последовательных снимков. Тесты на синтетических данных демонстрируют, что построенная траектория близка к истинной даже для длительных миссий. Средняя скорость работы навигационной системы составляет 3 кадра в секунду на процессоре среднего ценового диапазона.

Ключевые слова: АНПА, SLAM, особые точки, счисление пути, сопоставление изображений.

Цитирование: Камаев, А.Н. Визуальная навигация автономного необитаемого подводного аппарата с глобальным поиском связей между изображениями / А.Н. Камаев, Д.А. Карманов // Компьютерная оптика. - 2018. - Т. 42, № 3. - С. 457-467. - DOI: 10.18287/2412-6179-2018-42-3-457-467.

Введение

Использование фото- и видеоданных с бортовых камер позволяет наиболее точно позиционировать автономный необитаемый подводный аппарат (АНПА) вблизи дна. Существует большое количество методов построения навигационных систем (НС), опирающихся на визуальную информацию. Наиболее простой метод - это счисление пути, он использует найденные трансформации между соседними кадрами видеопоследовательности для вычисления положения АНПА [1]. Этому методу свойственны все недостатки инерциальной НС, но в ряде случаев он позволяет дольше сохранять приемлемую точность.

Наиболее распространённым подходом к визуальной навигации с одновременным построением карты (SLAM - Simultaneous Localization and Mapping) является использование расширенного фильтра Калма-на (EKF - Extended Kalman Filter), когда вектор состояния включает положение и ориентацию АНПА и положение ориентиров (landmarks) на карте [2, 3]. Новые ориентиры добавляются по мере получения новых изображений. При повторном посещении уже записанного ориентира неопределённость вектора параметров уменьшается (закрытие разрыва). Для закрытия разрыва используются точечные особенности на изображениях и эллипсоид поиска, ограничиваю-

щий количество потенциально соответствующих ориентиров. Так как матрица ковариаций вектора параметров при использовании EKF хранится в явном виде, эллипсоид поиска всегда можно оценить. Проблемой использования EKF является квадратичный рост сложности задачи по мере добавления новых ориентиров к вектору параметров. Эту проблему можно частично решить, разбив карту на локальные связанные зоны, каждая из которых строится с использованием EKF, тогда как сами зоны позиционируются путём решения задачи минимизации на графе [4].

Другим способом справиться с квадратично растущей сложностью EKF является использование расширенного фильтра информации (EIF - Extended Information Filter) [5]. В отличие от EKF, EIF оперирует матрицей информации, обратной к матрице ковариаций. Она является разреженной, что позволяет значительно сократить вычисления, необходимые для обновления вектора состояний и самой матрицы. В [6] показано, как можно оценить ковариации отдельных элементов вектора параметров для решения задачи закрытия разрывов, не вычисляя при этом матрицу, обратную к матрице информации.

Наиболее точным способом решения задачи навигации является подход, основанный на построении графа состояний АНПА, но он применяется значи-

тельно реже из-за сложности реализации в режиме реального времени. В [7] рассматривается оптимизация поиска связей между изображениями при использовании такого подхода. Оптимизация заключается в том, что для поиска отбираются изображения из некоторой окрестности предполагаемого положения только в том случае, если они имеют выраженный и достаточно уникальный текстурный рисунок. В [8] рассматривается подход, основанный на графе применительно к стереоизображениям.

Несмотря на то, что для решения задачи визуальной навигации АНПА могут быть применены самые разные подходы, их все объединяет один важный недостаток - для закрытия разрывов требуется сопоставлять текущее изображение со всеми попадающими в эллипсоид неопределённости положения. Для дна с хорошим текстурным покрытием и при условии большого пересечения последовательных снимков неопределённость положения невелика и для закрытия разрывов достаточно перебрать небольшое число изображений. При однотонной донной поверхности, малом перекрытии, пропуске кадров, длительных проходах без пересечения траектории неопределённость положения возрастает настолько, что задача закрытия разрывов становится нерешаемой. Чтобы иметь возможность закрывать разрывы независимо от знания о положении АНПА и неопределённости этого положения, в статье предлагается использовать глобальный поиск дескрипторов особых точек при сопоставлении изображений. Рассматривается система навигации по стереоизображениям, основанная на графе, использующая глобальный поиск связей между изображениями и способная работать в режиме реального времени.

Обозначения. Векторы и матрицы обозначены полужирным шрифтом: А - матрица, а, ат - векторы столбец и строка. Для обозначения поворота в М3 используется либо матрица поворота Я: К3^К3, либо вектор поворота геМ3, такой что |г| < п - угол, а г/|г| -ось поворота. Функции г(Я) и Я(г) отвечают за перевод матрицы поворота в вектор и обратно. Через diag(a), diag(a0, аь ...) обозначены диагональные матрицы с элементами а0, аь ... на диагонали.

1. Постановка задачи

На корпусе АНПА установлены две камеры, которые могут быть описаны с использованием модели камеры - обскуры (все нелинейные искажения должны быть устранены заранее). Будем считать, что матрицы внутренних параметров камер Щ и Щ известны.

Также известны матрицы поворота камер ^ , Я 2 и их положения ^ и 12 в системе координат АНПА. Известны ^ и г - положение и ориентация АНПА в глобальной системе координат в начальный момент времени. На аппарате работает инерциальная НС, позволяющая определить положение и ориентацию аппарата в момент получения к-го изображения: \\, г/ .

Инерциальная НС даёт накапливающуюся со временем ошибку. Обозначим через o't и о'г векторы средне-квадратических отклонений положений и ориентаций, вычисляемых инерциальной НС за 1 секунду.

Необходимо, используя визуальную информацию с камер, определять положение и ориентацию АНПА в момент получения к-го изображения tg, rg, к = 1, 2, ..., M- 1, где M - количество изображений, получаемых АНПА за миссию. В режиме фотообследования АНПА движется с небольшой скоростью, и для обеспечения достаточного перекрытия снимков требуемая частота съёмки составляет от 0,5 до 2 кадров в секунду. На участках, где удалённость или характер дна не позволяет использовать визуальные особенности, может быть использована инерциальная НС.

2. Построение облака точек

На изображениях, получаемых со стереокамер АНПА - стереопарах, выполняется корректировка освещения [9]. Рассмотрим работу с к-й стереопарой (при обозначении объектов, локальных для к-й стереопары, будем опускать индекс к). На изображениях стереопары извлекаются точечные особенности Xll = (Хщ, ХЦ2, 1)T, X2p = (X2p,1, X2p,2, 1)T и их SIFT-дескрипторы dn, d2p [10], где l = 1, 2, ..., L; p = 1, 2, ..., P; L и P - количество особенностей, обнаруженных на первом и втором изображении соответственно.

2.1. Модель связи

Соответствующие точки первого и второго изображения стереопары xj, x2 связаны эпиполярным ограничением x^^Fxj = 0, где F - фундаментальная матрица [11],

F = K2-TR2TR; [Rit12 - RjcTtc ]x K;-1.

Согласно эпиполярному ограничению точка на втором изображении x2, соответствующая точке x1 , лежит на линии, заданной коэффициентами Fx1.

Определим функцию расстояния от точки до прямой

d ((a,,)T, (Xo, X;,1)T ) =

= (a0x0 + alXj + a2)/sja, + a2 .

Точка первого изображения стереопары с номером i соответствует точке второго изображения стереопары с номером j, если выполнены следующие условия:

dv= d

id,i - d 2 j |< + mlb,

max

-^d-< d mi-mL+Щь,

mm(|d, - d2p|,|d,i - d2j|) j dmax 2Ь ^

p e (1,2, ...,P}\{j}, l e (1,2, ...,L}\{i},

где dmax, m1a, m1b, m2a, m2b - параметры алгоритма сопоставления. В параграфе 6 использованы следующие значения этих параметров: dmax = 2, m1a = 0,3, m1b = 0,5, m2a=0,8, m2b = 0,9. Пары индексов всех соответствую-

щих точек формируют множество P = {(i е {1, 2, ..., L}, j е {1, 2, ..., P}) | (1) л (2) л (3)}.

2.2. Данные облака Каждая стереопара позволяет сгенерировать облако точек - множество Bk. Каждая точка b е Bk описывается следующими свойствами: b.x = (b.x0, b.xb b.x2)T - положение в системе координат АНПА, b.C - матрица ковариаций вектора b.x, b.R - матрица, строки которой являются собственными векторами матрицы b.C, b.s = (b.s0, b.sj, b.s2)T - соответствующий вектор собственных чисел и b.d - SIFT-дескриптор.

Задача определения координат точки b.x по проекциям x1i, x2j , (i, j) е P и известным матрицам камеры K (R£T | -RkTtk), k = 1,2, называется триангуляцией и может быть решена различными способами [11]. Авторы для определения начального приближения использовали direct linear transformation (DLT) метод [11], а для уточнения результата - метод Ньютона. Матрица

b.C =

где J: = (JUp), J2 = (J2,lp), J\,p = (dxu,p /5b. xi), J2,lp = (dx2j,p /db. xl), l = 0, 1, 2, p = 1, 2, а Zx - матрица ковариаций вектора (xU1, xU2, x2j;1, x2j,2). Для тестирования алгоритма в параграфе 6 выберем £x = а2хE , где E -

единичная матрица, а стх = 2 пикселя - среднеквадрати-ческое отклонение измерений координат особых точек.

Дескриптор b.d рассчитывается как средний вектор между дескрипторами соответствующих точек на стереопаре: b.d = (d^+d2j) /2. Каждой паре сопоставленных точек (i,j) е P соответствует одна точка в множестве Bk. Таким образом, | P | = | Bk| - количество точек в k-м облаке.

2.3. Компактное представление данных

При большом количестве изображений в системе актуальным становится вопрос хранения данных облаков в оперативной памяти. Большая часть данных может храниться на диске и загружаться при необходимости, но для алгоритма поиска потенциально связанных облаков точек требуется присутствие в памяти дескрипторов b.d, b е Bk, k = 0, 1, ...,M- 1.

Эффективным способом уменьшения объема памяти, занимаемого дескрипторами b.d, является квантование. Компоненты дескрипторов предлагается квантовать по группам (product quantization) [12]. Приемлемая точность поиска связанных облаков (параграф 4.1) может быть обеспечена, если разбить вектор SIFT-дескриптора на 16 групп по 8 компонент в каждой. При таком разбиении достаточно словаря из 256 слов для каждой группы. Таким образом, 128-мерный SIFT-дескриптор b.d может быть представлен 16-мерным вектором b.d', каждый компонент которого представлен однобайтным индексом в словарь соответствующей группы. Операция построения словарей может быть выполнена с использованием алгоритма K-средних до выполнения миссии. При этом

для построения выбираются дескрипторы с изображений, типичных для акватории, в которой будет выполняться миссия. Подробно процесс квантования описан в работе [12].

3. Вычисление связи между облаками точек

Два облака точек Bk и B, могут содержать некоторое количество пар соответствующих точек (имеющих одинаковое положение в глобальной системе координат). Наличие таких точек позволяет рассчитать преобразование, связывающее точки из облака Bi с точками из облака Bk.

3.1. Определение точечных соответствий

Рассмотрим точки bk е Bk и bl е Bl, k Ф1. Близость этих точек характеризует функция асимметричного расстояния d(b, .d,bk.d') (ADC) [12]. Асимметричное расстояние, вычисляемое между полным дескриптором и его квантованной версией, согласно [12] точнее аппроксимирует истинное расстояние |b1.d - bk.d|, чем симметричное расстояние, вычисляемое между двумя квантованными дескрипторами.

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

d(b b) =_d(bi.d,bkd)_

rKk' min (d(b. .d, bk .d'),d(b1 .d, b.. .d'))'

bj е Bt \ [b,},bt е Bk \ [bk} .

Множество пар потенциально соответствующих точек bk е Bk и bi е Bi

B'„ = {(bk,b1) | d(b,.d,bk.d') < mlb лdr (bk,b) <m2b}.

3.2. Фильтрация ложных соответствий

Множество B'u может содержать ложные соответствия. Применение стандартной процедуры RANSAC [13] для фильтрации ложных соответствий затруднено тем, что для вычисления границ, в рамках которых точка может считаться инлайером для подобранной модели, требуются затратные по времени вычисления. Отфильтровать ложные соответствия можно опираясь на тот факт, что расстояние между любыми двумя точками одного облака должно быть равно расстоянию между соответствующими точками другого облака. Рассмотрим функцию - счётчик от четырёх точек (bk, bi), (ck, c,) е B'kl

I(bk ,ck ,bi ,ci) = Л (bk, c

1 A(bk, ck, b,, c,) < A(bk, ck, bi, ci ^ 0, A(bk,ck,bi,ci) > A(bk,ck,bi,ci), k, ck, bi, ci) = l! bkx - ckx | -| bix - cix I,

(

A(bk, ck, bi, ci) = 8

bt, cf

bkx - ckx \bk x - ck .x|

Л

+ 5

bi, cv

b 1 .x - c I .x

8(Ь,С,у ) = ЦЬ,у) + Цс,у),

где 1(Ь, у) - расстояние вдоль V от точки Ь до точки пересечения с границей эллипсоида ковариаций вектора Ь.х, построенного на расстоянии в 5 раз превышающем дисперсию Ь.х (предположив нормальное распределение ошибки в определении Ь.х, согласно правилу трёх сигм, выберем 5 = 32). Функцию 1(Ь, у), |у| = 1, можно определить, подставив точку 1(Ь, у)-Ь.Ку в уравнение эллипсоида и выразив

t(b,v) =

ф • b.s0 b.Sj b.s2 • g • diag (b.Rv)b.Rv g • diag(b.Rv)b.Rv

где g = (Ь.51 Ь.52, Ь.50 Ь.52, Ь.50 Ь.51). Используя функцию 1(Ьк, ск, Ь,, с), можно для каждой пары соответствующих точек определить количество других пар точек, расстояния до которых в обоих облаках одинаково:

П(Ь/, Ь,) = X I (Ь/, с/, Ь,, с), (Ь/, Ь,) е Б'и.

(ск .с )еБи

Введём множество

В/,(N) = {(Ь/,Ь,) е Б'|п(Ь/,Ь,) > Ж}, тогда множество истинных соответствий

БЫ = В/, (minN>|В/ (Ж)| Ж).

Рассмотрим зависимость |Б/| от | Б'и | для связанных облаков точек и несвязанных, где точечные соответствия определены случайным образом. Рис. 1 отображает зависимость, построенную для 20000 связанных и 20000 несвязанных облаков. Из рис. 1 видно, что эти множества могут быть разделены прямой линией.

\Вк1\\ '

0 10 20 30 40 50 60 70 80 90 \В'Н\

Рис. 1. Соотношение истинных и ложных соответствий для связанных и несвязанных облаков точек

Факт наличия связи можно установить, если выполняется условие Даь р1) = |БИ| - а1 - р1| Б'ы | > 0, где

а1 = 5 и р1 = 0,3 (такие значения позволяют отфильтровать ложные и неустойчивые связи). Надёжность связи опишем величиной

(

Vi = max

(

0,05, min

1,

f^P.)

f

//

x min

1,

B

K,

где К5 - минимальное количество точек, соответствующее надежной связи. В параграфе 6 используются следующие значения: а2 = 10 и р2 = 0,4 и К5 = 25.

3.3. Модель связи Для (Ь/, Ь,) е Б/1 справедливо соотношение

Ь/х = К Л х+1 «,

где ^ - вектор смещения между облаками, а Я« -матрица поворота. Чтобы найти К,, воспользуемся алгоритмом Кабша [14]. Для этого рассчитаем:

- 1 X Ьк.x, х, =7^ X Ь1х,

A = S (Ьк.x-x)(bl.x-x)T.

(Ьк A)eBu

Если A = UDVT - сингулярное разложение матрицы A, то матрица R« = U diag (1, 1, det (UVT))VT, а вектор ttt = xt - R^ xl. Хранить поворот будем в виде вектора поворота ги = r (R^).

3.4. Уточнение параметров связи

Для пары точек (Ьк, Ь) e Ba определим вектор

Ч ,Ь, (Гк/,t Ы ) =

= ((и )bl .x +t« - К .x)T , ( .x - b .x)T ) , где b .x - истинное положение точки bl. Допуская нормально распределённую погрешность в определении b^x и bl.x с соответствующими матрицами ковариаций Ьк.С и bl.C, запишем ошибку связи

EBtJ (r«,t Ы ) = ^ £bt ,b, (Гк/,t Ы ^Wbt ,b, Ч ,b, (Гк/,t Ы ) ,

Wbt A = diag(bt .С-1, b, .С-1).

Уточнить параметры связи ги и ta можно, минимизировав ошибку E B (rH, t ы) относительно параметров

rH, tH e К3 и координат точек Ь , .x e К3, (Ьк, Ь) e Ba. Такую задачу можно решить, используя метод bundle adjustment (BA) [11], выбрав в качестве начального приближения параметры связи, рассчитанные в подпара-графе 3.3 и точки Ь ,.x = Ь , x . Далее под параметрами

связи p ы = (r^, t будем понимать параметры, уточненные в результате решения задачи минимизации. На заключительной итерации BA рассчитывается матрица

- матрица ковариаций вектора рк.

Для повышения скорости вычислений и уменьшения машинной ошибки вместо матрицы в качестве веса связи будем использовать вектор w^ = (w^,), f Л '

wrn. = Vu min

w,.

№

,1

i = 0,1,...,5.

где Ери,й - диагональный элемент £Ри, Vkl определёна в подпараграфе 3.2, а wt - среднеквадратическое отклонение i-го параметра связи, при котором она может считаться надёжной. Возьмем w012 = 10-3 радиан

и W3 4 5 = 10-2 метров для компонент, соответствующих вектору поворота и вектору перемещения.

4. Поиск потенциально связанных облаков

Для поиска связей в режиме реального времени для каждого нового облака точек необходимо сформировать группу из M' е {10, 11, ..., 30} облаков, потенциально связанных с ним. Потенциально связанными будем считать облака, имеющие наибольшее количество точек с похожими дескрипторами.

4.1. Структура данных для хранения

дескрипторов точек

Для быстрого поиска похожих дескрипторов точек (задача поиска k ближайших соседей) будем использовать метрическое дерево [15] совместно с алгоритмом приближённого поиска best bin first [16, 17]. Алгоритм best bin first позволяет сократить время поиска, ограничив количество просматриваемых листовых узлов дерева до константы E. Квантование дескрипторов (подпараграф 2.3) позволяет ускорить поиск за счёт быстрого вычисления расстояния между дескрипторами. Рассмотрим, насколько эффективно работает структура данных (СД) при поиске среди 106 дескрипторов точек морского дна - рис. 2. Будем искать k = 10 ближайших соседей и оценивать количество

дескрипторов k , входящих в 10 ближайших, но не ненайденных при поиске. Параметр G определяет количество групп квантования дескриптора (подпараграф 2.3). Из рис. 2 видно, что k резко падает при переходе от 8 до 16 групп, а далее снижается медленнее. Значение G = 16 является хорошим компромиссом между объемом занимаемой памяти и точностью поиска. Увеличение E линейно увеличивает время, при этом значительное повышение точности поиска наблюдается в интервале значений [50; 100]. Параметр L определяет количество дескрипторов в листовых узлах дерева. Параметры E и L не влияют на потребляемую память, поэтому могут выбираться исходя из требуемого соотношения точности и времени поиска. В параграфе 6 используются значения E = 100, L = 20.

4.2. Обновление СД

Первое дерево поиска строится лишь тогда, когда в СД накопится достаточно большое количество дескрипторов, например, более 5000, до этого поиск выполняется полным перебором. Когда дерево построено, новые дескрипторы добавляются в него посредством записи в соответствующий листовой узел. Если количество дескрипторов в листовом узле превышает L, то он помечается не листовым и делится на два новых листовых узла гиперплоскостью, перпендикулярной направлению максимальной дисперсии дескрип-

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

п— ? - J Ю

F =40 ь -4

J Ь- 10 0

Е / JU -7ПП

V

Т =.1 б< V \

S V \ \

S ^— \

G 2 ¿4

S \ \ \

G \ \

ч.

о

время, л

0 0,5 1,0 1,5 2,0 2,5

Рис. 2. График зависимости количества ненайденных дескрипторов от времени при поиске 10 ближайших соседей в базе данных из 106 БШТ-дескрипторов

Для сохранения оптимальной глубины дерево периодически перестраивается. Запишем алгоритм обновления СД (АОСД). Через N обозначим количество дескрипторов в СД, через avg_d - усреднённую максимальную глубину посещаемых листьев по 1000 последних запросов.

Шаг 1: Если N < 5000, или новое дерево строится в настоящий момент, или avg Д - |_^2(М£) + 1] < 2, то закончить алгоритм.

Шаг 2: Присвоить N ^ N, построить новое дерево для N' дескрипторов, во время построения поисковые запросы обрабатывает старая СД.

Шаг 3: Пока N' Ф N добавлять в новое дерево дескриптор с номером N (предполагается нумерация с нуля) и увеличивать N : N ^ N + 1.

Шаг 4: Заблокировать возможность добавления дескрипторов в СД, повторить шаг 3, заменить прошлую СД новым деревом, разблокировать возможность добавления дескрипторов в СД.

АОСД выполняется в отдельном потоке и блокирует работу НС лишь на шаге 4. Для организации ал-

горитма удобно объединить квантованные дескрипторы всех облаков в непрерывном массиве данных. Тогда любой дескриптор в системе может быть найден по номеру /': , где I = 0, 1, ..., Ж- 1, а N - количество дескрипторов в системе. Дерево хранит в листовых узлах лишь индексы дескрипторов, а облака точек хранят указатели на группу дескрипторов облака в массиве. СД также хранит номера облаков, которым принадлежат дескрипторы - img(/).

4.3. Алгоритм поиска потенциально связанных облаков точек Пусть в навигационную систему добавлено М изображений с номерами т = 0, 1, ..., М- 1. Для нового изображения с номером М требуется найти М' < М изображений, потенциально с ним связанных. Введём алгоритм построения вспомогательного множества К^) для некоторого дескриптора d:

Шаг 1: Присвоить К^) ^ 0.

Шаг 2: Для всех d;, I е ОД) с {0,1,..., N -1}, где ОД) -подмножество индексов дескрипторов, близких к d, определяемое СД: если 3/ е К^), img(/) = img(k.n), то выполнять шаг 3а, иначе - шаг 3б. Шаг 3а: Если /Л1>( d;), то присвоить:

/.ё2 ^ / Д, /Д ^ (1 (d, dI'), /.п ^ /',

иначе, если /.ё2> (1 (d, di),

то присвоить /.ё2 ^ (1 (d, d'). Шаг 3б: Присвоить: ^ (1 (d, d,'), 1Л2 ^ 72, 1.п ^ I.

Если |К(Д)| < 10, то K(d) ^ K(d) и {,}, иначе:

р = argmax/еК № Ы^

если 1.ё1 < р.йх,

то К^ К(d)\{р} и {,}.

Алгоритм поиска потенциально связанных облаков точек (АППСОТ):

Шаг 1: Присвоить ут ^ 0.

Шаг 2: Для всех Ь е БМ и / е К(Ь^) присвоить

^img

+ wa (k.dl)wb

( k.d1 Л v k.d2 у

Шаг 3: Построить такое множество индексов I (M) с I = {0, 1, ..., M - 1}, что

min v. > max v, л 11(M)|= M'.

ieI

I(M) - множество индексов изображений, потенциально связанных с изображением с номером M. Задача весовых функций wa(x) и wb(x) - назначить больший вес дескрипторам, имеющим большие шансы соответствовать дескриптору запроса. В статье использованы следующие весовые функции:

Wa (x) = e

wb (x) = min

1,11 --

x - 0,7 0,3

5. Алгоритм навигации Основной операцией алгоритма навигации является поиск состояния АНПА pf = ((rf)T,(tf)T) в момент получения /-го изображения на основе обнаруженной связи pkl = (r¿,tT) (подпараграф 3.4), k< /:

pf = f (pf, p Uh где

f (pf, pkl) = (r (R(rkf )R(rw))T , (R(rkf )t„ + tf )T).

НС поддерживает список связей изображений L, в момент запуска L = 0. Запишем алгоритм, выполняемый НС при добавлении l-го изображения, / > 0 (алгоритм добавления изображения - АДИ):

Шаг 1: Вычислить Bl

(подпараграф 4.3, АППСОТ). Шаг 2: I'(l) 0, для k eI(/), если есть связь между облаками Bk и Bl: рассчитать pkl и wkl (4) (параграф 3), присвоить I'(l) ^ I'(l) u {k}. Шаг 3: Если / - 1 g I'(l), построить связь

p/-i,/ = (r/T-i,/, tT-1,/), где г/-1/ = r(R(r/I-i)T R(r/)), а t/-11 = R(r/-1 )T (tI -11-1), рассчитать wkl (4), взяв Vi-u =1 и =diag(A/-i,/о'rT, Д-i,/o'tT )2, где Д/_и -время между получением изображений l - 1 и l, присвоить I'(/) ^ Щ) u {/ - 1}. Шаг 4: Рассчитать pk = (r/kT,t\T) = f(pk,pk/), k eI'(/). Относительно параметров rf, tf e M3 решить задачу

X k diag(wM КT ^ min,

где г) = (r(R(rf )T R(r/ ))T, (tf -1] )T). Шаг 5: Рассчитать D¡ = maxkeI,(/} 11f -1k¡ |. Шаг 6: Дополнить список связанных изображений:

L ^ L u {(/, k) | k e I'(/)}. Шаг 7: Заменить дескрипторы точек облака B/ на квантованные (подпараграф 2.3) и поместить их в СД, выполнить АОСД (подпараграф 4.2).

Рассчитанная на шаге 5 АДИ величина D/ отслеживается процессом, работающим параллельно и отвечающим за корректировку траектории в случае, когда это необходимо согласно следующему алгоритму (алгоритм корректировки траектории - АКТ):

Шаг 1: Бездействовать, пока

Dj< Dmax, i = 1, 2, ...,M -1, где M - количество изображений в системе, а Dmax - допустимая погрешность связи. Шаг 2: Блокировать АДИ, присвоить M' M и L^ {(i,j) e L |j> k}, где k = min j,

D, >Dmax j)eL

разблокировать АДИ.

Шаг 3: Относительно параметров pf е М6,

i = к + 1, к+2, ..., M' - 1 решить задачу минимизации

X (i, У )er£iiTdiag(w « min

' r (R(rf )TR(rf )R(r,) f

£ =

j K)T

Шаг 4: Блокировать АДИ, применить решение задачи шага 3, пока M'фM: рассчитать pM' (шаг 4, АДИ), присвоить M': M' M' + 1. Разблокировать АДИ.

Задача минимизации на шаге 3 может быть сформулирована различным образом, в зависимости от того, какие данные навигационная система АНПА определяет точно (часто это может быть глубина, курс, крен и дифферент), а какие требуют уточнения. Подобная задача может быть эффективно решена, если учитывать разреженность матрицы системы. В работе [18] представлены несколько подходов к решению подобной задачи.

6. Результаты Рассмотрим работу НС на трёх наборах данных: набор из 3591 изображения с АНПА «Клавесин - 1Р» и два набора по 11000 и 10006 изображений виртуальной акватории.

Данные с АНПА «Клавесин - 1Р» обладают рядом особенностей: 3591 снимок покрывает лишь отдельные участки траектории, в общей сложности не более 40 %, последовательные снимки получены с интервалом времени в 2 секунды, что при средней скорости 0,7 м/с приводит к малому или отсутствующему перекрытию, дно на снимках вымощено однотипными панцирями ископаемых морских ежей - рис. 3.

Рис. 3. Пример снимков с АНПА «Клавесин - 1Р»

АНПА «Клавесин - 1Р» не обладает стереокамерой, поэтому облака точек для каждого изображения строятся на основе данных о расстоянии АНПА до дна в момент получения снимка и предположения о плоском дне в рамках получаемого снимка. На рис. 4а изображена траектория АНПА, полученная комплексирован-ной инерциальной и гидроакустической НС. В режиме постобработки была смоделирована работа предложенной НС - рис. 46. Связи между изображениями показаны черными линиями: сплошными - связи, полученные в результате сопоставления изображений, пунктирными - связи на основе данных инерциальной и гидроакустической НС. Рассмотрим вектор ошибки связи - вектор с началом в положении АНПА, предсказанном связью, и концом в положении, рассчитанном НС. Векторы ошибок показаны серыми линиями на рис. 4. Точность работы НС можно оценить с помощью длин векторов ошибок связей.

Рис. 4. Участки траектории АНПА «Клавесин-1Р»: полученной инерциальной НС с коррекцией по гидроакустическим маякам (а), полученной с использованием визуальной навигационной системы (б).

Размерность сетки - 10x10 метров

Рассмотрим связи, соединяющие непоследовательные изображения (т.е. связи, закрывающие разрывы). В общей сложности было найдено 3179 таких связей. Длины векторов ошибок этих связей варьировались от 8 • 10 - 4 до 0,96 м, средняя длина составила 0,07 м, а среднеквадратическое отклонение - 0,08 м. Если оценить траекторию АНПА, полученную без использования визуальной НС, с точки зрения этих же связей, то длины векторов их ошибок попадут в диапазон от 0,13 до 24,65 м, средняя длина составит 6,69 м, а среднеквадратическое отклонение - 4,18 м.

При построении траектории с помощью предложенной НС удалось достичь скорости 3 кадра в секунду на процессоре Intel Core Quad Q9550, что в шесть раз превысило скорость получения изображений в рассматриваемом примере. При этом разрешение снимков составляло 768 x 576 пикселей.

Для получения синтетических данных использовался алгоритм построения дна [19] и симулятор [20]. Стереоизображения с разрешением 640 х 480 и сте-реобазой 0,4 м были получены с помощью бортовой камеры виртуального АНПА - рис. 5. Синтетические данные позволяют сопоставить траекторию АНПА, рассчитываемую НС, с истинной траекторией, которая остаётся неизвестной при экспериментах в реальных акваториях.

Рис. 5. Пример снимков с бортовой камеры виртуального АНПА

В первом эксперименте в виртуальной акватории АНПА выполнял движение по окружности с радиусом 100 метров со смещающимся центром. На рис. 6 представлена траектория длиной 33 км, построенная НС при движении аппарата.

Рис. 6. Движение АНПА по окружности со смещающимся центром. Чёрной линией показана рассчитанная траектория АНПА, а серой - истинная.

Сетка имеет размерность 10х10 м

В общей сложности для 11000 изображений в ходе работы НС было обнаружено 94114 связей. При реализации АКТ использовались данные о глубине аппарата, поскольку стереобаза в 0,4 м не способна обеспечить приемлемую точность траектории по высоте. Из рис. 6 видно, что рассчитанная траектория (чёрная линия) близка к истинной (серая линия), а графики ошибки НС представлены на рис. 7а, б.

Рис. 7а показывает график зависимости ошибки НС от пройденного расстояния. Зависимость носит циклический характер, так как ошибка опускается к минимальным значениям каждый раз, когда АНПА проходит недалеко от начала координат. График на рис. 7б отображает зависимость ошибки НС от расстояния до начала координат. Данные усреднены по каждому метру расстояния.

Для второго теста в синтетической акватории была выбрана сложная траектория, чередующая зоны с полным покрытием дна и зоны прямолинейного движения - рис. 8. Траектория получена для 10006 изображений на основе 42985 рассчитанных связей. Как

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

12 15 18 21 24 27 31 км

м

б) 0 50 100 150 200 250 300 350 Рис. 7. Зависимость ошибки НС: от пройденного пути

Рис. 8. Траектория АНПА, чередующая зоны с полным покрытием и прямолинейные участки. Чёрной линией показана рассчитанная траектория, а серой - истинная. Сетка имеет размерность 100x100 м

м 60 50 40 30 20 10

0

27

3 6 9 12 15 18 21 24

Рис. 9. Зависимость ошибки НС от пройденного расстояния

Как в случае с реальными данными, в экспериментах в виртуальной акватории удалось достичь скорости обработки 3 кадра в секунду. При этом для поиска связей отбирались 20 потенциально связанных изображений (М' в подпараграфе 4.3). При моделировании работы инерциальной НС использовались следующие параметры: а1, = (0,15, 0,15, 0,15) м/с и а'г = (0,01, 0,01, 0,01) рад/с. Проведённые эксперименты показали, что глобальный поиск позволяет без труда находить связи между изображениями, даже если это слабоконтрастные, зашумлённые снимки морского дна с однообразными текстурными рисунками. При этом время поиска зависит от количества изображений лишь логарифмически, что позволяет пренебречь этой зависимостью для М < 100000 изображений, получаемых АНПА за миссию.

Заключение

Алгоритм визуальной навигации, описанный в статье, обладает рядом преимуществ по сравнению с алгоритмами навигации, использующими априорные знания о положении АНПА для поиска связей между изображениями. Самым важным преимуществом является устойчивость к участкам дна, не имеющим каких-либо визуальных особенностей. Даже если большая часть обследуемой территории представлена однородной поверхностью либо значительная часть траектории АНПА проходит слишком высоко для фотосъемки, небольшого количества визуально различимых участков будет достаточно для корректировки положения аппарата. Другим преимуществом системы является предсказуемое время работы, так как количество изображений, просматриваемых для установления связей, не зависит от точности, с которой положение АНПА известно в данный момент.

Предложенный алгоритм не требует высокопроизводительных систем для работы. Приемлемая частота обработки в 2-3 кадра в секунду может быть достигнута уже на 4-ядерных процессорах среднего ценового диапазона. Кроме двух параллельных потоков -потока АОСД и потока АКТ, параллельно могут выполняться следующие задачи потока АДИ: расчёт данных облаков, поиск похожих дескрипторов для каждого дескриптора запроса, поиск связей между парой изображений. Таким образом, увеличения скорости работы алгоритма при необходимости можно достичь, увеличив количество процессоров.

Предложенный алгоритм навигации можно рекомендовать к применению на АНПА, оснащённых сте-реокамерой и предназначенных для решения задачи фотообследования морского дна либо объектов, расположенных на этом дне. Глобальный поиск связей в равной степени может быть применен к НС, основанным на графе состояний АНПА, и к НС, использующим EKF или EIF.

Благодарности Работа выполнена при поддержке гранта РФФИ (№ 16-31-00187 мол_а).

Литература

1. Бобков, В.А. О некоторых алгоритмах визуальной навигации автономного необитаемого подводного аппарата с использованием стереоизображений / В.А. Бобков, С.В. Мельман, А.Ю. Толстоногов, А.Ф. Щербатюк // Подводные исследования и робототехника. - 2016. -№ 2(22). - С. 19-24.

2. Eustice, R.M. Visually augmented navigation for autonomous underwater vehicles / R.M. Eustice, O. Pizarro, H. Singh // IEEE Journal of Oceanic Engineering. - 2008. - Vol. 33, Issue 2. - P. 103-122. - DOI: 10.1109/J0E.2008.923547.

3. Salvi, J. Visual SLAM for 3D large-scale seabed acquisition employing underwater vehicles / J. Salvi, Y. Petilot, E. Battle // 2008 IEEE/RSJ International Conference on Intelligent Robots and Systems. - 2008. - P. 1011-1016. -DOI: 10.1109/IR0S.2008.4650627.

4. Aulinas J. Vision-based underwater SLAM for the SPA-RUS AUV / J. Aulinas, Y. Petillot, X. Llado, J. Salvi,

R. Garcia // Proceedings of the 10th International Conference on Computer and IT Applications in the Maritime Industries. - 2011. - P. 171-180.

5. Johnson-Roberson, M. Generation and visualization of large-scale three-dimensional reconstructions from underwater robotic surveys / M. Johnson-Roberson, O. Pizarro, S.B. Williams, I. Mahon // Journal of Field Robotics. -2010. - Vol. 27, Issue 1. - P. 21-51. - DOI: 10.1002/rob.20324.

6. Eustice, R.M. Visually mapping the RMS titanic: Conservative covariance estimates for SLAM information filters / R.M. Eustice, H. Singh J.J. Leonard, M.R. Walter // The International Journal of Robotics Research. - 2006. - Vol. 25, Issue 12. - P. 1223-1242. - DOI: 10.1177/0278364906072512.

7. Kim, A. Real-time visual SLAM for autonomous underwater HULL inspection using visual saliency / A. Kim, R.M. Eustice // IEEE Transaction on Robotics. - 2013. - Vol. 29, Issue 3. -P. 719-733. - DOI: 10.1109/TR0.2012.2235699.

8. Carraso, P.L.N. Stereo Graph-SLAM for autonomous underwater vehicles / P.L.N. Carraso, F. Bonin-Font, G.O Codina. - In: Intelligent Autonomous Systems 13: Proceedings of the 13th International Conference IAS-13 / ed. by E. Menegatti, N. Michael, K. Berns, H. Yamaguchi. - Heidelberg, New York, Dordrecht, London: Springer, 2014. -P. 351-360. - DOI: 10.1007/978-3-319-08338-4_26.

9. Камаев, А.Н. Относительное ориентирование снимков низкой контрастности, полученных с АНПА, на основе точечных особенностей / А.Н. Камаев // Учёные заметки ТОГУ. - 2011. - Т. 2, № 2. - С. 32-42.

10. Lowe, D.G. Distinctive image features from scale-invariant keypoints / D.G. Lowe // International Journal of Computer Vision. - 2004. - Vol. 60, Issue 2. - P. 91-110. - DOI: 10.1023/B:VISI.0000029664.99615.94.

11. Hartley, R. Multiple view geometry in computer vision / R. Hartley, A. Zisserman. - 2nd ed. - Cambridge: Cambridge University Press, 2004. - 672 p. - ISBN: 978-0521-54051-3.

12. Jegou, H. Product quantization for nearest neighbor search / H. Jegou, M. Douze, C. Schmid // IEEE Transactions on Pattern Analysis and Machine Intelligence. - 2011. - Vol. 33, Issue 1. - P. 117-128. - DOI: 10.1109/TPAMI.2010.57.

13. Fischler, M.A. Random sample consensus: A paradigm for model fitting with applications to image analysis and automated cartography / M.A. Fischler, R.C. Bolles // Communications of the ACM. - 1981. - Vol. 24, Issue 6. - P. 381395. - DOI: 10.1145/358669.358692.

14. Kabsch, W. A solution for the best rotation to relate two sets of vectors / W. Kabsch // Acta Crystallographica Section A. - 1976. - Vol. 32, Issue 5. - P. 922-923. - DOI: 10.1107/S0567739476001873.

15. Liu, T. An investigation of practical approximate nearest neighbor algorithms / T. Liu, A.W. Moore, A. Gray, K. Yang // Proceedings of the 17th International Conference on Neural Information Processing Systems. - 2004. - P. 825-832.

16. Beis, J. Shape indexing using approximate nearest-neighbour search in high-dimensional spaces / J. Beis, D.G. Lowe // Proceedings of the 1997 Conference on Computer Vision and Pattern Recognition (CVPR '97). - 1997. - P. 10001006. - DOI: 10.1109/CVPR.1997.609451.

17. Камаев, А.Н. Позиционирование изображений морского дна полученных с помощью АНПА / А.Н. Камаев // Подводные исследования и робототехника. - 2013. -№ 2(16). - С. 38-47.

18. Grisetti, G. A tutorial on graph-based SLAM / G. Grisetti, R. Kummerle, C. Stachniss, W. Burgard // IEEE Intelligent Transportation Systems Magazine. - 2010. - Vol. 2, Issue 4. - P. 31-43. - DOI: 10.1109/MITS.2010.939925.

19. Kamaev, A.N. Constructing and visualizing three-dimensional sea bottom models to test AUV machine vision systems / A.N. Kamaev, V.A. Sukhenko, D.A. Karmanov // Programming and Computer Software. - 2017. - Vol. 43, Issue 3. - P. 184-195. - DOI: 10.1134/S0361768817030070.

20. Kamaev, A.N. Synthetic data for AUV technical vision systems testing / A.N. Kamaev, S.I. Smagin, V.A. Sukhenko, D.A. Karmanov // CEUR Workshop Proceedings. - 2017. -Vol. 1839. - P. 126-140.

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

Камаев Александр Николаевич, кандидат технических наук, старший научный сотрудник Вычислительного центра ДВО РАН, лаборатория численных методов в математической физике, Хабаровск, Россия. Область научных интересов: компьютерное зрение, обработка изображений, компьютерная графика. E-mail: kamaev_an@mail.ru .

Карманов Дмитрий Александрович, старший программист Вычислительного центра ДВО РАН, лаборатория информационно-коммуникационных систем, Хабаровск, Россия. Область научных интересов: компьютерная графика, распределённые вычисления и системы искусственного интеллекта. E-mail: karmanno@gmail. com .

ГРНТИ: 28.23.15. Поступила в редакцию 7 февраля 2018 г.

VISUAL NAVIGATION OF AN AUTONOMOUS UNDERWATER VEHICLE BASED ON THE GLOBAL SEARCH OF IMAGE CORRESPONDENCES

A.N. Kamaev 1, D.A. Karmanov 1

1 Computing Center FEB RAS, Khabarovsk, Russia

Abstract

A task of autonomous underwater vehicle (AUV) navigation is considered in the paper. The images obtained from an onboard stereo camera are used to build point clouds attached to a particular AUV position. Quantized SIFT descriptors of points are stored in a metric tree to organize an effective search procedure using a best bin first approach. Correspondences for a new point cloud are searched in a compact group of point clouds that have the largest number of similar descriptors stored in the tree. The new point cloud can be positioned relative to the other clouds without any prior information about the AUV position and uncertainty of this position. This approach increases the reliability of the AUV navigation system and makes it insensitive to data losses, textureless seafloor regions and long passes without trajectory intersections. Several algorithms are described in the paper: an algorithm of point clouds computation, an algorithm for establishing point clouds correspondence, and an algorithm of building groups of potentially linked point clouds to speedup the global search of correspondences. The general navigation algorithm consisting of three parallel subroutines: image adding, search tree updating, and global optimization is also presented. The proposed navigation system is tested on real and synthetic data. Tests on real data showed that the trajectory can be built even for an image sequence with 60% data losses with successive images that have either small or zero overlap. Tests on synthetic data showed that the constructed trajectory is close to the true one even for long missions. The average speed of image processing by the proposed navigation system is about 3 frames per second with a middle-price desktop CPU.

Keywords: navigation, AUV, SLAM, feature points, dead reckoning, image matching.

Citation: Kamaev AN, Karmanov DA. Visual navigation of an autonomous underwater vehicle based on the global search of image correspondences. Computer Optics 2018; 42(3): 457-467. DOI: 10.18287/2412-6179-2018-42-3-457-467.

Acknowledgements: The work was partially funded by the Russian Foundation for Basic Research under the research project 16-31-00187 mol_a.

References

[1] Bobkov VA, Mel'man SV, Tolstonogov AJu, Shherbatjuk AF. On some algorithms for visual navigation of the autonomous underwater vehicle using stereo images [In Russian]. Underwater Investigations and Robotics 2016; 2(22): 19-24.

[2] Eustice RM, Pizarro O, Singh H. Visually augmented navigation for autonomous underwater vehicles. IEEE Journal of Oceanic Engeneering 2008; 33(2): 103-122. DOI: 10.1109/J0E.2008.923547.

[3] Salvi J, Petillot Y, Battle E. Visual SLAM for 3D large-scale seabed acqui-sition employing underwater vehicles. 2008 IEEE/RSJ International Conference on Intelligent Robots and Systems 2008: 1011-1016. DOI: 10.1109/IROS.2008.4650627.

[4] Aulinas J, Petillot Y, Llado X, Salvi J, Garcia R. Vision-based underwater SLAM for the SPA-RUS AUV. Proceedings of the 10th International Conference on Computer and IT Applications in the Maritime Industries 2011: 171-180.

[5] Johnson-Roberson M, Pizarro O, Williams SB, Mahon I. Generation and visualization of large-scale three-dimensional

reconstructions from underwater robotic surveys. Journal of Field Robotics 2010; 27(1): 21-51. DOI: 10.1002/rob.20324.

[6] Eustice RM, Singh H, Leonard JJ, Walter MR. Visually mapping the RMS titanic: Conservative covariance estimates for SLAM information filters. The International Journal of Robotics Research 2006; 25(12): 1223-1242. DOI: 10.1177/0278364906072512.

[7] Kim A, Eustice RM. Real-time visual SLAM for autonomous underwater HULL inspection using visual saliency. IEEE Transaction on Robotics 2013; 29(3): 719-733. DOI: 10.1109/TR0.2012.2235699.

[8] Carraso PLN, Bonin-Font F, Codina GO. Stereo GraphSLAM for autonomous underwater vehicles. In Book: Me-negatti E, Michael N, Berns K, Yamaguchi H, eds. Intelligent Autonomous Systems 13: Proceedings of the 13th International Conference IAS-13. Heidelberg, New York, Dordrecht, London: Springer; 2014: 351-360. DOI: 10.1007/978-3-319-08338-4_26.

[9] Kamaev AN. Relative features based orientation of low contrast images obtained from AUV [In Russian]. Scientists notes PNU 2011; 2(2): 32-42.

[10] Lowe DG. Distinctive image features from scale-invariant keypoints. International Journal of Computer Vision 2004; 60(2): 91-110. DOI: 10.1023/B:VISI.0000029664.99615.94.

[11] Hartley R, Zisserman A. Multiple view geometry in computer vision. 2nd ed. Cambridge: Cambridge University Press; 2004. DOI: 10.1023/B:VISI.0000029664.99615.94.

[12] Jegou H, Douze M, Schmid C. Product quantization for nearest neighbor search. IEEE Transactions on Pattern Analysis and Machine Intelligence 2011; 33(1): 117-128. DOI: 10.1109/TPAMI.2010.57.

[13] Fischler MA, Bolles RC. Random sample consensus: A paradigm for model fitting with applications to image analysis and automated cartography. Communications of the ACM 1981; 24(6): 381-395. DOI: 10.1145/358669.358692.

[14] Kabsch W. A solution for the best rotation to relate two sets of vectors. Acta Crystallographica Section A 1976; 32(5): 922-923. DOI: 10.1107/S0567739476001873.

[15] Liu T, Moore AW, Gray A, Yang K. An investigation of practical approximate nearest neighbor algorithms. Proceedings of the 17th International Conference on Neural Information Processing Systems 2004: 825-832.

[16] Beis J, Lowe DG. Shape indexing using approximate nearest-neighbour search in high-dimensional spaces. Proceedings of the 1997 Conference on Computer Vision and Pattern Recognition (CVPR '97) 1997: 1000-1006. DOI: 10.1109/CVPR.1997.609451.

[17] Kamaev AN. Computation of positions of sea floor images obtained from AUV [In Russian]. Underwater Investigations and Robotics 2013; 2(16): 38-47.

[18] Grisetti G, Kummerle R, Stachniss C, Burgard W. A Tutorial on graph-based SLAM. IEEE Intelligent Transportation Systems Magazine 2010; 2(4): 31-43. DOI: 10.1109/MITS.2010.939925.

[19] Kamaev AN, Sukhenko VA, Karmanov DA. Constructing and visualizing three-dimensional sea bottom models to test AUV machine vision systems. Program Comput Soft 2017; 43(3): 184-195. DOI: 10.1134/S0361768817030070.

[20] Kamaev AN, Smagin SI, Sukhenko VA, Karmanov DA. Synthetic data for AUV technical vision systems testing. CEUR Workshop Proceedings 2017; 1839: 126-140.

Author's information

Aleksandr Nikolaevich Kamaev, candidate of Technical Sciences, senior research fellow of the Computing Center FEB RAS, Numerical Methods in Mathematical Physics laboratory, Khabarovsk, Russia. His research interests are currently focused on computer vision, image processing and computer graphics. E-mail: kamaev an@mail.ru .

Dmitry Aleksandrovich Karmanov, senior software developer of the Computing Center FEB RAS, Information and Telecommunication Systems laboratory, Khabarovsk, Russia. His research interests are currently focused on computer graphics, distributed computations and artificial intelligence systems. E-mail: karmanno@gmail.com .

Received February 7, 2018.