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

Известия Тульского государственного университета Естественные науки. 2015. Вып. 3. С. 105-116

= Математика =

УДК 519.219.2

Анализ видов парно-потенциальных функций для байесовского восстановления

изображений

К. Т. Фам, А. В. Копылов

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

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

Введение

В настоящее время байесовская оценка является важным подходом по крайней мере в случаях фильтрации и восстановления изображений. В современных литературных источниках можно найти много методов шумоподавления и восстановления изображений с сохранением границ на основе этой оценки [1-9]. Основная идея этих методов заключается в поиске максимальной апостериорной (Maximum a posteriori, MAP) оценки истинных изображений с использованием марковского случайного поля (Markov random field, MRF) в рамках байесовского подхода [2, 3].

В рамках данного подхода, задача восстановления изображений может быть выражена как задача оценивания скрытой марковской компоненты двухкомпонентного MRF, роль наблюдаемой компоненты которого играет анализируемое изображение. Использование в качестве функции потерь, сингулярной функции, штрафующей сам факт отличия полученной оценки скрытого поля от «истинного» значения приводит к оптимизационной за-

даче поиска MAP оценки скрытой компоненты, относительно наблюдения [10-11]. Эквивалентная форма представления марковских случайных полей в виде гиббсовских случайных полей согласно теореме Хаммерсли-Клифорда [12], позволяет использовать для задания априорных вероятностных свойств скрытого марковского поля не переходные плотности распределения, а так называемые гиббсовские потенциалы на кликах [6].

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

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

В данной работе предлагаются альтернативные способы задания парных потенциалов. Изображение и результат его обработки рассматриваются как, соответственно, наблюдаемая и скрытая компонента двухкомпонентного случайного поля, а оператор оценивания скрытой компоненты на основе наблюдения строится исходя из принципа минимизации среднего риска ошибки для заданной функции потерь [17]. Разработанные процедуры основаны на этих способах задания априорных вероятностных свойств скрытого марковского поля при помощи функций специальных типов, определенных на кликах графа смежности [10-11, 18] и параметрической процедуре оптимизации на основе динамического программирования [19].

Байесовский подход для шумоподавления изображений

В рамках данного подхода задача обработки изображений может быть выражена как задача оценивания скрытой марковской компоненты X = = (xt,t € T) двухкомпонентного случайного поля T = (t = {tt,t2},tl = = l..Ni,t2 = I..N2), роль наблюдаемой компоненты Y = (yt,t € T) которого играет анализируемое изображение. Использование в качестве функции потерь, сингулярной функции, штрафующей сам факт отличия полученной оценки скрытого поля от истинного значения, приводит к оптимизационной задаче поиска MAP скрытой компоненты относительно наблюдения. Эквивалентная форма представления марковских случайных полей в виде гиббсовских случайных полей согласно теореме Хаммерсли-Клифорда [12] позволяет использовать для задания априорных вероятностных свойств скрытого марковского поля не переходные плотности распределения, а так называемые гиббсовские потенциалы на кликах [6]. Задача поиска максимума апостериорного распределения скрытой компоненты X с учетом данных

наблюдения Y может быть поставлена формально как задача минимизации действительнозначной парно-сепарабельной целевой функции (1), часто называемой функцией гиббсовской энергии [11, 16]:

ГX(Y) = arg min J(X\Y),

)J(X\Y) = E Mxt\Y)+ E Yt',t''X ,xt»). (1)

^ t&T (t',t'' )ec

Структура целевой функции (1) отражает свойство упорядоченности анализируемых данных и определяется при помощи неориентированного графа смежности G С T х T (рис. 1), имеющего вид решетки. Целевая функция (1) при этом представлена как сумма двух типов потенциальных функций на кликах, называемых узловыми функциями и функциями связи. В задаче сглаживания узловые функции фt(xt\Yt) играют роль штрафа на различие значений исходных данных Y и результирующей функции X, и обычно выбираются квадратичными. Функции связи Yt',t''(xt' ,xt») обычно называются парно-потенциальными функциями. Каждая функция связи Yt',t''(xt',xt»), налагает штраф на негладкость результата обработки на соответствующем ребре (t', t'') графа смежности G и может иметь различные формы представления. Так квадратичные функции связи соответствуют предположению о нормальном априорном распределении X.

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

Обобщенная процедура динамического программирования

Для сингулярной функции потерь [17], штрафующей лишь сам факт несовпадения оценки с истиной реализацией скрытого процесса, оценка основывается на максимизации апостериорной плотности распределения скрытого процесса и тем самым сводится к оптимизационной задаче, эффективно решаемо на основе принципа динамического программирования. Вычислительная процедура при этом состоит в рекуррентной декомпозиции исходной задачи оптимизации функции многих переменных на последовательность элементарных подзадач оптимизации функции лишь одной переменной. Элементарные функции одной переменной, подлежащие оптимизации на каждом шаге процедуры минимизации целевой функции, играют ту же роль, что и функции Беллмана в классическом динамическом программировании. Мы будем называть их расширенными функциями Беллмана. Численная

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

Основной принцип процедуры динамического программирования заключается в рекуррентной декомпозиции исходной задачи минимизации функции многих переменных в последовательности частных задач, каждая из которых состоит в минимизации функции только одной переменной. При анализе данных такая процедура выполняет оптимизацию целевой функции (1) за два прохода — в прямом и обратном направлении [18-20].

На прямом ходе t = l, 2,.., N соответствии с прямым рекуррентным соотношением определяются функции Беллмана (2) [10-11]:

Jt(xt) = ^t(xt) +min [yt(xt-i,xt) + Jt-i(xt-i)] . (2)

xt-1

Глобальный минимум функции Беллмана от последней переменной (3) minXN Jn (xn) совпадает с глобальным минимумом полной целевой функции по всем переменным min J (X), поэтому оптимальное значение последней переменной (3) является минимумом соответствующей функции Беллмана:

xN = argmin JN (xN). (3)

xN

Остальные переменные могут быть найдены путём применения обратного рекуррентного соотношения (4) [10-11]:

xt-i(xt) = arg min [jt(x—i,xt) + Jt-i(xt-i)] , t = N - l,.., 2, l. (4)

xt-1

Способы задания парных потенциалов на кликах в виде функции двух соседних переменных

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

В данной работе представлены собой альтернативные способы для шумоподавления изображений в рамках байесовского подхода с использованием других видов парно-потенциальных функций (табл. 1), таких как функция минимума набора квадратичных функций (Multi-quadratic function, MQF)

[16, 18], функция Хубера (Huber function, HF) [14], полуфункция Хубера (semi-Huber, SH) [15], обобщенная гауссовская функция (generalized gaussian function, GGF) [7], функция Блейка-Зиссермана (Blake-Zisserman, BZ) [8], функция Бешага (Besag function, BF) [5].

Таблица 1

Виды парно-потенциальных функций

Название Вид функции Yt(xt-i, xt)

MQF umin [Ai(xt - x—i)2,..., \b(xt — xt-1 )2]

HF ug(x), ((xt — xt-i)2, если x ^ Д g(X) [2Д \xt — xt-1\— Д2, если x> Д

SH u ? 4(%^1)2 — !)

GGF u \xt — xt-i\p , 1 ^ p < 2

BZ u min [(xt — xt-1)2, Д2]

BF u \xt — Xt-i\

Функция Блейка-Зиссермана является частным случаем функции минимума набора квадратичных функций в случае Ь = 2 и вторая функция имеет вид постоянной функции. Функция Бешага является частным случаем обобщенной гауссовской функции при р = 1.

В данной работе узловые функции фг(хг выбираются квадратичными фг(%г\!г) = (жг — уг)2. Функции связи выбраны как минимум из конечного набора квадратичных функций (в случаях видов MQF и BZ) и, то функции Беллмана (2) на каждом шаге динамического программирования также представлены как минимум из конечного набора квадратичных функций, что приводит к так называемой мультиквадратичной процедуре динамического программирования. Это процедура была предложена в работах [16, 18].

Для остальных случаев видов функций связи (HF, БЫ, GGF, BF) функции Беллмана (2) выражаются конечным набором значений и приводит к так называемой дискретной процедуре динамического программирования [10].

На рис.1 представлены примеры видов функций (табл. 1) на интервале оси [-4,4].

Экспериментальные результаты

В данной работе проводилось тестирование на стандартных изображениях в градациях серого цвета с 8 битами на пиксель (рис. 2) при различных уровнях аддитивного белого гауссовского шума со стандартным отклонением и = 5... 30. В экспериментах сравнивались результаты шумоподавления изображений процедуры динамического программирования с использова-

Рис. 1. Примеры график парно-потенциальных функций при u =1, Д = 1 (MQF и BZ), Д = 1(HF), в = 1 (SM), p = 1.2 (GGF)

нием других видов парно-потенциальных функций, которые представлены (в табл. 1) по скорости и точности.

Для оценки результатов шумоподавления использованы среднее значения индекса структурного сходства (MSSIM-Mean Structure Similarity Index) [20] и пиковое отношение сигнала к шуму (PSNR - peak-to-signal ratio). MSSIM определяется как среднее значение SSIM:

(2^p hQ + ci )(2ap,Q + c2)

SSIM (P, Q) = 2 2 2

V (ß2P + ß2Q + c1)(a2P + a2Q + C2)

PSNR (дБ) определяется как:

( 2552 ^

10

PSNR = 101og4 MSm ,

(5)

(6)

M N

где MSE = M1N ¿2 E (P) - Q(i,j))2; ßP, ßQ — средние значения изоб-i=i j=i

ражений; ар,uq — стандартные отклонения изображений; upq — ковари-ация изображений P(i,j) и Q(i,j); c1 =6.5025, c2=58.5225; P(i,j) и Q(i,j) обозначим значения пикселей исходного изображения и реконструированного или шумного изображения соответственно.

В экспериментах узловые функции выбраны квадратичной формой фг1г2(xt1t2) = (xt1t2 - уЫ2)2. Узловые функции фг1г2(Xt1t2) = (Xt1t2 - уЫ2)2 являются одинаковыми как для горизонтального Yh(xtl,t2-i,xtl,t2), так и для вертикального Yv(xtl-1,t2,xtl,t2) направления смежности Для функции связи MQF выбраны значения L = 2, L = 3. При L = 2, если вторая функция имеет вид постоянной функции, то получим форму функции BZ, p = 1.2. В случае функции GGF получены BF, когда p = 1.

В табл. 2 представлены парно-потенциальные функции для шумоподавления изображений процедуры динамического программирования. Экспериментальные параметры выбраны: Х^ = XV = 0.2, й = 0.5Д2 для функции MQF, Т = 1 для функции ИР, в = 1 для функции БЫ, р = 1.2 для функции СОР. Значения параметров пь, щ каждого вида парно-потенциальной функции выбраны потому, что результат измерения РБКИ имеет максимальное значение.

На рис.2 представлены исходные изображения (а, в) и их зашумленные (б, д) для экспериментальных исследований с соответственно белым гауссов-ским шумом со стандартным отклонением а = 10 и а = 20.

Таблица 2

Парно-потенциальные функции для экспериментальных исследований

Название

мдр ик ш1п[(же1,е2-1 - хч,г2)2, Хн(хч,г2-1 - хн,г2 )2 + й, Д2] пу ш1п[(х41-1,42 - Хг1М)2, Ху(хН-1М - х^2)2 + й, Д2]

ИР 7к = инЯн, Ъ = иуду, Г(xí1,í2-l - х^,г2)2 если х < Т {2Т М-1 - \- Т2 если х>Т' \(хг1-1,г21 - хг1,г2 )2 если х < Д ик ^ I 2 [2Д \хН-1М - хг1,г2 \- Д если х> Д

БЫ и^ + ^^ ^^ )2 1), Пув2 У1+ )2 1)

ССР ик \х41,42-1 - хи,12 \р , Пу \х41-1,42 - хг1,12 \р, 1 < р < 2

BZ ик Ш1п[(х41,42-1 - )2, Д2], Пу Ш1п[(х41-1,42 - )2, Д2]

ВР ик \хг1,г2-1 - хг1,г2 \ ,иу \хг1-1,г2 - хг1,г2 \ ,

На рис. 3, 4 представлены примеры результатов шумоподавления процедур динамического программирования с использованием различных видов парно-потенциальных функций при случае зашумленного изображения. В табл. 3-4 представлено сравнение результатов измерения РБКИ шумоподавления изображений (рис. 2.б, 2.д) процедуры динамического программирования с использованием различных видов функции.

Наилучшие результаты шумоподавления изображений представлены полужирными шрифтами (табл. 3-4). В табл. 5 представлено сравнение времени выполнения шумоподавления изображений процедур динамического программирования с использованием различных видов парно-потенциальных функций.

Результаты из табл. 3-5 показывают, что мультиквадратичная процедура динамического программирования позволяет получить высокую эффектив-

Рис. 2. Исходные изображения (а, в) и их зашумленные (б, д) для экспериментальных исследований. Зашумленные изображения (б, д) с соответственно белым гауссовским шумом со стандартным отклонением

а = 10 и а = 20

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

Рис. 4. Результаты шумоподавления изображения (рис. 1,д) процедуры динамического программирования с различными видами парно-потенциальных функций

Таблица 3

Сравнение результатов измерения РБМИ шумоподавления изображений процедуры динамического программирования с использованием различных

видов функции связи

а 5 10 20 30 5 10 20 30

РБМЕ^п 34.12 28.10 22.08 18.56 34.22 28.31 22.23 18.61

Вид функции Изображение (рис. 2.а) 128x128 Изображение (рис. 2.б) 256x256 РЯМКоредпее

мд 36.98 32.18 27.88 25.29 38.56 34.70 31.20 29.44 32.03

НЕ 36.46 32.05 27.71 25.29 38.67 34.58 31.06 29.39 31.90

БЫ 36.54 32.09 27.80 25.35 38.61 34.51 31.10 29.54 31.94

СОЕ 36.16 31.69 27.34 25.07 38.52 34.63 31.15 29.43 31.75

BZ 36.90 32.17 27.76 25.28 38.49 34.67 31.17 29.40 31.98

ВЕ 36.48 32.06 27.77 25.33 38.52 34.73 30.80 29.03 31.84

Таблица 4

Сравнение результатов измерения М881М шумоподавления изображений процедуры динамического программирования с использованием различных

видов функции связи

а 5 10 20 30 5 10 20 30

34.12 28.10 22.08 18.56 34.22 28.31 22.23 18.61

Вид функции Изображение (рис. 2.а) 128x128 Изображение (рис. 2.б) 256х256 ИЯЯШсреднее

мд 0.9635 0.9241 0.8458 0.7991 0.9570 0.9217 0.8539 0.8115 0.8846

ИР 0.9619 0.9128 0.8201 0.7572 0.9625 0.9011 0.8181 0.7949 0.8661

яы 0.9647 0.9178 0.8428 0.7600 0.9590 0.8977 0.8189 0.7896 0.8688

ООР 0.9559 0.8977 0.7747 0.7225 0.9572 0.8995 0.8282 0.7804 0.852

Б2 0.9637 0.9238 0.8348 0.7989 0.9567 0.9219 0.8537 0.8071 0.8826

ВР 0.9619 0.9128 0.8201 0.7527 0.9625 0.9067 0.8038 0.7877 0.8635

Таблица 5

Сравнение времени шумоподавления изображений процедуры динамического программирования с использованием различных видов

функции

Время выполнения шумоподавления изображений (в секундах)

Изображение а мд BZ НЕ БЫ СОЕ ВЕ

(рис. 2.а) 128x128 10 3.1 2.2 24.78 30.48 25.80 26.50

(рис. 2.б) 256x256 20 13.2 8.47 136.19 138.75 139.85 164.30

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

Заключение

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

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

1. Xiaojuan G, Li G. A new method for parameter estimation of edge-preserving regularization in image restoration // Journal of Computational and Applied Mathematics. 2009 V. 225. P. 478-486.

2. Stan Z. Li Markov Random Field Modeling in Image Analysis. Springer: Advances in Pattern Recognition Series, 2009. 371 p.

3. Geman S, German D. Stochastic Relaxation, Gibbs Distributions, and the Bayesian Restoration of Images // Pattern Analysis and Machine Intelligence. IEEE Transactions. 1984. PAMI-6(6). P. 721-741.

4. A method for approximating the density of maximum likelihood and maximum a posteriori estimates under a Gaussian noise model / C.K. Abbey [et al.] // Med. Im. Anal. 1998. V.2(4). P. 395-403.

5. Besag J. On the statistical analysis of dirty pictures // Journal of the Royal Statistical Society (Series B). 1986. V. 48. P. 259-302.

6. Geman S, German D. Stochastic Relaxation, Gibbs Distributions, and the Bayesian Restoration of Images // Pattern Analysis and Machine Intelligence. IEEE Transactions. 1984. P. 721-741.

7. Bouman C, Sauer K. A generalized Gaussian image model for edge-preserving MAP estimation // Image Processing. IEEE Transactions. 1993. V. 2(3). P. 296-310.

8. Blake A., Zisserman A. Visual Reconstruction. Cambridge, Massachusetts: MIT Press, 1987. 232 p. Bouman C., Sauer K. A generalized Gaussian image model for edge-preserving MAP estimation // Image Processing. IEEE Transactions. 1993. V. 2(3). P. 296-310.

9. Besag J. Towards Bayesian image analysis // Journal of Applied Statistics. 1989. V. 16(3). P. 395-407.

10. Kopylov A.V. Dynamic programming procedures for image analysis // Proceedings of the Eight IASTED International Conference Intelligent systems and control. Cambridge: ACTA Press. 2005. P. 404-409.

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

12. Hammersley J. M, Clifford P.E. Markov random fields on finite graphs and lattices. 1971.

13. Nikolova M. Analysis of the recovery of edges in images and signals by minimizing non-convex regularized least-squares // SIAM J. Multiscale Model. Sim. 2005. V. 4. № 3. P. 960-991.

14. Stevenson R., Delp E. Fitting curves with discontinuities // Proc. of the first international workshop on robust computer vision. 1990. P. 127-136.

15. De la Rosa J. I., Fleury G. Bootstrap methods for a measurement estimation problem // IEEE T. Instrumentation and Measurement. 2004. V. 55. № 3. P. 820-827.

16. Pham C.T. and Kopylov A.V. Multi-quadratic dynamic rogramming procedure of edge-preserving denoising for medical images // Int. Arch. Photogramm. Remote Sens. Spatial Inf. Sci. 2015. P. 101-106.

17. Моттль В.В., Мучник И.Б. Скрытые марковские модели в структурном анализе сигналов. М.: ФИЗМАТЛИТ, 1999. 352 с.

18. Фам К.Т., Копылов А.В. Мультиквадратичная процедура динамического программирования для восстановления изображений с сохранением локальных особенностей // Известия Тульского государственного университета. Естественные науки. 2014. Вып. 4. С. 143-158.

19. Bellman R. Dynamic Programming. Princeton, N.J.: Princeton University Press, 1957.

20. Image quality assessment: From error visibility to structural similarity / Z. Wang [et al.] // IEEE Transactions on Image Processing. 2004. V. 13 (4). P. 600-612.

Фам Конг Тханг (pacotha@gmail.com), аспирант, кафедра информационной безопасности, Тульский государственный университет.

Копылов Андрей Валериевич (And.Kopylov@gmail.com), к.т.н., доцент, кафедра информационной безопасности, Тульский государственный университет.

Analysis of pairwise potential functions for bayesian image

restoration

CongThang Pham, A.V. Kopylov

Abstract. The results of the comparative analysis of pairwise potential functions for image denoising and restoration on the Bayesian approach. Considered computationally efficient methods of edge preserving image

reconstruction with preservation based on dynamic programming procedure. The procedures were developed on these kinds of a priori probability properties of a hidden Markov field with the help of special types of functions, which defined on the cliques of the adjacency graph and parametric optimization procedure based on dynamic programming. A comparative study of image denoising methods in Bayesian framework for speed and accuracy.

Keywords: dynamic programming, Bayesian framework, Markov random fields, edge preserving smoothing, pairwise-separable optimization, potential functions.

Pham Cong Thang (pacotha@gmail.com), postgraduate student, Department of Information Security, Tula State University, Tula.

Kopylov Andrey (And.Kopylov@gmail.com), PhDoc, Associated Professor, Department of Information Security, Tula State University, Tula.

Поступила 12.06.2015