图像处理中重要的一换就是图像增强，而图像增强分为空域处理和频域处理，这里只讨论频域处理
大家知道频域处理方法，大致的步骤是：图像傅立叶变换到频域，经过传递函数即滤波器，再傅立叶逆变换，就可得到处理过的图像。
而其中关键的步骤就是滤波器的选择，作为低通滤波来说，可以有多种选择，基本来讲大致有理想低通滤波，巴特沃斯低通滤波，指数低通滤波，梯形低通滤波等等。
低通滤波的作用就是使图像平滑，滤掉高频的成分，把细节滤掉，结果会使图像模糊，去掉噪声
大多数的vc数字图像处理书中都有傅立叶变换编程的介绍，但是具体到滤波器的实现就很少了，所以我通过研究，实现了滤波器编程，下面把主要代码贴出，供大家参考
/****************************************************
BLPF()
参数：
hDIB为输入的DIB句柄
order为阶数
stopfre为截止频率
返回值：
成功为TRUE；失败为FALSE
说明：
本函数实现DIB位图的快速傅立叶变换,并对傅立叶变换图进行
直方图均衡化
****************************************************/
BOOL BLPF(HDIB hDIB,int order,double stopfre)
{
if (hDIB == NULL)
return FALSE;
// start wait cursor
WaitCursorBegin();
HDIB hDib = NULL;
HDIB hNewDib = NULL;
// we only convolute 24bpp DIB, so first convert DIB to 24bpp
WORD wBitCount = DIBBitCount(hDIB);
if (wBitCount != 24)
{
hNewDib = ConvertDIBFormat(hDIB, 24, NULL);
hDib = CopyHandle(hNewDib);
}
else
{
hNewDib = CopyHandle(hDIB);
hDib = CopyHandle(hDIB);
}
if (hNewDib == NULL && hDib == NULL)
{
WaitCursorEnd();
return FALSE;
}
// process!
LPBYTE lpSrcDIB = (LPBYTE)GlobalLock(hDib);
LPBYTE lpDIB = (LPBYTE)GlobalLock(hNewDib);
LPBYTE lpInput = FindDIBBits(lpSrcDIB);
LPBYTE lpOutput = FindDIBBits(lpDIB);
int nWidth = DIBWidth(lpSrcDIB);
int nHeight = DIBHeight(lpSrcDIB);
int w=1,h=1,wp=0,hp=0;
while(w*2<=nWidth)
{
w*=2;
wp++;
}
while(h*2<=nHeight)
{
h*=2;
hp++;
}
int x,y;
BYTE *lpPoints=new BYTE[nWidth*nHeight];
GetPoints(nWidth,nHeight,lpInput,lpPoints);
COMPLEX *TD=new COMPLEX[w*h];
COMPLEX *FD=new COMPLEX[w*h];
COMPLEX *FTD=new COMPLEX[w*h];
COMPLEX *TTD=new COMPLEX[w*h];
for(y=0;y<h;y++)
{
for(x=0;x<w;x++)
{
TD[x+w*y].re=Point(x,y);
TD[x+w*y].im=0;
}
}
/////FFT()
for(y=0;y<h;y++)
{
FFT(&TD[w*y],&FD[w*y],wp);
}
for(y=0;y<h;y++)
{
for(x=0;x<w;x++)
{
TD[y+h*x]=FD[x+w*y];
// TD[x+w*y]=FD[x*h+y];
}
}
for(x=0;x<w;x++)
{
FFT(&TD[x*h],&FD[x*h],hp);
}
memset(lpPoints,0,nWidth*nHeight);
double m,d,H,mid1;
double mid2;
int i;
//////转置
for(y=0;y<h;y++)
{
for(x=0;x<w;x++)
{
FTD[x+w*y]=FD[y+x*h];
}
}
//////移到中心
for(y=0;y<h;y++)
{
for(x=0;x<w;x++)
{
FD[(x<w/2?x+w/2:x-w/2)+w*(y<h/2?y+h/2:y-h/2)]=FTD[x+w*y];
}
}
//////BLPF
for(y=0;y<h;y++)
{
for(x=0;x<w;x++)
{
mid1=(x-w/2)*(x-w/2)+(y-h/2)*(y-h/2);
d=sqrt(mid1);
for(i=1;i<=order;i++)
{
mid2=1.0;
mid2*=(d/stopfre);
mid2*=(d/stopfre);
}
// mid2=(d/stopfre)*(d/stopfre)*(d/stopfre)*(d/stopfre);
H=1/(1+mid2);
FD[x+w*y].re*=H;
FD[x+w*y].im*=H;
}
}
//////移回原位
for(y=0;y<h;y++)
{
for(x=0;x<w;x++)
{
FTD[(x<w/2?x+w/2:x-w/2)+w*(y<h/2?y+h/2:y-h/2)]=FD[x+w*y];
}
}
//////IDFT()
for(y=0;y<h;y++)
{
IFFT(&FTD[w*y],&TD[w*y],wp);
}
for(y=0;y<h;y++)
{
for(x=0;x<w;x++)
{
FTD[y+h*x]=TD[x+w*y];
// TD[x+w*y]=FD[x*h+y];
}
}
for(x=0;x<w;x++)
{
IFFT(&FTD[x*h],&TD[x*h],hp);
}
//////转置
for(y=0;y<h;y++)
{
for(x=0;x<w;x++)
{
TTD[x+w*y]=TD[y+x*h];
}
}
//////转成图像
for(y=0;y<h;y++)
{
for(x=0;x<w;x++)
{
m=TTD[x+w*y].re;
Point(x,y)=(BYTE)(m);
}
}
delete TD;
delete FD;
delete FTD;
delete TTD;
PutPoints(nWidth,nHeight,lpOutput,lpPoints);
delete lpPoints;
// recover
DWORD dwSize = GlobalSize(hDib);
memcpy(lpSrcDIB, lpDIB, dwSize);
GlobalUnlock(hDib);
GlobalUnlock(hNewDib);
if (wBitCount != 24)
{
hNewDib = ConvertDIBFormat(hDib, wBitCount, NULL);
lpSrcDIB = (LPBYTE)GlobalLock(hDIB);
lpDIB = (LPBYTE)GlobalLock(hNewDib);
dwSize = GlobalSize(hNewDib);
memcpy(lpSrcDIB, lpDIB, dwSize);
GlobalUnlock(hDIB);
GlobalUnlock(hNewDib);
}
else
{
lpSrcDIB = (LPBYTE)GlobalLock(hDIB);
lpDIB = (LPBYTE)GlobalLock(hDib);
dwSize = GlobalSize(hDib);
memcpy(lpSrcDIB, lpDIB, dwSize);
GlobalUnlock(hDIB);
GlobalUnlock(hDib);
}
// cleanup
GlobalFree(hDib);
GlobalFree(hNewDib);
// return
WaitCursorEnd();
return TRUE;
}
其中调用了快速傅立叶变换函数FFT()，和反变换函数IFFT()
/****************************************************
FFT()
参数：
TD为时域值
FD为频域值
power为2的幂数
返回值：
无
说明：
本函数实现快速傅立叶变换
****************************************************/
void FFT(COMPLEX * TD, COMPLEX * FD, int power)
{
int count;
int i,j,k,bfsize,p;
double angle;
COMPLEX *W,*X1,*X2,*X;
/*计算傅立叶变换点数*/
count=1<<power;
/*分配运算所需存储器*/
W=(COMPLEX *)malloc(sizeof(COMPLEX)*count/2);
X1=(COMPLEX *)malloc(sizeof(COMPLEX)*count);
X2=(COMPLEX *)malloc(sizeof(COMPLEX)*count);
/*计算加权系数*/
for(i=0;i<count/2;i++)
{
angle=-i*PI*2/count;
W[i].re=cos(angle);
W[i].im=sin(angle);
}
/*将时域点写入存储器*/
memcpy(X1,TD,sizeof(COMPLEX)*count);
/*蝶形运算*/
for(k=0;k<power;k++)
{
for(j=0;j<1<<k;j++)
{
bfsize=1<<(power-k);
for(i=0;i<bfsize/2;i++)
{
p=j*bfsize;
X2[i+p]=Add(X1[i+p],X1[i+p+bfsize/2]);
X2[i+p+bfsize/2]=Mul(Sub(X1[i+p],X1[i+p+bfsize/2]),W[i*(1<<k)]);
}
}
X=X1;
X1=X2;
X2=X;
}
/*重新排序*/
for(j=0;j<count;j++)
{
p=0;
for(i=0;i<power;i++)
{
if (j&(1<<i)) p+=1<<(power-i-1);
}
FD[j]=X1[p];
}
/*释放存储器*/
free(W);
free(X1);
free(X2);
}
/****************************************************
IFFT()
参数：
FD为频域值
TD为时域值
power为2的幂数
返回值：
无
说明：
本函数利用快速傅立叶变换实现傅立叶反变换
****************************************************/
void IFFT(COMPLEX * FD, COMPLEX * TD, int power)
{
int i, count;
COMPLEX *x;
/*计算傅立叶反变换点数*/
count=1<<power;
/*分配运算所需存储器*/
x=(COMPLEX *)malloc(sizeof(COMPLEX)*count);
/*将频域点写入存储器*/
memcpy(x,FD,sizeof(COMPLEX)*count);
/*求频域点的共轭*/
for(i=0;i<count;i++)
x[i].im = -x[i].im;
/*调用FFT*/
FFT(x, TD, power);
/*求时域点的共轭*/
for(i=0;i<count;i++)
{
TD[i].re /= count;
TD[i].im = -TD[i].im / count;
}
/*释放存储器*/
free(x);
}
这个程序代码是基于周长发的《精通visual c++图像处理编程》代码开发出来的，大家可以参考该书的各种变量定义等等。
希望大家指正，可以留言讨论哦
参考书目：
《精通visual c++图像处理编程（第三版）》周长发 电子工业出版社
《数字图像处理（第二版）》冈萨雷斯 电子工业出版社
《visual c++数字图像处理（第二版）》何斌 人民邮电出版社
《VC++图像处理程序设计（第二版）》杨淑英 清华大学出版社
《C++程序设计教程》钱能 清华大学出版社
