摘要：本文以VC++ 6.0為編程工具，講述了采取逆濾波和維納濾波兩種圖像恢復算法對退化圖像的恢復實現過程。

　　引言

　　圖像恢復技術是圖像處理領域一類重要的處理技術，與圖像增強等其他基本圖像處理技術類似，該技術也是以獲取視覺質量得到某種程度改善為目的的，所不同的是圖像恢復過程需要根據指定的圖像退化模型來完成，根據這個退化模型對在某種情況下退化或惡化了的退化圖像進行恢復，以獲取到原始的、未經過退化的原始圖像。換句話說，圖像恢復的處理過程實際是對退化圖像品質的提升，并通過圖像品質的提升來達到圖像在視覺上的改善。本文以VC++作為開發工具，講述了對退化圖像進行逆濾波和維納濾波處理算法。

　　逆濾波處理

　　對圖像進行恢復處理通常需要根據一定的圖像退化模型來進行，一個簡單的通用圖像退化模型可將圖像的退化過程模型化為一個作用在原始圖像f(x,y)上的退化系統H，作用結果與一個加性噪聲n(x,y)的聯合作用導致產生出了退化圖像g(x,y)，表現為數學形式為g(x,y)=H[f(x,y)]+n(x,y)。根據上述退化系統H可以從給定的退化圖像g(x,y)得到原始圖像f(x,y)的一個近似結果。逆濾波處理就是其中一種無約束恢復的圖像恢復技術，其恢復過程的數學形式可表示為F (u,v)=G(u,v)/H(u,v) (u,v=0,1,…,M-1),其中F(u,v)和G(u,v)分別為圖像f(x,y)和g(x,y)的頻域變換,H(u,v)可看作是一個濾波函數。由于圖像在退化過程中存在噪聲的干擾，因此通常情況下的濾波器往往不是正好的1/H(u,v),而是關于u和v的某個非線形的恢復轉移函數M(u,v)。經過以上的分析，圖像的退化和恢復過程（模型）大致可用下圖來表示：


　　一種簡便的恢復方法是在選取恢復轉移函數M(u,v) 時，如果u2+v2≤w2,則取值1/H(u,v)，否則為1。這樣處理雖然簡單，但是恢復后的圖像往往存在較明顯的振鈴現象，通常為了消除振鈴現象，以H(u,v)的值作為判據，如不大于d(0
　　由于恢復過程需要在頻域進行，因此需要通過二維傅立葉變換將圖像由空域變換到頻域。二維的傅立葉變換較一維傅立葉變換要復雜的多，一般采取連續2次運用一維離散快速傅立葉變換的方法來實現，即先沿f(x,y)的每一個x對y求變換再乘以N得到F(x,v)，完成第一步變換。然后再將得到的F(x,v)沿f(x,v)的每一個v對x求變換即可得到f(x,y)的最終變換F(u,v)，這兩步的數學表達式如下：

F(x,v)=N*[(1/N)* f(x,y)exp[-j2πvy/N]] (v=0,1,……,N-1) F(u,v)=(1/N)* F(x,v)exp[-j2πux/N] (u,v=0,1,……,N-1)

　　類似也可以得出二維離散傅立葉變換逆變換用一維變換計算的表達式：

F(x,v)= F(u,v)exp[j2πux/N] (x,y=0,1,……,N-1) f(x,y)=(1/N)* F(x,v)exp[j2πvy/N]] (y=0,1,……,N-1)

　　在分布進行一維傅立葉變換時，多采用"蝴蝶圖"的快速算法（詳見信號處理方面資料），其核心算法如下：

int N=(int)pow(2,M); file://N：序列長度（2的整數次冪） ReverseOrder(A,N); file://對空間序列進行倒序 for(int i=1;i<=M;i++){ 　int b=(int)pow(2,(i-1)); 　for(int j=0;j<=(b-1);j++) { 　　float p=(float)(pow(2,(M-i))*j*2.0*PI/(float)N); 　　for(int k=j;k<=(N-1);){ 　　　float tr=(float)(A[k+b].Re*cos(p)+A[k+b].Im*sin(p)); file://計算復數運算A*U 　　　float ti=(float)(A[k+b].Im*cos(p)-A[k+b].Re*sin(p)); 　　　A[k+b].Re=A[k].Re-tr; file://復數運算A-tr 　　　A[k+b].Im=A[k].Im-ti; 　　　A[k].Re+=tr; file://復數運算A+tr 　　　A[k].Im+=ti; 　　　k+=b*2; 　　} 　} }

　　傅立葉逆變換的同傅立葉變換比較相似，只是在計算exp[j2πvy/N]時同正變換有符號的區別，以及在計算完成后逆變換需要將值除以N，因此不難寫出一維傅立葉逆變換的實現代碼。在進行二維傅立葉變換將圖像由空域變換到頻域之前，首先需要通過補0的手段將點數非2的整數次冪的非正方型網格采樣構造為一個長寬均為2的整數次冪的正方型網格：

int WM=(int)(log(W)/log(2)+1.0f); file://計算圖像寬應為2的多少次冪 int HM=(int)(log(H)/log(2)+1.0f); file://計算圖像高應為2的多少次冪 WM=HM=max(WM,HM); file://取二者大值 int WN=(int)pow(2,WM); file://構造網格寬度 int HN=(int)pow(2,HM); file://構造網格高度 for{int i=0;i;for(int j=0;j 　　if(i 　　　U[i*WN*3+j].Re=D[i*W*3+j]; file://D為圖像序列 　　U[i*WN*3+j].Im=0.0f; 　　}else file://缺位補0 　　U[i*WN*3+j].Re=U[i*WN*3+j].Im=0.0f; 　} }

　　預處理完畢后，可對構造網格的每一列分別進行一維快速傅立葉變換，并將結果存放在原位置，結果乘以N，完成第一步的轉換，求得F(x,v)：

for(i=0;i 　for(int j=0;j 　　UH[j].Re=U[j*WN*3+i].Re; 　　UH[j].Im=U[j*WN*3+i].Im; 　} 　DFT_FFT(UH,HM); file://對UH進行快速離散傅立葉變換 　for(j=0;j 　　U[j*WN*3+i].Re=HN*UH[j].Re; file://N=HN 　　U[j*WN*3+i].Im=HN*UH[j].Im; 　} }

　　隨即對構造網格的每一行進行傅立葉變換，得到最終的變換結果F(u,v)：

for(i=0;i 　for(int k=0;k<3;k++){ file://對24位位圖的R、G、B三分量均各自進行變換 　　for(int j=0;j 　　　UW[j].Re=U[i*WN*3+j*3+k].Re; 　　　UW[j].Im=U[i*WN*3+j*3+k].Im; 　　} 　　DFT_FFT(UW,WM); file://對UW序列進行快速離散傅立葉變換 　　for(j=0;j 　　　U[i*WN*3+j*3+k].Re=UW[j].Re; 　　　U[i*WN*3+j*3+k].Im=UW[j].Im; 　　} 　} }

　　至于二維傅立葉逆變換則基本上是上述過程的逆過程，在此就不再贅述。根據逆濾波圖像恢復的設計方案，先通過前面的二維傅立葉變換將退化圖像g(x,y)從空域變換到頻域得到G(u,v)，然后在頻域經過恢復轉移函數M(u,v)的恢復處理并經過二維傅立葉逆變換將結果由頻域轉換回空域，就可得到經過恢復處理的近似原始圖像：

…… dsp.DFT_2D_FFT(m_cpBuffer+54,m_nWidth,m_nHeight,U); file://進行二維傅立葉變換 for(int i=0;i 　for(int j=0;j 　　int k=(int)(j/3); 　　D1=(float)sqrt(i*i+k*k); 　　H=1.0f/(1+(D1/D0)*(D1/D0)); file://H(u,v)=1/(1+(u2+v2)/D02)) 　　if(H>0.45f){ file://閥值 d取0.45 　　　U[i*3*WN+j].Re/=H; file://在頻域與M(u,v)相乘 　　　U[i*3*WN+j].Im/=H; 　　}else{ 　　　U[i*3*WN+j].Re*=0.6f; file://如未超過閥值則M(u,v)取常數k=0.6 　　　U[i*3*WN+j].Im*=0.6f; 　　} 　} } dsp.DFT_2D_IFFT(m_cpBuffer+54,m_nWidth,m_nHeight,U); file://進行傅立葉逆變換

　　這里的逆濾波處理算法采用的是經過改進的恢復轉移函數M(u,v)，因此恢復后的圖像不會出現振鈴現象。以標準檢測圖像Lina為處理對象應用以上恢復處理算法，效果如下圖所示。其中間圖像為未經過改進的簡單算法，在胳膊和臉部存在較明顯的振鈴現象，而采取了改進措施的圖像則沒有任何振鈴現象出現，圖像得到了較好的恢復。

維納濾波處理

　　維納（Wiener）濾波是對退化圖像進行恢復處理的另一種常用算法，是一種有約束的恢復處理方法，其采用的維納濾波器是一種最小均方誤差濾波器，其數學形式比較復雜：

F(u,v)=[(1/H(u,v))*(|H(u,v)|2)/(|H(u,v)|2+s*[Sn(u,v)/Sf(u,v)])]*G(u,v)

　　當s為1時，上式就是普通的維納濾波；如果s為變量，則為參數維納濾波，如果沒有噪聲干擾，即Sn(u,v)=0時，上式實際就是前面的逆濾波。從其數學形式可以看出：維納濾波比逆濾波在對噪聲的處理方面要強一些。以上只是理論上的數學形式，在進行實際處理時，往往不知道噪聲函數Sn(u,v)和Sf(u,v)的分布情況，因此在實際應用時多用下式進行近似處理：

F(u,v)=[(1/H(u,v))* (|H(u,v)|2)/(|H(u,v)|2+K)]*G(u,v)

　　其中K是一個預先設定的常數。由此可以寫出維納濾波的實現代碼：

…… float K=0.05f; file://預先設定常數K dsp.DFT_2D_FFT(m_cpBuffer+54,m_nWidth,m_nHeight,U); file://轉換到頻域 for(int i=0;i 　for(int j=0;j 　　int k=(int)(j/3); 　　D1=(float)sqrt(i*i+k*k); 　　float H=1.0f/(1+(D1/D0)*(D1/D0));//H(u,v)= 1/(1+(u2+v2)/D02)) 　　U[i*3*WN+j].Re=(U[i*3*WN+j].Re*H)/(H*H+K); file://維納濾波 　　U[i*3*WN+j].Im=(U[i*3*WN+j].Im*H)/(H*H+K); 　} } dsp.DFT_2D_IFFT(m_cpBuffer+54,m_nWidth,m_nHeight,U);//返回到空域

　　對經過退化的Lina圖像應用維納濾波處理，可得到如右圖所示的恢復效果圖。由于維納濾波在進行恢復時對噪聲進行了處理，因此其恢復效果要比逆濾波要好，尤其是退化圖像的噪聲干擾較強時效果更為明顯。


　　小結

　　本文對比較常用的兩種圖像恢復算法逆濾波和維納濾波的實現過程作了較為詳細的講述，通過對圖像質量較低的退化圖像應用上述算法可以使圖像質量得到一定程度的改善，在視覺上可以得到較好的改觀。類似的圖像恢復算法還有有約束最小平方恢復算法等多種，應視具體情況靈活選擇合適的算法以獲取最佳的恢復效果。本文所述程序在Windows 98下，由Microsoft Visual C++ 6.0編譯通過。