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

ШБН 2079-3316 ПРОГРАММНЫЕ СИСТЕМЫ: ТЕОРИЯ И ПРИЛОЖЕНИЯ № 4(8), 2011, с. 3-15 УДК 517.977

А.А.Ардентов, Ю.Л.Сачков

Антропоморфное восстановление поврежденных изображений на основе методов субримановой

геометрии

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

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

1. Введение

Задача восстановления поврежденных или скрытых изображений является одной из актуальных проблем компьютерной графики. Для ее решения разработан ряд методов, многие из которых используют технику вариационного исчисления и оптимального управления [1—3]. Данная работа основана на положениях одного из новых направлений нейрофизиологии — нейрогеометрии зрения [4,5], а также на недавних результатах по субримановой геометрии [6-8]. На основе результатов этих исследований разработан алгоритм для восстановления серии изофот на поврежденных штриховых изображениях.

Работа поддержана Госконтрактом № 07.514.11.4033 по ФЦП «Исследования и разработки по приоритетным направлениям развития научно-технологического комплекса России на 2007—2013 годы», а также Российским фондом фундаментальных исследований, проект 09-01-00246-а.

© А.А.Ардентов, Ю.Л.Сачков, 2011 © Программные системы: теория и приложения, 2011

2. Постановка задачи восстановления изображения и метод ее решения

Задача восстановления изображения может быть формализована следующим образом. Дана прямоугольная область D С R2, взаимно непересекающиеся подобласти 0i, ...,0n С D, и функция f : D\(UN 1 Oi) ^ [0,1]. Требуется восстановить функцию f в областях

0i,...,0n. Здесь D — область исходного изображения, 0i —подобласти с поврежденным изображением; функция f задает доступное наблюдателю изображение (для полутонового изображения яркость, а для штрихового изображения это функция, линии уровня которой совпадают с кривыми, составляющими изображение). Восстановление изображения в областях 0i должно быть антропоморфным, т.е. естественным для глаза человека.

Согласно результатам нейрогеометрии зрения [4,5], человеческий мозг восстанавливает изофоты (x(t), y(t)), t € [tо, ti], на поврежденных изображениях на основе вариационного принципа

(1) f X2 + у2 + а2 в2 dt ^ min, а > 0,

Jta

где 9 = arctg(y/X). Таким образом, минимизируется интегральный компромисс между линейной (x, ) и угловой скоростями восстанавливающей кривой. Параметр а задает масштаб на изображении, и должен подбираться в зависимости от изображения. В данной работе вариационный принцип (1) положен в основу восстановления поврежденных изображений.

При математическом моделировании использовались следующие допущения:

• изображение может быть представлено как портрет линий уровня гладкой функции f : D ^ [0,1],

• функция не имеет критических точек в областях повреждения 0i, ..., 0N С D,

• известна информация о точках пересечения линий уровня функции f с границами областей 0i.

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

Замена переменных

(2) x = ах, у = ау, в = в

преобразует задачу (1) к случаю а = 1. Поэтому при математическом рассмотрении задачи будем считать а = 1, а в алгоритме восстановления изображений параметр а будет подбираться с целью устранения точек возврата на восстанавливающих кривых. Описанный вариационный метод восстановления изображений рассмотрен в работе [9]. В данной работе представлены следующие новые результаты:

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

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

в) детально описан алгоритм восстановления изофот.

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

Задача минимизации функционала f \Jх2 + у2 + в2 dt для кривых (x(t),y(t),6(t)) с закрепленными концами может быть сформулирована как следующая задача оптимального управления:

(3) х = u1 cos в, у = и1 sin в, в = и2,

(4) q = (х, у, в) Є М = М2 у х Sfj, и = (и\, U2) Є R2,

(5) 9(0) = qo = (0, 0, 0), q(t\) = q\ = (х1,у1,в1),

Пространство состояний М = М2 у х Sg естественным образом отождествляется с группой SE(2) движений двумерной плоскости, сохраняющих ориентацию. Эта группа представляется матрицами 3 х 3 следующим образом:

{/ cos в — sin в х \ Ї

I sin в cos в у І | в Є S1 = М/(2-їїХ), х,у Є М > .

Задача (3)-(6) естественно переформулируется как левоинвариантная субриманова задача на группе Ли SE(2) [10,11]. Рассмотрим неин-тегрируемую левоинвариантную субриманову структуру ранга 2 на SE(2), т. е. неинтегрируемое левоинвариантное распределение Д ранга 2 на SE(2) вместе с левоинвариантым скалярным произведением

Jи\ + м2 dt ^ min .

(•, •} на Д. Легко показать, что такая структура единственна, с точностью до скалярного множителя в скалярном произведении. Рассмотрим следующую модель такой субримановой структуры:

(где Ец есть матрица 3 х 3 с единственным единичным элементом в г-ой строке и j-ом столбце, и нулевыми остальными элементами) и соответствующую задачу оптимального управления:

д = иіХі(д) + «2^2(9), Я Є М = БЕ(2), и = (мі,И2) Є К2,

^(0) = до = И, д(іі) = ді,

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

Дч = врап(Хі(д),Х2(д)), (Хі,Хі} = Зіі, і, і = 1, 2,

^і(я) = 4Е13, Х2(д) = д(Е2і — Е12), д Є ЯЕ(2),

X

Рис. 1. Постановка задачи (3)-(6)

В работе [6] показано, что для любой конечной точки д\ € БЕ(2) задача (3)-(6) имеет оптимальное решение.

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

В этом разделе мы опишем некоторые результаты работ [6-8], которые позволяют свести задачу (3)-(6) к решению систем уравнений в функциях Якоби.

Из принципа максимума Понтрягина следует, что экстремальные траектории в задаче (3)-(6) параметризуются точками фазового цилиндра С = (25'1) х Rc математического маятника

(7) 7 = с, с = — sin

Семейство параметризованных длиной дуги нормальных экстремальных кривых в задаче (3) -(6) описывается экспоненциальным отображением

Exp : N ^ М, N = С х R+,

Exp(^) = Exp(A, t) = q(t), v = (A, t) = (7, c, t) G N.

с2

Уравнение маятника (7) имеет интеграл энергии Е = —— cos 7 G

[—1, +то). Рассмотрим следующее разбиение цилиндра С на взаимно

непересекающиеся инвариантные подмножества фазового цилиндра маятника:

5

(8) С = U Си

г= 1

CÍ = = {А є С Е є (-1,1)},

С2 = = {А є С Е є (1, +то)},

Сз = = {А є С Е = 1, с = 0},

04 = = {А є с Е = -1} = {(7, с) є С | 7 = 2пп, с = 0}, п Є N,

С5 = = {А є с Е = 1, с = 0} = {(7, с) Є С | 7 = я + 2пп, с = 0}.

В работе [6] были введены эллиптические координаты (р, к) в области Сі и С2 и Сз цилиндра С, где к — перепараметризованная энергия, а ¡р — время движения маятника (7). В эллиптических координатах поток маятника (7) выпрямляется: ф =1, к = 0, и получена следующая параметризация экстремальных траекторий. Далее используются функции Якоби аіп(<р,к), сп(<р,к), вп(р,к), Ап(р,к), Е(<р,к); кроме того, К (к) есть полный эллиптический интеграл первого рода [13].

Если А = (р, к) € Ci, то ipt = ^ + t и cos = cn ^ cn ipt + sn ip sn pt, sin 6t = si (sn p cn pt — cn p sn <fit),

(9) 6t = si (am p — am pt) (mod 2^),

(10) xt = (si/fc)[cn ^(dn ^ — dn ft) +sn ^(t + E(y>) — E(y>t))],

(11) yt = (1/^)[sn ^(dn p — dn ft) — cn y>(i + E(y>) — E(y>t))], s1 = sgncos(7/2).

В области удобно использовать координату ф = p/k, t^t = pt/k = 'ф + t/k. Если А € С2, то

(12) cos 6t = к2 sn ф sn фt + dn ф dn ф-t,

(13) sin 6t = k(sn ф dn ф-t, — dn ф sn ф^,

(14) xt = s2fc[dn ^(cn ф — cn ф^ + sn ф(Ъ/к + E(^) — E(^t))],

(15) yt = S2[k2 snф(cnф — cnф^ — dnф(Ь/к + E(^) — E(^t))],

S2 = sgn c.

Если A € U¿=3^i, то экстремальные траектории параметризуются элементарными функциями [6].

Рассмотрим следующее разбиение пространства М = SE(2) = Re у х , определяемое значениями функций Ri = у cos 2 — х sin R2 = х cos 2 + y sin 2:

M = {q € M | R1(q)R2(q) sin в = 0}, M = Llf=1M¿,

M' = {g € M | R1(q)R2(q) sin0 = 0},

где каждое из множеств Mi определяется знаками функций sin в, Ri, R2, описанными в Таблице 1.

Mi M1 М2 м3 MA M5 M6 M7 Ms

sgn(sin в) - - - - + + + +

sgn(ñi) + + - - - - + +

sgn(ñ2) + - - + + - - +

Таблица 1. Определение областей Мі

В работах [7,8] было получено глобальное описание времени разреза

tcutW = supjti > 0 | qs оптимальна при s G [0,i 1]}, A G С,

вдоль экстремальных траекторий. В работе [8] были определены множества

М = М \{qo}, N = {(A,i) G N 11 < tcut(A)}.

Так как для любого qi G M существует оптимальное управление, то отображение Exp : IV ^ М сюрьективно; так как оно имеет кратные точки (Максвелла), то отображение не является инъективным. Далее, в работах [7,8] были определены такие плотные подмножества N С N и их разбиение IV = L\8=iDi, что экспоненциальное отображение имеет следующую глобальную структуру:

(16) Exp : IV ^ М и все Exp : Di ^ Mi —диффеоморфизмы.

Поэтому поиск оптимальной кривой в задаче (3)-(6) для конечной точки qi G Mi сводится к решению системы уравнений Exp(A, t) = qi, имеющей единственный корень (\,t) G Di.

Итак, в случае общего положения qi G М решение задачи (3)-(6) сведено к решению систем уравнений

xt(X) = xi, yt(X) = yi, 6t(\) = 6i, (X,t) G Di,

где функции xt, yt, Ot, —компоненты оптимальных траекторий qt, описанные формулами (9)—(15), а (xi,yi, 9i) = qi —конечная точка искомой траектории.

5. Алгоритм восстановления изображения

В этом разделе описывается алгоритм GlobalSolve восстановления изображения на основе вариационного принципа (1), базирующийся на результатах, описанных в предыдущем разделе.

Алгоритм GlobalSolve состоит из следующих подалгоритмов:

• RestoreDomain восстанавливает поврежденные изофоты в подобласти Ok, к = 1, ..., N, см. подраздел 5.4.

• RemoveCusp удаляет точки возврата у оптимальной траектории (x(t),y(t)) при помощи соответствующего подбора значения параметра а (1), см. подраздел. 5.3.

• FindRoot по граничным условиям qi = (xi,yi,9i) (5) вычисляет параметры, определяющие соответствующую поврежденную изофоту (x(t),y(t)), см. подраздел 5.2.

• Solver решает систему из трех алгебраических уравнений с тремя неизвестными

(17) Exp(^)= qi, qi € Mi, v € Di,

где v € Di П Cj, j € {1, 2}, см. подраздел. 5.1.

5.1. Алгоритм Solver решения системы уравнений (17)

Вход: qi € Mi, initial € {true, false}, j € {1,2}, i € {1,..., 4}, v € Di Pi Cj.

Выход: v € Di P Cj.

Выполняемые действия: Алгоритм вычисляет значение корня v € Di P Cj системы (17) с заданной точностью е. Переменная v является одновременно входным параметром (начальное приближение для корня при условии, что initial = true) и выходным параметром (искомый корень).

Константы алгоритма: maxiteration, maxiterrnd € N, £ > 0.

Шаги алгоритма:

(1) Если initial = true, то выполняется переход к шагу 3.

(2) В подобласти Di P Cj случайно выбирается начальное приближение V.

(3) Из начальной точки v запускается итерационный алгоритм поиска приближенного значения корня системы (17). Если число итераций этого алгоритма превышает maxiteration, то выполняется переход к шагу 2. Если v € Di P Cj, то выполняется переход к шагу 2. Если общее число итераций алгоритма Solver превышает maxiterrnd, то алгоритм завершается (корень не найден).

(4) Если | Exp(^) — qil < £, то алгоритм завершается и возвращается значение корня v. Иначе выполняется переход к шагу 2.

5.2. Алгоритм FindRoot вычисления оптимальной траектории задачи (3)—(6)

Вход: qi € М, initial € {true, false}, j € {1, 2}, v € Cj.

Выход: j € {1, 2}, v € Cj.

Выполняемые действия: Алгоритм численно находит корень v системы (17) и номер j такой области Cj, что v € Cj. Параметр а (2) приравнивается 1. Переменные v и j являются одновременно входными параметрами (начальное приближение корня и номер области, содержащей этот корень при условии initial = true) и выходными параметрами (искомый корень).

Шаги алгоритма:

(1) С помощью Таблицы 1 вычисляется номер г такой области Mi, что qi € Mi.

(2) Если initial = true, то выполняется переход к шагу (j + 2), т. е. 3 или 4.

(3) Запускается алгоритм Solver с параметрами qi, initial, i, j = 1, v. Если корень найден, то алгоритм завершается успешно, иначе initial := false.

(4) Запускается алгоритм Solver с параметрами qi, initial, i, j = 2, v. Если корень найден, то алгоритм завершается успешно, иначе initial := false и выполняется переход к шагу 3.

В итоге для корня V = (A,ti) системы (17) вычисляется оптимальная траектория q(t) = Exp(A,í), t € [0, íi], задачи (3)—(6) . Однако оптимальная траектория задачи (3)—(6) может иметь точки возврата (см. рис. 1), которые непригодны при восстановления изофот поврежденного изображения. Такие точки возврата в большинстве случаев могут быть удалены подбором подходящего значения параметра а (2), см. подраздел 5.3.

5.3. Алгоритм RemoveCusp удаления точек возврата на изофотах

Вход: qi = (xi,yi,0i) € М. Выход: а > 0, j € {1, 2}, V € Cj.

Выполняемые действия: Алгоритм находит такое значение а > 0, что соответствующая оптимальная траектория (x(t), y(t)) для задачи (3)—(5) с функционалом качества

(18) Iа = j \jv2 + а2 и\ dt ^ min

не имеет точек возврата. С помощью замены координат (2) новая задача (3)—(5), (18) сводится к изначальной (3)—(6).

Константы алгоритма: Да = 0.1,

06^ - 19l -ЭТ| +°.? 1ж1+10*.2| + 02 ainit = 0.6------------,------------+ 0.2.

0.05 +^Jx2 +у'2

Шаги алгоритма:

(1) Задаются начальные значения а := alnlt, initial := false.

(2) (xa, ya) := (a.xi,ayi).

(3) Запускается алгоритм FindRoot с параметрами (xa, уа,вi),

initial, , .

(4) Если корень v, вычисленный с помощью FindRoot, соответствует кривой (x(t), y(t)) без точек возврата (т.е. X2(t) + y2(t) = 0), то алгоритм завершается успешно.

(5) initial := true, а := а + Да и выполняется переход к шагу 2.

Примеры восстановленных изображений с точками возврата и без точек возврата приведены на рис. 2 и 3 соответственно. Эти изображения получены с помощью алгоритмов FindRoot и RemoveCusp.

В этом алгоритме используется свойство диффеоморфности (16) отображения Exp: близкие точки прообраза экспоненциального отображения переходят под действием Exp в близкие точки его образа. Это свойство используется также далее в алгоритме RestoreDomain, см. подраздел 5.4.

5.4. Алгоритм RestoreDomain восстановления одной поврежденной области

Вход: к Є {1, .. ., N}.

Выход: текстовый файл outputs<k>.

Выполняемые действия: Номер к соответствует входному файлу inputs<k>, который содержит координаты конечных точек всех

Рис. 2. Кривые с Рис. 3. Кривые

точками возврата без точек возврата

изофот в поврежденной области Ok С D. Алгоритм находит все параметры для восстановления этих изофот и записывает их в файл outputs<k>.

Шаги алгоритма:

(1) initial := false.

(2) Из файла inputs<k> считываются входные данные с координатами точек qi = ) для текущей изофоты. Алгоритм за-

вершается, если файл пустой.

(3) Если initial = true, то а0 := а + sign(a — ainit)Аа и запускается алгоритм FindRoot с параметрами (а0Х\,а0у\,0\), initial, v, j. Если корень v, вычисленный с помощью алгоритма FindRoot, соответствует кривой без точек возврата, то искомый корень найден, выполняется переход к шагу 5. Иначе полагается а0 := «init+sign(ainit —а)Аа и запускается алгоритм FindRoot с параметрами (а0х\, а0у\,в\), initial, v, j. Если корень v, найденный с помощью FindRoot, соответствует кривой без точек возврата, то искомый корень найден, выполняется переход к шагу 5.

(4) Запускается алгоритм RemoveCusp с входным значением q\.

(5) Параметры а, и, j записываются в файл outputs<k>, выполняется переход к шагу 2.

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

6. Заключение

В работе описан метод восстановления изофот штрихового изображения с помощью оптимальных кривых в левоинвариантной суб-римановой задаче на группе движений плоскости SE(2). Метод и алгоритм, представленные в этой статье, были реализованы в параллельном программном комплексе OptimalInpainting. Описание этого программного комплекса и демонстрация эффективности подхода будут приведены в последующей работе.

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

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

[1] Chan T. F. Kang S. H. Shen J. Euler’s elastica and curvature based inpainting // SIAM Journal of Applied Math. 2002. 63, no. 2, p. 564—592 f1

[2] Citti G. Sarti A. A cortical based model of perceptual completion in the roto-translation space // J. Math. Imaging Vision, 2006. 24, no. 3, p. 307-326 f

[3] Duits R. Franken E. M. Left-invariant parabolic Evolutions on SE(2) and Contour Enhancement via Invertible Orientation Scores. Part I: Linear Left-invariant Diffusion Equations on SE(2) // Quarterly of Applied Mathematics, June 2008. 68, no. 2, p.255-292 f1

[4] Petitot J. The neurogeometry of pinwheels as a sub-Riemannian contact structure // J. Physiology —Paris, 2003, no. 97, p. 265-309 f1, 2

[5] Petitot J. Neurogeometrie de la vision. Modeles mathematiques et physiques des architectures fonctionelles : Editions de l’Ecole Polytechnique, 2008. f 1, 2

[6] Moiseev I. Sachkov Yu. L. Maxwell strata in sub-Riemannian problem on the group of motions of a plane // ESAIM: Control, Optimisation and Calculus of Variations, 29 July 2010. 16, p. 380-399 f1, 3, 4, 4, 4

[7] Sachkov Yu.L. Conjugate and cut time in the sub-Riemannian problem on the group of motions of a plane // ESAIM: Control, Optimisation and Calculus of Variations, 2010. 16, p. 1018-1039 f4

[8] Sachkov Yu.L. Cut locus and optimal synthesis in the sub-Riemannian problem on the group of motions of a plane // ESAIM: Control, Optimisation and Calculus of Variations, 2011, no. 17, p. 293-321 f1, 4, 4

[9] Сачков Ю. Л., Ардентов А. А., Маштаков А. П. Параллельный алгоритм и программа восстановления изофот для поврежденных изображений // Программные системы: теория и приложения, 2010. 1, № 1, с. 3—20 ^2

[10] Аграчев А. А., Сачков Ю. Л. Геометрическая теория управления. М. : Физ-матлит, 2005. ^3

[11] Montgomery R. A Tour of Subriemannian Geometries, Their Geodesics and Applications : American Mathematical Society, 2002. ^3

[12] Laumond J.P. Nonholonomic motion planning for mobile robots. Lecture notes in Control and Information Sciences, Vol. 256 : Springer, 1998. ^3

[13] Whittaker E.T. Watson G. N. A Course of Modern Analysis. An introduction to the general theory of infinite processes and of analytic functions; with an account of principal transcendental functions. Cambridge : Cambridge University Press, 1996. t4

A. A. Ardentov, Yu.L. Sachkov. Antropomorphic recovery of corrupted images via methods of sub-Riemannian geometry.

Abstract. According to neurogeometry of vision, human brain recovers isophotes of corrupted images via optimal curves for a certain optimal control problem—sub-Riemannian problem on the group of motions of the lane. The paper describes an algorithm for antropo-morphic inpainting of corrupted binary images on the basis of optimal synthesis for this problem.

Key Words and Phrases: image inpainting, optimal control, sub-Riemannian geometry.

Образец ссылки на статью:

А.А.Ардентов, Ю.Л.Сачков. Антропоморфное восстановление поврежденных изображений на основе методов субримановой геометрии // Программные системы: теория и приложения : электрон. научн. журн. 2011. № 4(8), с. 3-15. URL: http://psta.psiras.ru/read/psta2011_4_3-15. pdf