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

图像处理与傅里叶变换 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 ),( 1 0 1 0 M x N y N vy M ux iyxf MN vuF )2(2exp),( 1 ),( 1 0 1 0 M u N v N vy M ux ivuF MN yxf ),(),(),( 22 vuIvuRvuF ),(/ ),(),(vuRvuIarctgvu ),(),(),(),( 22 2 vuIvuRvuFvuE 1.2 快速傅里叶变化快速傅里叶变化 可分离性的优点是二维的傅立叶变换或逆变换由两个连续的一 维傅立叶变换变换来实现，对于一个影像 f(x,y)，可以先沿着其每一 列求一维傅立叶变换，再对其每一行再求一维变换 正变化正变化 逆变换 由于二维的傅立叶变换具有可分离性，故只讨论一维快速傅立叶由于二维的傅立叶变换具有可分离性，故只讨论一维快速傅立叶 变换。变换。 正变换 逆变换 由于计算机进行运算的时间主要取决于所用的乘法的次数。 按照上式进行一维离散由空间域向频率域傅立叶变换时，对于 N 个 1 0 1 0 1 0 1 0 )(2exp),( 1 )(2exp),( 1 )(2exp),( 1 ),( N v N u N u N v N vy ivuF NN ux ivuF N N vyux ivuF NN yxf 1 0 2 exp)( 1 )( N x N uxi xf N uF 1 0 1 0 1 0 1 0 )(2exp),( 1 )(2exp),( 1 )(2exp),( 1 ),( N y N x N x N y N vy iyxf NN ux iyxf N N vyux iyxf NN vuF 1 0 2 exp)( 1 )( N u N uxi uF N xf F（u)值，中的每一个都要进行 N 次运算，运算时间与运算时间与 N2成正比。成正比。 1965 年库里-图基（ Cooly-Tudey)提出将运算操作降到 Nlog2N 数量级的算法，即 N 可以分解为一些较小整数的乘积，当 N 为 2 的 幂（即 N=2P，其中 P 是整数时） ，效率最高，实现起来也最简单。 这就是快速傅立叶变换。 1.3 关于基图像（频率矩形）关于基图像（频率矩形） 由二维离散傅里叶反变换式。可知，由于由二维离散傅里叶反变换式。可知，由于 u 和和 v 均有均有 0,1,N-1 的的 N 个可能的取值，所以个可能的取值，所以 f(x,y)由由 N2个频率分量组成，所以每个频率分个频率分量组成，所以每个频率分 量都与一个特定的量都与一个特定的(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 1 2 , N vNuN j N vuN j N vuN j N vNu j N vu j N vu j N vNu j N vu j N vu j N f vu 证明证明 周期性与共轭对称性周期性与共轭对称性 F(u,v)=F(u+mM,v+Nn) F(u,v)=F(-u,-v) 也就是说，图也就是说，图 3.7 的频谱图的频谱图(a)和和(b)实质上是函数实质上是函数的的 ( 1)( + )(,) 傅里叶频谱图。傅里叶频谱图。 1.3 运用 1.3.1 频率滤波 在频域中，频率越大说明原始信号变化速度越快；频率越小说明 原始信号越平缓。当频率为 0 时，表示直流信号，没有变化。因此， 频率的大小反应了信号的变化快慢。高频分量解释信号的突变部分， 而低频分量决定信号的整体形象。 在图像处理中，频域反应了图像在空域灰度变化剧烈程度，也就 是图像灰度的变化速度，也就是图像的梯度大小。对图像而言，图 像的边缘部分是突变部分，变化较快，因此反应在频域上是高频分 )2/, 2/() 1(),( )( NvMuFyxf yx 量；图像的噪声大部分情况下是高频部分；图像平缓变化部分则为 低频分量。也就是说，傅立叶变换提供另外一个角度来观察图像， 可以将图像从灰度分布转化到频率分布上来观察图像的特征。书面 一点说就是，傅里叶变换提供了一条从空域到频率自由转换的途径。 对图像处理而言，以下概念非常的重要： 图像高频分量：图像突变部分；在某些情况下指图像边缘信息， 某些情况下指噪声，更多是两者的混合； 低频分量：图像变化平缓的部分，也就是图像轮廓信息 高通滤波器：让图像使低频分量抑制，高频分量通过 低通滤波器：与高通相反，让图像使高频分量抑制，低频分量 通过 带通滤波器：使图像在某一部分的频率信息通过，其他过低或 过高都抑制 模板运算与卷积定理 在时域内做模板运算，实际上就是对图像进行卷积。模板运算 是图像处理一个很重要的处理过程，很多图像处理过程，比如增强/ 去噪（这两个分不清楚） ，边缘检测中普遍用到。根据卷积定理，时 域卷积等价与频域乘积。因此，在时域内对图像做模板运算就等效 于在频域内对图像做滤波处理。 比如说一个均值模板，其频域响应为一个低通滤波器；在时域 内对图像作均值滤波就等效于在频域内对图像用均值模板的频域响 应对图像的频域响应作一个低通滤波。 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 #include #include #include /* * /src IPL_DEPTH_8U /dst IPL_DEPTH_64F * */ /傅里叶正变换 void fft2(IplImage *src, IplImage *dst) /实部、虚部 IplImage *image_Re = 0, *image_Im = 0, *Fourier = 0; / int i, j; image_Re = cvCreateImage(cvGetSize(src), IPL_DEPTH_64F, 1); /实部 /Imaginary part image_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, image_Re, 1, 0); / Imaginary part (zeros) cvZero(image_Im); / Join real and imaginary parts and stock them in Fourier image cvMerge(image_Re, image_Im, 0, 0, Fourier); / Application of the forward Fourier transform cvDFT(Fourier, dst, CV_DXT_FORWARD); cvReleaseImage( cvReleaseImage( cvReleaseImage( /* * /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 part image_Im = cvCreateImage(cvGetSize(src), IPL_DEPTH_64F, 1); cvSplit( src, image_Re, image_Im, 0, 0 ); /具体原理见冈萨雷斯数字图像处理 p123 / Compute the magnitude of the spectrum Mag = sqrt(Re2 + Im2) /计算傅里叶谱 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 ); /对数变换以增强灰度级细节(这种变换使以窄带低灰度输入图像值映射 /一宽带输出值，具体可见冈萨雷斯数字图像处理 p62) / Compute log(1 + Mag); cvAddS( image_Re, cvScalar(1.0), image_Re ); / 1 + Mag cvLog( image_Re, image_Re ); / log(1 + Mag) /Rearrange the quadrants of Fourier image so that the origin is at the image center nRow = src-height; nCol = src-width; cy = nRow/2; / image center cx = nCol/2; /CV_IMAGE_ELEM 为 OpenCV 定义的宏，用来读取图像的像素值，这一部分就是进行中心 变换 for( j = 0; j depth,src-nChannels); ImageDst = cvCreateImage(cvGetSize(src),src-depth,src-nChannels); fft2(src,Fourier); /傅里叶变换 fft2shift(Fourier, Image); /中心化 cvDFT(Fourier,dst,CV_DXT_INV_SCALE);/实现傅里叶逆变换，并