1.图像的傅里叶频谱的意义
之前的博文其实已经归纳过这方面的内容了。我们常用的图像平滑处理,其实就是一个低通滤波,一定程度上去除高频信号,可以使得图像变得柔和(也就是平滑)。但是,在去除周期性噪声时候,空间域内的滤波(卷积)就不是那么好操作了。所以,这里时候,无论是理解起来方便,还是其他原因,都需要在频域内进行滤波。
详细的叙述还是在下面的博文里面啦!!!!
2. 傅里叶频谱的计算
这部分的内容,主要就是使用OpenCV自带的函数
void cvDFT( const CvArr* src, CvArr* dst, int flags, int nonzero_rows=0 )
去求取图像的傅里叶变换。这里,输出结果CvArr* dst由两个通道组成,分别代表了实部与虚部。我们再根据如下算式,就可以得到傅里叶频谱了。
|F(u,v)|=R2(u,v)+F2(u,v)‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾√2
我自己也参考了很多人的代码,然后实现的代码如下。
IplImage* fft2(IplImage* image_input)
{
int dftWidth = getOptimalDFTSize(image_input->width);
int dftHeight = getOptimalDFTSize(image_input->height);
IplImage* image_padded = cvCreateImage(cvSize(dftWidth,dftHeight),
IPL_DEPTH_8U,
1);
cvCopyMakeBorder( image_input, image_padded, cvPoint(0,0), IPL_BORDER_CONSTANT,cvScalarAll(0));
IplImage *image_Re =0 , *image_Im = 0, *image_Fourier = 0;
image_Re = cvCreateImage(cvSize(dftWidth,dftHeight),IPL_DEPTH_64F,1);
image_Im = cvCreateImage(cvSize(dftWidth,dftHeight),IPL_DEPTH_64F,1);
image_Fourier = cvCreateImage(cvSize(dftWidth,dftHeight),IPL_DEPTH_64F,2);
cvConvertScale(image_padded,image_Re);
cvZero(image_Im);
cvMerge(image_Re,image_Im,0,0,image_Fourier);
cvDFT(image_Fourier,image_Fourier,CV_DXT_FORWARD);
cvSplit(image_Fourier,image_Re,image_Im,0,0);
cvPow(image_Re,image_Re,2.0);
cvPow(image_Im,image_Im,2.0);
cvAdd(image_Re,image_Im,image_Re);
cvPow(image_Re,image_Re,0.5);
cvAddS(image_Re,cvScalar(1),image_Re );
cvLog (image_Re,image_Re);
IplImage *Fourier = cvCreateImage(cvSize(dftWidth,dftHeight),IPL_DEPTH_64F,1);
cvZero(image_Fourier);
int cx = image_Re->width/2;
int cy = image_Re->height/2;
cvSetImageROI(image_Re,cvRect( 0, 0,cx,cy));
cvSetImageROI( Fourier,cvRect(cx,cy,cx,cy));
cvAddWeighted(image_Re,1,Fourier,0,0,Fourier);
cvSetImageROI(image_Re,cvRect(cx,cy,cx,cy));
cvSetImageROI( Fourier,cvRect( 0, 0,cx,cy));
cvAddWeighted(image_Re,1,Fourier,0,0,Fourier);
cvSetImageROI(image_Re,cvRect(cx, 0,cx,cy));
cvSetImageROI( Fourier,cvRect( 0,cy,cx,cy));
cvAddWeighted(image_Re,1,Fourier,0,0,Fourier);
cvSetImageROI(image_Re,cvRect( 0,cy,cx,cy));
cvSetImageROI( Fourier,cvRect(cx, 0,cx,cy));
cvAddWeighted(image_Re,1,Fourier,0,0,Fourier);
cvResetImageROI(image_Re);
cvResetImageROI( Fourier);
cvNormalize(Fourier,Fourier,1,0,CV_C,NULL);
return(Fourier);
}
- 1
- 2
- 3
- 4
- 5
- 6
- 7
- 8
- 9
- 10
- 11
- 12
- 13
- 14
- 15
- 16
- 17
- 18
- 19
- 20
- 21
- 22
- 23
- 24
- 25
- 26
- 27
- 28
- 29
- 30
- 31
- 32
- 33
- 34
- 35
- 36
- 37
- 38
- 39
- 40
- 41
- 42
- 43
- 44
- 45
- 46
- 47
- 48
- 49
- 50
- 51
- 52
- 53
- 54
- 55
- 56
- 57
- 58
- 59
- 60
- 61
- 62
- 63
- 64
- 65
- 66
- 67
- 68
- 69
- 70
- 71
- 72
- 73
- 74
- 75
- 76
- 77
- 78
- 1
- 2
- 3
- 4
- 5
- 6
- 7
- 8
- 9
- 10
- 11
- 12
- 13
- 14
- 15
- 16
- 17
- 18
- 19
- 20
- 21
- 22
- 23
- 24
- 25
- 26
- 27
- 28
- 29
- 30
- 31
- 32
- 33
- 34
- 35
- 36
- 37
- 38
- 39
- 40
- 41
- 42
- 43
- 44
- 45
- 46
- 47
- 48
- 49
- 50
- 51
- 52
- 53
- 54
- 55
- 56
- 57
- 58
- 59
- 60
- 61
- 62
- 63
- 64
- 65
- 66
- 67
- 68
- 69
- 70
- 71
- 72
- 73
- 74
- 75
- 76
- 77
- 78
从这里开始,还是简单的分析一下代码吧。
int dftWidth = getOptimalDFTSize(image_input->width);
int dftHeight = getOptimalDFTSize(image_input->height);
IplImage* image_padded = cvCreateImage(cvSize(dftWidth,dftHeight),
IPL_DEPTH_8U,
1);
cvCopyMakeBorder( image_input, image_padded, cvPoint(0,0), IPL_BORDER_CONSTANT,cvScalarAll(0));
这里参考了文献[2]中的说法,在尺寸数为2,3,5的倍数的场合,计算的速度是最快的。所以使用函数getOptimalDFTSize()
来寻找最匹配的尺寸,然后再同伙cvCopyMakeBorder()
进行多余部分的填充,这里选的配置是将图放在从点(0,0)开始的位置,其余不足的地方,用0进行填充。
IplImage *image_Re =0 , *image_Im = 0, *image_Fourier = 0;
image_Re = cvCreateImage(cvSize(dftWidth,dftHeight),IPL_DEPTH_64F,1);
image_Im = cvCreateImage(cvSize(dftWidth,dftHeight),IPL_DEPTH_64F,1);
image_Fourier = cvCreateImage(cvSize(dftWidth,dftHeight),IPL_DEPTH_64F,2);
//image_Re <-
cvConvertScale(image_padded,image_Re);
//image_Im <-
cvZero(image_Im);
//image_Fourier[0] <-
//image_Fourier[1] <-
cvMerge(image_Re,image_Im,0,0,image_Fourier);
cvDFT(image_Fourier,image_Fourier,CV_DXT_FORWARD);
//image_Fourier[0]
//image_Fourier[1]
cvSplit(image_Fourier,image_Re,image_Im,0,0);
- 1
- 2
- 3
- 4
- 5
- 6
- 7
- 8
- 9
- 10
- 11
- 12
- 13
- 14
- 15
- 16
- 17
- 18
- 19
- 1
- 2
- 3
- 4
- 5
- 6
- 7
- 8
- 9
- 10
- 11
- 12
- 13
- 14
- 15
- 16
- 17
- 18
- 19
其实这里的很好理解的,将填充到最适尺寸的图像赋值给image_Re,将image_Im赋值为0。让后将这两层图复制到image_Fourier的两个通道里,然后使用函数cvDFT()
进行傅里叶变换。得到结果还是存在于image_Fourier的两个通道里,分别代表实部与虚部,然后通过cvSplit()
将其抽出到image_Re与image_Im里。
cvPow(image_Re,image_Re,2.0);
cvPow(image_Im,image_Im,2.0);
cvAdd(image_Re,image_Im,image_Re);
cvPow(image_Re,image_Re,0.5);
cvAddS(image_Re,cvScalar(1),image_Re );
cvLog (image_Re,image_Re);
以上代码,实现了以下计算。
|F(u,v)|=R2(u,v)+F2(u,v)‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾√2
还有就是进行了一个对数变换,这个也没的说,看傅里叶频谱的标配操作。
// |-----|-----| |-----|-----|
// | 1 | 3 | | 4 | 2 |
// |-----|-----| ---> |-----|-----|
// | 2 | 4 | | 3 | 1 |
// |-----|-----| |-----|-----|
IplImage *Fourier = cvCreateImage(cvSize(dftWidth,dftHeight),IPL_DEPTH_64F,1);
cvZero(image_Fourier);
int cx = image_Re->width/2;
int cy = image_Re->height/2;
cvSetImageROI(image_Re,cvRect( 0, 0,cx,cy)); // 1
cvSetImageROI( Fourier,cvRect(cx,cy,cx,cy)); // 4
cvAddWeighted(image_Re,1,Fourier,0,0,Fourier);
cvSetImageROI(image_Re,cvRect(cx,cy,cx,cy)); // 4
cvSetImageROI( Fourier,cvRect( 0, 0,cx,cy)); // 1
cvAddWeighted(image_Re,1,Fourier,0,0,Fourier);
cvSetImageROI(image_Re,cvRect(cx, 0,cx,cy)); // 3
cvSetImageROI( Fourier,cvRect( 0,cy,cx,cy)); // 2
cvAddWeighted(image_Re,1,Fourier,0,0,Fourier);
cvSetImageROI(image_Re,cvRect( 0,cy,cx,cy)); // 2
cvSetImageROI( Fourier,cvRect(cx, 0,cx,cy)); // 3
cvAddWeighted(image_Re,1,Fourier,0,0,Fourier);
cvResetImageROI(image_Re);
cvResetImageROI( Fourier);
cvNormalize(Fourier,Fourier,1,0,CV_C,NULL);
return(Fourier);
- 1
- 2
- 3
- 4
- 5
- 6
- 7
- 8
- 9
- 10
- 11
- 12
- 13
- 14
- 15
- 16
- 17
- 18
- 19
- 20
- 21
- 22
- 23
- 24
- 25
- 26
- 27
- 28
- 29
- 30
- 31
- 32
- 33
- 34
- 1
- 2
- 3
- 4
- 5
- 6
- 7
- 8
- 9
- 10
- 11
- 12
- 13
- 14
- 15
- 16
- 17
- 18
- 19
- 20
- 21
- 22
- 23
- 24
- 25
- 26
- 27
- 28
- 29
- 30
- 31
- 32
- 33
- 34
其实重头戏在这里,这里需要一个交换操作。至于为何所求得的傅里叶频谱为什么需要交换的原因是,这个代码求得的结果其实是范围
[0,2π]
内的傅里叶变换。而我们所需要的是
[−π,π]
内的结果。详细的原因,还是在我以前的博客里有说明。
[数字图像处理]频域滤波(1)–基础与低通滤波器
这里,我使用了ROI操作与cvAddWeighted()
函数进行了实现。其运行的结果如下所示。
恩,基本可以看出来,直流分量也被我移动到了中心,以上代码实现了傅里叶频谱的计算与显示。
3. 不用交换操作的代码
使用MATLAB去求取尺寸为
M×N
图像
f
的傅里叶频谱时候,通常会用fft2(f,2*M,2*N)
。使用此函数求得的福利叶变换,其实还是
[0,2π]
范围内的傅里叶变换。要想使得傅里叶频谱的DC分量位于中心的话,其实还是需要一些别的操作的。冈萨雷斯的《数字图像处理》一书中,对于这个问题,是利用了傅里叶变换的评议特性,即
f(x,y)=f(x,y)×(−1)x+y
然后再对函数
f(x,y)
进行福利叶变换,所得到结果,就是所需要的是
[−π,π]
内傅里叶变换的结果。具体的原理与推导,还是参看冈萨雷斯的《数字图像处理》英文原本的p258页左右,中文译本的p148左右。当然,嫌麻烦的话,还可以看我的博文,博文中也简明推导了平移特性。
[数字图像处理]频域滤波(1)–基础与低通滤波器
为此,实现的代码变为了如下形式。
IplImage* fft2_New(IplImage* image_input)
{
int dftWidth = getOptimalDFTSize(image_input->width);
int dftHeight = getOptimalDFTSize(image_input->height);
cout<< " Width" << image_input->width << " " << dftWidth << "\n";
cout<< "Height" << image_input->height << " " << dftHeight << "\n";
IplImage* image_padded = cvCreateImage(cvSize(dftWidth,dftHeight),
IPL_DEPTH_8U,
1);
cvCopyMakeBorder( image_input, image_padded, cvPoint(0,0), IPL_BORDER_CONSTANT,cvScalarAll(0));
IplImage *image_Re =0 , *image_Im = 0, *image_Fourier = 0;
image_Re = cvCreateImage(cvSize(dftWidth,dftHeight),IPL_DEPTH_64F,1);
image_Im = cvCreateImage(cvSize(dftWidth,dftHeight),IPL_DEPTH_64F,1);
image_Fourier = cvCreateImage(cvSize(dftWidth,dftHeight),IPL_DEPTH_64F,2);
double pixel;
for(int y=0;y<image_padded->height;y++)
{
for(int x=0;x<image_padded->width;x++)
{
pixel = cvGetReal2D(image_padded,x,y);
pixel = ((x+y)%2 == 0)?(pixel):((-1)*pixel);
cvSetReal2D(image_Re,x,y,pixel);
}
}
cvZero(image_Im);
cvMerge(image_Re,image_Im,0,0,image_Fourier);
cvDFT(image_Fourier,image_Fourier,CV_DXT_FORWARD);
cvSplit(image_Fourier,image_Re,image_Im,0,0);
cvPow(image_Re,image_Re,2.0);
cvPow(image_Im,image_Im,2.0);
cvAdd(image_Re,image_Im,image_Re);
cvPow(image_Re,image_Re,0.5);
cvAddS(image_Re,cvScalar(1),image_Re );
cvLog (image_Re,image_Re);
cvNormalize(image_Re,image_Re,1,0,CV_C,NULL);
return(image_Re);
}
- 1
- 2
- 3
- 4
- 5
- 6
- 7
- 8
- 9
- 10
- 11
- 12
- 13
- 14
- 15
- 16
- 17
- 18
- 19
- 20
- 21
- 22
- 23
- 24
- 25
- 26
- 27
- 28
- 29
- 30
- 31
- 32
- 33
- 34
- 35
- 36
- 37
- 38
- 39
- 40
- 41
- 42
- 43
- 44
- 45
- 46
- 47
- 48
- 49
- 50
- 51
- 52
- 53
- 54
- 55
- 56
- 57
- 1
- 2
- 3
- 4
- 5
- 6
- 7
- 8
- 9
- 10
- 11
- 12
- 13
- 14
- 15
- 16
- 17
- 18
- 19
- 20
- 21
- 22
- 23
- 24
- 25
- 26
- 27
- 28
- 29
- 30
- 31
- 32
- 33
- 34
- 35
- 36
- 37
- 38
- 39
- 40
- 41
- 42
- 43
- 44
- 45
- 46
- 47
- 48
- 49
- 50
- 51
- 52
- 53
- 54
- 55
- 56
- 57
在这里,由于考虑到计算的原因,我将
f(x,y)=f(x,y)×(−1)x+y
这个计算,用下面代码去进行了实现。
for(int y=0;y<image_padded->height;y++)
{
for(int x=0;x<image_padded->width;x++)
{
pixel = cvGetReal2D(image_padded,x,y);
pixel = ((x+y)%2 == 0)?(pixel):((-1)*pixel);
cvSetReal2D(image_Re,x,y,pixel);
}
}
其实也就相当于,
x+y
是偶数的话
f(x,y)
不变,
x+y
是奇数的话
f(x,y)
变为负。用这样一个简单的判断去实现。其运行的结果,如下所示。
原图的傅里叶频谱
使用
11×11
高斯滤波器后的傅里叶频谱
从实验结果看来,可以看出以下两点
- 对于原图而言两个函数的结果基本一致,两个函数都得到了正确的结果。
- 使用平滑处理后,频谱的高频成分明显变小,对于空间域的图像而言,图像变得模糊了。