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

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

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

УДК 517.97, 519.63

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

Н. Х. Т. Данг, С. Д. Двоенко

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

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

Введение

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

В настоящее время разработано много методов устранения шумов для случаев, когда тип шума известен. Например, метод полной вариации [1-13] является известным и эффективным методом.

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

модель ROF предназначена только для устранения гауссовского шума [12, 13].

Другим популярным шумом является пуассоновский шум. Например, такой шум появляется в рентгеновских снимках. Модель ROF не сможет эффективно устранить такой шум. Поэтому в работе Ли [14] была построена другая модель, известная как модифицированная модель ROF.

Оба типа шумов (гауссовский и пуассоновский) распространены, но их комбинация также важна [15]. Она появляется на биомедицинских изображениях, например, на изображениях электронной микроскопии [16, 17].

В данной статье предлагается совместно применить модель ROF и модифицированную модель ROF для устранения комбинации шумов.

В экспериментах были использованы реальные изображения с натуральным неизвестным шумом. Предполагалось, что это смесь гауссовского и пуас-соновского шумов. Качество обработки сравнивалось с другими методами устранения шума: модель ROF, модифицированная модель ROF, медианный фильтр [18], фильтр Винера [19], регуляризация Бельтрами (Beltrami) [20]. Для сравнения качества изображений до и после их обработки был использован известный критерий BRISQUE (Blind/Referenceless Image Spatial QUality Evaluator) [21].

1. Модель устранения шума

Пусть в пространстве R2 задана ограниченная область Q С R2. Назовем функции u(x, y) Е R2 и v(x, y) Е R2, соответственно идеальным (без шума) и реальным (зашумлённым) изображениями, где (x,y) Е О.

Если функция u гладкая, то ее полная вариация имеет вид

VT[u] = J \Vu\dxdy, п

где Vu = (ux,uy) — градиент, ux = du/dx, uy = du/dy, \Vu\ = ^uX + u2. В

этой работе мы считаем, что полная вариация функции u ограничена Vt[u] < < ж.

Согласно работам [1, 2, 12, 13], гладкость изображений можно характеризовать их полной вариацией. Полная вариация зашумлённых изображений всегда больше полной вариации соответствующих гладких изображений. Для решения задачи Vt[u] ^ min необходимо ввести ограничение на

вариацию гауссовского шума:

j(v - ufdxdy =

п

В этих условиях модель ROF для устранения гауссовского шума на изображении имеет вид [12]

u* = arg min | J \Vu\dxdy + Л J (v — u)2dxdy \п п

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

Для устранения пуассоновского шума на основе модели ROF была предложена другая модель [14]. Такая модель получается при решении задачи Vt[u] — min с ограничением:

/эдф))«* = /(„ — v ln(u))dxdy = const

п п

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

u* = arg min I I \ Vu\dxdy + ß / (u — v ln(u))dxdy

и

где ß > 0 - коэффициент регуляризации. Такая модель известна как модифицированная модель ROF для устранения пуассоновского шума.

Чтобы построить модель устранения смешанного шума, будем также решать задачу устранения шума, основанную на свойстве гладкости полной вариации Vt[u] — min.

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

1 v:>)"xdy =const' (1)

п

где p(v \u) - условная вероятность наблюдения реального изображения v при заданном идеальном изображении u.

Рассмотрим гауссовский шум. Его плотность распределения с дисперсией

2

и определяется как:

pi(v\u) = exp ^—/(aV2n).

Плотность пуассоновского шума

р2(у|и) = ехр(—и) и"/у\.

Обратим внимание, что значения функций яркости изображения и и V - это целые числа (например, для восьмибитового изображения интервал яркости определяется значениями от 0 до 255).

Для устранения комбинации гауссовского и пуассоновского шумов рассмотрим следующую линейную комбинацию:

\п(р(у\п)) = Х\ \п(р1(у\п)) + Л2 1п(р2Им)),

где Л1 > 0, Л2 > 0, Л1 + Л2 = 1.

Согласно (1), получим задачу устранения шума с ограничениями:

и* = а^шт/ \ Vu\dxdy, и п

/ (2^"(^ — и)2 + Л2(и — V 1п(и))) dxdy = к , п а

где к — постоянное значение.

Сведем эту задачу к задаче безусловной оптимизации с использованием функционала Лагранжа

Ь(и,т) = у \Vu\dxdy + т ( j(V — и)2dxdy+ Л2 ^(и — V 1n(u))dxdy — к

п V п п

чтобы найти решение:

(и*,т * ) = ащшт Ь(и,т), (2)

и, т

где т > 0 — множитель Лагранжа.

В данной модели, если Л1 =0 и Л2 = 1, то при в = тЛ2 = т будет получена модифицированная модель ИОР для устранения пуассоновского шума. Если Л2 = 0 и Л1 = 1, то при Л = тЛ 1 /(2а2) = т/(2а2) будет получена модель ИОР для устранения гауссовского шума. Если Л 1 > 0, Л2 > 0, то будет получена модель для устранения смеси гауссовского и пуассоновского шумов.

2. Дискретная модель устранения шума

Для решения задачи (2) можно применить метод множителей Лагранжа [22, 23, 24]. Его применение в данной задаче оказывается не совсем удобным, поэтому в данной работе использовано уравнение Эйлера-Лагранжа [24].

Пусть функция f ^,у) определена в ограниченной области О С И2 и непрерывно дифференцируема до второго порядка по x и у при ^,у) £ О. Пусть Ру, ¡, ¡х, ¡у) - выпуклый функционал, где ¡х = дf/дx, ¡у = д//ду. Решение задачи оптимизации

У Ру, ¡, ¡х, ¡у)dxdy ^ шт

ху

п

удовлетворяет уравнению Эйлера-Лагранжа:

дд ^.У, Цх,^) — д^^х^^¡¡х^у) — д^^у^^¡¡х^у) = 0

где Еу = дЕ/д/, Еи = дЕ/д/х, Еу = дЕ/д/у.

Тогда решение задачи (2) удовлетворяет следующему уравнению Эйлера-Лагранжа:

А1 л V)_ д | их д | иу

а2 У и 2 и ^дх\ + иу) + иу

где ц = 1/т.

Представим уравнение (3) в следующем виде:

Х1 , , . , иххиУ — 2ихиу иху + иХиУУ „

(V — и) — (1 — и ) + Ц-(их + и2)3/2- = 0' (4)

где

д2и д2и д ди

ихх - м 2 ' иуу - м 2 ' иху - П П - иух-

дх2 ду2 дх \ду)

Для получения дискретной модели (4) добавим искусственный параметр времени и = и(х,у,г). Уравнение (4) соответствует следующему уравнению диффузии:

ди X1, х V, иххиУ — 2ихиуиху + и1иуу

и = я7 = ~2 ^ — и) — Х2(1 — 7)+ Ц-( 2 , 2)3/2-• (5)

дг а2 и (их + иу )3/2

Пусть размер изображения N х N2. Тогда дискретная форма уравнения (5) имеет вид

икг+ = ик3 + (^ — и%) — Х2(1 — ) + цфк3) , (6)

V а и1] )

г = ^хх(икгз )(Уу ))2 +

Фг3 ((Ух(ик3 ))2 + (Уу (и* ))2)3/2 +

+ —2Ух(ик3 )Уу (и% )Уху «•) + (Ух(ик3 ))2Ууу )

((Ух(ик3 ))2 + (Уу (и* ))2)3/2 '

и / к _ и / к и / к _ и / к

V (иг.) = _V (ик-) = и^+1_

Ух(иг3) = 2Дх > Уу (игз) = 2Ду ,

ик к к

У хх(игЗ) —

иг+1,з 2игз + иг-1, хх^г3)- (Дх)2

XI (ик ) = и1,з + 1 2ик + икз-1

УууУ и'з) - (Ду)2

Vху (ик ) =

<7 /к _ II0 _ II0 I II0

иг+1 1+1 иг+1 1— 1 и4—1 1+1 "г

4AxAy

и01 = и11; = ; ий) = и01; и0^2+1 = и>0^2 ; г = 1,...,М; ; = 1,...,^;

к = 0,1,..., К; Ax = Ау = 1; 0 < { < 1, где К — достаточно большое число, К = 500.

3. Параметры модели устранения шума

Процедура (6) может быть использована для устранения шума на изображениях, если значения параметров Л1, Л2, у, а заданы. Часто на практике эти параметры неизвестны, и их нужно оценить. Тогда параметры Л1, Л2, у в процессе (6) нужно представить как Л°, Л°, ук на каждой итерации к. В новой процедуре такие параметры будут вычисляться на каждом шаге итерации.

3.1. Оптимальные параметры Л1 и Л2

Пусть и,т являются решением задачи (2). Тогда мы получим условие дЬ(и,т)/ди = 0. Данное условие позволяет вычислить оптимальные параметры линейной комбинации шумов Л1, Л2:

/ (1 — и ^У

л 1__п_ Л2_1_Л1

/(V — u)dxdy + /(1 — и)dxdy,

пп

Дискретная форма для вычисления параметров имеет вид

N1 N2

£ £(1 — ик)

\к _ г=11=1 гз \к 1 хк

Л1 = N1 N2 ик ' Л2 = 1 — Л1'

£ £ (^ + 1 — ^)

г=11=1 ч

где к = 0, 1, . . . , К.

3.2. Оптимальный параметр у

Для поиска оптимального параметра сглаживания у умножим (3) на (V — — и) и проинтегрируем по частям по области О. В итоге получим выражение для оптимального параметра у:

/ (— % (V — и)2 — Л2^ )dxdy

п

уп

/ (^ и2 + и2 — ^У

Его дискретная форма имеет вид

^ (-£ (% - ик)2 - А| Ч4;

к _ г=13=1 "

¡1 _

N1 N2

£ £ пк

г=13=1

где

п% ))2 + (V,(и3))2 -

Ух{и% ) + VУ (и% )Vу (-ц )

"Ц) *х\игз) I ^уК^Ц! V У\иг]

У 3 3 ' 3))2 + (Vy(икз))2

^J(Vx(uk3 ))2 + (V, (икз))

и Iк _ ик ик _ ик

V (ик) _ иг+1'3_и'г-1'3 V (ик) _ иг'3+1_з3-1

Ух(из) _ 2Дж ' Уу(из) _ 2Ду

о —укл, . , . , 1 - -к

V (-к) _ -±13_—3 ^ V (-к) _ -31_-г) _ 2Дж - Уу(-гз) _ 2Ду ,

ик3 _ икз ; uN1+l'j _ ; ик0 _ ик1, ик^2+1 _ uk'N2;

-оз _ -13; —N1+1,3 _ —N1 з; -го _ —1; —'N2+1 _ —,N2; г _ 1,N1, ^ _ 1, ...,N2; к _0,1,..,К; Дж _ Ду _ 1.

3.3. Оптимальный параметр а

Для вычисления параметра а здесь использован метод Иммеркера [25]:

/П/2 N1 N2

а - 2) £* А|. (7)

г=1 3=1

/ 1 -2 1 \

где Л _ -2 4 -2 - маска изображения. Оператор * - это оператор

V 1 -2 1 У

конволюции (свертки):

иг3 * Л _ иг-1'3-1Лзз + иг,з-1Л32 + щ+1,3-1Лз1 + Щ-1,3Л23+

+иг3 Л22 + иг+1'3 Л21 + иг-^+^з + иг^Л^ + иг+1'3+1Лп,

где г _ 1,..., .1; j _ 1,..., .2; иг3- _ 0, если г _ 0, или j _ 0, или г _ + 1, или j _ N2 + 1. Параметр а вычисляется сразу же на первой итерации.

4. Поиск начального решения

Очевидно, что в локальной итерационной процедуре (6) результат обработки в общем случае зависит от начальных значений параметров А1, А„, Л0.

Если сначала задать параметры Л0, Л0, /0, то их неудачные значения определят не очень хорошие оценки Uj, а через них также не очень хорошие оценки параметров распределений.

Случайный выбор параметров Л0, ЛЦ, /0 также неприемлем, так как фактически вносит дополнительный шум в изображение.

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

5. Оценка качества изображений

Оценка качества реальных зашумленных изображений, для которых неизвестно, как они выглядят без шума (отсутствуют т.н. опорные изображения) - это трудная теоретическая задача. В настоящий момент разработано несколько методов для решения такой задачи. В данной статье используется критерий BRISQUE, предложенный Анишом (Anish) в работе [21]. Чем больше значение критерия QBRISQUE, тем лучше качество изображения.

Оценка качества изображения на основе метода BRISQUE заключается в выполнении следующих шагов:

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

2) построение специального вектора, компоненты которого получены из коэффициентов НСС;

3) применение методов регрессии для оценки значения критерия качества изображения.

6. Эксперимент

Предложенная модель была протестирована на реальных изображениях с неизвестным естественным шумом. Например, здесь показан результат обработки изображения «aphid» (рис. 1а) из базы данных Университета Альберты (Канада) [26]. Данное изображение получено на электронном микроскопе. На рис. 1б показан увеличенный фрагмент.

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

При обработке предполагалось, что шум на изображении представляет собой смесь гауссовского и пуассоновского шумов в неизвестной пропорции. Поэтому в ходе обработки были автоматически определены следующие значения параметров: Л1 = 0.4042, Л2 = 0.5958, / = 0.5238, а = 8.3884. Значение Qbrisque показано в табл. 1.

в Д

Рис. 1. Устранение шума на реальном изображении: а) - реальное изображение; б) - увеличенный фрагмент реального изображения; в) - изображение после устранения шума; д)- увеличенный фрагмент изображения после устранения шума

Таблица 1

Сравнение качества методов устранения шума на реальном изображении

по критерию qbrisque

Шум ROF Модиф. ROF Медиана фильтр Винера фильтр Бельтрами метод Предложенный метод

20.0455 29.0058 27.5682 30.2719 23.4842 31.9801 35.7603

Заключение

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

Качество результата устранения шума зависит от значений коэффициентов линейной комбинации Х\, и коэффициента сглаживания

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

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

1. Chan T.F., Shen J. Image processing and analysis: Variational, PDE, Wavelet and stochastic methods. SIAM, 2005. 400 p.

2. Level set and PDE based reconstruction methods in imaging / M. Burger [et al.]. Springer, 2008. 319 p.

3. An introduction to total variation for image analysis / A. Chambolle [et al.] // Theoretical foundations and numerical methods for sparse recovery. 2009. V. 9. P. 263-340.

4. Xu J., Feng X., Hao Y. A coupled variational model for image denoising using a duality strategy and split Bregman // Multidimensional systems and signal processing. 2014. V. 25. P. 83-94.

5. Rankovic N, Tuba M . Improved adaptive median filter for denoising ultrasound images // Advances in computer science, WSEAS ECC'12. 2012. P. 169-174.

6. Lysaker M, Tai X . Iterative image restoration combining total variation minimization and a second-order functional // International journal of computer vision. 2006. V. 66. P. 5-18.

7. Li F, Shen C, Pi L .A new diffusion-based variational model for image denoising and segmentation // Journal mathematical imaging and vision. 2006. V. 26. No. 1-2. P. 115-125.

8. Noise reduction with low dose CT data based on a modified ROF model / Y. Zhu [et al.] // Optics express. 2012. V. 20. No. 16. P. 17987-18004.

9. Tran M.P., Peteri R., Bergounioux M. Denoising 3D medical images using a second order variational model and wavelet shrinkage // Image analysis and recognition. 2012. V. 7325. P. 138-145.

10. Getreuer P. Rudin- Osher-Fatemi total variation denoising using split Bregman // IPOL 2012. URL: http://www.ipol.im/pub/art/2012/g-tvd/ (дата обращения: 23.07.2015).

11. Caselles V., Chambolle A., Novaga M. Handbook of mathematical methods in imaging. Springer, 2011. 1607 p.

12. Rudin L.I., Osher S, Fatemi E. Nonlinear total variation based noise removal algorithms // Physica D. 1992. V. 60. P. 259-268.

13. Chen K. Introduction to variational image processing models and application // International Journal of computer mathematics. 2013. V. 90. No. 1. P. 1-8.

14. Le T, Chartrand R., Asaki T.J. A variational approach to reconstructing images corrupted by Poisson noise // Journal of mathematical imaging and vision. 2007. V. 27. No. 3. P. 257-263.

15. Luisier F, Blu T, Unser M. Image denoising in mixed Poisson-Gaussian noise // IEEE transaction on Image processing. 2011. V. 20. No. 3. P. 696-708.

16. Poisson-Gaussian noise parameter estimation in fluorescence microscopy imaging / A. Jezierska [et al.] // IEEE International Symposium on Biomedical Imaging 9th. 2012. P. 1663-1666.

17. An EM approach for Poisson-Gaussian noise modeling / A. Jezierska [et al.] // EUSIPCO 19th. 2011. V. 62. No. 1. P. 13-30.

18. Wang C., Li T. An improved adaptive median filter for Image denoising // ICCEE. 2012. V. 53. No. 2.64. P. 393-398.

19. Abe C., Shimamura T. Iterative Edge-Preserving adaptive Wiener filter for image denoising // ICCEE. 2012. V. 4. No. 4. P. 503-506.

20. Zosso D, Bustin A. A Primal-Dual Projected Gradient Algorithm for Efficient Beltrami Regularization. Computer Vision and Image Understanding, 2014. http://www.math.ucla.edu/~zosso/ (дата обращения: 23.07.2015).

21. Mittal A., Moorthy A.K., Bovik A.C. No-Reference Image Quality Assessment in the Spatial Domain // IEEE Transactions on Image Processing. 2012. V. 21. No. 12. P. 4695-4708.

22. Zeidler E. Nonlinear functional analysis and its applications: Variational methods and optimization. Springer, 1985. 662 p.

23. Rubinov A., Yang X. Applied Optimization: Lagrange-type functions in constrained non-convex optimization. Springer, 2003. 286 p.

24. Gill P.E., Murray W. Numerical methods for constrained optimization. Academic Press Inc., 1974. 283 p.

25. Immerker J . Fast noise variance estimation // Computer vision and image understanding. 1996. V. 64. No. 2. P. 300-302.

26. База данных биомедицинских изображений университета Алберта (Alberta University), Канада: http://www.ualberta.ca/~mingchen/index.htm (дата обращения: 23.07.2015).

Данг Нгок Хоанг Тхань ([email protected]), аспирант, кафедра информационной безопасности, Институт прикладной математики и компьютерных наук, Тульский государственный университет.

Двоенко Сергей Данилович ([email protected]), д.ф.-м.н., профессор, кафедра информационной безопасности, Институт прикладной математики и компьютерных наук, Тульский государственный университет.

On restoration of an image intensity function based on the total

variation method

N.H.T. Dang, S.D. Dvoenko

Abstract. we consider the restoration problem of intensity function of an original noiseless image based on the total variation method. Today, there are many modern devices to create digital images. These devices have optical defects to create images. Hence, the quality of images depends on the quality of optical sensors. Because of the limit of technology, such defects are always included, for example, as noise. Noise reduces image quality and result of image processing. In order to reduce noises effectively, many special methods are used. For electronic microscope images the noise includes Gaussian noise, Poisson noise and also their linear combination. The solution of the restoration problem of image intensity function allows us to remove those noises more effectively.

Keywords: restoration problem of intensity function, total variation, ROF model, Gaussian noise, Poisson noise, image processing, biomedical image, Euler-Lagrange equation.

Dang Ngoc Hoang Thanh ([email protected]), postgraduate student, department of Information security, Institute of Applied Mathematics and Computer Sciences, Tula State University.

Dvoenko Sergey ([email protected]), doctor of physical and mathematical sciences, professor, department of Information security, Institute of Applied Mathematics and Computer Sciences, Tula State University.

nocmynuAa 18.09.2015