Устранение шума на изображении при помощи обратной свертки

Открыть пример в Engee

Применим алгоритмы устранения шума и размытия на изображении при помощи библиотеки Deconvolution.

Описание задачи

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

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

Вот библиотеки, которые нам понадобятся:

Pkg.add( ["Images", "TestImages", "Deconvolution", "FFTW", "ZernikePolynomials"] )

Обратная свертка Винера

Загрузим подопытное изображение и добавим к нему шум. Заменив функцию textimage("image") на load("image.png") вы сможете открыть любое изображение с диска.

using Images, TestImages, Deconvolution, FFTW

img = channelview(testimage("cameraman"))

# Создадим размывающий фильтр в частотной области
x = hcat(ntuple(x -> collect((1:512) .- 257), 512)...)
k = 0.001
blurring_ft = @. exp(-k*(x ^ 2 + x ^ 2)^(5//6))

# а еще у нас будет аддитивный шум
noise = rand(Float64, size(img))

# Преобразование Фурье позволяет получить картинку с шумом
blurred_img_ft = fftshift(blurring_ft) .* fft(img) .+ fft(noise)

# Преобразуем изображение из частотной в пространственную область (восстановим из спектра)
blurred_img = real(ifft(blurred_img_ft))

# Также получим фильтр в пространственной области
blurring = ifft(fftshift(blurring_ft));

Теперь выполним операцию обратной свертки и выведем три изображения, склеенных в одно:

• исходное изображение,
• изобржение с добавленным шумом,
• отфильтрованное изображение.

polished = wiener(blurred_img, img, noise, blurring)

[ Gray.( img ) Gray.(blurred_img) Gray.( polished ) ]

Мы видим исходное изображение, затем размытое изображение с добавочным шумом, и восстановленное изборажение.

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

img2 = channelview( testimage("livingroom") ) # Загрузим другое изображение
noise2 = rand(Float64, size(img)) # Другая модель шума с 

# Уберем шум с пролого изображения при помощи обратной свертки
polished2 = wiener( blurred_img, img2, noise2, blurring );

[ Gray.( img2 ) Gray.( polished2 ) ]

Изображение слева дало алгоритму представление о том, как должен выглядеть спектр очищенного изображения. А поскольку мы не моделировали размытие при создании фильтра noise2, результат получился менее четким.

Обратная свертка Richardson-Lucy

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

using Images, TestImages, Deconvolution, FFTW, ZernikePolynomials

img = channelview( testimage("cameraman") )

# Модели аберрации линзы
blurring = evaluatezernike( LinRange(-16,16,512), OSA.([12, 4, 0]), [1.0, -1.0, 2.0] )
blurring = fftshift(blurring)
blurring = blurring ./ sum(blurring)

blurred_img = fft(img) .* fft(blurring) |> ifft |> real

@time restored_img = lucy(blurred_img, blurring, iterations=1000)

[ Gray.(img)  Gray.( fftshift(blurring .* 255));
  Gray.(blurred_img) Gray.(restored_img) ]

 55.638180 seconds (36.03 k allocations: 35.187 GiB, 0.79% gc time)

Заключение

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