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

УПРАВЛЕНИЕ, ВЫЧИСЛИТЕЛЬНАЯ ТЕХНИКА И ИНФОРМАЦИОННЫЕ ТЕХНОЛОГИИ

УДК 004.932

АДАПТИВНАЯ ПРОЦЕДУРА ОБОБЩЕННОГО СГЛАЖИВАНИЯ ИЗОБРАЖЕНИЙ НА ОСНОВЕ СТАТИСТИЧЕСКОГО ПОДХОДА

И. А. Грачева, А.В. Копылов, О.В. Красоткина

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

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

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

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

Большинство популярных методов сглаживания изображений были разработаны на основе уравнений в частных производных (Partial differential equations, PDE) и вариационных моделей. Например, нелинейная диффузия (Anisotropic diffusions, AD) [1], а также методы регуляризации, основанные на полной вариации (Total variation, ТУ) [2], являются наиболее популярными и широко используются как нелинейные методы фильтрации обработки сигналов и изображений. В целом, сглаженное изображение постепенно приближается к исходному варианту, однако является более гладким и простым в некотором смысле. Эти методы достаточно эффективны при сглаживании изображений с сохранением границ, однако они реализуются как итерационный процесс, требующий значительных вычислительных затрат. Хорошей альтернативой алгоритмам, основанным на полной вариации, является двусторонний фильтр, который был впервые назван алгоритмом Томази и Мандучи [3] на основе работы [4], [5], а позже модифицирован и улучшен в [6]. Его формулировка проста, и данный неитерационный метод достигает удовлетворительных результатов только за один проход. Однако непосредственное выполнение двустороннего фильтра, как известно, медленно. Хотя несколько методов [7] - [10] предлагают ускорить оценку двустороннего фильтра, его быстрая реализация все еще является сложной проблемой. Так же к наиболее известным методам сглаживания изображений с сохранением границ относятся метод Haar-Fisz [11] на основе Haar преобразования, Platelet алгоритм [12], разработанный Виллеттом и Новаком. Пожалуй, одним из самых концептуально привлекательных и популярных является байесовский подход к восстановлению изображений. В частности в работе [13] предложен BLS-GSM алгоритм для удаления части шума на основе многомерной Байессовской оценки результата методом наименьших квадратов (Bayesian least-squares, BLS), а в работе [14] предложен подход Байессовской оценки интенсивности изображения, которая считается кусочно - гладкой. Чаще всего в работах данного типа задача сглаживания сводится к оптимизационной задаче поиска максимума апостериорной вероятности скрытой компоненты двух-компонентного случайного поля, используемого в качестве модели изображения и желаемого результата обработки [15]. Существующие параметрические беспереборные процедуры поиска максимума апостериорного распределения, как правило, позволяют учитывать лишь весьма узкий класс априорных предположений о решаемой прикладной задаче, оставаясь в рамках линейной нормальной модели скрытого случайного поля. В ряде работ, например [16], рассматриваются и негауссовские априорные распределения скрытой компоненты, однако процедуры минимизации для таких задач остаются чрезвычайно трудоемкими.

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

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

Гамма - нормальная модель скрытого случайного поля. В рамках Байессовского подхода, задача восстановления изображений может быть выражена как задача оценивания скрытой марковской компоненты X = (х,г = 1,...,М) Т={,=,2):^ = 1,...,N1,= 1,...,М2} двухкомпонентного случайного поля, роль наблюдаемой компоненты которого играет анализируемое изображение У = (у, ,е Т).

Использование в качестве функции потерь, сингулярной функции, штрафующей сам факт отличия полученной оценки скрытого поля X = (х1, 1 е Т) от «истинного» значения приводит к оптимизационной задаче

поиска максимума апостериорного распределения скрытой компоненты, относительно наблюдения.

Вероятностные свойства двухкомпонентного случайного поля (X, У) полностью определяются совместной условной плотностью вероятности Ф(У | X,5) исходной функции У = (у, ^ = 1,..., М) по отношению к вторичному массиву данных X = (х(, , = 1,..., М), и априорным распределением ¥( X | Л, 5) скрытой компоненты X = (х(, , = 1,..., М).

Условную плотность вероятности Ф(У | X,5) примем нормальной:

где Е(е2) = 5 дисперсия шума в наблюдении, которая считается неизвестной.

Так же, как в [5, 6], мы априори имеем марковский случайный процесс, который начинается с нормально распределенного первого значения Е(х1) = 0, с шумом Е(ХХТ) = 0. Но, в отличие от [17, 18], мы не предполагаем, что общая дисперсия независимых компонент остается одинаковой со временем, а изменяется Е(££Т) = . Неизвестные дисперсии (г, ^ = 1,..., М) считаются пропорциональными дисперсии шума в наблюдении г = 15 с коэффициентами пропорциональности 1, действующими как факторы неизвестного мгновенного изменения во времени вектора X = (х(, , = 1,..., М), которые так же неизвестны.

В соответствии с этим предположением, априорное совместное распределение оценок последовательности X = (х(, t = 1,..., М) условно нормально по отношению к изменяющейся во времени последовательности факторов Л = (1, t=2,..., М). Таким образом, мы пришли к априорной совместной плотности распределения вида:

1 1 м

Ф(У\X ,5) = 5^2 ехр(-^ У,- X )2),

(1)

где V - граф соседства элементов изображения, имеющий вид решетки.

Наконец, предположим, что обратные факторы 1/1 должны быть априори независимы и одинаково гамма - распределены g(1/1 | a, J) - (1/1 )a-1 exp(-J(1 /1)) на положительной полуоси 1t > 0. Математическое ожидание и дисперсия гамма - распределения вычисляются как отношения a/J и a/J соответственно. Априорное распределение плотности всей последовательности изменяющихся во времени факторов имеет вид:

N f N i a f N i

G(L | a, J) = nr(1t|a,J) - П"Т exp ~¿LT

V t=2 Tt J

t=2t j

= exp

N

N 1

-(a-1)2 ln Tt -J T

t=2 t=2 Tt

(3)

Мы можем представить параметры а и г? с помощью новых параметров X и ^

a = — 2

d

1

Л

1 + —

v mj

+1

,J=

T

2dm'

(4)

предполагая тем самым параметрическое семейство гамма - распределении случайных величин 1/1

2т+1 с 1 ^

g(1/T |d,T,m) - (1/T)2dm exp

с математическим ожиданием и дисперсией

(1/ T)

^(1 / T) =(1+d)m+1, Уаг(1 /1) = 2dm

2dm

(1+d)m+1

(5)

1 __г. 12 .

С точки зрения этой параметризации, независимое априорное распределение каждого мгновенно обработанного фактора 1/1 почти полностью сконцентрированы вокруг математического ожидания 1/1, если т®0. При т®^ коэффициент 1/1 имеет тенденцию к почти равномерному распределению.

Таким образом, мы пришли к априорному распределению вида:

G(L | S,1, m) = exp

1

N í

z

11

1-+ — ln1

1 1 t

(6)

28т =2

и также, полностью определили совместное априорное нормальное гамма -распределение обоих последовательностей X = (х, [ = 1,..., V) и л = (1,1 = 2,..., V):

н( х, л 18,1,т) = ^ (х | л, 8)в(л 18,1,т).

В сочетании с условной плотностью (2), это дает основу для байесовского оценивания вектора X = (х, ? = 1,..., V).

Байессовское оценивание скрытой компоненты. Совместное апостериорное распределение скрытых последовательностей, а именно

=2

1

1

массива вторичных данных и их мгновенно изменяющихся факторов, полностью определяется (1) и (2), и, с точки зрения исходных параметров (сс,$), в силу (3):

Р( Л-, А|У Л**> - х\кз)ама,ажг\хл-

Л ¥( X А ',^)С(А '| а,$)Ф(У | X \8)бХ' б А' Байесовская оценка (X, А) является точкой максимума в числителе

(X, Л | d, a, J) = arg max[ln Ф( Y | X, d) + ln Y(X | Л, d) + ln G(L | a, J)] =

X ,Л

f

аг§ тах

X ,А

1 N 1 N 11

-^Z(y-x)2-ZindA-28 z J_(x-xr)2-

28 = 2 t=2 28 e,ceV Ae

N 1 ^

(a- 1)Z in A-JZt

t =2

t=2 At J

или

(X, Л | d, a, J) = arg min J( X, Л | Y, d, a, J),

X ,Л N

J (X, Л | Y,d,a,J) = Z (7t - x )2 +

t=1

Z \A- [(X - x■ )2 + 2dj] + d(2a - 1)ln A [.

^и "еУ I1'

Замена новыми параметрами (4) делает байесовские оценки, независимыми от дисперсии шума в наблюдениях 5:

(X, Л | A,m) = arg min J( X, Л | Y, A, m),

X ,Л

J(X,Л | Y,A,m) = Z(7t-x)2 +

t=1

z [(X-xO2+A/m]+(1+1/m)lnA[.

(7)

t \t "eV IAt'

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

Алгоритм сохранения границ. Условно оптимальные изменяющиеся во времени факторы А( XI, т) - [1 (XI, т), I - 2,..., Щ определяются

независимо друг от друга:

Л( X, A, m) = arg min ЛЛ | X, A, m):

Л

-^8 \-i[( xf,- xr)2 +a/m]+(1+1/m)lnAt \ = o.

dlt I Ät

Нулевые условия для производных, за исключением тривиальных решений 1 ® ¥, приводят к равенствам

1

и, следовательно,

1 [(*г)2 +1/= (1+1/т),

(1/1)( хе- х( „)2+1/ т

1 (х ,1,т)=1

1+1/т

(9)

Подстановка (9) в (7) дает эквивалентную форму, которая позволяет избежать немедленного нахождения самих факторов:

(X, Л11, т) = агБ ш1й Д X, Л | У, 1, т),

X ,Л

Л (X, ЛI У,1,т)=£ (у - X )2 +

f=l

Х|(1+1/1)( х- * -)2+1/т

(10)

^^у 1 1+1/ т

Эта квадратичная в непосредственной близости от нулевой точки функция (х(,-х(„)2 = 0, и остается таковой практически на всей числовой оси, если значение ^ достаточно мало (рис.1).

(1 + 1///) 1л

1 + 1/ и

-10

V ■ ! 1 : ■ и=0.001

и^.1 //

..... £ \ Jigf.fi

- 1 \ / " 1 !

Рис. 1. "Эффект насыщения " из за штрафа при достаточно большом значении параметра ц и фиксированном значении 1=1

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

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

395

метров ^ и X, мы применяем итерационную процедуру Гаусса-Зейделя для обеих групп переменных X = (х(, t = 1,..., V) и Л = (1*, t=2,..., V), при начальных значениях Л0 = (10 = 1, ^ = 2,..., V).

На каждой итерации текущая последовательность Лк = (11=1, г = 1,..., V) в соответствии с (7) позволяет вычислить вектор X = (х*, t = 1,..., V), причем минимизация целевой функции дает новое при-б)лижение к оценке последовательности взаимосогласованных оценок Xк = (), * = 1,..., V):

Xк = (X)1,1 = 1,..., V) = а^шт J(X,Лк | У,1,т)

к = (х , I = 1,..., V) = агБ

X

(11)

агБ тп ■] I (- х)2 + I -1 (хс - X-)2 [.

X [ *=1 * и "е К 1

Такая оптимизация имеет линейную вычислительную сложность по отношению к размеру входного массива изображения и может быть легко решена с помощью фильтра Калмана-Бьюси [18, 19]. После нахождения оценок X1' = (хк, *= 1,..., V), находят очередное приближение к оценкам факторов Лк = (1к=1, ^ = 1,..., V) в соответствии с (9):

Г = (х -^ +Л/т, *= 2,..., V, (12)

1+1/т

которые, в соответствии с (8), дают решение задачи условной оптимизации

Лк+1 = агБ шт J( Xк, Л | У, 1, т) =

Л

V

агБ шт

V 1 I

I (у* - х )2 + ) - )')2|.

л [ (=1 * и "ек Л*'

Сама структура итерационной процедуры (11) - (12) обеспечивает удовлетворение целевой функции неравенству

J( Xк+, Лк+11 у, 1, т) < J( хк, Лк | у, 1, т),

которое, как это показано ранее, имеет место быть только в стационарной модели.

Остается только указать способ выбора значений структурных параметров X - сглаживание вектора X = (х(, * = 1,..., V) и ^ - коэффициент нестабильности вектора X = (х(, * = 1,..., V).

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

Результаты. Для измерения уровня искажений изображения наиболее часто используется значение Р8КК - пиковое отношение сигнала к

396

шуму (англ. peak signal-to-noise ratio). PSNR является соотношением между максимумом возможного значения изображения и мощностью шума, искажающего значения изображения. Проще всего его определить через среднеквадратичную ошибку, которое для двух монохромных изображений € и x размера m*n, одно из которых считается зашумленным приближением другого, вычисляется так:

PSNR = 10log

10

^ max(x2) ^

I € - x |2

В таблице 1 приведено сравнение алгоритмов сглаживающих изображение по Р8КК при наложении на исходное изображение Гауссовского шума. А в таблице 2 приведено сравнение времени работы алгоритмов.

Таблица 1

Сравнение алгоритмов по РБМЯ

Изображение Cameraman 256x256

исходного изображения 24.08 21.07 18.05 16.29 13.28 10.27 6.29 3.28

Иааг-Р1в7[11] 29.73 27.98 26.35 25.45 23.70 22.55 20.77 19.77

БЬ8-08М[12] 30.85 29.13 27.54 26.56 24.63 22.50 19.07 14.44

РМеВДЗ] 30.54 29.10 27.66 26.80 25.14 23.56 21.72 20.57

Наш метод 30.77 29.03 27.53 26.67 25.07 23.45 21.55 20.38

Изображение Barbara 512x512

исходного изображения 24.00 20.99 17.98 16.22 13.21 10.20 6.22 3.21

Иааг-Р1в7[11] 28.36 26.59 25.06 24.28 23.09 22.18 21.24 20.47

БЬ8-08М[12] 31.24 29.40 27.55 26.46 24.65 22.61 19.90 14.97

РМеВДЗ] 28.62 26.53 24.59 23.74 22.96 22.32 21.47 20.65

Наш метод 28.98 27.76 26.11 25.40 23.67 22.35 21.37 20.54

Таким образом, из табл. 1 вид, что предложенный в данной

397

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

Таблица 2

Сравнение времени работы алгоритмов

Алгоритм/Изображение Cameraman 256x256 Barbara 512x512

PSNR Время (с) PSNR Время (с)

Иааг-Р1в7[11] 29.73 1.3 28.36 1.3

BLS-GSM[12] 30.85 4.6 31.24 4.6

РМе^Э] 31.03 1112 29.55 891

Наш метод 30.77 0.5 28.98 0.8

На рис. 2 представлен график зависимости времени работы алгоритма от размера исходного изображении.

Рис. 2. График зависимости времени работы алгоритма от размера изображения

На рис. 3-4 представлены результаты работы алгоритма на разных изображениях при значениях параметров сглаживания изображения X и жесткости сглаживания ^ таких, чтобы Р8КК имело максимально возможное значение.

а б в

Рис. 3. Результат работы алгоритма на примере изображения Cameraman 256x256: а - исходное изображение; б - зашумленное изображение; в - сглаженное изображение

в

а б

Рис. 4. Результат работы алгоритма на примере изображения Barbara 512x512: а - исходное изображение; б - зашумленное изображение; в - сглаженное изображение

Как видно из представленных результатов, алгоритм выдает хорошие результаты как при наложении на исходное изображение Гауссовско-го шума (рис. 4, б), так и Пуассоновского шума (рис. 3, б). Таким образом, алгоритм является не только быстрым по скорости обработки исходного изображения, но и универсальным. Минимальное время работы предложенного в данной статье алгоритма, по сравнению с другими алгоритмами, позволяет с его помощью обрабатывать изображения высокого разрешения за приемлемое время, что в современном мире является актуальной задачей, так как с ростом разрешения растет и размер требующего обработки изображения.

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

1. P. Perona and J. Malik, "Scale-space and edge detection using anisotropic diffusion," IEEE Trans. Pattern Anal. Mach. Intell. Vol. 12. No. 7. 1990. P. 629-639.

2. L. Rudin, S. Osher, and E. Fatemi "Nonlinear total variation based noise removal algorithms," Phys. D. Vol. 60. 1992. P. 259-268.

3. C. Tomasi and R. Manduchi, "Bilateral filtering for gray and color images," in Proc. 6th Int. Conf. Comput. Vis., Bombay, India, Jan. 1998. P. 839-846.

4. V. Aurich and J. Weule, "Non-linear Gaussian filters performing edge preserving diffusion," in Proc. DAGM Symp., Bielefeld, Germany, 1995. P. 538-545.

5. S. M. Smith and J. M. Brady, "SUSAN—A new approach to low level image processing," Int. J. Comput. Vis. Vol. 23. No. 1. 1997. P. 45-78.

6. M. Elad, "On the origin of the bilateral filter and ways to improve it," IEEE Trans. Image Process. Vol. 11. No. 10. 2002. P. 1141-1151.

7. F. Durand and J. Dorsey, "Fast bilateral filtering for the display of high- dynamic range images," ACMTrans. Graph. Vol. 21. No. 3. 2002. P. 257266.

8. S. Paris and F. Durand, "A fast approximation of the bilateral filter using signal processing approach," in Proc. Eur. Conf. Comput. Vis., 2006. P. 568-580.

9. F. Porikli, "Constant time O(1) bilateral filtering," in Proc. IEEE Conf. Comput. Vis. Pattern Recognit., Anchorage, AK, Jun. 2008. P. 1-8.

10. Q. Yang, K. H. Tan, and N. Ahuja, "Real-time O(1) bilateral filtering," in Proc. IEEE Conf. Comput. Vis. Pattern Recognit.,Miami,FL, Jun. 2009.

11. P. Fryzlewicz and G. P. Nason, "A Haar-Fisz algorithm for Poisson intensity estimation," J. Comput. Graph. Stat. Vol. 13. No. 3. 2004. P. 621-638.

12. R.M.Willett and R.D. Nowak, "Multiscale Poisson intensity and density estimation," IEEE Trans. Inf. Theory. Vol. 53. No. 9. 2007. P. 3171-3187.

13. J. Portilla, V. Strela, M. J. Wainwright, and E. P. Simoncelli, "Image denoising using scale mixtures of Gaussians in the wavelet domain," IEEE Trans. Image Process. Vol. 12. No. 11. 2003. P. 1338-1351.

14. Raj A., Singh G., Zabih R. Kressler, B. Wang, Y. Schuff N., Weiner M.: Bayesian parallel imaging with edge-preserving priors. Magn. Reson. Med. 57. 2007. P. 8-21.

15. Humphrey D., Taubman D. A filtering approach to edge preserving MAP estimation of images. IEEE Trans. Image Process. 20, 2011. P. 1234-48.

16. Bouman C., Sauer K. A generalized Gaussian image model for edge-preserving MAP estimation. Image Process. IEEE Trans. 2, 1993. P. 296-310.

17. И.А. Грачева, А.В. Копылов. Адаптивный параметрический алгоритм обработки изображений. Известия ТулГУ. Технические науки.

Вып. 9. Ч. 2. Тула: Изд-во ТулГУ, 2013. С. 61-67.

18. M. Markov, V. Mottl, I. Muchnik. Principles of Nonstationary regression Estimation: A New Approach to Dynamic Multi-factor Models in Finance. DIMACS Technical Report 2004-47, Rutgers University, USA, 2004.

19. M. Markov, I. Muchnik, V. Mottl, O. Krasotkina. Dynamic analysis of hedge funds. Proceedings of the 8th IASTED International Conference on Financial Engineering and Applications. MIT, Cambridge, Massachusetts, USA, October 9-11, 2006.

Грачева Инесса Александровна, gia1509@mail.ru, Россия, Тула, Тульский государственный университет,

Копылов Андрей Валериевич, канд. техн. наук, доц., and.kopylov@gmail.com, Россия, Тула, Тульский государственный университет,

Красоткина Ольга Вячеславовна, канд. ф.-м. наук, доц., O. V.Krasitkina@yandex.ru, Россия, Тула, Тульский государственный университет

ADAPTIVE PROCEDURE FOR GENERALIZED SMOOTHING OF IMAGES ON THEBASIS

OF STA TISTICAL APPROA CH

I.A.Gracheva, A.V.Kopylov, O.V.Krasotkina

In this paper, we consider a generalization of the Bayesian approach to the image analysis, which involves automatic adjustment of the changing smoothness parameters to the observable data. Such generalized approach allows us to preserve substantial local image features, in particular edges of objects. Proposed gamma-normal model of the image and the expected result of processing, considered as a two-component random field, allowed to develop parametric adaptive algorithm for image smoothing.

Key words: image processing, hidden Markov models, dynamic programming, antialiasing.

Gracheva Inessa Aleksandrovna, gia1509@mail. ru, Russia, Tula, Tula State University,

Kopylov Andrey Valerievich, candidate of technical science, docent, and.kopylov@gmail. com, Russia, Tula, Tula State University,

Krasotkina Olga Vjacheslavovna, candidate of mathematical science, docent, O. V.Krasitkina@yandex.ru, Russia, Tula, Tula State University