{"id":180,"date":"2007-11-02T07:00:25","date_gmt":"2007-11-02T11:00:25","guid":{"rendered":"https:\/\/blogs.mathworks.com\/steve\/2007\/11\/02\/image-deblurring-wiener-filter\/"},"modified":"2019-10-23T15:25:29","modified_gmt":"2019-10-23T19:25:29","slug":"image-deblurring-wiener-filter","status":"publish","type":"post","link":"https:\/\/blogs.mathworks.com\/steve\/2007\/11\/02\/image-deblurring-wiener-filter\/","title":{"rendered":"Image deblurring – Wiener filter"},"content":{"rendered":"
\r\n \r\n

I'd like to welcome back guest blogger Stan Reeves, professor of Electrical and Computer Engineering at Auburn University.\r\n Stan will be writing a few blogs here about image deblurring.<\/i><\/p>\r\n <\/introduction>\r\n

In my last blog<\/a>, I looked at image deblurring using an inverse filter and some variations. The inverse filter does a terrible job due to\r\n the fact that it divides in the frequency domain by numbers that are very small, which amplifies any observation noise in\r\n the image. In this blog, I'll look at a better approach, based on the Wiener filter.\r\n <\/p>\r\n

I will continue to assume a shift-invariant blur model with independent additive noise, which is given by<\/p>\r\n

y(m,n) = h(m,n)*x(m,n) + u(m,n)<\/p>\r\n

where * is 2-D convolution, h(m,n) is the point-spread function (PSF), f(m,n) is the original image, and u(m,n) is noise.<\/p>\r\n

The Wiener filter can be understood better in the frequency domain. Suppose we want to design a frequency-domain filter G(k,l)\r\n so that the restored image is given by\r\n <\/p>\r\n

<\/p>\r\n

We can choose G(k,l) so that we minimize<\/p>\r\n

<\/p>\r\n

E[] is the expected value of the expression. We assume that both the noise and the signal are random processes and are independent\r\n of one another. The minimizer of this expression is given by\r\n <\/p>\r\n

<\/p>\r\n

<\/p>\r\n

<\/p>\r\n

This filter gives the minimum mean-square error estimate of X(k,l). Now that sounds really great until we realize that we\r\n must supply the signal and noise power spectra (or equivalently the autocorrelation functions of the signal and noise).\r\n <\/p>\r\n

The noise power spectrum is fairly easy. We usually assume white noise, which makes the noise power spectrum a constant.\r\n But how do we determine what this constant is? If\r\n <\/p>\r\n

<\/p>\r\n

then the noise power spectrum is given by<\/p>\r\n

<\/p>\r\n

for an MxN image. This noise variance may be known based on knowledge of the image acquisition process or may be estimated\r\n from the local variance of a smooth region of the image.\r\n <\/p>\r\n

The signal power spectrum is a little more challenging in principle, since it is not flat. However, we have two factors working\r\n in our favor: 1) most images have fairly similar power spectra, and 2) the Wiener filter is insensitive to small variations\r\n in the signal power spectrum.\r\n <\/p>\r\n

Consider two very different images -- the cameraman and house.<\/p>

cam = im2double(imread('cameraman.tif'<\/span>));\r\nhouse_url = 'https:\/\/blogs.mathworks.com\/images\/steve\/180\/house.tif'<\/span>;\r\nhouse = im2double(imread(house_url));\r\n\r\nsubplot(1,2,1)\r\nimshow(cam)\r\ncolormap(gray(256))\r\ntitle('cameraman'<\/span>)\r\nsubplot(1,2,2)\r\nimshow(house)\r\ntitle('house'<\/span>)<\/pre> 
We show the actual (log) spectrum for these two images.<\/p>
cf = abs(fft2(cam)).^2;\r\nhf = abs(fft2(house)).^2;\r\nsubplot(1,2,1)\r\nsurf([-127:128]\/128,[-127:128]\/128,log(fftshift(cf)+1e-6))\r\nshading interp<\/span>, colormap gray<\/span>\r\ntitle('cameraman power spectrum'<\/span>)\r\nsubplot(1,2,2)\r\nsurf([-127:128]\/128,[-127:128]\/128,log(fftshift(hf)+1e-6))\r\nshading interp<\/span>, colormap gray<\/span>\r\ntitle('house power spectrum'<\/span>)<\/pre> 
Let's see what happens if we restore the cameraman with the actual power spectrum.<\/p>
h = fspecial('disk'<\/span>,4);\r\ncam_blur = imfilter(cam,h,'circular'<\/span>);\r\n% 40 dB PSNR<\/span>\r\nsigma_u = 10^(-40\/20)*abs(1-0);\r\ncam_blur_noise = cam_blur + sigma_u*randn(size(cam_blur));\r\n\r\ncam_wnr = deconvwnr(cam_blur_noise,h,numel(cam)*sigma_u^2.\/cf);\r\nsubplot(1,1,1)\r\nimshow(cam_wnr)\r\ncolormap(gray(256))\r\ntitle('restored image with exact spectrum'<\/span>)<\/pre> 
For comparison purposes, we restore the cameraman using the power spectrum obtained from the house image.<\/p>
cam_wnr2 = deconvwnr(cam_blur_noise,h,numel(cam)*sigma_u^2.\/hf);\r\nimshow(cam_wnr2)\r\ncolormap(gray(256))\r\ntitle('restored image with house spectrum'<\/span>)<\/pre> 
Visually, the two are very similar in quality.  In terms of mean-square error (MSE), the former is better (lower), as the\r\n      theory predicts:\r\n   <\/p>
format short<\/span> e<\/span>\r\nmse1 = mean((cam(:)-cam_wnr(:)).^2)\r\nmse2 = mean((cam(:)-cam_wnr2(:)).^2)<\/pre>
\r\nmse1 =\r\n\r\n  2.9905e-003\r\n\r\n\r\nmse2 =\r\n\r\n  3.9191e-003\r\n\r\n<\/pre>
In both cases, the spectrum is concentrated in the low frequencies and falls away at higher frequencies.  A smoothed version\r\n      of the spectra would look even more similar.  Rather than using the power spectrum from a specific image, one can either average\r\n      a large number of images or use a simple model of the power spectrum or autocorrelation function.\r\n   <\/p>\r\n   
A common model for the image autocorrelation function is<\/p>\r\n   
 <\/p>\r\n   
 <\/p>\r\n   
 <\/p>\r\n   
We calculate these parameters for the cameraman, averaging over horizontal and vertical shifts to get a single parameter for\r\n      the correlation coefficient.  Then we calculate an autocorrelation function.\r\n   <\/p>
sigma2_x = var(cam(:))\r\nmean_x = mean(cam(:))\r\ncam_r = circshift(cam,[1 0]);\r\ncam_c = circshift(cam,[0 1]);\r\nrho_mat = corrcoef([cam(:); cam(:)],[cam_r(:); cam_c(:)])\r\nrho = rho_mat(1,2);\r\n[rr,cc] = ndgrid([-128:127],[-128:127]);\r\nr_x = sigma2_x*rho.^sqrt(rr.^2+cc.^2) + mean_x^2;\r\n\r\nsurf([-128:127],[-128:127],r_x)\r\naxis tight<\/span>\r\nshading interp<\/span>, camlight, colormap jet<\/span>\r\ntitle('image autocorrelation model approximation'<\/span>)<\/pre>
\r\nsigma2_x =\r\n\r\n  5.9769e-002\r\n\r\n\r\nmean_x =\r\n\r\n  4.6559e-001\r\n\r\n\r\nrho_mat =\r\n\r\n  1.0000e+000  9.4492e-001\r\n  9.4492e-001  1.0000e+000\r\n\r\n<\/pre> 
From this we calculate another restored image:<\/p>
cam_wnr3 = deconvwnr(cam_blur_noise,h,sigma_u^2,r_x);\r\nimshow(cam_wnr3)\r\ncolormap(gray(256))\r\ntitle('restored image using correlation model'<\/span>)\r\nmse3 = mean((cam(:)-cam_wnr3(:)).^2)<\/pre>
\r\nmse3 =\r\n\r\n  3.5255e-003\r\n\r\n<\/pre> 
The restored image from this method is better than the restored image using the house spectrum but not quite as good as the\r\n      one using the exact cameraman spectrum.  All in all, we can see that the exact noise-to-signal spectrum isn't necessary to\r\n      yield good results.\r\n   <\/p>\r\n   
- by Stan Reeves, Department of Electrical and Computer Engineering, Auburn University<\/i><\/p>