Python下opencv使用筆記（十）（圖像頻域濾波與傅里葉變換）

前面以前介紹過空間域濾波，空間域濾波就是用各種模板直接與圖像進行卷積運算，實現對圖像的處理，這個方案直接對圖像空間操作，操作簡單。所以也是空間域濾波。

頻域濾波說究竟終於可能是和空間域濾波實現相同的功能，比方實現圖像的輪廓提取，在空間域濾波中我們使用一個拉普拉斯模板就能夠提取，而在頻域內，我們使用一個高通濾波模板（由於輪廓在頻域內屬於高頻信號）。能夠實現輪廓的提取，后面也會把拉普拉斯模板頻域化。會發現拉普拉斯事實上在頻域來講就是一個高通濾波器。

既然是頻域濾波就涉及到把圖像首先變到頻域內。那么把圖像變到頻域內的方法就是傅里葉變換。關於傅里葉變換，感覺真是個偉大的發明。尤其是其在信號領域的應用，對於傅里葉變換的理解。要是剛接觸這個東西，想要理解還真是非常的困難。除非你的數學功底特別好。

這里推薦一個非常好的通俗易懂的博客(待會再去o-_-)：

傅里葉分析之掐死教程（完整版）

傅里葉的原理表明，不論什么連續測量的時序或信號，都能夠表示為不同頻率的正弦波信號的無限疊加。利用傅立葉變換算法直接測量原始信號。以累加方式來計算該信號中不同正弦波信號的頻率、振幅和相位就能夠表示原始信號。這里借用上述博客的一個圖：
這個圖就是把時域圖像（大概是方波）變成了一系列的正弦波的線性疊加。其等價關系能夠表示為：

f(原始信號)=A1sin(w1x+ϕ1)+A2sin(w2x+ϕ2)+...

那么 w1,w2,... 能夠看成是頻率的變化（一般覺得就是從1,2,…n定死了）。全部的A就是相應頻率下的振幅，全部的 ϕ 就是相應頻率下的相位。那么對於任一個信號，假設都覺得頻率w是從1,2,3…一直添加的話。那么每一個信號就僅僅由一組振幅與一組 ϕ 來決定，他們的不同決定了終於信號的不同。

再來理解下什么是振幅，振幅就是各個頻率下的信號的決定程度有多大。假設某個頻率的振幅越大，那么它對原始信號的的重要性越大。像上圖，當然是w=1的時候振幅最大，說明它對總的信號影響最多（去掉w=1的信號，原始信號講嚴重變形）。越往后面，也就是越高頻。振幅逐漸減小，那么他們的作用就越小。而他們對於總體信號又有什么影響呢？既然越小，那就是影響小，所以事實上去掉。原始信號也基本上不變，他們影響就在於對原始信號的細節上的表現。比方原始信號上的邊邊角角，偶爾有個小凸起凹槽什么的，這些小細節部分都是靠這些個影響不大的高頻信號來表現出來的。

深入推廣一下。這就非常好理解為什么圖像的高頻信號事實上表現出來的就是圖像的邊緣輪廓、噪聲等等這些細節的東西了，而低頻信號，表現的卻是圖像上整塊整塊灰度大概一樣的區域了（這些個區域又稱為直流分量區域）。

再來理解下什么是相位。相位表示事實上表面相應頻率下的正弦分量偏離原點的程度，再借用下上述博客中的一個圖。把分量示意圖放大了：

上圖看到，假設各個頻率的分量相位都是0的話。那么每一個正弦分量的最大值（在頻率軸附近的那個最大值）都會落在頻率軸為0上，然而上述圖並非這樣。

在說簡單一點。比方原始信號上有個凹槽，正好是由某一頻率的分量疊加出來的，那么假設這個頻率的相位變大一點或者變小一點的話，帶來的影響就會使得這個凹槽向左或者向右移動一下。也就是說，相位的作用就是精確定位到信號上一點的位置的。

好了。有了上述的概念，再來看看圖像的傅里葉變換。上述舉得樣例是一維信號的傅里葉變換。而且信號是連續的。我們知道圖像是二維離散的。連續與離散都能夠用傅里葉進行變換，那么二維信號無非就是在x方向與y方向都進行一次一維的傅里葉變換得到，這么看來，能夠想象，它的頻率構成就是一個網格矩陣了，橫軸從w=1到n，縱軸也是這樣。

全部圖像的頻率構成都覺得是這種，那么不同的就是一幅圖的振幅與相位了（振幅與相位此時相同是一個網格矩陣）。也就是說你在opencv或者matlab下對圖像進行傅里葉變換后事實上是能夠得到圖像的振幅圖與相位圖的，而想把圖像從頻域空間恢復到時域空間，必需要同一時候有圖像的振幅圖與相位圖才干夠。缺少一個就恢復的不完整（后面會實驗看看）。
以下看看二維傅里葉變換，一個圖像為M*N的圖像f(x,y)進過離散傅里葉變換得到F(u,v)，那么一般的公式為：

F(u,v)=∑x=0M−1∑y=0N−1f(x,y)e−j2π(ux/M+vy/N)

它的反變換就是

f(x,y)=1MN∑x=0M−1∑y=0N−1F(u,v)ej2π(ux/M+vy/N)

反變換就能夠實現將頻域圖像恢復到時域圖像。

對正變換分析，當u=v=0時，那么：

F(0,0)=∑x=0M−1∑y=0N−1f(x,y)=f¯(x,y)

看看這個式子發現它是什么？f(x,y)可是時域里面圖像的灰度。非常明顯，這個東西事實上是整個圖像的灰度求平均了（這里說的都是灰度圖像），當圖像進行傅里葉變換以后，你去把F(0,0)與原始圖像平均灰度去比較看看是不是一樣的。那么F(0,0)在頻域內稱為直流分量，其它的全部F稱為交流分量，直流分量能夠看到是在0處獲得的。所以非常明顯是存在於低頻分量下的。

說了這么多。來上個實驗看看究竟什么是傅里葉變換吧。

在python+opencv下想實現圖像傅里葉變換有兩種途徑，一種採用numpy包能夠實現，還有opencv自帶的能夠實現。當中numpy帶的使用方便。直觀易懂。

import cv2
import numpy as np
import matplotlib.pyplot as plt

img = cv2.imread('flower.jpg',0) #直接讀為灰度圖像
f = np.fft.fft2(img)
fshift = np.fft.fftshift(f)
#取絕對值：將復數變化成實數
#取對數的目的為了將數據變化到較小的范圍（比方0-255）
s1 = np.log(np.abs(f))
s2 = np.log(np.abs(fshift))
plt.subplot(121),plt.imshow(s1,'gray'),plt.title('original')
plt.subplot(122),plt.imshow(s2,'gray'),plt.title('center')

注意的是，上圖事實上並沒有什么含義，顯示出來的能夠看成是頻域后圖像的振幅信息，並沒有相位信息，圖像的相位 ϕ=atan(實部虛部) ，numpy包中自帶一個angle函數能夠直接依據復數的實部與虛部求出角度（默認出來的角度是弧度）。像上述程序出來的f與fshift都是復數，就能夠直接angle函數一下，比方試試並把相應的相位圖像顯示出來：

import cv2
import numpy as np
import matplotlib.pyplot as plt

img = cv2.imread('flower.jpg',0) #直接讀為灰度圖像
f = np.fft.fft2(img)
fshift = np.fft.fftshift(f)
#取絕對值：將復數變化成實數
#取對數的目的為了將數據變化到較小的范圍（比方0-255）
ph_f = np.angle(f)
ph_fshift = np.angle(fshift)

plt.subplot(121),plt.imshow(ph_f,'gray'),plt.title('original')
plt.subplot(122),plt.imshow(ph_fshift,'gray'),plt.title('center')

這就是圖像上每一個像素點相應的相位圖，事實上是毫無規律的，理解就是偏移的角度。

Ok再來說說程序中為什么要有一個np.fft.fftshift(f)中心化操作，整個圖像是在傅里葉變換的一個周期內完畢的，將其看成橫縱兩個方向的一維傅里葉變換，在每一個方向上都會有高頻信號和低頻信號，那么傅里葉變換將低頻信號放在了邊緣。高頻信號放在了中間。然而一副圖像，非常明顯的低頻信號多而明顯，所以將低頻信號採用一種方法移到中間，在時域上就是對f乘以（-1）^(M+N)，換到頻域里面就是位置的移到了。

圖像變換到頻域后就能夠進行操作了，眼下接觸到的頻域操作似乎也就是一些濾波操作，如同空域里面的濾波操作一樣，只是原理不一樣了，后面再說說一些頻域濾波方法。好了一旦操作完，得到的數據還是頻域數據。那么怎樣將其變換到時域呢？這里就是傅里葉反變換了。公式表示就如同前面那樣。這個頻域變換到時域的操作就是逆向傅里葉變換再走一遍（比方先反中心化。在逆變換）。

一個實比例如以下：

import cv2
import numpy as np
import matplotlib.pyplot as plt

img = cv2.imread('flower.jpg',0) #直接讀為灰度圖像
f = np.fft.fft2(img)
fshift = np.fft.fftshift(f)
#取絕對值：將復數變化成實數
#取對數的目的為了將數據變化到0-255
s1 = np.log(np.abs(fshift))
plt.subplot(131),plt.imshow(img,'gray'),plt.title('original')
plt.subplot(132),plt.imshow(s1,'gray'),plt.title('center')
# 逆變換
f1shift = np.fft.ifftshift(fshift)
img_back = np.fft.ifft2(f1shift)
#出來的是復數，無法顯示
img_back = np.abs(img_back)
plt.subplot(133),plt.imshow(img_back,'gray'),plt.title('img back')

能夠看到恢復的一模一樣。

我們說，恢復一個頻域圖像需要圖像的振幅以及相位。而一個復數也正好包括這些，振幅就是實部虛部的平方和開方，相位就是 atan(實部虛部) 。前面說過。那么如今假設我們僅僅用一副圖像的振幅或者相位來將頻域內圖像恢復到時域會怎么樣呢？以下給出僅僅用振幅、僅僅用相位以及兩者在聯合起來來恢復的程序：

import cv2
import numpy as np
import matplotlib.pyplot as plt

img = cv2.imread('flower.jpg',0) #直接讀為灰度圖像
f = np.fft.fft2(img)
fshift = np.fft.fftshift(f)
#取絕對值：將復數變化成實數
#取對數的目的為了將數據變化到0-255
s1 = np.log(np.abs(fshift))
plt.subplot(221),plt.imshow(img,'gray'),plt.title('original')
plt.xticks([]),plt.yticks([])
#---------------------------------------------
# 逆變換--取絕對值就是振幅
f1shift = np.fft.ifftshift(np.abs(fshift))
img_back = np.fft.ifft2(f1shift)
#出來的是復數，無法顯示
img_back = np.abs(img_back)
#調整大小范圍便於顯示
img_back = (img_back-np.amin(img_back))/(np.amax(img_back)-np.amin(img_back))
plt.subplot(222),plt.imshow(img_back,'gray'),plt.title('only Amplitude')
plt.xticks([]),plt.yticks([])
#---------------------------------------------
# 逆變換--取相位
f2shift = np.fft.ifftshift(np.angle(fshift))
img_back = np.fft.ifft2(f2shift)
#出來的是復數，無法顯示
img_back = np.abs(img_back)
#調整大小范圍便於顯示
img_back = (img_back-np.amin(img_back))/(np.amax(img_back)-np.amin(img_back))
plt.subplot(223),plt.imshow(img_back,'gray'),plt.title('only phase')
plt.xticks([]),plt.yticks([])
#---------------------------------------------
# 逆變換--將兩者合成看看
s1 = np.abs(fshift) #取振幅
s1_angle = np.angle(fshift) #取相位
s1_real = s1*np.cos(s1_angle) #取實部
s1_imag = s1*np.sin(s1_angle) #取虛部
s2 = np.zeros(img.shape,dtype=complex) 
s2.real = np.array(s1_real) #又一次賦值給s2
s2.imag = np.array(s1_imag)

f2shift = np.fft.ifftshift(s2) #對新的進行逆變換
img_back = np.fft.ifft2(f2shift)
#出來的是復數，無法顯示
img_back = np.abs(img_back)
#調整大小范圍便於顯示
img_back = (img_back-np.amin(img_back))/(np.amax(img_back)-np.amin(img_back))
plt.subplot(224),plt.imshow(img_back,'gray'),plt.title('another way')
plt.xticks([]),plt.yticks([])

能夠看到，僅僅振幅的恢復圖啥也不是。僅僅的相位圖還有那么點意思，當然也是啥也不是。最后是把振幅與相位分別作為頻域內復數的實部和虛部。得到的恢復圖才與原來的一樣。

基於此，我們來做一個有趣的實驗，假設有兩幅圖像，將這兩幅圖像進行傅里葉變換到頻域，然后把用一個圖像的振幅做為振幅，用另一幅圖像的相位作為相位生成一副新的圖像，那么，這個圖像會怎么樣呢？你覺得生成的圖像會更像取振幅的那副還是取相位的那副呢？來看看吧：

import cv2
import numpy as np
import matplotlib.pyplot as plt

img_flower = cv2.imread('flower.jpg',0) #直接讀為灰度圖像
img_man = cv2.imread('woman.jpg',0) #直接讀為灰度圖像
plt.subplot(221),plt.imshow(img_flower,'gray'),plt.title('origial1')
plt.xticks([]),plt.yticks([])
plt.subplot(222),plt.imshow(img_man,'gray'),plt.title('origial_2')
plt.xticks([]),plt.yticks([])
#--------------------------------
f1 = np.fft.fft2(img_flower)
f1shift = np.fft.fftshift(f1)
f1_A = np.abs(f1shift) #取振幅
f1_P = np.angle(f1shift) #取相位
#--------------------------------
f2 = np.fft.fft2(img_man)
f2shift = np.fft.fftshift(f2)
f2_A = np.abs(f2shift) #取振幅
f2_P = np.angle(f2shift) #取相位
#---圖1的振幅--圖2的相位--------------------
img_new1_f = np.zeros(img_flower.shape,dtype=complex) 
img1_real = f1_A*np.cos(f2_P) #取實部
img1_imag = f1_A*np.sin(f2_P) #取虛部
img_new1_f.real = np.array(img1_real) 
img_new1_f.imag = np.array(img1_imag) 
f3shift = np.fft.ifftshift(img_new1_f) #對新的進行逆變換
img_new1 = np.fft.ifft2(f3shift)
#出來的是復數，無法顯示
img_new1 = np.abs(img_new1)
#調整大小范圍便於顯示
img_new1 = (img_new1-np.amin(img_new1))/(np.amax(img_new1)-np.amin(img_new1))
plt.subplot(223),plt.imshow(img_new1,'gray'),plt.title('another way')
plt.xticks([]),plt.yticks([])
#---圖2的振幅--圖1的相位--------------------
img_new2_f = np.zeros(img_flower.shape,dtype=complex) 
img2_real = f2_A*np.cos(f1_P) #取實部
img2_imag = f2_A*np.sin(f1_P) #取虛部
img_new2_f.real = np.array(img2_real) 
img_new2_f.imag = np.array(img2_imag) 
f4shift = np.fft.ifftshift(img_new2_f) #對新的進行逆變換
img_new2 = np.fft.ifft2(f4shift)
#出來的是復數。無法顯示
img_new2 = np.abs(img_new2)
#調整大小范圍便於顯示
img_new2 = (img_new2-np.amin(img_new2))/(np.amax(img_new2)-np.amin(img_new2))
plt.subplot(224),plt.imshow(img_new2,'gray'),plt.title('another way')
plt.xticks([]),plt.yticks([])

這就是合成的圖像，是不是非常有趣。圖像3是圖1的振幅加圖2的相位，圖像4是圖1的相位加上圖2的振幅。非常明顯的能夠看到，新圖像占用誰的相位就越像誰，為什么會這樣？非常easy，能夠理解振幅只是描寫敘述圖像灰度的亮度。占用誰的振幅只是使得結果哪些部分偏亮或者暗而已。而圖像是個什么樣子是由它的相位決定的。相位描寫敘述的是一個方向，方向正確了。那么終於的結果離你的目的就不遠了。

可想而知，方向對於一件事物是多么的重要。大自然的規律尚且如此。更別說做人做事了，找准並相信一個方向。慢慢的走下去吧，總有一天會看到成果的，哪怕你的振幅不正確。走的慢一點，可是終於也能走的像模像樣的，不會說到了人生的最后。回頭一看，再來感嘆哎呀這是什么玩意，那就非常悲哀了。

好了。扯回來我們的問題，關於傅里葉變換下的圖像恢復基本上就是這了。可是我們發現。似乎我們僅僅說了開頭（怎么變換）和結尾（怎么變回去）。那么中間我們要做的東西才是傅里葉變換的目的—頻域下的各種變化操作。

這里主要介紹頻域下的濾波–低通濾波器。高通濾波器，帶通帶阻濾波器。

我們知道，圖像在變換加移動中心后。從中間到外面，頻率上依次是從低頻到高頻的。那么我們假設把中間規定一小部分去掉，是不是相對於把低頻信號去掉了呢？這也就是相當於進行了高通濾波。

這個濾波模板畫出來可能就是這種：

黑色為0，白色為1，把這個模板去和圖像進過傅里葉變換的頻域矩陣去與（相乘）一下就實現了高通濾波。比方以下：

import cv2
import numpy as np
import matplotlib.pyplot as plt

img_man = cv2.imread('woman.jpg',0) #直接讀為灰度圖像
plt.subplot(121),plt.imshow(img_man,'gray'),plt.title('origial')
plt.xticks([]),plt.yticks([])
#--------------------------------
rows,cols = img_man.shape
mask = np.ones(img_man.shape,np.uint8)
mask[rows/2-30:rows/2+30,cols/2-30:cols/2+30] = 0
#--------------------------------
f1 = np.fft.fft2(img_man)
f1shift = np.fft.fftshift(f1)
f1shift = f1shift*mask
f2shift = np.fft.ifftshift(f1shift) #對新的進行逆變換
img_new = np.fft.ifft2(f2shift)
#出來的是復數。無法顯示
img_new = np.abs(img_new)
#調整大小范圍便於顯示
img_new = (img_new-np.amin(img_new))/(np.amax(img_new)-np.amin(img_new))
plt.subplot(122),plt.imshow(img_new,'gray'),plt.title('Highpass')
plt.xticks([]),plt.yticks([])

能夠看到。高通濾波器有利於提取圖像的輪廓，這里我們從原理上分析一下為什么，圖像的輪廓或者邊緣或者一些噪聲處，灰度變化劇烈，那么在把它們經過傅里葉變換后。就會變成高頻信號（我們知道高頻時捕捉細節的），所以在把圖像低頻信號濾掉以后剩下的自然就是輪廓了。

反過來我們來看看空間域濾波中的拉普拉斯模板。我們知道這個模板是這種

M=⎡⎣⎢0−10−14−10−10⎤⎦⎥

如今我們把這個模板進行傅里葉變換到頻域看看：

import numpy as np
import matplotlib.pyplot as plt

laplacian = np.array([[0,1,0],[1,-4,1],[0,1,0]])
f = np.fft.fft2(laplacian)
f1shift = np.fft.fftshift(f)
img = np.log(np.abs(f1shift))
plt.imshow(img,'gray')

能夠看到，這個模板的頻域變換下四周特別亮。也就是是個高通濾波器。

其它空間域下的模板都能夠轉到頻域下來看看。

以下再來試試低通濾波器什么效果，構造個低通濾波器也非常easy，把上述模板中的1變成0,0變成1。事實上就是低通了：

import cv2
import numpy as np
import matplotlib.pyplot as plt

img_man = cv2.imread('woman.jpg',0) #直接讀為灰度圖像
plt.subplot(121),plt.imshow(img_man,'gray'),plt.title('origial')
plt.xticks([]),plt.yticks([])
#--------------------------------
rows,cols = img_man.shape
mask = np.zeros(img_man.shape,np.uint8)
mask[rows/2-20:rows/2+20,cols/2-20:cols/2+20] = 1
#--------------------------------
f1 = np.fft.fft2(img_man)
f1shift = np.fft.fftshift(f1)
f1shift = f1shift*mask
f2shift = np.fft.ifftshift(f1shift) #對新的進行逆變換
img_new = np.fft.ifft2(f2shift)
#出來的是復數。無法顯示
img_new = np.abs(img_new)
#調整大小范圍便於顯示
img_new = (img_new-np.amin(img_new))/(np.amax(img_new)-np.amin(img_new))
plt.subplot(122),plt.imshow(img_new,'gray'),plt.title('Highpass')
plt.xticks([]),plt.yticks([])

能夠看到低通濾波后圖像除了輪廓模糊了外，基本上沒什么變化，圖像的大部分信息基本上都保持了。從原理上來看，圖像的主要信息都集中在低頻上，所以低通濾波器的效果是這樣也是能夠理解的。上述的高通、低通濾波器的構造有0,1構成的理想濾波器，也是最簡單的濾波器，另一些其它的濾波器。比方說高斯濾波器，butterworth濾波器等等。以下是參考岡薩雷斯書上的圖：

還有就是把高通低通都結合一部分到一個模板上形成的帶通濾波器，這在一些場合會非常實用，還是以理想的帶通濾波器實驗下：

import cv2
import numpy as np
import matplotlib.pyplot as plt

img_man = cv2.imread('woman.jpg',0) #直接讀為灰度圖像
plt.subplot(121),plt.imshow(img_man,'gray'),plt.title('origial')
plt.xticks([]),plt.yticks([])
#--------------------------------
rows,cols = img_man.shape
mask1 = np.ones(img_man.shape,np.uint8)
mask1[rows/2-8:rows/2+8,cols/2-8:cols/2+8] = 0
mask2 = np.zeros(img_man.shape,np.uint8)
mask2[rows/2-80:rows/2+80,cols/2-80:cols/2+80] = 1
mask = mask1*mask2
#--------------------------------
f1 = np.fft.fft2(img_man)
f1shift = np.fft.fftshift(f1)
f1shift = f1shift*mask
f2shift = np.fft.ifftshift(f1shift) #對新的進行逆變換
img_new = np.fft.ifft2(f2shift)
#出來的是復數。無法顯示
img_new = np.abs(img_new)
#調整大小范圍便於顯示
img_new = (img_new-np.amin(img_new))/(np.amax(img_new)-np.amin(img_new))
plt.subplot(122),plt.imshow(img_new,'gray')
plt.xticks([]),plt.yticks([])

這就是帶通的效果。既能夠保留一部分低頻。也能夠保留一部分高頻。至於保留多少。怎么保留就視問題的不同而不同了。

當然頻域濾波的應用絕不僅限於此，很多其它的應用還有待探索學習。