I'm testing Wiener filter in MATLAB to restore a blurred image. I was using conv2() to blur the original image. If I use the 'full' option for conv2(), everything works well. But when I change to 'same' or 'valid', suddenly a lot of artifacts appeared in the restored image and Wiener filter failed. Please see below for blurred image, restoration from 'full' convolution, restoration from 'same' convolution.
Here's my implementation of the Wiener filter:
% load image
img = rgb2gray(imread('cameraman.jpg'));
[W, H] = size(img);
dim = 300;
img_fft = fft2(img,dim,dim);
% create blur kernel
kernel = ones(5) / 25;
kernel_fft = fft2(kernel,dim,dim);
% The option here makes huge difference 'same'/'full'/'valid'
img_blur = conv2(img,kernel,'same');
img_blur_fft = fft2(img_blur,dim,dim);
% Wiener filtering
k = 1e-5;
kernel_fft_conj = conj(kernel_fft);
img_wiener_freq = kernel_fft_conj .* img_blur_fft ./ (kernel_fft .* kernel_fft_conj + k);
img_wiener_ifft = ifft2(img_wiener_freq);
img_wiener_ifft = img_wiener_ifft(1:W,1:H);
As in real life the blurred image never has the form from circular or full convolution, how can I properly implement Wiener filter so it doesn't depends on the border of the image?
You have chosen a very small value for the regularization parameter k
. This works fine in the 'full'
case because there is no noise and the input exactly matches the expected. However, in the 'same'
case the input does not exactly match. If you were to use a circular convolution (e.g. by multiplication in the Fourier domain) then you would also get an exact result.
The regularization parameter exists to prevent small deviations from blowing up and corrupting the whole output image.
I got reasonable results using k = 1e-2
(1e-1
leaves the image blurry, 1e-3
still shows a lot of artefacts throughout the image, but it is possible to further fine-tune this value, I have not put more effort in than that).
This is the code I used, there are some important differences:
img = imread('cameraman.tif');
% create blur kernel
kernel = ones(5) / 25;
% The option here makes huge difference 'same'/'full'/'valid'
img_blur = conv2(img,kernel,'same');
img_blur_fft = fft2(img_blur);
% Wiener filtering
kernel_fft = padarray(kernel,size(img)-size(kernel),0,'post');
kernel_fft = circshift(kernel_fft,-floor(size(kernel)/2));
kernel_fft = fft2(kernel_fft);
k = 1e-2;
kernel_fft_conj = conj(kernel_fft);
img_wiener_freq = kernel_fft_conj .* img_blur_fft ./ (kernel_fft .* kernel_fft_conj + k);
img_wiener_ifft = ifft2(img_wiener_freq);
Notice that I do not use the size parameters to the fft2
and ifft2
functions. It is always best to not pad images with zeros in this case.
The kernel
image was padded in a very specific way. The FFT assumes that the origin is in the top-left pixel of the input. fft2(kernel,dim,dim)
leads to padding the kernel to the right and bottom with zeros, but this leaves the kernel shifted with respect the FFT's origin. This shift causes the deconvolve image to be shifted as well (by only 2 pixels, it's hard to notice, but look at img_wiener_ifft-double(img)
for OP's code and this code to see this shift).