图像处理与傅里叶变换原理与运用

图像处理与傅里叶变换1背景傅里叶变换是一个非常复杂的理论，我们在图像处理中集中关注于其傅里叶离散变换离散傅立叶变换(Discrete Fourier Transform) 。

1.1离散傅立叶变换图象是由灰度（RGB ）组成的二维离散数据矩阵，则对它进行傅立叶变换是离散的傅立叶变换。

对图像数据f(x,y)（x=0,1,… ,Ｍ-1； y=0,1,… ,N-1）。

则其离散傅立叶变换定义可表示为：式中，u=0,1,…, M-1;v= 0,1,…, N-1 其逆变换为式中，x=0,1,…, M-1;y= 0,1,…, N-1在图象处理中，一般总是选择方形数据，即Ｍ＝Ｎ影像f(x,y)的振幅谱或傅立叶频谱： 相位谱：能量谱（功率谱） )1(2exp ),(1),(101∑∑-=-=⎥⎦⎤⎢⎣⎡⎪⎭⎫ ⎝⎛+-=M x N y N vy M uxi y x f MNv u F π)2(2exp ),(1),(101∑∑-=-=⎥⎦⎤⎢⎣⎡⎪⎭⎫ ⎝⎛+=M u N v N vy M uxi v u F MNy x f π),(),(),(22v u I v u R v u F +=[]),(/),(),(v u R v u I arctg v u =ϕ),(),(),(),(222v u I v u R v u F v u E +==1.2快速傅里叶变化可分离性的优点是二维的傅立叶变换或逆变换由两个连续的一维傅立叶变换变换来实现，对于一个影像f(x,y)，可以先沿着其每一列求一维傅立叶变换，再对其每一行再求一维变换正变化逆变换由于二维的傅立叶变换具有可分离性，故只讨论一维快速傅立叶变换。

正变换 逆变换由于计算机进行运算的时间主要取决于所用的乘法的次数。

按照上式进行一维离散由空间域向频率域傅立叶变换时，对于N 个F∑∑∑∑-=-=-=-=⎥⎦⎤⎢⎣⎡⨯⎥⎦⎤⎢⎣⎡=⎥⎦⎤⎢⎣⎡+=110101)(2exp ),(1)(2exp ),(1)(2exp ),(1),(N v N u N u N v N vy i v u F NN ux i v u F N N vy ux i v u F NNy x f πππ∑-=⎥⎦⎤⎢⎣⎡-=12exp )(1)(N x N ux i x f Nu F π∑∑∑∑-=-=-=-=⎥⎦⎤⎢⎣⎡-⨯⎥⎦⎤⎢⎣⎡-=⎥⎦⎤⎢⎣⎡+-=11101)(2exp ),(1)(2exp ),(1)(2exp ),(1),(N y N x N x N y N vy i y x f NN ux i y x f NN vy ux i y x f NNv u F πππ∑-=⎥⎦⎤⎢⎣⎡=12exp )(1)(N u N ux i u F Nx f π（u)值，中的每一个都要进行N 次运算，运算时间与N 2成正比。

1965年库里-图基（ Cooly-Tudey)提出将运算操作降到Nlog 2N 数量级的算法，即N 可以分解为一些较小整数的乘积，当N 为2的幂（即N=2P ，其中P 是整数时），效率最高，实现起来也最简单。

这就是快速傅立叶变换。

1.3关于基图像（频率矩形）由二维离散傅里叶反变换式。

可知，由于u 和v 均有0,1,…,N-1的N 个可能的取值，所以f(x,y)由N 2个频率分量组成，所以每个频率分量都与一个特定的(u,v)值相对应；且对于某个特定的(u,v)值来说，当(x,y)取遍所有可能的值(x=0，1，…，N-1；y=0，1，…，N-1）时，就可得到对应于该特定的(u,v)值的一幅基图像。

基图像可表示为。

所以，一幅图像的灰度平均值可由DFT 在原点处的值求得⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡-+-+-+--+++-+++=)])1()1((2exp[)]1)1((2exp[)]0)1((2exp[)])1(1(2exp[)]11(2exp[)]01(2exp[)])1(0(2exp[)]10(2exp[)]00(2exp[12,N v N u N j N vu N j N v u N j Nv N u j N v u j N v u j NvN u j N v u j N v u j N f vu πππππππππ证明周期性与共轭对称性F(u,v)=F(u+mM,v+Nn) F(u,v)=F(-u,-v)也就是说，图3.7的频谱图(a)和(b)实质上是函数的傅里叶频谱图。

1.3运用1.3.1频率滤波在频域中，频率越大说明原始信号变化速度越快；频率越小说明原始信号越平缓。

当频率为0时，表示直流信号，没有变化。

因此，频率的大小反应了信号的变化快慢。

高频分量解释信号的突变部分，而低频分量决定信号的整体形象。

在图像处理中，频域反应了图像在空域灰度变化剧烈程度，也就是图像灰度的变化速度，也就是图像的梯度大小。

对图像而言，图像的边缘部分是突变部分，变化较快，因此反应在频域上是高频分量；)2/,2/()1(),()(N v M u F y x f y x --⇔-⋅+图像的噪声大部分情况下是高频部分；图像平缓变化部分则为低频分量。

也就是说，傅立叶变换提供另外一个角度来观察图像，可以将图像从灰度分布转化到频率分布上来观察图像的特征。

书面一点说就是，傅里叶变换提供了一条从空域到频率自由转换的途径。

对图像处理而言，以下概念非常的重要：图像高频分量：图像突变部分；在某些情况下指图像边缘信息，某些情况下指噪声，更多是两者的混合；低频分量：图像变化平缓的部分，也就是图像轮廓信息高通滤波器：让图像使低频分量抑制，高频分量通过低通滤波器：与高通相反，让图像使高频分量抑制，低频分量通过带通滤波器：使图像在某一部分的频率信息通过，其他过低或过高都抑制模板运算与卷积定理在时域内做模板运算，实际上就是对图像进行卷积。

模板运算是图像处理一个很重要的处理过程，很多图像处理过程，比如增强/去噪（这两个分不清楚），边缘检测中普遍用到。

根据卷积定理，时域卷积等价与频域乘积。

因此，在时域内对图像做模板运算就等效于在频域内对图像做滤波处理。

比如说一个均值模板，其频域响应为一个低通滤波器；在时域内对图像作均值滤波就等效于在频域内对图像用均值模板的频域响应对图像的频域响应作一个低通滤波。

1.3.2图像去噪图像去噪就是压制图像的噪音部分。

因此，如果噪音是高频额，从频域的角度来看，就是需要用一个低通滤波器对图像进行处理。

通过低通滤波器可以抑制图像的高频分量。

但是这种情况下常常会造成边缘信息的抑制。

常见的去噪模板有均值模板，高斯模板等。

这两种滤波器都是在局部区域抑制图像的高频分量，模糊图像边缘的同时也抑制了噪声。

还有一种非线性滤波-中值滤波器。

中值滤波器对脉冲型噪声有很好的去掉。

因为脉冲点都是突变的点，排序以后输出中值，那么那些最大点和最小点就可以去掉了。

中值滤波对高斯噪音效果较差。

椒盐噪声：对于椒盐采用中值滤波可以很好的去除。

用均值也可以取得一定的效果，但是会引起边缘的模糊。

高斯白噪声：白噪音在整个频域的都有分布，好像比较困难。

冈萨雷斯版图像处理P185：算术均值滤波器和几何均值滤波器（尤其是后者）更适合于处理高斯或者均匀的随机噪声。

谐波均值滤波器更适合于处理脉冲噪声。

1.3.3图像增强有时候感觉图像增强与图像去噪是一对矛盾的过程，图像增强经常是需要增强图像的边缘，以获得更好的显示效果，这就需要增加图像的高频分量。

而图像去噪是为了消除图像的噪音，也就是需要抑制高频分量。

有时候这两个又是指类似的事情。

比如说，消除噪音的同时图像的显示效果显著的提升了，那么，这时候就是同样的意思了。

常见的图像增强方法有对比度拉伸，直方图均衡化，图像锐化等。

前面两个是在空域进行基于像素点的变换，后面一个是在频域处理。

我理解的锐化就是直接在图像上加上图像高通滤波后的分量，也就是图像的边缘效果。

对比度拉伸和直方图均衡化都是为了提高图像的对比度，也就是使图像看起来差异更明显一些，我想，经过这样的处理以后，图像也应该增强了图像的高频分量，使得图像的细节上差异更大。

同时也引入了一些噪音。

1.4实现在MATLAB中F=imread(filename);F=fft2(f,P,Q);%完成FFT变换FC=fftshift(F):%实现居中S=abs(F(或Fc));%取得傅里叶频谱f=real(ifft2(F));%实现傅里叶逆变换基于OPENCV库#include <stdio.h>#include <cv.h>#include <c xcore.h>#include <highgui.h>/**************************************************************** **********//src IPL_DEPTH_8U//dst IPL_DEPTH_64F***************************************************************** *********///傅里叶正变换void fft2(IplImage *src, IplImage *dst){ //实部、虚部IplImage *image_Re = 0, *im age_Im = 0, *Fourier = 0;// int i, j;im age_Re = cvCreateIm age(cvGetSize(src), IPL_DEPTH_64F, 1); //实部//Im aginary partim age_Im = cvCreateImage(cvGetSize(src), IPL_DEPTH_64F, 1); //虚部//2 channels (image_Re, image_Im)Fourier = cvCreateImage(cvGetSize(src), IPL_DEPTH_64F, 2);// Real part conversion from u8 to 64f (double)cvConvertScale(src, im age_Re, 1, 0);// Im aginary part (zeros)cvZero(im age_Im);// Join real and imaginary parts and stock them in Fourier imagecvMerge(im age_Re, image_Im, 0, 0, Fourier);// Application of the forward Fourier transformcvDFT(Fourier, dst, CV_DXT_FORWARD);cvReleaseImage(&image_Re);cvReleaseImage(&image_Im);cvReleaseImage(&Fourier);}/**************************************************************** **********//src IPL_DEPTH_64F//dst IPL_DEPTH_8U***************************************************************** *********/void fft2shift(IplImage *src, IplImage *dst){IplImage *image_Re = 0, *image_Im = 0;int nRow, nCol, i, j, cy, cx;double scale, shift, tmp13, tmp24;image_Re = cvCreateImage(cvGetSize(src), IPL_DEPTH_64F, 1);//Imaginary partimage_Im = cvCreateImage(cvGetSize(src), IPL_DEPTH_64F, 1);cvSplit( src, im age_Re, image_Im, 0, 0 );//具体原理见冈萨雷斯数字图像处理p123// Com pute the m agnitude of the spectrum Mag = sqrt(Re^2 + Im^2)//计算傅里叶谱cvPow( im age_Re, image_Re, 2.0);cvPow( im age_Im, image_Im, 2.0);cvAdd( im age_Re, image_Im, image_Re);cvPow( im age_Re, image_Re, 0.5 );//对数变换以增强灰度级细节(这种变换使以窄带低灰度输入图像值映射//一宽带输出值，具体可见冈萨雷斯数字图像处理p62)// Com pute log(1 + Mag);cvAddS( im age_Re, cvScalar(1.0), image_Re ); // 1 + MagcvLog( im age_Re, image_Re ); // log(1 + Mag)//Rearrange the quadrants of Fourier image so that the origin is at the im age centernRow = src->height;nCol = src->width;cy = nRow/2; // image centercx = nCol/2;//CV_IMAGE_ELEM为OpenCV定义的宏，用来读取图像的像素值，这一部分就是进行中心变换for( j = 0; j < cy; j++ ){for( i = 0; i < cx; i++ ){//中心化，将整体份成四块进行对角交换tmp13 = CV_IMAGE_ELEM( image_Re, double, j, i);CV_IMAGE_ELEM( image_Re, double, j, i) = CV_IMAGE_ELEM(im age_Re, double, j+cy, i+c x);CV_IMAGE_ELEM( image_Re, double, j+cy, i+cx) = t m p13;tmp24 = CV_IMAGE_ELEM( image_Re, double, j, i+cx);CV_IMAGE_ELEM( image_Re, double, j, i+c x) =CV_IMAGE_ELEM( image_Re, double, j+cy, i);CV_IMAGE_ELEM( image_Re, double, j+cy, i) = t m p24;}}//归一化处理将矩阵的元素值归一为[0,255]//[(f(x,y)-minVal)/(maxVal-minVal)]*255double minVal = 0, m axVal = 0;// Localize minimum and m aximum valuescvMinMaxLoc( im age_Re, &minVal, &maxVal );// Normalize image (0 - 255) to be observed as an u8 imagescale = 255/(maxVal - minVal);shift = -minVal * scale;cvConvertScale(image_Re, dst, scale, shift);cvReleaseImage(&image_Re);cvReleaseImage(&image_Im);}/**************************************************************** *******/int m ain(){IplImage *src; //源图像IplImage *Fourier; //傅里叶系数IplImage *dst ;IplImage *ImageRe;IplImage *ImageIm;IplImage *Image;IplImage *ImageDst;double m,M;double scale;double shift;src = cvLoadImage("D:\\main.jpg",0); //加载源图像，第二个参数表示将输入的图片转为单信道Fourier = cvCreateImage(cvGetSize(src),IPL_DEPT H_64F,2);dst = cvCreateIm age(cvGetSize(src),IPL_DEPTH_64F,2);ImageRe = cvCreateImage(cvGetSize(src),IPL_DEPT H_64F,1);ImageIm = cvCreateIm age(cvGetSize(src),IPL_DEPTH_64F,1);Image = cvCreateImage(cvGetSize(src),src->depth,src->nChannels);ImageDst = cvCreateIm age(cvGetSize(src),src->depth,src->nChannels);fft2(src,Fourier); //傅里叶变换fft2shift(Fourier, Image); //中心化cvDFT(Fourier,dst,CV_DXT_INV_SCALE);//实现傅里叶逆变换，并对结果进行缩放cvSplit(dst,Im ageRe,ImageIm,0,0);cvNam edWindow("源图像",0);cvShowIm age("源图像",src);//对数组每个元素平方并存储在第二个参数中cvPow(Im ageRe,Im ageRe,2);cvPow(Im ageIm,ImageIm,2);cvAdd(Im ageRe,ImageIm,ImageRe,NULL);cvPow(Im ageRe,Im ageRe,0.5);cvMinMaxLoc(Im ageRe,&m,&M,NULL,NULL);scale = 255/(M - m);shift = -m * scale;//将shift加在ImageRe各元素按比例缩放的结果上，存储为ImageDst cvConvertScale(Im ageRe,ImageDst,scale,shift);cvNam edWindow("傅里叶谱",0);cvShowIm age("傅里叶谱",Image);cvNam edWindow("傅里叶逆变换",0);cvShowIm age("傅里叶逆变换",ImageDst);//释放图像cvWaitKey(10000); cvReleaseImage(&src); cvReleaseImage(&Image); cvReleaseImage(&ImageIm); cvReleaseImage(&ImageRe); cvReleaseImage(&Fourier); cvReleaseImage(&dst); cvReleaseImage(&ImageDst); return 0;}。