СПОСОБ УСТРАНЕНИЯ ШУМА НА РАСТРОВЫХ ИЗОБРАЖЕНИЯХ НА ОСНОВЕ ПОЛНОЙ ВАРИАЦИИ
METHOD OF REMOVAL OF NOISE ON RASTER IMAGES BASED ON FULL VARIATION

Данное изобретение относится к области обработки изображений и может быть использовано для устранения шума на растровых изображениях в научно-технической сфере, медицине, биологии и т.д.
В современных исследованиях растровые изображения являются важным типом информации. В растровых изображениях наблюдаются различные эффекты, в том числе и шум. Шум снижает качество изображений и качество их обработки. Поэтому проблема устранения шума на цифровых изображениях по-прежнему остается актуальной. В настоящее время разработано много способов устранения шума на растровых изображениях (Chan T.F., Shen J. Image processing and analysis: Variational, PDE, Wavelet, and stochastic methods. - SIAM, 2005, 400 c; Burger M., Mennuci A.C.G, Osher S., Rumpf M. Level set and PDE based reconstruction methods in imaging. - Springer, 2008, 319 c).
Близким к заявленному изобретению является способ устранения шума, предложенный авторами: Рудиным, Ошером, Фатеми (L.I. Rudin, S. Osher, Е. Fatemi. Nonlinear total variation based noise removal algorithms. Physica D, 1992, Vol. 60, 259-268). Ими было предложено использовать полную вариацию функции яркости растрового изображения для устранения шума:
,
2
2
x
y
где u(х, у)∈R - идеальное (без шума) и ν(x, y)∈R - реальное (зашумленное) изображения, , ∇u=(u, u), , , , λ>0 - множитель Лагранжа,  - полная вариация.
Недостаток данного способа заключается в эффективном устранении только гауссовского шума, а для других типов шумов этот способ оказывается не очень эффективным.
Также близким к заявленному изобретению является способ Ганнана Юаня и соавторов (Yuan Gannan, Han Zifa, Zhang Jie, Dong Jing, Zhao Yuxin, Li Tao, Song Chengye, Li Qianf, Guo Ruiliang - патент КНР CN103198455 A, опубликован 10.07.2013). Данный способ заключается в выполнении следующих шагов:
σ
2
1. Определить фильтр Гаусса зашумленного изображения Н=ν*G, где  - гауссовское ядро, σ - дисперсия распределения Гаусса.
2. Определить четыре матрицы совпадений градаций серого (Gray-Level Co-Occurrence Matrix - GLCM) при обработке фильтром Гаусса по разным направлениям (0°, 45°, 90°, 135°). Четыре матрицы GLCM определяются на основе соседних пар пикселей, определенных заданным направлением и расстоянием. Так, если расстояние равно одному пикселю, а направление - горизонтальное (90°), то для определения сопряженности берутся все пары пикселей в окне, расположенные рядом по горизонтали. Аналогично определяются пары для остальных направлений.
3. Вычислить контрастное изображение как среднее значение четырех матриц GLCM. Модель устранения шума построена на основе минимизации функционала энергии полной вариации от функции яркости зашумленного изображения и функции яркости контрастного изображения.
4. Построить уравнение Эйлера-Лагранжа и решить его итерационным методом.
Недостаток данного способа заключается в фильтрации только гауссовского шума. Данный способ также не рассчитан на эффективное устранение других типов шумов.
Еще одно изобретение, которое является близким к заявленному изобретению, - это способ устранения шума, предложенный авторами: Ли, Чартрандом и Асаки (Т. Le, R. Chartrand, T.J. Asaki. A variational approach to reconstructing images corrupted by Poisson noise. Journal of mathematical imaging and vision, 2007, Vol. 27, Is. 3, 257-263):
,
где β>0 - коэффициент регуляризации. Недостаток данного способа заключается в устранении только пуассоновского шума, а для других типов шумов этот способ также оказывается не очень эффективным.
Все вышеуказанные способы заключаются в том, что в ходе процедуры устранения шума использовано характерное свойство полной вариации функции яркости: полная вариация зашумленных изображений всегда больше полной вариации соответствующих гладких изображений. Сущность данных способов устранения шума на основе полной вариации функции яркости заключается в выполнении следующих действий:
1. Ввести полную вариацию функции яркости зашумленного изображения. Пусть в пространстве  задана ограниченная область . Если функция и является гладкой, то ее полная вариация имеет вид: .
T
2. Определить ограничение для решения задачи оптимизации V[u]→min.
3. Свести задачу условной оптимизации (с ограничением) к задаче безусловной оптимизации (без ограничения).
4. Построить уравнение Эйлера-Лагранжа задачи оптимизации без ограничения.
5. Построить численную схему для решения уравнения Эйлера-Лагранжа.
6. Оценить качество устранения шума через оценки качества восстановленного изображения.
Основной недостаток вышеприведенных способов заключается в отсутствии гарантии того, что будут эффективно устранены распространенные типы шумов, так как эти способы предназначены только для устранения гауссовского шума или только для устранения пуассоновского шума. Рассмотренные способы неэффективно устраняют другие типы шумов и особенно смесь гауссовского и пуассоновского шумов. Смесь шумов этих двух типов часто рассматривается в медицине, биологии, космонавтике и т.д.
Задачей, на решение которой направлено заявляемое изобретение, является обеспечение гарантии эффективного устранения, как смеси шумов, так и отдельно гауссовского и пуассоновского шумов.
Данная задача решается за счет того, что в заявленном изобретении решается задача минимизации полной вариации функции яркости зашумленного изображения с ограничением на интенсивность шума. Для заявленного способа вычисляются параметры, определяющие оптимальное решение.
Способ устранения шума на основе полной вариации характеризуется тем, что восстанавливается оригинальное изображение и для заданного зашумленного изображения ν следующим образом:
1. вычисляется полная вариация функции яркости ;
2. формулируется задача минимизации полной вариации функции яркости с ограничением на интенсивность шума;
3. полученная задача условной оптимизации с ограничением сводится к задаче безусловной оптимизации без ограничения в виде функционала Лагранжа;
4. строится уравнение Эйлера-Лагранжа для решения задачи безусловной оптимизации;
5. строится численная схема решения данного уравнения: ищутся оптимальные входные параметры, сравниваются отклонения значений функций яркости двух изображений на двух последовательных шагах итерации с заданной точностью ε, если Δu≤ε, то решение найдено, если Δu>ε, то процесс итерации продолжается.
Структурная схема для способа устранения шума (фиг. 1) содержит: входной блок 1, на который поступает зашумленное изображение. Выход блока 1 соединен с входом блока подготовки 2, в котором вычисляется полная вариация функции яркости зашумленного изображения и определяется ограничение на интенсивность шума. Выход блока подготовки 2 является входом блока преобразования 3, выход которого соединен с входом блока вычислений и оптимизации 4. В блоке преобразования 3 задача оптимизации с ограничением сводится к задаче без ограничения и строится уравнение Эйлера-Лагранжа. Численная схема определяется в блоке 4. Выход блока 4 соединен со входом блока корректировки 5 с обратной связью. На основе сравнения отклонений функций яркости двух изображений на двух последовательных шагах итерации с заданной точностью, которые определяются в блоке 4, проверяется условие остановки. Выход блока 5 соединен со вторым входом блока 4. Выход блока 5 является результатом устранения шума в исходном изображении.
Техническим результатом является эффективное устранение как гауссовского и пуассоновского шумов в отдельности, так и линейной комбинации этих шумов.
Сущность изобретения поясняется чертежами, на которых изображено:
Фиг. 1 - структурная схема для устранения шума;
Фиг. 2 - сравнение эффективности использования нескольких способов устранения шума на примере реальных изображений с искусственными шумами.
Способ устранения шума на растровых изображениях осуществляется следующим образом.
1. Взять зашумленное изображение ν.
T
2. Вычислить полную вариацию функции яркости зашумленного изображения ν и сформулировать задачу оптимизации V[u]→min.
3. Определить ограничения. Предполагается, что для заданного изображения вариация шума постоянна (пуассоновский шум не изменяется, а гауссовский шум зависит только от его вариации):
где p(ν|u) - условная вероятность наблюдения реального изображения ν при заданном идеальном изображении u.
Рассмотрим гауссовский шум. Плотность распределения определяется как:
.
2
ν
Плотность пуассоновского шума определяется как: p(ν|u)=exp(-u)u/ν!.
Обратим внимание, что значения функций яркости изображения u и ν - это целые числа (например, для восьмибитового изображения интервал яркости определяется значениями от 0 до 255).
Для устранения смеси гауссовского и пуассоновского шумов рассмотрим их линейную комбинацию:
1
l
2
2
1
2
1
2
ln(p(ν|u))=λln(p(ν|u))+λln(p(ν|u)), где λ>0, λ>0, λ+λ=1. Согласно (1), получим задачу устранения шума с ограничением:
,
где κ - постоянное значение.
4. Свести задачу условной оптимизации с ограничением к задаче безусловной оптимизации с использованием функционала Лагранжа:
,
чтобы найти решение:
где τ>0 - множитель Лагранжа.
1
2
2
2
1
1
1
2
2
2
В данной модели, если λ=0 и λ=1, то при β=τλ=τ будет получена модифицированная модель ROF для устранения пуассоновского шума. Если λ=0 и λ=1, то при λ=τλ/(2σ)=τ/(2σ) будет получена модель ROF для устранения гауссовского шума. Если λ>0, λ>0, то будет получена модель для устранения смеси гауссовского и пуассоновского шумов.
5. Построить численную схему для модели устранения смеси шумов. Для решения задачи (2) можно применить метод множителей Лагранжа. Но здесь используется уравнение Эйлера-Лагранжа.
2
Пусть функция ƒ(x, y) определена в ограниченной области Ω⊂R и непрерывно дифференцируема до второго порядка по х и у при (х, у)∈Ω.
x
y
Пусть F(x, y, ƒ, ƒ, ƒ) - выпуклый функционал, где , . Решение задачи оптимизации  удовлетворяет уравнению Эйлера-Лагранжа:
,
где , , .
Тогда решение задачи (2) удовлетворяет следующему уравнению Эйлера-Лагранжа:
где μ=1/τ.
Представим уравнение (3) в следующем виде:
где .
Для получения дискретной модели (4) добавим искусственный параметр времени u=u(x, y, t). Уравнение (4) соответствует следующему уравнению диффузии:
Тогда дискретная форма уравнения (5) имеет вид:
,
, , ,
, ,
1
2;
; ; ; ; i=1, …, N; j=1, …, N
k=0, 1, …, K; Δx=Δу=1; 0<ξ<1, K - достаточно большое число, K=500.
1
2
1
2
k
6. Найти оптимальные параметры модели устранения смеси шумов. Процедура (6) может быть использована для устранения шума на изображениях, если значения параметров λ, λ, μ, σ заданы. Часто на практике эти параметры неизвестны, и их нужно оценить. Тогда параметры λ, λ, μ в процессе (6) нужно представить как , , μ на каждой итерации k. В новой процедуре такие параметры будут вычисляться на каждом шаге итерации.
1
2
Пусть (u, τ) является решением задачи (2). Тогда мы получим условие . Данное условие позволяет вычислить оптимальные параметры линейной комбинации шумов λ, λ:
2
1.
, λ=1-λ
Дискретная форма для вычисления этих параметров имеет вид:
, , k=0,1, …, K.
Для поиска оптимального параметра μ умножим (3) на (ν-u) и проинтегрируем по частям по всей области Ω. В итоге получим формулу для поиска оптимального параметра μ:
.
Дискретная форма вычисления этого параметра имеет вид:
,
где
,
, , ,
, ; ; ; ;
0j
1j
ν=ν; ; ;
k=0, 1, …, K; Δx=Δy=1.
Для вычисления параметра о здесь использован способ Иммеркера (J. Immerker, Fast noise variance estimation. Computer vision and image understanding, 1996, Vol. 64, No. 2, 300-302):
1
2
где  - маска изображения, N и N определяют размер изображения.
Оператор* - это оператор конволюции, выполняющий преобразование:
l
2
ij
1
2
где i=1, …, N; j=1, …, N; u=0, если i=0, или j=0, или i=N+1, или j=N+1.
Параметр σ вычисляется сразу же на первой итерации.
7. Оценить качество изображений и качество устранения шума. Для оценки качества изображений после устранения шума использованы критерии PSNR, MSE и SSIM:
,
, ,
где , , ,
, ,
1
1
2
2
1
2
2
2
С=(KL), С=(KL); K<<1; K<<1.
1
2
-6
8
Например, K=K=10, L=2-1=255 - яркость 8-битового серого изображения.
PSNR
PSNR
Чем больше Q, тем лучше качество изображения. Если значение Q лежит в интервале между 20 и 25, то качество изображения приемлемо, например, для беспроводной передачи.
MSE
MSE
MSE
MSE
PSNR
Значение Q использовано для оценки различия между двумя изображениями, где Q - среднеквадратичная ошибка. Чем меньше значение Q, тем лучше результат восстановления. Значение Q прямо связано со значением Q.
SSIM
SSIM
Значение Q использовано для оценки качества изображения с помощью сравнения сходства двух изображений. Его значение лежит в интервале [-1, 1]. Чем больше значение Q, тем лучше качество изображения.
0
0
0
ij
8. Оценить зависимость результата устранения шума от начального решения. Очевидно, что в локальной итерационной процедуре (6) результат в общем случае зависит от начальных значений параметров , , μ. Если сначала задать параметры , , μ, то неудачные значения определят не очень хорошие оценки u, а через них - оценки параметров распределений. Случайный выбор параметров , , μ также неприемлем, т.к., фактически, вносит дополнительный шум в изображение.
0
0
Очевидно, что начальные значения параметров , , μ должны быть, по возможности, достаточно близки к тем значениям, которые будут найдены. Поэтому оценим параметры , , μ как средние по соседним пикселям изображения, используя, например, метод Иммеркера.
Рассмотрим пример реализации способа. Предложенная модель была протестирована на реальных изображениях. В качестве примера показан результат для изображения черепа человека размером 300×300 пикселей (фиг. 2а). На остальных изображениях (фиг. 2б-2з) показан увеличенный фрагмент этого изображения.
1
2
1
2
Для получения зашумленного изображения сначала был добавлен гауссовский шум (фиг. 2в), потом - пуассоновский шум (фиг. 2д). На фиг. 2ж показано зашумленное изображение для смеси двух шумов с параметрами λ=0,8, λ=0,2. Параметры линейной комбинации λ и λ были определены следующим образом.
2
ij
l
2
Сначала был рассмотрен пуассоновский шум с плотностью распределения p(ν|u) и вариацией , относительно среднего u в каждом пикселе с координатами (i, j) i=1, …, N; j=1, …, N.
(2)
(2)
ij
ij
Функция яркости такого изображения обозначена как ν. Ее значения должны находиться в интервале от 0 до 255. Если значение выходит за этот диапазон, то оно не изменяется ν=u. На данном изображении оказалось всего пять таких значений (0,0056%). Общая вариация пуассоновского шума определена как среднее .
1
2
(1)
(1)
(1)
(2)
(1)
(2)
Далее, нами было принято, что вариация гауссовского шума σ=40,2412 будет в четыре раза больше σ. Функция яркости такого изображения обозначена как ν. Как и раньше, значения функции яркости ν также должны быть в интервале от 0 до 255. В этом случае оказалось, что 5780 (6,42%) пикселей со значением яркости  находятся за пределами этого диапазона. Результирующее изображение (фиг. 1ж) образовано двумя зашумленными изображениями в пропорции 0,5 для ν и 0,5 для ν. Это означает, что ν=0,5ν+0,5ν. Следовательно:
.
1
2
В итоге получим, что коэффициенты линейной комбинации имеют значения, соответственно, λ=4/5=0,8 и λ=1/5=0,2.
PSNR
MSE
SSIM
Для зашумленного изображения критерии качества имеют значения, соответственно, Q=21,4168, Q=427,9526 и Q=0,4246.
В табл. 1-3 показаны результаты устранения шума на данном изображении для случаев заранее заданных и автоматически определенных параметров.
PSNR
SSIM
Заметим, что в этом случае значение Q после устранения шума для заранее заданных (идеальных) параметров лучше, чем для автоматически полученных оценок, хотя для значения Q иногда наблюдается и обратное.
Для создания начального изображения был использован оператор конволюции (7). В табл. 4 показана зависимость восстановленного результата от начального решения, где:
(а) начальные параметры , μ=1;
(б) начальные параметры , μ=1;
0
(в) начальное решение μ - это случайная матрица заданного размера;
0
0
(г) начальное решение μ найдено как усреднение соседних пикселей u=ν*А оператором конволюции, где .
В табл. 4 показано, что наилучший результат устранения комбинированного шума соответствует случаю (г) выбора начального решения по критериям PSNR и MSE.
1
2
Таким образом, здесь предложен способ устранения смеси гауссовского и пуассоновского шумов на основе известного вариационного подхода. Качество результата устранения шума зависит от значений коэффициентов линейной комбинации λ и λ. Их значения должны быть заранее заданы, или они должны быть автоматически определены, что важно в случае обработки реальных изображений. Для реальных изображений данный способ с автоматически определяемыми параметрами дает результат, близкий по качеству к идеальному, когда правильные значения параметров заданы заранее.
Данный способ может быть использован для устранения отдельно как гауссовского, так и пуассоновского шумов. В этом случае качество обработки практически не уступает качеству способов, специально предназначенных для устранения шума только одного вида.
Табл. 1. Сравнение качества способов устранения шума на реальном изображении для смеси шумов
Табл. 2. Сравнение качества способов устранения шума на реальном изображении с гауссовским шумом
Табл. 3. Сравнение качества способов устранения шума на реальном изображении с пуассоновским шумом
Табл. 4. Зависимость результата устранения комбинированного шума от начального решения