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

Известия Тульского государственного университета Естественные науки. 2012. Вып. 1. С. 119-129

= ИНФОРМАТИКА

УДК 004.932

Алгоритм поиска квазиоптимальной разметки для обработки изображений с построчным комбинированием переменных

П. А. Мельников, А. В. Копылов

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

Ключевые слова: обработка изображений, оптимизационный подход, (min, +) задачи разметки, динамическое программирование, агрегирование переменных.

Введение

Несмотря на большое разнообразие задач и множество существующих алгоритмов низкоуровневой обработки изображений, оказывается возможным выделить достаточно большие подклассы, допускающие единую математическую постановку. В большинстве случаев исходное изображение представляется как числовой массив, У = (уt,t £ Т), который обычно принимает значения на непрерывной или дискретной оси уровня яркости Уt £ У и определен на дискретном множестве элементов Т = {Ч = (¿1,£2) : 1 ^ ¿1 ^ N1, 1 ^ ¿2 ^ N2}. Результатом анализа данных является вторичный массив X = (xt, t £ Т) определенный на том же множестве элементов t £ Т и принимающий значения xt из множества Х^ специфичного для каждой конкретной задачи. Декартово произведение множеств Xt для всех t £ Т определяет множество X всех возможных значений вторичного массива данных.

При использовании оптимизационного подхода, вторичный массив играет роль аргумента целевой функции .](X\У), оценивающей несоответствие между каждой допустимой версией результата X £ X и исходным массивом данных У. Важно лишь так выбрать класс целевых функций, чтобы было

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

X (Y ) = argmin J (X\Y). (1)

X ex

Алгоритм минимизации целевых функций из выбранного класса и будет играть роль универсального обобщенного алгоритма обработки изображений.

Будем учитывать, что каждый элемент xt, t Е T массива данных, естественным образом связан с рядом соседей. Удобно выражать такое отношение соседства в виде неориентированного графа G, который понимается как множество пар соседних элементов G С T х T.

В качестве универсального носителя всей используемой в дальнейшем информации наблюдения об исходном массиве данных Y в окрестности каждой его точки t е T удобно использовать локальные оценочные функции ^t(xt) = ^t(xt\Y), каждая из которых определена на множестве значений целевой переменной и принимает тем меньшие значения, чем более правдоподобной представляется предположение, что xt как раз и есть искомое значение этой переменной. Функция ^t(xt) соответствует одной целевой переменной xt и, следовательно, одной вершине (узлу) графа соседства G элементов изображения, поэтому будем называть данные функции узловыми.

Для выражения априорных (модельных) предположений о результате обработки используются функции, определяемые модельными представлениями о предпочтительных сочетаниях значений пар соседних целевых переменных. Каждая такая функция задается как функция соответствующих двух соседних переменных Yt’,t"(xt',xtn), (t', t") Е G, возрастающая

при увеличении взаимного рассогласования их значений в том смысле, который диктуется имеющимися модельными представлениями о желаемом результате обработки. Будем называть в дальнейшем функции Yt',t'' (xt' ,xt») функциями связи, т.к. они отождествляются с ребрами графа соседства (t', t'') Е G.

Таким образом, для заданного изображения Y решение о вторичном массиве X в целом следует искать минимизируя как узловые функции, так и функции связи. Естественно принять их сумму

J (X\Y) = ^ ■0t(xt|Yt)+ ^ Yt',t''(xt', xt'') =

teT (t',t'')eG

= Y^ Mxt)+ ¿2 Yt',t'' (xt', xt'') = J (X) (2)

teT (t',t'')eG

в качестве комбинированного критерия принятия решения (1). Заметим, что целевая функция (2), подлежащая минимизации, представляет собой сумму

некоторого числа элементарных функций, каждая из которых зависит не более чем от двух переменных, причем структура смежности переменных задается графом соседства элементов массива. Мы будем называть такие функции парно-сепарабельными с графом смежности переменных G. Число переменных такой целевой функции равно числу вершин графа, т.е. числу элементов массива \Т\.

Для решения оптимизационной задачи (2) в общем виде используются процедуры итерационного случайного поиска, такие как стохастическая релаксация и моделируемый отжиг (Simulated Annealing) [1]. Медленная сходимость такого рода методов значительно ограничивает класс практических задач, для решения которых они могут быть использованы.

В случае конечного множества значений целевых переменных x Е X = {1,...,K}, оптимизационная проблема (2) носит название задачи разметки, известной также как (max, +) или (min, +) задачи, которая для произвольного графа соседства G является NP-полной. Однако, существуют частные случаи, для которых возможно построение достаточно эффективных итерационных алгоритмов [2]: субмодулярные задачи,

сводящиеся к поиску оптимального сечения графа [3]; задачи в которых граф соседства не имеет циклов (является деревом), решение которых основано на принципе динамического программирования [4]; задачи «эквивалентные тривиальным», включающие и первые два типа, решение которых сводится к специфической задаче линейного программирования [5].

Вычислительная эффективность процедуры динамического программирования способствовала созданию приближенных алгоритмов на ее основе для разрешения оптимизационной проблемы (2) с решеткой в качестве графа смежности переменных, среди которых Loopy Belief Propagation (LBP) [6] и Tree-Reweighed Message Passing (TRW) [7]. Данные алгоритмы занимают в настоящее время лидирующие позиции по точности решения и скорости работы, однако у них имеются и существенные недостатки: 1) итерационная процедура, лежащая в их основе, не позволяет точно спрогнозировать время поиска решения; 2) нет гарантии, что решение вообще будет найдено, процедура может зациклиться и не дать приемлемого результата, ведь алгоритмические схемы, лежащие в ее основе, были разработаны для ациклических графов.

Построчное комбинирование целевых переменных

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

Пусть Xtl = (xtlt2 ,¿2 = 1, ■ ■ ■, N2) — совокупности всех переменных, соответствующих узлам решетки в ¿i-ой строке. Тогда все строки решетки естественным образом образуют вертикальную цепочку групповых

переменных Х\,Хъ,... ,Х^1 с парными связями между соседними переменными (рис. 1).

Рис. 1. Последовательная смежность групп переменных в строках

решетки

Запишем критерий (2) для решетчатого графа смежности переменных в виде

3(х1,1 , Х1,2, • • • , Х N1 N2 ) ^ ' ^tlt2 (Х^1^2 ) +

-^2 = 1

+ ^ 2 *2=2 ^ (xtl, t2 —1, Xtlt2 ) + ^ ' *2 = 1 ^ (Х^-1^,2 ,Xtlt2 )

= 3 (Х1, ...,Х^) = (Xtl) + Е^=2 Гtl-1,*1 №1-1,^) (3)

где (Xtl) = Е^2=1 ^*1*3 (хМ2 ) + Е¿2=2 1”(ХН^2-1,Х^2) - узловая

функция ¿1-ой строки, Г^-М1 (Xtl-l,Xtl) = Е ¿=1 7/(xtl-l,t2 ,х^2 ) —

функция связи между ¿1-ой и (¿1 + 1)-ой строками. Очевидно, что

представление критерия в форме (3) абсолютно эквивалентно представлению (2), но использование в последнем переменных Xl, X2, . . . , XN1, образующих цепочечный граф смежности, позволяет построить алгоритм оптимизации данного критерия, основанный на принципе динамического программирования. Сущность процедуры заключается в рекуррентной декомпозиции исходной задачи оптимизации по всем переменным на

последовательность частных задач оптимизации промежуточных функций только одной переменной.

Рассмотрим кратко основные моменты классической процедуры

динамического программирования. Пусть 3^(Xl,...,Xtl) — частичная целевая функция, в точности повторяющая структуру полной функции с тем единственным различием, что суммирование ведется лишь по аргументам, индексы которых не превосходят текущего значения ¿1:

(Xl,...,Xtl) = У \ ^s(Xs) + У*\Гs-l,s(Xs-l,Xs), ¿1 = 1,...,^1.

.¿—/5=1 ^=2

Зафиксируем значение любой переменной Xt1 и мысленно проведем глобальную минимизацию частичной функции (4) по всем переменным, индекс которых меньше t\. Мы получим функцию одной переменной

Btl (Xtl )= min Jt-1 (Xi,...,Xtl), Bi(Xi) = Ф1Х1),

Xi,. . . , Xt1-i

называемую функцией Беллмана по имени создателя динамического программирования.

Функции Беллмана связаны рекуррентным соотношением [9]:

Bt1 (Xt1 )=^t1 (Xt1)+ min [rt1-i,t1 (Xt1-i,Xt1 )+Bt1-i(Xt1-i)] , ti=2,...,^i.

Xt1-i

(5)

Очевидно, что глобальный минимум функции Беллмана от последней переменной min Bn1 (Xn1 ) совпадает с глобальным минимумом полной целевой функции по всем переменным min J(Xi,... ,Xn1 ), поэтому значение последней переменной, минимизирующее соответствующую функцию Беллмана

XN1 = arg min BN1 (Xn), (6)

Xn1

является оптимальным значением этой переменной в составе искомой оптимальной совокупности X.

Остальные элементы искомого решения (Xi,..., X^-i) могут быть найдены путем применения обратного рекуррентного соотношения, представляющего собой обращенную форму прямого рекуррентного соотношения (5):

Xt1-i(Xt1 ) = argmin rt1 -i,t1 (Xt1-i,-Xt1)+ Bt1-i(Xt1-i) . (7) Xt1-1 L J

Таким образом, каковы бы ни были узловые функции и функции связи в составе парно-сепарабельной целевой функции J(X) с цепочечной смежностью переменных, алгоритм динамического программирования находит точку ее глобального минимума, выполняя при этом заранее известное конечное число операций, пропорциональное числу строк решетки

Ni.

Однако в задачах оптимизации (5)-(7) процедуры динамического программирования для дискретного случая xt1t2 Е Xt112 = {1,...,K} требуется выполнить полный перебор всех KN2 комбинаций значений переменных в строке Xt1 = (xt1i,..., xt1N2) Е Пt=i n2 Xt1t2. Избежать полного перебора можно было, если бы функции Беллмана были парно-сепарабельными функциями с последовательной смежностью N2 переменных Xt1 = (xt1i,..., xt1N2). В этом случае, для решения каждой промежуточной задачи оптимизации можно было бы использовать обычную последовательную процедуру динамического программирования, вычислительная сложность которой линейна относительно числа

переменных в строке N2. Обеспечить линейную вычислительную сложность можно лишь заменив истинные функции Беллмана их парно-сепарабельными аппроксимациями, сохраняющими некоторые

основные черты исходных функций, и отказавшись, тем самым, от точного решения задачи оптимизации.

Заметим, что первая функция Беллмана

В1(Х1) = Ф1(Х1) = ^2=1 ^U2 (x1t2 ) +YlN=2 Y'{Xl,t2-1,XU2 )

парно-сепарабельна по определению.

Предположим, что очередная функция Беллмана для (¿1 — 1)-ой строки решетки парно-сепарабельна

Bti — 1(Xti — 1) = Bti-1(xti-1, t2, ¿2 = 1, ■ ■ ■, N2) =

EN2 ^-^,N2

t2_1 ntl-1, t2 (xti — 1, t2 ) + / j t2=2 ^tl-1, t2 (xti-1, t2 — 1, Xti-1, t2 )

тогда функция Беллмана для следующей ¿1-ой строки в соответствии с правилами рекуррентного пересчета (5) будет иметь вид:

Bti (Xti) = (Xti) + min { r(Xti—1,Xti) + Bt(i—Xti—1)} =

Xti-i < J

N2 N2

^ ^ ^tit2 (xtit2 ) + Y (xti, t2 — 1,Xtit2 ) +

t2_1 t2_2

( N2 N2

£ Y (xti — 1, t2 ,Xtit2 ) + ''У ' nti — 1, t2 (xti — 1, t2 ) +

t2=i,. ,N2 l t2_1 t2_1

N2

+ У 2 №ti — 1, t2 (xti — 1, t2 — 1, Xti — 1, t2 ) f ■ t2_2

(8)

Введем отдельное обозначение для операции минимизации, присутствующей в формуле (8)

N2 N2

*ti (Xti )= mint \ У 1 Y(xti — 1,t2 ,Xtit2 )+2^t 1 nti — 1,t2 (xti — 1,t2) +

xti-i, t2 > K. *—^t2_1 *—^t2_1

t2=i,. . . ,N2

+ У 2¿2_2 ^ti — 1,t2 (xti — 1,t2 — 1,xti — 1,t2 ) j1 ,

тогда Би (Хи) = Ф^ (Х^)+ (Хи). Здесь функция Ф^ (Х^) обладает

цепочечной парной сепарабельностью по определению, однако функция Е*1 (Х*1) = Ft1 (ж*1*2, ¿2 = N2) уже не обладает этим свойством,

следовательно, и функция Б^ (Х^) не будет являться парно-сепарабельной.

В работе [8] предлагается эвристическое разрешение описанной ситуации путем замены функции (Х*1) на функцию (Х^), обладающую

свойством цепочечной парной сепарабельности и имеющую одинаковые с исходной Fti (Xti), во-первых, так называемые маргинальные функции каждой переменной

Ft-i_t2 (xtit2 )=Ftit2 (xtit2 ) = min Fti (xtis, s=1, ■ ■ ■ ,N2) , xM2, ¿2=1, ■■■,N2,

xtis: s_ t2

(9)

а также маргинальные функции пар смежных переменных в строке '^ti, (t2 —1, t2)(xti, t2 — 1, Xti t2 )= Fti, (t2 — 1, t2)(xti, t2 —1, Xtit2 );

Fti, (t2 — 1, t2) (xti, t2 — 1, xtit2 ) min Fti (xtis, s 1, ■ ■ ■ , N2) , xti , t2 — 1,

xtis'- s_ t2

xtit2, ¿2 = 2,■■■,N2■ (10)

Нетрудно доказать, что для выполнения требований (9) и (10) достаточно принять

N2 — 1 N2

Fti (xtit2 , ¿2 ^■■■,N2) — У ' Ftit2 (xtit2) + У ' Fti, (t2 — 1,t2)(xti,t2 — 1,xtit2)■

t2 _2 t2 _2

Последовательность маргинальных функций Ftit2 (xtit2), ¿2 = 1,---,N2, и Fti, (t2—1,t2 )(xti,t2—1, xtit2), ¿2 = 2, ■■■,N2, легко вычисляется в каждой строке ¿1 в виде последовательности таблиц, состоящих, соответственно, из K и K2 элементов, через уже найденные ранее функции nti—1, s(xti—1, s), s = 1,^^ N2, и ßti—1, s(xti—1, s—1,xti—1, s), S = 2, ■■■,N2, для предыдущей строки ¿1 — 1 и функции связи J/(xti — 1, s, xtis), s = 1,■■■,N2, с помощью процедуры динамического программирования, имеющей линейную вычислительную сложность относительно числа переменных в строке N2.

Таким образом, очередная приближенная функция Беллмана определяется по правилу:

N2 N2

Bti (Xti ) Bti (xtit2 , ¿2 1, ■ ■ ■ , N2) ^ 2 Vtit2 (xtit2 ) + ^ 2 №it2 (xti, t2 — 1,xtit2 ) ,

t2 _1 t2 _2

где

Vti t2 (xtit2 ) = ^tit2 (xtit2 ), ¿2 = 1,N2,

Vtit2 (xtit2 ) *0tit2 (xti t2 ) Ftit2 (xtit2 ) , ¿2 2, ■ ■ ■ , N2 1,

ßtit2 (xti, t2 — 1, xtit2 ) T (xti, t2 — 1, xtit2 )+F?ti, (t2 — 1, t2) (xti, t2 — 1, xtit2 ) , ¿2 2, ■ ■ ■ , N2 ■

Расчет следует начинать с исходных присвоений

n1t2 (x1, t2 —1,x1t2 ) = Ф1, t2 — 1(x1, t2 — l), ¿2 = 1, ■ ■ ■ , N2,

№1t2 (x1, t2 — 1, x1t2 ) 1 (x1,t2 — 1,x1t2 ), ¿2 ^, ■ ■ ■ , N2 ■

Оптимальные значения переменных Xtl = (xt1t2, ¿2 = !,•••,N2) определяются путем вычисления выражений (6) и (7) с помощью алгоритма динамического программирования.

Эксперименты

Экспериментальная проверка точности разработанного алгоритма проводилась на задачах текстурной сегментации, для которых известно точное решение (200x200 с двумя типами областей (рис. 2, а) и 300x300 с тремя типами областей (рис. 2, б)). Сравнение быстродействия и точности работы было выполнено среди алгоритмов, использующих древовидные структуры для оптимизации парно-сепарабельных целевых функций с графом соседства в виде решетки, а именно: последовательный алгоритм Tree-reweighted message passing (TWRS) [7] и алгоритм циклического распространении доверия — Max Product Belief Propagation (MaxProdBP). Алгоритмы поиска минимального сечения графа не применимы к задачам данного типа.

На рис. 3 показаны результаты сравнения алгоритмов на примере решения двух задач текстурной сегментации. Как видно из рисунка, для данных задач алгоритм на основе агрегирования переменных нашел точное решение за наименьшее время.

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

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

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

___и ребра имеют нулевое значение функции связи

Рис. 4. Пример «сложной» (шш,+) задачи

Все четыре возможные значения целевых переменных изображены маленькими окружностями в узлах графа. Узловые функции ф; (х;) равны нулю для всех XI Е {1, 2, 3, 4}. Значения соответствующих функций связи 11';'' (х; ,х;") представлены ребрам графа. Легко заметить, что такой

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

Заключение

Представленная в работе версия алгоритма на основе построчного комбинирования целевых переменных отличается от оригинальной [8] применением процедуры динамического программирования «вперед и обратно» (классическая схема) вместо процедуры «вперед и навстречу»

[10]. Использованная версия процедуры динамического программирования в сочетании с аппроксимациями функций Беллмана позволила решать задачи, которые алгоритмы на основе древовидной декомпозиции и алгоритмы поиска минимального сечения графа решить не в состоянии. Для других классов задач данный алгоритм находит достаточно точное решение за приемлемое время. Основной недостаток алгоритма на основе построчного комбинирования целевых переменных — высокое требование к объему оперативной памяти.

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

1. Geman S., Geman D. Stochastic Relaxation, Gibbs Distributions, and the Bayesian Restoration of Images // IEEE Trans. on PAMI. 1984. V.6. №6. P.721-741.

2. Schlesinger M.I., Flach B. Some solvable subclasses of structural recognition problems // Proceedings of Czech Pattern Recognition Workshop, Praha. 2000. P.55-62.

3. Boykov Y, Veksler O., Zabih R. Fast approximate energy minimization via graph cuts // PAMI. 2001. V.23. №11. P.1222-1239.

4. Optimization techniques on pixel neighborhood graphs for image processing / V. Mottl [et al.] // Graph-Based Representations in Pattern Recognition. Computing, Supplement 12. Wien: Springer-Verlag, 1998. P.135-145.

5. Шлезингер М.И., Гигиняк В.В. Решение (max,+)-задач структурного распознавания с помощью их эквивалентных преобразований // УСиМ. 2007. №1. С.3-15.

6. Felzenszwalb P.F., Huttenlocher D.P. Efficient belief propagation for early vision // Conference version from IEEE CVPR. 2004. V.1. P.261-268.

7. Kolmogorov V. Convergent tree-reweighted message passing for energy minimization // Pattern Analysis and Machine Intelligence. 2006. V.28. № 10. P.1568-1583.

8. Dmitriev D.A., Mottl V.V., Kopylov A.V. Algorithms of Approximate Pairwise Separable Optimization for Image Processing // Pattern Recognition and Image Analysis. 2003. V.13. №1. P.90-94.

9. Muchnik I.B., Mottl V.V. Bellman Functions on Trees for Segmentation, Generalized Smoothing, Matching and Multi-Alignment in Massive Data Sets. DIMACS Technical Report 8-15 / Rutgers University, USA. 1998. P.63.

10. Алгоритм динамического программирования для анализа нестационарных сигналов / А.А. Костин [и др.] // Журнал вычисл. матем. и матем. физики. 2004. Т.44. №1. С.70-86.

Мельников Петр Александрович (petrmelnikov@mail.ru), ведущий инженер-программист, отдел информационных систем, ОАО «АК «Туламашзавод», Тула.

Копылов Андрей Валерьевич (kopylov@uic.tula.ru), к.т.н., доцент, кафедра автоматики и телемеханики, Тульский государственный университет.

An algorithm of the quasioptimal labeling search for image processing with rowwise aggregation of variables

P. A. Melnikov, A. V. Kopylov

Abstract. Many problems of image processing leads to the specific discrete optimization problems, known as (min, +) labeling problems, that are NP-hard in general formulation. This paper describes the algorithm that allows to solve some special cases of (min, +) problems, unsolvable by methods based on acyclic decomposition of the image lattice.

Keywords: image processing, optimization approach, (min, +) labeling problems, dynamic programming, row-wise aggregation of objective variables.

Melnikov Peter (petrmelnikov@mail.ru), lead part-programming engineer, department of information systems, OAO «AK «Tulamashzavod», Tula.

Kopylov Andrey (kopylov@uic.tula.ru), candidate of technical sciences, associate professor, department of automatics and telemechanics, Tula State University.

Поступила 25.12.2011