Removing noise from an image using inverse convolution¶
Let's apply algorithms for noise and blur removal on image using the library Deconvolution.
Problem description¶
Deconvolution or inverse convolution (deconv) is an operation, inverse to convolution, that allows us to retrieve the original signal (or reconstruct the image) under a known or estimated filter model. In our scenario, when we will be doing image enhancement, by filter we mean describing the process that introduced interference into our image (blurring, blurring, noise, etc.).
The characteristics of this process are usually difficult to know in advance, and the inverse convolution operation itself is very sensitive to errors in the filter coefficients. This should be taken into account when creating a procedure for estimating the parameters of the filter model.
Here are the libraries we will need:
Pkg.add( ["Images", "TestImages", "Deconvolution", "FFTW", "ZernikePolynomials"] )
Wiener's inverse convolution¶
Let's load a test image and add noise to it. By replacing the function textimage("image") with load("image.png") you can open any image from the disc.
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));
Now let's perform the reverse convolution operation and output three images glued into one:
- the original image,
- the image with added noise,
- the filtered image.
polished = wiener(blurred_img, img, noise, blurring)
[ Gray.( img ) Gray.(blurred_img) Gray.( polished ) ]
*We see the original image, then the blurred image with added noise, and the reconstructed image.
Wiener convolution also works when we do not have the original image and filter model (that is, in a typical realistic scenario). Inverse convolution can be performed with only a rough understanding of the power spectra of the original signal and the filter.
img2 = channelview( testimage("livingroom") ) # Загрузим другое изображение
noise2 = rand(Float64, size(img)) # Другая модель шума с
# Уберем шум с пролого изображения при помощи обратной свертки
polished2 = wiener( blurred_img, img2, noise2, blurring );
[ Gray.( img2 ) Gray.( polished2 ) ]
The image on the left gave the algorithm an idea of what the spectrum of the cleaned image should look like. And since we did not simulate blurring when creating the filter noise2, the result is less clear.
Richardson-Lucy inverse convolution¶
As a second example, we use the function lucy, which implements the Richardson-Lucy inverse convolution method over an image that has been blurred using an aberration model.
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) ]
Conclusion¶
We have successfully tested a couple of algorithms for image restoration and got good results.