« 形态学运算小结标准MFC在OpenGL »

快速傅立叶变换

傅立叶变换,在图形学上貌似应用不怎么多,海水波浪的运动算法方面需要吧,然后不怎么见了。而在图像处理上,频域处理这个大领域,傅立叶变换可是算法基础口牙~——ZwqXin.com

其实啊,这个傅立叶的实现同是为了作业,但就在形态学处理之前。之所以迟迟不记录在博,主要还是因为自己对其原理还是云里雾里的。嘛~有时候写多了日志就有种“什么东西要说就得说明白啊”的错觉。其实嘛,本博客也不是什么传道授业解惑的地方,只是写着写着就有点忘记初衷了。

废话太多了。其实我只想说,找傅立叶原理算法的同学,接下来这篇文章没啥太大价值的。姑且当作一个被傅立叶实现搞晕了一两天的傻子的自言自语,恩,这里是我的记录簿。

本文来源于 ZwqXin (http://www.zwqxin.com/), 转载请注明
      原文地址:http://www.zwqxin.com/archives/image-processing/dif-fast-fourier.html

最初被傅立叶缠上应该是上自动控制原理的时候吧,印象大概就是时域下的东西转去频率域吧,然后扔出一个复杂的公式。都知道,连续波都可以由正弦搅成。那么,代表任意连续波的函数由正弦函数搅成如何呢?正弦波状态特征是频率,振幅,相位吧,那么,一段时间内的波的函数形式分成几种不同周期(不同频率)的正弦表示的“有机组合”,就是所谓的时域 - 频域转换了。这很有用吧,在很多领域,包括自动控制.....不过其实我不太清楚具体起来是什么用,包括自动控制领域。

而在图像上的应用,其实也就是空间域 - 频率域。一张500*400的灰度图像,不妨假想为一张有500*400个方格的地形图(啥叫3D中毒...),某坐标处的灰度自然就是地形在该点的高度了(当然,没可能就像一张地形灰度图那样能形成比较连续的地形)。图像是按像素记数的,每个相邻像素的灰度或多或少都有个差值,形成一个梯度——关键就在于这个“或多或少”上。现在,忽略相邻像素灰度间的插值等问题,就有一个地形出来了,地形当然有陡坡有缓坡(正是由梯度造成)——这就是我们要处理的空间波了。傅立叶变换功效在于,变换结果是频率(梯度)——经过变换所得的是一堆按灰度大小聚集的新像素(总的像素数量不变,但意义改变了,反映的是梯度而不是原来的“灰度”,而且与原来的空间位置无关)。

该看什么呢?就看“聚集的情况”(如,图中的白色区域和黑色区域,注意,它们没有明显的分界,因此只能直观去感受~)对于中部高亮区域,反映了与相邻像素梯度比较大的像素量;周围的黑色区域,反映了与相邻像素梯度比较小的像素量。这么说,就是图像中的高频部分和低频部分了——图像中的边缘处就是最直观的“高频部分”了。至于要分离出高频和低频有什么用,那就是频率域处理的事情了。傅立叶变换作用就只是生成一张这样的图以把高低频部分分开。

图像中离散的像素可以认为是二维波的采样点。因此对每行进行一维离散傅立叶变换,再在此基础上对每列进行一维离散傅立叶变换,则能得出两个纬度结合的傅立叶频率图(如上)。离散傅立叶变换的计算形式对计算机而言太残酷了(前面说了,式子复杂得很),因此就有人发明出快速傅立叶变换这种东西。

最容易理解的是“基-2”的快速傅立叶,因为它的采样点尺寸是2的N次方(N> 0),能直接执行蝶形算法。一般图像长宽都不一定满足——因此比较受限。补或截都是迫不得已的方法。N为合数的 FFT 算法 (混合基)可以不受限,不过复杂死了,我不想去深究。对“基-2”的快速傅立叶,又分为“时域抽取基-2”和“频域抽取基-2”,其实差不多,不过中间的步骤顺序有点差异而已(前者还要倒码之类的)。

我用的是“频域抽取基-2快速傅立叶”(DIF-FFT)。算法推导,诸君在网络上找找吧,反正结论就是一个符合蝶形算法的形式。

  1.  
  2. //  DIT – FFT [频域抽取基-2FFT] 蝶形算法:
  3. double angle;
  4.     for(i  =  0;  i  <  Num / 2;  i++)    
  5.     {    
  6.       angle  =  -i  *  PI  *  2  /  Num;    
  7.       if(!InverseTranform)
  8.           W[i]  =  complex <double>  (cos(angle),  sin(angle)); 
  9.       else
  10.           W[i]  =  complex <double>  (cos(angle),  -sin(angle)); 
  11.      } 
  12.  
  13.     int Apartsize = 0, halfPartsize = 0, Partoffset = 0;
  14.     for(k = 0; k < r; k++)
  15.     {
  16.         for(j = 0; j < (1<<k); j++)
  17.         {
  18.             Apartsize = 1 << (r - k);
  19.             halfPartsize = Apartsize / 2;
  20.  
  21.             for(i = 0; i < halfPartsize; i++)
  22.             {
  23.                 Partoffset = j * Apartsize;
  24.                 X2[i + Partoffset] 
  25.                     = X1[i + Partoffset] + X1[i + Partoffset + halfPartsize];
  26.  
  27.                 X2[i + Partoffset + halfPartsize] 
  28.                     = (X1[i + Partoffset] - X1[i + Partoffset + halfPartsize]) * W[i * (1<<k)];
  29.             }
  30.         }
  31.         X = X1;//
  32.         X1 =X2;
  33.         X2 = X;//
  34.     }
  35.  
  36. //到二维,横竖方向各转一次:
  37.     for(j = 0; j <Exheight; ++j)
  38.         DifFastFourierTransform1D(&input[j * Exwidth], &output[j * Exwidth], Exwtimes);
  39.  
  40.     for(j = 0; j <Exheight; ++j)
  41.       for(i = 0; i <Exwidth; ++i)
  42.           input[i + j * Exwidth] = output[j + i * Exheight];
  43.  
  44.      for(i = 0; i <Exwidth; ++i)
  45.           DifFastFourierTransform1D(&input[i * Exheight], &output[i * Exheight], Exhtimes);

另外,对彩色图片,一般是针对每个颜色通道都做一个快速傅立叶变换哦。

再另外,其实傅立叶变换除了可以生成频率域图,还能生成相位图哦。只是输出的时候输出结果的振幅就可了,相位图能反映出梯度点在整张图片上的方位呢~

最后一句话:反变换跟正变换差不多的。不多得结合频率域图和相位图的信息哦,才能还原。

程序截图:
[若看不到本博客图片,可参照此帖,简单的复原法:[显示本站所有图片] ]

 快速傅立叶变换[DifFFTProcess]
 快速傅立叶 www.zwqxin.com
(原始图)
快速傅立叶 www.zwqxin.com
(Fourier)

快速傅立叶 www.zwqxin.com
(彩色图,原图)
快速傅立叶 www.zwqxin.com
      快速傅立叶 www.zwqxin.com
[对非2幂大小的某图像,按需(见对话框)取其各个角域。图上为变换(上)和反变换]
快速傅立叶 www.zwqxin.com
[按需察看彩色图的不同通道下的频域图,图左为上图的蓝色通道;图右:上图的相位图]

本文来源于 ZwqXin (http://www.zwqxin.com/), 转载请注明
      原文地址:http://www.zwqxin.com/archives/image-processing/dif-fast-fourier.html

  • quote 1.fan
  • 文章不是一般的精彩! 看你两段写的内容 比学校半学期讲的 都受用
  • 2009-11-8 23:13:33 回复该留言

发表评论:

◎欢迎参与讨论,请在这里发表您的看法、交流您的观点。

IE下本页面显示有问题?

→点击地址栏右侧【兼容视图】←

日历

Search

网站分类

最新评论及回复

最近发表

Powered By Z-Blog 1.8 Walle Build 100427

Copyright 2008-2013 ZwqXin. All Rights Reserved. Theme edited from ipati.