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

Известия Тульского государственного университета Естественные науки. 2014. Вып. 4. С. 143-157 Прикладная математика и информатика

УДК 004.932

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

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

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

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

Введение

Предварительная обработка изображений является необходимым этапом в работе практически любой системы компьютерного зрения, играющим важную роль в повышении качества дальнейшего анализа и интерпретации данных. Одной из основных проблем, возникающих в процессе предварительной обработки и восстановления изображений, в особенности полученных в плохих условиях регистрации, является сохранение или восстановление исходной локальной структуры границ и резких перепадов яркости, имеющихся на изображении. Важность данной проблемы для построения систем компьютерного зрения привела к появлению большого количества научных исследований [1-3] по этой тематике. В целом рассматриваемые в литературе подходы можно разделить на локальные и глобальные методы обработки с сохранением границ. К первой группе относятся такие методы, как анизотропная диффузия (anisotropic Diffusion) [2, 4, 5], основанная на решении дифференциального уравнения диффузии в частных производных, служащего моделью

преобразования изображения, билатеральная фильтрация (bilateral filtering) [6-8], использующая два фильтрующих ядра: пространственное ядро для учета пространственной зависимости и ядро интенсивности для учета различий в интенсивности между центральным пикселем и его соседями, наведенная фильтрация (guided filtering) [3, 9, 10] использующая при фильтрации данные о структуре некоторого другого, «управляющего» изображения. Ко второй группе относятся фильтрация на основе вейвлет-преобразования (Wavelet-Transform) [11-13] и байесовская фильтрация (Bayes Filtering) на основе вероятностных графовых моделей [14-16].

Байесовский подход к принятию решений является, пожалуй, одним из самых концептуально привлекательных и популярных подходов к восстановлению изображений. В рамках данного подхода задача восстановления изображений может быть выражена как задача оценивания скрытой марковской компоненты двухкомпонентного случайного поля (Markov random fields-MRFs), роль наблюдаемой компоненты которого играет анализируемое изображение. Использование в качестве функции потерь, сингулярной функции, штрафующей сам факт отличия полученной оценки скрытого поля от «истинного» значения приводит к оптимизационной задаче поиска максимума апостериорного распределения (MAP) скрытой компоненты, относительно наблюдения [15, 17-21]. Эквивалентная форма представления марковских случайных полей в виде гиббсовских случайных полей согласно теореме Хаммерсли-Клифорда [22], позволяет использовать для задания априорных вероятностных свойств скрытого марковского поля не переходные плотности распределения, а так называемые гиббсовские потенциалы на кликах [23].

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

Подходы, использующие негауссовские MRFs c невыпуклыми потенциальными функциями для устранения этого эффекта, были введены в работах [24-29], но поскольку такие функции являются невыпуклыми, глобальная оптимизация в оценке MAP не может быть точно вычислена, и обычно является разрывной функцией данных. Однако такие функции имеют значительные преимущества: сохранение границ и общий класс потенциальных функций.

В данной работе предлагается задавать парные потенциалы в виде функции двух соседних переменных [19, 23, 24], являющейся минимумом некоторого набора квадратичных функций вместо единственной квадратичной функции. Частным случаем такой функции является широко распространенная функция Блейка и Зиссермана («просевшая рессора») [26, 27]. Разработанная процедура анализа изображений, основанная на новых моделях скрытой компоненты, позволяет существенно расширить класс решаемых прикладных задач, учесть наличие неоднородностей и разрывов в исходных данных, сохранив при этом высокую вычислительную эффективность процедур динамического программирования [30] и фильтрации-интерполяции Калмана-Бьюси [31].

Параметрическая процедура восстановления сигналов с сохранением границ

В данной работе под задачей восстановления изображения будем понимать задачу преобразования исходных данных У = (у t € Т), заданных на дискретном растре, в другую функцию X = (хь,£ € Т), определенную на той же совокупности аргументов £ € Т и найденную в предположении, что ее значение изменяется в основном плавно лишь с относительно небольшим числом резких скачков. Формально задача восстановления изображений может быть выражена как задача оценивания скрытой марковской компоненты двухкомпонентного случайного поля, роль наблюдаемой компоненты которого играет анализируемое изображение. Использование в качестве функции потерь сингулярной функции, штрафующей сам факт отличия полученной оценки скрытого поля от «истинного» значения, приводит к оптимизационной задаче поиска максимума апостериорного распределения скрытой компоненты, относительно наблюдения. Эквивалентная форма представления марковских случайных полей в виде гиббсовских случайных полей позволяет использовать для задания априорных вероятностных свойств скрытого марковского поля не переходные плотности распределения, а так называемые Гиббсовские потенциалы на кликах [23]. В этом случае такая задача может быть сведена к задаче минимизации действительнозначной целевой функции следующего вида [30]:

3

г&т

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

сумма двух типов функций, называемых узловыми функциями и функциями связи. Узловые функции ^(х ) выбираются в зависимости от конкретной решаемой задачи таким образом, что каждая из них принимает тем большее значение, чем более очевидно расхождение между гипотезой о том, что хявляется именно тем локальным значением, которое мы ищем, и соответствующей окрестностью массива данных. Каждая функция связи 7£/ ¿п (х£/ ) налагает штраф на негладкость результата обработки на соответствующем ребре (£') графа соседства С [30, 32].

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

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

а б

Рис. 1. Граф С соседства на множестве элементов изображения: а — прямоугольная решетка; б — дерево

Парно-сепарабельная целевая функция (1) допускает построение эффективных неитерационных процедур глобальной оптимизации, в том случае, если граф соседства ее переменных представляет не имеет циклов, то есть представляет собой дерево [30, 33].

Процедура оптимизации при этом состоит в рекуррентной декомпозиции исходной задачи оптимизации функции |Т| переменных, где |Т| означает количество элементов в множестве Т, на последовательность |Т | элементарных подзадач оптимизации функции лишь одной переменной. Элементарные функции одной переменной 3 (х^), подлежащие оптимизации на каждом шаге процедуры минимизации целевой функции с древовидной

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

Фундаментальное свойство функции Беллмана [32]

3

будет называться прямое рекуррентное соотношение. Обращенная форма этого отношения

Хь-1

будет называться обратное рекуррентное соотношение.

В случае непрерывных переменных, например, при X = (хь, £ € Т), хь € € Ка, численная реализация процедуры динамического программирования возможна лишь если существует конечно-параметрическое семейство функций 3(х,а), замкнутое относительно узловых функций фь(хь) и функций связи 7ь',ь''(хь/,хЬ") в том смысле, что функции Беллмана 3ь(хь) на каждом шаге принадлежат этому семейству. В этом случае, процедура оптимизации состоит в рекуррентном пересчете параметроваь, которые полностью представляют функции Беллмана 3ь(х) = (х, аь).

Пусть функции связи в критерии (1) имеют следующую форму:

Ъ',г''

Такая узловая функция обладает необходимыми свойствами, для сохранения резких изменений в анализируемых данных [19]. Узловые функции, также выбираются в квадратичной форме (хь) = (уь — — хь)т (уь — хь). В этом случае можно построить параметрическую процедуру динамического программирования для минимизации критерия в виде (1).

Рассмотрим случай, когда граф смежности С имеет вид цепи (рис. 2), так что источник данных У = (уь, £ € Т) и результат обработки X = (хь, £ € Т) — сигналы и £ = 1, 2,..., N.

Рис. 2. Граф отношения соседства в виде цепи

Можно доказать, что если функция Беллмана в узле £ — 1:

3-1(хь-1) = ш1п[3ь(.-)1(хь-1), 3-1 (хь-1),..., 3-1-1 )(хь-1)], то в следующем шаге функция Беллмана будет

3(X) = шт[3(1)(х), 3(2)(х),..., (X)],

где Lt = Lt-i + 1, и в соответствии с

ко.

Jt(i)(xt) = + min [(xt - Xt-i)T Ut(xt - xt-i) + Ji(-)1(xi-i)j,

xt-i€X

i = 1, 2,

,Lt - 1,

Jt(Lt)(xt) = ^t(xt) + min [5 + Jt-1 (xt-i)].

xt-iex

Легко видеть, что каждая функция Jt(i)(xt) имеет параметрическое представление:

Jt(i)(xt) = (xt - x«)TR(i)(xt - x«) + dti), i = 1,2,..., Lt.

Обратное рекуррентное соотношение (3) принимает вид

xti)i(xt) = arg min <j min f(i)(xt-i, xt), ft(2)(xt-i, xt),...,

xt-iex

(2) /

ft(Lt-l)(xt-i ,xt),ft(Lt)(xt-i ,xt)]}

где ft(i)(xt-i,xt) = (xt - xt-i)'rUt(xt - xt-i) + J-^xt-i), i = 1,2,... ,Lt-i

ft(Lt )(xt-i) = 5 + Jt-i(xt-i).

Пусть Jt

(2) /

t(l)i

Если для любого xt € R Jt(i) ^ min

Jt

, то

Jt

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

min

Jt(i)(xt),...,Jt(L)(xt) =min[Jt(i)(xt),...,Jt(M)(xt)], M < L. (5

j-(L)/

(M),

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

сохранением границ

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

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

[32].

Сначала такая одномерная процедура применяется ко всем горизонтальным строкам £ 1 = 1,...,Ж 1 независимо для каждого ¿2 = 1,...^2.

Введем здесь дополнительное обозначение <Л1)Ь2 (хь1,ь2) для так называемых маргинальных функций узлов, каждый из которых показывает, как полная целевая функция зависит от значения одной переменной на соответствующем узле (£1 ,£2), если другие узловые переменные принимают оптимальные значения в соответствии хь1,ь2. Маргинальные функции Беллмана могут быть вычислены в результате начального горизонтального прохода по каждой из строк изображения.

Полученные маргинальные функции узлов (хь1,ь2) должны быть

сохранены в памяти. Затем процедура применяется ко всем столбцам £ 2 = 1,..., N2 изображения независимо для каждого £ 1 = 1,..., N1 с единственным изменением: соответствующие маргинальные функции узлов <Л1)Ь2 (хь1,ь2) принимаются вместо узловых функций ^^ (хь1,ь2) -оптимальные значения переменных хЬ1 ,ь2 сохраняются в качестве результатов.

Пусть узловые функции имеют вид ^^(х^2) = (х^2 — у^2)2, а функции связи (4) принимают следующую форму:

- функции связи по горизонтали:

Yh(xibi2-i,xtbi2) = uh min [(xtl,i2-i - xtl,i2)2 , Д2] ,

- функции связи по вертикали:

Yv(xti-i,t2,xtbt2) = Uv min [(xti-i,t2 - Xti,t2)2 , Д2] .

Пусть граф соседства изображения имеет вид (рис. 1), тогда целевая функция (1) запишется так:

N1 N2 N1 N2

3(Х ) = (Шь*2 - Х«1,*2 )2 + ^2 (Ж*1,*2-1,Ж*1,*2 ) +

¿1 = 1 ¿2 = 1 ¿1 = 1 ¿2 =2

N1 N2 у 7

+ ^ У] 7у (х*1-1 ,¿2 , ,¿2 ) • ¿1=2 ¿2 = 1

Здесь используем одинаковые функции связи в горизонтальных направлениях ¿1 и ¿2 в силу их одинакового физического смысла.

Предположим, что на шаге (¿ь Ь2 — 1), функция Беллмана <ЛЬ*2-1(^,^1) представлена в виде:

3*1,«2-1(х*1,*2-1) = ™П [311,)*2-1(х*1,*2-1)>3(12,)42-1(х*1,*2-1)> • • • '

Т№«1,«2-1)

^¿1,*2-1 ,¿2 1)

где 31)42-1(х*1,*2-1) = 9(1)42-1(х*1,*2-1 — х(Ь42-1)2 + ^8*2-1' ^ 0' * =

= 1, • • • ,

Тогда, в соответствии с прямым рекуррентным соотношением (2) процедуры динамического программирования, функция Беллмана (ж*1,*2) на следующем шаге (Ь1,Ь2) будет

,¿2 ,¿2 ) = т1п [^к ,¿2 )' ,¿2 )' • • • ' ''2 ) ,¿2 ^ '

3й2 ,¿2 ) = ,¿2 — ,¿2 )2 + Ш1П и  (^1,*2 — ,¿2 — 1) 2 + Ль*2-1 ,¿2 — 1) ^Й^2 ) ,¿2 ) = ,¿2 — ,¿2 )2 + т1П [и  Д2 + 31,*2-1(ж*1 ,¿2-1 ^ •

ХЬ1 ,¿2 — 1

Это означает, что каждая функция Беллмана имеет конечно-параметрическое представление. Примечательно также, что на каждом шаге число квадратичных функций, которые необходимы для представления функции Беллмана <Льг2), вообще говоря, увеличивается не более чем на единицу.

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

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

\2 , т

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

В соответствии с (5) и (7) рассмотрим квадратичные функции как точки в

(  Т

трехмерном пространстве с координатами ж^ , г = 1,..., М.

иу

Можно заметить, что согласно виду критерия (6), оцененное значение скрытого случайного поля не может превышать максимальное значение наблюдения Ь = утах и быть меньше минимального значения наблюдения а = ут;п. Поэтому расстояние между квадратичными функциями можно вычислить как сумму квадратов разности по ограниченной значениями [а, Ь] области определения:

ь

ЭД(ж),/ (ж)) = / [¿(ж) - / (ж)]2 ^ =

а

ь

= J [&  (ж - ж^)2 + ^ - ^  (ж - ж^-)2 - ^] ,

а

•ВД(ж), £ (ж)) = С^2 + С2з3 + С3  (2  + 8283) + 648183 + 658182

(Ь5 - а5) . 2(Ь3 - а3) . 2 ,2, . 4 ,4,

С1 = ---—1, С2 = (Ь - а), С3 = -, С4 = 2(а2 - Ь2); С5 = (а4 - Ь4),

5 3

81 = фж* - ^ ж^-, 82 = Зг - , 83 = ^ ж2 - ^ ж2 + ^ - ^, г = г,; = 1,...,М.

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

Использование алгоритма к-средних для редукции количества квадратичных функций позволяет получить эффективную реализацию процедуры обработки изображения. Экспериментальные исследования показывают, что для подавляющего большинства массивов исходных данных двух или трех квадратных функций, то есть к = 2 или к = 3, как правило, достаточно полностью отражать каждую функцию Беллмана в критерии

(7).

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

На рис. 3 представлены исходное простое модельное изображение и результат его сглаживания. Можно видеть, что в результате сглаживания границы изображения сохранились.

б г е

Рис. 3. Исходное изображение (а) и его деталь (б); результат обработки изображения гауссовской модели (в) и его деталь (г), результат обработки изображения разработанного алгоритма (д) и его деталь (е) при u = 30, А = 10

Для оценки эффективности шумоподавления изображений, такого как восстановление шумного изображения, проводилось тестирование на множестве изображений в градациях серого цвета с 8 битами на пиксель при различных уровнях аддитивного белого Гауссовского шума со стандартным отклонением: а = 10,15,20,25,30. Чтобы оценить производительность предложенного метода, проведено его сравнение с некоторыми существующими методами [12, 35, 36], имеющими наилучшие характеристики. Для оценки качества будем использовать пиковое отношение сигнала к шуму (PSNR - peak signal-to-noise ratio), которое определяется так:

/ 2552 \

PSNR =10ьЦ mSB) дБ'

M-1 N-1

где MSE = Y, (P(i, j) - Q(i, j))2 — среднеквадратичная

i=0 j=0

ошибка (Mean Square Error) между исходным изображением P(i,j) и восстановленным изображением Q(i, j) размера M x N.

На рис. 4 представлены результаты сравнения работы методов: Stein's Unbiased Risk Estimator for a linear expansion of threshold

[12], Probability shrinkage (ProbShrink) [35], the Bayesian least squares with Gaussians Scale Mixture (BLS-GSM) [36] и предлагаемого метода. Результаты измерения PSNR представлены в таблице.

Результаты экспериментов показывают, что разработанный алгоритм имеет лучшее качество восстановления изображения по сравнению с другими методами. Среднее значениеР£^Дсреднее = 30.37 превышает данное значение для других сравниваемых методов примерно на один процент.

Д е

Рис. 4. а - исходное изображение без шумов, б - его шумовая версия

при а = 20, PSNR = 22.11, в - шумоподавление с использованием SURE-LET при PSNR = 29.33, г - шумоподавление с использованием BLS-GSM при PSNR = 29.07, д - шумоподавление с использованием ProbShrink при PSNR = 28.85, е - шумоподавление с использованием нашей MQDP при PSNR = 29.56

Сравнение результатов обработки изображений

пиковое отношение сигнала к шуму - PSNR

Изображение и размера а ProbShrink [35] BLS-GSM [36] SURE-LET [12] MQDP

Peppers (256 х 256) 10 32.68 32.86 33.18 33.14

House (256 х 256) 10 33.84 34.26 34.29 34.67

Al (512 х 512) 10 34.58 34.83 34.90 35.15

GoldHill (512 х 512) 10 32.30 32.61 32.69 33.09

Peppers (256 х 256) 15 30.41 30.62 30.91 32.21

House (256 х 256) 15 31.74 32.23 32.32 32.29

Al (512 х 512) 15 32.64 32.93 32.97 32.94

GoldHill (512 х 512) 15 30.35 30.68 30.76 31.08

Peppers (256 х 256) 20 28.85 29.07 29.33 29.56

House (256 х 256) 20 30.29 30.79 30.93 31.23

Al(512 х 512) 20 31.28 31.58 31.64 31.98

GoldHill (512 х 512) 20 29.07 29.41 29.52 29.91

Peppers (256 х 256) 25 27.67 27.90 28.12 28.45

House (256 х 256) 25 29.20 29.65 29.86 30.12

Al(512 х 512) 25 30.08 30.53 30.64 31.06

GoldHill (512 х 512) 25 28.13 28.47 28.60 28.54

Peppers (256 х 256) 30 26.70 26.97 27.13 27.42

House (256 х 256) 30 28.35 28.72 28.98 29.41

Al(512 х 512) 30 29.32 29.68 29.84 30.22

GoldHill (512 х 512) 30 27.43 27.73 27.89 28.16

Peppers (256 х 256) 50 23.85 24.40 24.43 24.59

House (256 х 256) 50 25.99 26.15 26.58 27.32

Al(512 х 512) 50 27.18 27.35 27.61 29.98

GoldHill (512 х 512) 50 25.62 25.73 26.06 26.24

Средняя оценка P ^среднее 29.48 29.8 29.97 30.37

Заключение

На основе процедуры динамического программирования обработки сигналов [33], разработана мульти-квадратичная процедура динамического программирования обработки двухмерного изображения. С применением алгоритма к-средних выполнена редукция количества квадратичных функций и получена высокая вычислительная эффективность процедур динамического программирования для задачи обработки изображения. Среднее значение РБКЯ показывает, что разработанный алгоритм имеет

более высокое качество восстановления изображения по сравнению с

другими методами.

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

1. Gijbels I., Lambert A., Qiu P. Edge-preserving image denoising and estimation of discontinuous surfaces // Pattern Analysis and Machine Intelligence. IEEE Transactions. 2006. V. 28(7). P. 1075-1087.

2. Anisotropic diffusion-based detail-preserving smoothing for image restoration / M.C. Shin [et al.] // Image Processing (ICIP), 2010 17th IEEE International Conference, 2010. P. 4145-4148.

3. He K., Sun J. and Tang X. Guided image filtering // Proceedings of the 11th European conference on Computer vision: Part I (ECCV'10). Berlin, Heidelberg: Springer-Verlag, 2010. P. 1-14.

4. Perona P. and Malik J. Scale space and edge detection using anisotropic diffusion // IEEE Transactions on Pattern Analysis and Machine Intelligence, 1990. V. 12. P. 629-639.

5. A time-dependent anisotropic diffusion image smoothing method / X. Yu [et al.] // Proceedings of 2nd International Conference on Intelligent Control and Information Processing. 2011. P. 859-862.

6. Tomasi C. and Manduchi R. Bilateral filtering for gray and color images // Proceedings of Sixth International Conference on Computer Vision. 1998. P. 839-846.

7. Constant Time Joint Bilateral Filtering Using Joint Integral Histograms / K. Zhang [et al.] // IEEE Transactions on Image Processing, 2012. V. 21(9). P. 4309-4314.

8. Malik K. and Smolka B. Improved Bilateral Filtering Scheme For Noise Removal in Color Images // Proceedings of the International Conference on Informatics and Applications. 2012. P. 118-130.

9. Pham C.C., Ha S.V.U., Jeon J.W. Adaptive Guided Image Filtering for Sharpness Enhancement and Noise Reduction // Advances in Image and Video Technology. 2012. P. 323-334.

10. Feng Z., Liang L. Contrast Enhancement of Mammographic Images Using Guided Image Filtering // Advances in Image and Graphics Technologies, Communications in Computer and Information Science. 2013. P. 300-306.

11. Gnanadurai D. and Sadasivam V. An Efficient Adaptive Thresholding Technique for Wavelet Based Image Denoising // World Academy of Science, Engineering and Technology. 2008. V. 2: P. 922-927.

12. Luisier FBlu T., Unser M.A. New SURE Approach to Image Denoising: Interscale Orthonormal Wavelet Thresholding // IEEE transactions on image processing. 2007. V. 16(3). P. 593-606.

13. Saeedi J., Moradi M.H., Faez K. A new wavelet-based fuzzy single and multi-channel image denoising // Image and Vision Computing, 2010. V. 28(12). P. 1611-1623.

14. German D., Reynolds G. Constrained restoration and the recovery of distributions// IEEE Transactions on Pattern Analysis and Machine Intelligence. 1992. V. 14. P. 367-383.

15. 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.

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

17. Ameeramol P.M, Greeshma T.R. Bayesian MAP Model for Edge Preserving Image Restoration: A Survey // International Journal of Computer Applications. 2012. V. 1. P. 14-18.

18. Bayesian parallel imaging with edge-preserving priors / A. Raj [et al.]. Magn Reson Medicine. 2007. V. 57(1). P. 8-21.

19. A Comparative Study of Energy Minimization Methods for Markov Random Fields / R. Szeliski [et al.] // Pattern Analysis and Machine Intelligence. IEEE Transactions. 2006. V. 30. № 6. P. 1068-1080.

20. Jeffs B.D., Gardiner A.H. Markov random field image prior models for Map reconstruction of magnetoencephalogram images // Signals, Systems & Computers, 1998. Conference Record of the Thirty-Second Asilomar Conference. 1998. V. 1. P. 314-318.

21. Fessler J.A. Mean and variance of implicitly defined biased estimators (such as penalized maximum likelihood): Applications to tomography // IEEE Tr. Im. Proc. 1996. V. 5(3). P. 493-506.

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

23. 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.

24. 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.

25. German S., McClure D. Bayesian images analysis: An application to single photon emission tomography // Proc. Statist. Comput. sect. Amer. Stat. Assoc. Washington: DC, 1985. P. 12-18.

26. Blake A. Comparison of the efficiency of deterministic and stochastic algorithms for visual reconstruction // IEEE Trans. on Pattern Analysis and Machine Intelligence, 1989. V. 11(1). P. 2-30.

27. Blake A., Zisserman A. Visual Reconstruction. Cambridge, Massachusetts: MIT Press, 1987. 232 p.

28. German S., McClure D. Statistical methods for tomographic image reconstruction // Bull. Int. Stat. Inst. 1987. LII-4. P. 5-21.

29. Hebert T., Leahy R. A generalized EM algorithm for 3-D Bayesian reconstruction from Poisson data using Gibbs priors // IEEE Trans. on Medical Imaging. 1989. V. 8(2). P. 194-202.

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

31. Kalman R.E., Bucy R.S. New Results in Linear Filtering and Prediction Theory // Journal of Basic Engineering. 1961. V. 83. P. 95-108.

32. Kopylov A.V. Parametric dynamic programming procedures for edge preserving in smoothing of signals and images // Pattern recognition and image analysis. 2005. V. 15. P. 227-229.

33. A Signal Processing Algorithm Based on Parametric Dynamic Programming / A. Kopylov [et al.] // Lecture Notes in Computer Science. 2010. V. 6134. P. 280-286.

34. MacQueen J. Some methods for classification and analysis of multivariate observations // Proceedings of 5-th berkeley symposium on mathematical statistics and probability. Berkeley: University of California press, 1967. P. 281-297.

35. Pizurica A., Philips W. Estimating the probability of the presence of a signal of interest in multiresolution single- and multiband image denoising // IEEE transactions on image processing. 2006. V. 15(3). P. 645-665.

36. Image Denoising Using Scale Mixtures of Gaussians in the Wavelet Domain / J. Portilla [et al.] // IEEE transactions on image processing. 2003. V. 12(11). P. 1338-1351.

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

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

Multi - Quadratic Dynamic Programming Procedure for Image Restoration Preserving Local Features

C. T. Pham, A. V. Kopylov

Abstract. An approach to image restoration with edge preserving feature is developed on the basis of dynamic programming multi-quadratic procedure. The procedure of image analysis, based on the new data models and clustering, significantly expands the class of applied problems, and can take into account the presence of heterogeneities and discontinuities in the source data, while retaining high computational efficiency of the dynamic programming procedure and Kalman-Bucy filter-interpolator.

Keywords: multi-quadratic procedure, dynamic programming, filtering -interpolation Kalman-Bucy, edges-preserving restoration.

Pham Cong Thang (pacotha@gmail.com), postgraduate student, department of information security, Tula State University.

Kopylov Andrey (And.Kopylov@gmail.com), candidate of technical sciences, associate professor, department of information security, Tula State University.

Поступила 12.11.2014