在【影象處理中的數學原理】專欄(該專欄中的文章已經結集出版,書名為《影象處理中的數學修煉》)之前的一篇文章中,我們曾經討論過一種“自適應影象降噪濾波器的設計與實現”。彼時,也曾經提過其中運用了維納濾波器的一些方法,但我們並未深入討論關於維納濾波的更多內容。本文作為這個系列中的一個續篇,繼續來深入研究著名的維納濾波,特別是其背後的數學原理。這也涉及到了限制性影象復原和非限制性影象復原的一些話題。

歡迎關注白馬負金羈的部落格 http://blog.csdn.net/baimafujinji,為保證公式、圖表得以正確顯示,強烈建議你從該地址上檢視原版博文。本部落格主要關注方向包括:數字影象處理、演算法設計與分析、資料結構、機器學習、資料探勘、統計分析方法、自然語言處理。

影象復原的逆濾波方法

如果H是一個線性移不變系統,那麼在時域中給出的退化過程可由如下公式給出:

g(x,y)=h(x,y)f(x,y)+η(x,y).其中,h(x,y)是退化函式在時域下的表示,運算子表示時域卷積。由卷積定理可知,時域上的卷積等同於頻域上的乘積,所以上式在頻域中的表述如下:
G(u,v)=H(u,v)F(u,v)+N(u,v)
其中的大寫字母項是之前公式裡對應項的傅立葉變換。退化函式H通常是指模糊、抖動等影響。
如果認為噪聲的影響很小,那麼上面的式子就可以寫成G(u,v)=H(u,v)F(u,v),假設H(u,v)是可逆的,那便可以得出F(u,v)=G(u,v)/H(u,v)。但實際上噪聲是在所難免的,因而只能設法求出F(u,v)的估計值F̂ (u,v),此時用H(u,v)F(u,v)+N(u,v)來替換G(u,v),則有:
F̂ (u,v)=H(u,v)F(u,v)+N(u,v)H(u,v)=F(u,v)+N(u,v)H(u,v)

這也就是逆濾波(inverse filter)的基本原理。從這個公式可以看出:首先,如果希望復原後的F̂ (u,v)F(u,v)儘量接近,就需要min||N(u,v)||;其次,如果噪聲比較大,特別是噪聲的影響大過退化函式的影響,那麼復原的效果就會很差。

不妨在MATLAB中做一些簡單的實驗。在後面我們實現了一個用於影象復原的維納濾波函式mydeconvwr(I, PSF, NSR),其中I是待處理的退化影象,PSF是退化函式(以矩陣形式給出),當NSR=0時,這個函式就變成了一個逆濾波器。下面來實驗一下,不加任何噪聲的情況下呼叫mydeconvwr函式來對影象進行復原。

I = im2double(imread('cameraman.tif'));
LEN = 21; THETA = 11;
PSF = fspecial('motion', LEN, THETA);
blurred = imfilter(I, PSF, 'conv', 'circular');

result1 = mydeconvwnr(blurred, PSF, 0);

subplot(1,3,1),imshow(I),title('original image');
subplot(1,3,2),imshow(blurred),title('blurred image');
subplot(1,3,3),imshow(result1),title('restored image');

執行上述程式碼,所得之結果如下圖所示。可見在不引入噪聲的情況下,逆濾波的影象復原效果非常好。


下面我們試著向退化的影象中加入一些高斯噪聲,然後再用逆濾波來對影象進行復原,實驗程式碼如下。
noise_mean = 0;
noise_var = 0.0001;
blurred_noisy = imnoise(blurred, 'gaussian', noise_mean, noise_var);
estimated_nsr = 0;
result2 = mydeconvwnr(blurred_noisy, PSF, estimated_nsr);
subplot(1,2,1),imshow(blurred_noisy ),title('blurred image with noise');
subplot(1,2,2),imshow(result2),title('resotred image without considering noise');

執行上述程式碼,所得之結果如下圖所示。可見在引入噪聲的情況下,逆濾波的影象復原效果非常不理想。

在MATLAB中實現維納濾波

維納濾波的基本公式是

F̂ (u,v)=[H(u,v)|H(u,v)|2+PN(u,v)/PS(u,v)]G(u,v)=[1H(u,v)|H(u,v)|2|H(u,v)|2+PN(u,v)/PS(u,v)]G(u,v)
其中H(u,v)表示退化函式,H(u,v)2=H(u,v)H(u,v)H(u,v)表示H(u,v)的複共軛;PN(u,v)=|N(u,v)|2表示噪聲的功率譜;PS(u,v)=|F(u相關文章