图像的傅里叶变换
from: http://www.techug.com/image-vary?ref=myread
傅里葉變換(Fourier Transform)是非常重要的數(shù)學(xué)分析工具,同時也是一種非常重要的信號處理方法。我記得本科課程電路原理中有提到過,但由于計(jì)算過于復(fù)雜,好像是超出考試范圍了,所以并沒有深入學(xué)習(xí)。最近實(shí)驗(yàn)中需要對圖像進(jìn)行濾波處理,文獻(xiàn)中提到的方法通常是經(jīng)過傅里葉變換之后對頻域進(jìn)行過濾,將圖像中的低頻信息與高頻信息區(qū)分開來。
理解傅里葉變換對非數(shù)學(xué)專業(yè)的人來說比較難以理解的原因主要有兩方面。首先,由于涉及到比較復(fù)雜繁瑣的數(shù)學(xué)操作,看到下面的這兩個公式,一般人可能當(dāng)時就蒙了:
f?(t)=∫∞?∞f(x)e?2πixtdx;
f(x)=∫∞?∞f?(t)e2πitxdt;
另一方面的原因則是由于變換過程比較抽象,很難從直覺上去把握在傅里葉變換過程中到底發(fā)生了什么。關(guān)于傅里葉變換的科普文章,有一篇是知乎上傳閱較廣的《傅里葉分析之掐死教程》,作者用了盡可能少的數(shù)學(xué)公式和圖形分解來解決這兩方面的問題。在此之前,國外有個專門科普數(shù)學(xué)概念的網(wǎng)站(Better Explained)也寫了一篇類似的科普文章,但更徹底的是,全文都不涉及到任何數(shù)學(xué)公式和推斷,完全用英語來向讀者解釋傅里葉變換過程。第一篇文章以音樂和樂譜為例,第二篇作者用“奶昔”作為例子,我想了半天終于找到一個更通俗的例子:
我們來想象一下,假設(shè)這塊面餅的厚度是N層面條,每一層都是由一根彎曲成正弦曲線形狀的面條排列而成,有些面條波浪較大,也就是排列較為稀疏,而密有些排列較密集。我們把這N層面條擠壓到一起,就得到上圖這一塊雜亂無序、世間獨(dú)一無二的面餅。傅里葉變換所做的事就是把上面的過程反過來,我們可以從一塊完整的面餅得到最初的N層面條。如果有些人的口味比較特殊,喜歡波浪大又稀疏的面,于是我們就將排列太緊湊(高頻面)剔除之后再重新壓制成一塊新的面餅,這就是我們最終想要的濾波(Filter)的過程。
帶著對面餅的想象,我們來看一種更為抽象、優(yōu)雅的描述(from Wikipedia):
一般的波形或者說信號(Signal)都是基于時間尺度上的采樣結(jié)果,因此也稱為時域(Time Domain),而上面泡面的例子和我們將要處理的圖像信號則是基于空間尺度上的采樣,但好像并沒有“空域(Space Domain)”這一說,畢竟我們對空間的感知仍然依賴于時間。不過在空間尺度上我們可以更直觀地認(rèn)為信號是靜止,例如下面這張圖像(灰度圖),其實(shí)是由250×250個像素點(diǎn)組成,每個像素點(diǎn)的灰度值([0,255])就是基于像素坐標(biāo)的空間采樣的結(jié)果:
右邊的3D Fourier就是一塊長相奇怪的面餅。
實(shí)踐
對傅里葉變換有了大概的了解之后可以先動手嘗試一下,來更加直觀地感受一下(實(shí)際上完全可以在不理解的情況下,直接上手)。這里用到的是OpenCV + Numpy,實(shí)際上OpenCV和Numpy都提供了快速傅里葉變換(FFT)算法:
import cv2 as cv import numpy as np from matplotlib import pyplot as pltimg = cv.imread('Joseph_Fourier_250.jpg', 0) f = np.fft.fft2(img) # 快速傅里葉變換算法得到頻率分布 fshift = np.fft.fftshift(f) # 默認(rèn)結(jié)果中心點(diǎn)位置是在左上角,轉(zhuǎn)移到中間位置fimg = np.log(np.abs(fshift)) # fft 結(jié)果是復(fù)數(shù),求絕對值結(jié)果才是振幅# 展示結(jié)果 plt.subplot(121), plt.imshow(img, 'gray'), plt.title('Original Fourier') plt.subplot(122), plt.imshow(fimg, 'gray'), plt.title('Fourier Fourier') plt.show()右邊圖就是頻率分布圖譜,其中越靠近中心的位置頻率越低,越亮(灰度值越高)的位置代表該頻率的信號振幅越大。fft的結(jié)果是復(fù)數(shù)形式,保留了圖像的全部信息,但去絕對值得到的頻譜圖只表現(xiàn)了振幅而沒有體現(xiàn)相位。
回想一下高中時候?qū)W過的三角函數(shù):
f(x)=Asin(ωx+φ)=Asin(2πfx+φ)
一個正弦波是由下面三個參數(shù)決定的:
- 角速度(頻率)ω=2πf ;
- 振幅A;
- 相位φ。
除了上面這個公式之外,還可以用另外一種形式來(唯一地)表示一個正弦波(from BetterExplained):
即:
cos(x)+isin(x)?a+ib
所以說,fft的復(fù)數(shù)結(jié)果保留了正弦波成分的所有信息,但頻譜圖只展現(xiàn)了頻率和振幅的分布。因此可以根據(jù)fft的結(jié)果還原原始圖像,但是我們做傅里葉變換的目的并不是為了觀察圖像的頻率分布(至少不是最終目的),更多情況下是為了對頻率進(jìn)行過濾。過濾的方法一般有三種:低通(Low-pass)、高通(High-pass)、帶通(Band-pass)。所謂低通就是保留圖像中的低頻成分,過濾高頻成分,可以把過濾器想象成一張漁網(wǎng),根據(jù)上文對頻譜圖的解讀,想要低通過濾器,就是將高頻區(qū)域的信號全部拉黑,而低頻區(qū)域全部保留:
img = cv.imread('Joseph_Fourier_250.jpg', 0) f = np.fft.fft2(img) # 快速傅里葉變換算法得到頻率分布 fshift = np.fft.fftshift(f) # 默認(rèn)結(jié)果中心點(diǎn)位置是在左上角,轉(zhuǎn)移到中間位置lpButterMask = butterWorthLPF(img.shape[:2], 24, 2) hpButterMask = butterWorthHPF(img.shape[:2], 36, 2)lpFshift = fshift * lpButterMask maskedInvf = np.fft.ifft2(np.fft.ifftshift(lpFshift)) lpfImg = np.abs(maskedInvf)hpFshift = fshift * hpButterMask maskedInvf = np.fft.ifft2(np.fft.ifftshift(hpFshift)) hpfImg = np.abs(maskedInvf)plt.subplot(221), plt.imshow(lpButterMask, 'gray'), plt.title('Butterworth LPF') plt.subplot(222), plt.imshow(lpfImg, 'gray'), plt.title('LPF Image') plt.subplot(223), plt.imshow(hpButterMask, 'gray'), plt.title('HPF Spectrum') plt.subplot(224), plt.imshow(hpfImg, 'gray'), plt.title('HPF Image')plt.show()很顯然,濾波器的選擇也是很重要的,這里用到的是 Butterworth 濾波器,有興趣的可以自己實(shí)現(xiàn)一下:P
總結(jié)
- 上一篇: 关于x86、i386、i486、i586
- 下一篇: Ubuntu 16.04更新软件提示需要