Python下opencv使用筆記(圖像頻域濾波與傅里葉變換)
本文轉載自 https://blog.csdn.net/on2way/article/details/46981825
首先謝謝原創博主了,這篇文章對我幫助很大,記錄下方便再次閱讀。
Python下opencv使用筆記(圖像頻域濾波與傅里葉變換)
前面曾經介紹過空間域濾波,空間域濾波就是用各種模板直接與圖像進行卷積運算,實現對圖像的處理,這種方法直接對圖像空間操作,操作簡單,所以也是空間域濾波。
頻域濾波說到底最終可能是和空間域濾波實現相同的功能,比如實現圖像的輪廓提取,在空間域濾波中我們使用一個拉普拉斯模板就可以提取,而在頻域內,我們使用一個高通濾波模板(因為輪廓在頻域內屬於高頻信號),可以實現輪廓的提取,后面也會把拉普拉斯模板頻域化,會發現拉普拉斯其實在頻域來講就是一個高通濾波器。
既然是頻域濾波就涉及到把圖像首先變到頻域內,那么把圖像變到頻域內的方法就是傅里葉變換。關於傅里葉變換,感覺真是個偉大的發明,尤其是其在信號領域的應用,對於傅里葉變換的理解,要是剛接觸這個東西,想要理解還真是非常的困難,除非你的數學功底特別好。這里推薦一個非常好的通俗易懂的博客(待會再去o-_-):
傅里葉分析之掐死教程(完整版)(o(╯□╰)o,鏈接文章已經找不到了)
傅里葉的原理表明,任何連續測量的時序或信號,都可以表示為不同頻率的正弦波信號的無限疊加。利用傅立葉變換算法直接測量原始信號,以累加方式來計算該信號中不同正弦波信號的頻率、振幅和相位就可以表示原始信號。這里借用上述博客的一個圖:
那么可以看成是頻率的變化(一般認為就是從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),那么一般的公式為:

反變換就可以實現將頻域圖像恢復到時域圖像。對正變換分析,當u=v=0時,那么:
說了這么多,來上個實驗看看到底什么是傅里葉變換吧。在python+opencv下想實現圖像傅里葉變換有兩種途徑,一種采用numpy包可以實現,還有opencv自帶的可以實現,其中numpy帶的使用方便,直觀易懂。
-
1 import cv2
-
2 import numpy as np
-
3 import matplotlib.pyplot as plt
-
4
-
5 img = cv2.imread('flower.jpg',0) #直接讀為灰度圖像
-
6 f = np.fft.fft2(img)
-
7 fshift = np.fft.fftshift(f)
-
8 #取絕對值:將復數變化成實數
-
9 #取對數的目的為了將數據變化到較小的范圍(比如0-255)
-
10 s1 = np.log(np.abs(f))
-
11 s2 = np.log(np.abs(fshift))
-
12 plt.subplot(121),plt.imshow(s1,'gray'),plt.title('original')
-
13 plt.subplot(122),plt.imshow(s2,'gray'),plt.title('center')
注意的是,上圖其實並沒有什么含義,顯示出來的可以看成是頻域后圖像的振幅信息,並沒有相位信息,圖像的相位
ϕ=atan(實部虛部),numpy包中自帶一個angle函數可以直接根據復數的實部與虛部求出角度(默認出來的角度是弧度)。像上述程序出來的f與fshift都是復數,就可以直接angle函數一下,比如試試並把對應的相位圖像顯示出來:
-
1 import cv2
-
2 import numpy as np
-
3 import matplotlib.pyplot as plt
-
4
-
5 img = cv2.imread('flower.jpg',0) #直接讀為灰度圖像
-
6 f = np.fft.fft2(img)
-
7 fshift = np.fft.fftshift(f)
-
8 #取絕對值:將復數變化成實數
-
9 #取對數的目的為了將數據變化到較小的范圍(比如0-255)
-
10 ph_f = np.angle(f)
-
11 ph_fshift = np.angle(fshift)
-
12
-
13 plt.subplot(121),plt.imshow(ph_f,'gray'),plt.title('original')
-
14 plt.subplot(122),plt.imshow(ph_fshift,'gray'),plt.title('center')
這就是圖像上每個像素點對應的相位圖,其實是毫無規律的,理解就是偏移的角度。
Ok再來說說程序中為什么要有一個np.fft.fftshift(f)中心化操作,整個圖像是在傅里葉變換的一個周期內完成的,將其看成橫縱兩個方向的一維傅里葉變換,在每個方向上都會有高頻信號和低頻信號,那么傅里葉變換將低頻信號放在了邊緣,高頻信號放在了中間,然而一副圖像,很明顯的低頻信號多而明顯,所以將低頻信號采用一種方法移到中間,在時域上就是對f乘以(-1)^(M+N),換到頻域里面就是位置的移到了。
圖像變換到頻域后就可以進行操作了,目前接觸到的頻域操作似乎也就是一些濾波操作,如同空域里面的濾波操作一樣,不過原理不一樣了,后面再說說一些頻域濾波方法。好了一旦操作完,得到的數據還是頻域數據,那么如何將其變換到時域呢?這里就是傅里葉反變換了,公式表示就如同前面那樣。這個頻域變換到時域的操作就是逆向傅里葉變換再走一遍(比如先反中心化,在逆變換)。一個實例如下:
-
1 import cv2
-
2 import numpy as np
-
3 import matplotlib.pyplot as plt
-
4
-
5 img = cv2.imread('flower.jpg',0) #直接讀為灰度圖像
-
6 f = np.fft.fft2(img)
-
7 fshift = np.fft.fftshift(f)
-
8 #取絕對值:將復數變化成實數
-
9 #取對數的目的為了將數據變化到0-255
-
10 s1 = np.log(np.abs(fshift))
-
11 plt.subplot(131),plt.imshow(img,'gray'),plt.title('original')
-
12 plt.subplot(132),plt.imshow(s1,'gray'),plt.title('center')
-
13 # 逆變換
-
14 f1shift = np.fft.ifftshift(fshift)
-
15 img_back = np.fft.ifft2(f1shift)
-
16 #出來的是復數,無法顯示
-
17 img_back = np.abs(img_back)
-
18 plt.subplot(133),plt.imshow(img_back,'gray'),plt.title('img back')
可以看到恢復的一模一樣。
我們說,恢復一個頻域圖像需要圖像的振幅以及相位,而一個復數也正好包含這些,振幅就是實部虛部的平方和開方,相位就是
atan(實部虛部),前面說過。那么現在假設我們只用一副圖像的振幅或者相位來將頻域內圖像恢復到時域會怎么樣呢?下面給出只用振幅、只用相位以及兩者在聯合起來來恢復的程序:
-
1 import cv2
-
2 import numpy as np
-
3 import matplotlib.pyplot as plt
-
4
-
5 img = cv2.imread('flower.jpg',0) #直接讀為灰度圖像
-
6 f = np.fft.fft2(img)
-
7 fshift = np.fft.fftshift(f)
-
8 #取絕對值:將復數變化成實數
-
9 #取對數的目的為了將數據變化到0-255
-
10 s1 = np.log(np.abs(fshift))
-
11 plt.subplot(221),plt.imshow(img,'gray'),plt.title('original')
-
12 plt.xticks([]),plt.yticks([])
-
13 #---------------------------------------------
-
14 # 逆變換--取絕對值就是振幅
-
15 f1shift = np.fft.ifftshift(np.abs(fshift))
-
16 img_back = np.fft.ifft2(f1shift)
-
17 #出來的是復數,無法顯示
-
18 img_back = np.abs(img_back)
-
19#調整大小范圍便於顯示
-
20 img_back = (img_back-np.amin(img_back))/(np.amax(img_back)-np.amin(img_back))
-
21 plt.subplot(222),plt.imshow(img_back,'gray'),plt.title('only Amplitude')
-
22 plt.xticks([]),plt.yticks([])
-
23 #---------------------------------------------
-
24 # 逆變換--取相位
-
25 f2shift = np.fft.ifftshift(np.angle(fshift))
-
26 img_back = np.fft.ifft2(f2shift)
-
27 #出來的是復數,無法顯示
-
28 img_back = np.abs(img_back)
-
29 #調整大小范圍便於顯示
-
30 img_back = (img_back-np.amin(img_back))/(np.amax(img_back)-np.amin(img_back))
-
31 plt.subplot(223),plt.imshow(img_back,'gray'),plt.title('only phase')
-
32 plt.xticks([]),plt.yticks([])
-
33 #---------------------------------------------
-
34 # 逆變換--將兩者合成看看
-
35 s1 = np.abs(fshift) #取振幅
-
36 s1_angle = np.angle(fshift) #取相位
-
37 s1_real = s1*np.cos(s1_angle) #取實部
-
38 s1_imag = s1*np.sin(s1_angle) #取虛部
-
39 s2 = np.zeros(img.shape,dtype=complex)
-
40 s2.real = np.array(s1_real) #重新賦值給s2
-
41 s2.imag = np.array(s1_imag)
-
42
-
43 f2shift = np.fft.ifftshift(s2) #對新的進行逆變換
-
44 img_back = np.fft.ifft2(f2shift)
-
45 #出來的是復數,無法顯示
-
46 img_back = np.abs(img_back)
-
47 #調整大小范圍便於顯示
-
48 img_back = (img_back-np.amin(img_back))/(np.amax(img_back)-np.amin(img_back))
-
49 plt.subplot(224),plt.imshow(img_back,'gray'),plt.title('another way')
-
50 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的振幅。很明顯的可以看到,新圖像占用誰的相位就越像誰,為什么會這樣?很簡單,可以理解振幅不過描述圖像灰度的亮度,占用誰的振幅不過使得結果哪些部分偏亮或者暗而已,而圖像是個什么樣子是由它的相位決定的。相位描述的是一個方向,方向正確了,那么最終的結果離你的目的就不遠了。可想而知,方向對於一件事物是多么的重要,大自然的規律尚且如此,更別說做人做事了,找准並相信一個方向,慢慢的走下去吧,總有一天會看到成果的,哪怕你的振幅不對,走的慢一點,但是最終也能走的像模像樣的,不會說到了人生的最后,回頭一看,再來感嘆哎呀這是什么玩意,那就很悲哀了。
好了,扯回來我們的問題,關於傅里葉變換下的圖像恢復基本上就是這了。但是我們發現,似乎我們只說了開頭(怎么變換)和結尾(怎么變回去),那么中間我們要做的東西才是傅里葉變換的目的—頻域下的各種變化操作。
這里主要介紹頻域下的濾波–低通濾波器,高通濾波器,帶通帶阻濾波器。
我們知道,圖像在變換加移動中心后,從中間到外面,頻率上依次是從低頻到高頻的,那么我們如果把中間規定一小部分去掉,是不是相對於把低頻信號去掉了呢?這也就是相當於進行了高通濾波。這個濾波模板畫出來可能就是這樣的:
黑色為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⎤⎦⎥
現在我們把這個模板進行傅里葉變換到頻域看看:
-
1 import numpy as np
-
2 import matplotlib.pyplot as plt
-
3
-
4 laplacian = np.array([[0,1,0],[1,-4,1],[0,1,0]])
-
5 f = np.fft.fft2(laplacian)
-
6 f1shift = np.fft.fftshift(f)
-
7 img = np.log(np.abs(f1shift))
-
8 plt.imshow(img,'gray')
可以看到,這個模板的頻域變換下四周特別亮,也就是是個高通濾波器。其它空間域下的模板都可以轉到頻域下來看看。
下面再來試試低通濾波器什么效果,構造個低通濾波器也很簡單,把上述模板中的1變成0,0變成1,其實就是低通了:
-
1 import cv2
-
2 import numpy as np
-
3 import matplotlib.pyplot as plt
-
4
-
5 img_man = cv2.imread('woman.jpg',0) #直接讀為灰度圖像
-
6 plt.subplot(121),plt.imshow(img_man,'gray'),plt.title('origial')
-
7 plt.xticks([]),plt.yticks([])
-
8 #--------------------------------
-
9 rows,cols = img_man.shape
-
10 mask = np.zeros(img_man.shape,np.uint8)
-
11 mask[rows/2-20:rows/2+20,cols/2-20:cols/2+20] = 1
-
12 #--------------------------------
-
13 f1 = np.fft.fft2(img_man)
-
14 f1shift = np.fft.fftshift(f1)
-
15 f1shift = f1shift*mask
-
16 f2shift = np.fft.ifftshift(f1shift) #對新的進行逆變換
-
17 img_new = np.fft.ifft2(f2shift)
-
18 #出來的是復數,無法顯示
-
19 img_new = np.abs(img_new)
-
20 #調整大小范圍便於顯示
-
21 img_new = (img_new-np.amin(img_new))/(np.amax(img_new)-np.amin(img_new))
-
22 plt.subplot(122),plt.imshow(img_new,'gray'),plt.title('Highpass')
-
23 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([])
這就是帶通的效果,既可以保留一部分低頻,也可以保留一部分高頻,至於保留多少,怎么保留就視問題的不同而不同了。
當然頻域濾波的應用絕不僅限於此,更多的應用還有待探索學習。