麻辣GIS微信平台

更多 GIS 干货

微信关注不错过

IDL实现傅里叶(傅立叶)变换

关于傅里叶变换

这个网上的描述实在太多,原理上就是把空域的图片转换成频域。可以看这张维基百科的示意图

实现公式

IDL实现

;Author:尹全超
;2014-01-17

PRO DFTSELF

  ;读取图像
  file =  DIALOG_PICKFILE(/READ, FILTER = '*.bmp')
  binary_img = READ_BMP(file)
  img01 = image(binary_img,LAYOUT=[2,2,1],title='Original')

  ;傅里叶变换
  N_col=N_elements(binary_img[*,0])
  M_row=N_elements(binary_img[0,*])
  M=M_row
  N=N_col
  n_dft=complexarr(N,M)
  n_dft_tp=complexarr(N,M)
  dft_img=complexarr(N,M)
  ;列(一维DFT)
  FOR v=0,N-1,1 DO BEGIN 
    FOR x=0,M-1,1 DO BEGIN
      fxv=complex(0.0d,0.0d)
      FOR y=0,N-1,1 DO BEGIN
        fxv=fxv+binary_img[y,x]*complex(cos(2*!pi*v*y/N),-1*sin(2*!pi*v*y/N))
      ENDFOR
      n_dft_tp[v,x]=fxv
    ENDFOR  
  ENDFOR
  ;行(一维DFT)
  FOR u=0,M-1,1 DO BEGIN 
    FOR y=0,N-1,1 DO BEGIN
      fxv=complex(0.0d,0.0d)
      FOR x=0,M-1,1 DO BEGIN
        fxv=fxv+n_dft_tp[y,x]*complex(cos(2*!pi*u*x/M),-1*sin(2*!pi*u*x/M))
      ENDFOR
      n_dft[y,u]=fxv
    ENDFOR  
  ENDFOR
  n_dft=n_dft/(M*N)

  powerSpectrum = ABS(n_dft)^2
  scaledPowerSpect = ALOG10(powerSpectrum)
  img02 = IMAGE(scaledPowerSpect,/CURRENT,LAYOUT=[2,2,2],title='MyFFT')

  ;自带FFT
  ffTransform = FFT(binary_img)
  powerSpectrum = ABS(ffTransform)^2
  scaledPowerSpect = ALOG10(powerSpectrum)
  img03 = IMAGE(scaledPowerSpect,/CURRENT,LAYOUT=[2,2,3],title='FFT')

  ;反傅立叶变换
  ;列(一维DFT)
  FOR v=0,N-1,1 DO BEGIN 
    FOR x=0,M-1,1 DO BEGIN
      fxv=complex(0.0d,0.0d)
      FOR y=0,N-1,1 DO BEGIN
        fxv=fxv+n_dft[y,x]*complex(cos(2*!pi*v*y/N),sin(2*!pi*v*y/N))
      ENDFOR
      n_dft_tp[v,x]=fxv
    ENDFOR  
  ENDFOR
  ;行(一维DFT)
  FOR u=0,M-1,1 DO BEGIN 
    FOR y=0,N-1,1 DO BEGIN
      fxv=complex(0.0d,0.0d)
      FOR x=0,M-1,1 DO BEGIN
        fxv=fxv+n_dft_tp[y,x]*complex(cos(2*!pi*u*x/M),sin(2*!pi*u*x/M))
      ENDFOR
      dft_img[y,u]=fxv
    ENDFOR  
  ENDFOR
  ;dft_img=dft_img/(M*N)

  img04 = IMAGE(dft_img,/CURRENT,LAYOUT=[2,2,4],title='DFT to IMAGE')

END

这个代码实现了DFT的变换方法,将图像做傅立叶变换,并和IDL自带FFT函数做了一个对比,然后再做了一个反变换。

结果如图

缺陷

目前反傅立叶变换的时候,会产生色差,具体的bug没有找到。如果你找到了相关的方法,欢迎留言,或者给本站投稿。

相关阅读

麻辣GIS-Sailor

作者:

GIS爱好者,学GIS,更爱玩GIS。

声明

1.本文所分享的所有需要用户下载使用的内容(包括但不限于软件、数据、图片)来自于网络或者麻辣GIS粉丝自行分享,版权归该下载资源的合法拥有者所有,如有侵权请第一时间联系本站删除。

2.下载内容仅限个人学习使用,请切勿用作商用等其他用途,否则后果自负。

手机阅读
公众号关注
知识星球
手机阅读
麻辣GIS微信公众号关注
最新GIS干货
关注麻辣GIS知识星球
私享圈子

已有 9 条评论

  1. 麻辣GIS-Cook
    1#
    Cook  · 2015-03-24 15:20

  2. 麻辣GIS-wwww
    2#
    wwww  · 2018-11-12 14:29

    有一句语句是错误的,但我不知道怎么改,想问一下

      1. 麻辣GIS-Sailor
        Sailor  · 2018-11-12 16:25

        哪句?

          1. 麻辣GIS-wwww
            wwww  · 2018-11-12 16:44

            fxv=fxv+binary_img[y,x]complex(cos(2!pivy/N),-1sin(2!pivy/N))

              1. 麻辣GIS-Sailor
                Sailor  · 2018-11-13 09:35

                我更新了一下,你再试试

                  1. 麻辣GIS-yusw
                    yusw  · 2019-03-20 12:33

                    你好,请问更新的版本在哪

                  2. 麻辣GIS-Sailor
                    Sailor  · 2019-03-21 13:10

                    就是文章里,目前是改过之后的。

  3. 麻辣GIS-wwww
    3#
    wwww  · 2018-11-13 10:28

    可以了,感谢!但是为什么会出现浮点型的错误。

留言板(小编看到第一时间回复)