VC＋＋实现对退化图像的恢复

摘要：本文以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为处理对象应用以上恢复处理算法，效果如下图所示。其中间图像为未经过改进的简单算法，在胳膊和脸部存在较明显的振铃现象，而采取了改进措施的图像则没有任何振铃现象出现，图像得到了较好的恢复。