圖像處理與傅里葉變換原理與運(yùn)用

1、圖像處理與傅里葉變換1背景傅里葉變換是一個(gè)非常復(fù)雜的理論，我們?cè)趫D像處理中集中關(guān)注于其傅里葉離散變換離散傅立葉變換 (Discrete Fourier Transform)。1.1離散傅立葉變換圖象是由灰度(RGB組成的二維離散數(shù)據(jù)矩陣，則對(duì)它進(jìn)行傅 立葉變換是離散的傅立葉變換。對(duì)圖像數(shù)據(jù)f(x,y) (x=0,1,M -1; y=0,1,N-1)。則其離散傅立葉變換定義可表示為：1M JN J'uxvy-zzf (x, y)expi 2兀十 “MNXy =0<MN丿一F(u,v)二(1)式中，u=0,1,M-1;v= 0,1,N-1其逆變換為f (x, y)1M -4N -4-

2、''uxvyzIF (u, v) exp |i2兀+MNu mv=0<MN丿-式中，x=0,1，,M-1;y= 0,1,N-1在圖象處理中，一般總是選擇方形數(shù)據(jù)，即M = N影像f(x,y)的振幅譜或傅立葉頻譜：F (u,v) = Jr2 (u, v) + I 2(u, v)相位譜：(u, v)二 arctg 11 (u, v) / R(u, v)丨12快速傅里葉變化可分離性的優(yōu)點(diǎn)是二維的傅立葉變換或逆變換由兩個(gè)連續(xù)的一維傅立葉變換變換來(lái)實(shí)現(xiàn)，對(duì)于一個(gè)影像f(x,y)，可以先沿著其每一列求一維傅立葉變換，再對(duì)其每一行再求一維變換正變化F (u,v)1NNN J N JE

3、Ex -0 y -0_ux + vyf(x,y)exPi2 二需)二'"exp "(巴)丄二(x,y)exp "(空)N x出-N N y乂-N逆變換/ NN1f (x, y)：NN u =o v =oux vy N -11ZN u出-ux 1F (u, v) exp |i2兀()1N 一- I .F (u, v) exp i2：()ILN N 二1zN v=0- vy 1 F(u,v)exp2q由于二維的傅立葉變換具有可分離性，故只討論一維快速傅立葉變換。正變換F(u)二1 N 一1 zN xn一 -i 2 兀 ux f(x)exp-逆變換f (x)二/

4、 N -11zN u =0_ i ux 1 FT 丁由于計(jì)算機(jī)進(jìn)行運(yùn)算的時(shí)間主要取決于所用的乘法的次數(shù)。按照上式進(jìn)行一維離散由空間域向頻率域傅立葉變換時(shí)，對(duì)于N個(gè)F（u）值，中的每一個(gè)都要進(jìn)行 N次運(yùn)算，運(yùn)算時(shí)間與N2成正比1965年庫(kù)里-圖基（Cooly-Tudey提出將運(yùn)算操作降到 NIogzN數(shù)量級(jí)的算法，即N可以分解為一些較小整數(shù)的乘積，當(dāng)N為2的幕（即N=2P，其中P是整數(shù)時(shí)），效率最高，實(shí)現(xiàn)起來(lái)也最簡(jiǎn)單。這就 是快速傅立葉變換。13關(guān)于基圖像（頻率矩形）由二維離散傅里葉反變換式?？芍?，由于u和v均有0,1,N-1的N個(gè)可能的取值，所以f（x,y）由N2個(gè)頻率分量組成，所以每個(gè)頻率分

5、量 都與一個(gè)特定的（u,v）值相對(duì)應(yīng)；且對(duì)于某個(gè)特定的（u,v值來(lái)說(shuō)，當(dāng)（x,y） 取遍所有可能的值（x=0, 1,，N-1; y=0, 1,，N-1 ）時(shí)，就可得 到對(duì)應(yīng)于該特定的（u,v值的一幅基圖像?；鶊D像可表示為。u,vexpexpOu +0v j2:（）N1u +0v j2：（）N0u +1 vexp j2：（）N1u +1vexp j2:（）N0u +（ N -1）v exp j 2 二（）N1u +（N 1）v exp j2:（）N（N 1）u +0vexp j2Jl（）N（N 1）u +1v exp j2：（）N（N 1）u + （N 1）vexp j 2 二（）N所以，一幅圖

6、像的灰度平均值可由 DFT在原點(diǎn)處的值求得十 3.2.4圖像的傅里葉頻譜特性分析 F圖像傅里葉頻譜關(guān)于(M/2, N/2)的對(duì)稱性設(shè)是一幅大小為MXN的圖像，根據(jù)離散傅立 葉變換的周期性公式(3. 40): F(",v)=F(u+mMmN) 有=|F(3)|=|F(+M* + E(3.45)再根據(jù)離散傅立葉變換的共輒對(duì)稱性式(3. 42):(3. 46)|F(",u) |=| F(-w-v) | 就可得：| F(s)冃F(M仏Nt)|3.2.4圖像傅里葉頻譜關(guān)于(M/2, N/2)的對(duì)稱性根據(jù)(3. 46),對(duì)于u二0：當(dāng) v二0 時(shí):円Q0)冃|當(dāng)心 時(shí)：卩(0,1)冃F

7、(M,N-T) | 當(dāng)v二2時(shí)：|尸(0,2» 斗 F(M,N-2)IIIIII當(dāng) V二N/2 時(shí):|F(0,N二 |F(M,N/2)|0M/2MN?2 NA(M/2. N/2ic(M/2. N)(亂 N/2)(亂 N)4324圖像的傅里葉頻譜特性分析 尹石亦贏頻譜關(guān)于(M/2, N/2)的對(duì)稱性同理，對(duì)于廠0：當(dāng) g)時(shí)：F(OF(M,N)0N72NA I BM/2(M/2. N/2(M/C I Dk»*1(MP N/2)<(札N)u當(dāng) uJ 時(shí)：|aO)|=|F(A/-lA9|當(dāng) u 二2 時(shí)：|f(2?0)| =| F(M-2.N) 當(dāng)二/2 時(shí):喊頻譜圖A區(qū)與

8、D區(qū)和B區(qū)與C區(qū)由此可得:關(guān)于坐標(biāo)(M/ZN/2)對(duì)稱。證明 周期性與共軛對(duì)稱性F(u,v)=F(u+mM,v+Nn) F(u,v)=F(-u,-v)1 JV-IjV-1F(it + mN、mN)-NN j=o j=oCq迦沁)旳®(哼)exp- i27T(mx+nyj因?yàn)?exp-z2(/?7A'+n>')為單位值&3 2.4圖像的傅里葉頻譜特性分析像傅里葉頻譜特性及其頻譜圖3.4(b/圖3-5b)的原頻譜圖一原點(diǎn)在(60)時(shí)的頻譜圖ffl3 7»頻譜團(tuán)何和原點(diǎn)平移到(M/2t N/勿后的頻譜圖f (x, y) (-1)以 y) = F (u

9、 - M /2, v - N / 2)也就是說(shuō)，圖3.7的頻譜圖(a)和(b)實(shí)質(zhì)上是函數(shù)行| y的傅里葉頻譜圖。1.3運(yùn)用1.3.1頻率濾波在頻域中，頻率越大說(shuō)明原始信號(hào)變化速度越快；頻率越小說(shuō)明原始信號(hào)越平緩。當(dāng)頻率為0時(shí)，表示直流信號(hào)，沒(méi)有變化。因此， 頻率的大小反應(yīng)了信號(hào)的變化快慢。高頻分量解釋信號(hào)的突變部分， 而低頻分量決定信號(hào)的整體形象。在圖像處理中，頻域反應(yīng)了圖像在空域灰度變化劇烈程度，也就是圖像灰度的變化速度，也就是圖像的梯度大小。對(duì)圖像而言，圖像的邊緣部分是突變部分，變化較快，因此反應(yīng)在頻域上是高頻分量；圖像的噪聲大部分情況下是高頻部分；圖像平緩變化部分則為低頻分 量。也就是

10、說(shuō)，傅立葉變換提供另外一個(gè)角度來(lái)觀察圖像，可以將圖 像從灰度分布轉(zhuǎn)化到頻率分布上來(lái)觀察圖像的特征。書(shū)面一點(diǎn)說(shuō)就 是，傅里葉變換提供了一條從空域到頻率自由轉(zhuǎn)換的途徑。對(duì)圖像處理而言，以下概念非常的重要：圖像高頻分量：圖像突變部分；在某些情況下指圖像邊緣信息，某些情況下指噪聲，更多是兩者的混合；低頻分量：圖像變化平緩的部分，也就是圖像輪廓信息高通濾波器：讓圖像使低頻分量抑制，高頻分量通過(guò)低通濾波器：與高通相反，讓圖像使高頻分量抑制，低頻分量 通過(guò)帶通濾波器：使圖像在某一部分的頻率信息通過(guò)，其他過(guò)低或過(guò)高都抑制模板運(yùn)算與卷積定理在時(shí)域內(nèi)做模板運(yùn)算，實(shí)際上就是對(duì)圖像進(jìn)行卷積。模板運(yùn)算是 圖像處理一個(gè)很

11、重要的處理過(guò)程，很多圖像處理過(guò)程，比如增強(qiáng)/去噪（這兩個(gè)分不清楚），邊緣檢測(cè)中普遍用到。根據(jù)卷積定理，時(shí)域 卷積等價(jià)與頻域乘積。因此，在時(shí)域內(nèi)對(duì)圖像做模板運(yùn)算就等效于在 頻域內(nèi)對(duì)圖像做濾波處理。比如說(shuō)一個(gè)均值模板，其頻域響應(yīng)為一個(gè)低通濾波器； 在時(shí)域內(nèi) 對(duì)圖像作均值濾波就等效于在頻域內(nèi)對(duì)圖像用均值模板的頻域響應(yīng) 對(duì)圖像的頻域響應(yīng)作一個(gè)低通濾波1.3.2 圖像去噪圖像去噪就是壓制圖像的噪音部分。因此，如果噪音是高頻額， 從頻域的角度來(lái)看， 就是需要用一個(gè)低通濾波器對(duì)圖像進(jìn)行處理。 通 過(guò)低通濾波器可以抑制圖像的高頻分量。 但是這種情況下常常會(huì)造成 邊緣信息的抑制。常見(jiàn)的去噪模板有均值模板，高斯模

12、板等。這兩種 濾波器都是在局部區(qū)域抑制圖像的高頻分量， 模糊圖像邊緣的同時(shí)也 抑制了噪聲。還有一種非線性濾波 - 中值濾波器。中值濾波器對(duì)脈沖 型噪聲有很好的去掉。 因?yàn)槊}沖點(diǎn)都是突變的點(diǎn)， 排序以后輸出中值， 那么那些最大點(diǎn)和最小點(diǎn)就可以去掉了。 中值濾波對(duì)高斯噪音效果較 差。椒鹽噪聲： 對(duì)于椒鹽采用中值濾波可以很好的去除。 用均值也可以取 得一定的效果，但是會(huì)引起邊緣的模糊。高斯白噪聲：白噪音在整個(gè)頻域的都有分布，好像比較困難。岡薩雷斯版圖像處理 P185:算術(shù)均值濾波器和幾何均值濾波器（尤 其是后者） 更適合于處理高斯或者均勻的隨機(jī)噪聲。 諧波均值濾波器 更適合于處理脈沖噪聲。1.3.3

13、 圖像增強(qiáng)有時(shí)候感覺(jué)圖像增強(qiáng)與圖像去噪是一對(duì)矛盾的過(guò)程， 圖像增強(qiáng)經(jīng) 常是需要增強(qiáng)圖像的邊緣， 以獲得更好的顯示效果， 這就需要增加圖 像的高頻分量。 而圖像去噪是為了消除圖像的噪音， 也就是需要抑制 高頻分量。有時(shí)候這兩個(gè)又是指類似的事情。比如說(shuō)，消除噪音的同 時(shí)圖像的顯示效果顯著的提升了，那么，這時(shí)候就是同樣的意思了。常見(jiàn)的圖像增強(qiáng)方法有對(duì)比度拉伸， 直方圖均衡化，圖像銳化等。前 面兩個(gè)是在空域進(jìn)行基于像素點(diǎn)的變換， 后面一個(gè)是在頻域處理。我 理解的銳化就是直接在圖像上加上圖像高通濾波后的分量，也就是圖像的邊緣效果。對(duì)比度拉伸和直方圖均衡化都是為了提高圖像的對(duì)比 度，也就是使圖像看起來(lái)差異

14、更明顯一些， 我想，經(jīng)過(guò)這樣的處理以 后，圖像也應(yīng)該增強(qiáng)了圖像的高頻分量，使得圖像的細(xì)節(jié)上差異更大。 同時(shí)也引入了一些噪音。1.4實(shí)現(xiàn)在MATLAB中F=imread(file name);F二fft2(f,P,Q);儘成 FFT變換FC二fftshift(F):%實(shí)現(xiàn)居中S=abs(F或 Fc);%取得傅里葉頻譜 f=real(ifft2(F);%實(shí)現(xiàn)傅里葉逆變換基于OPENCV庫(kù)#i nclude <stdio.h>#in clude < cv.h>#in clude < cxcore.h>#in clude <highgui.h>/*src

15、 IPL_DEPTH_8U /dst IPL DEPTH 64F*/ 傅里葉正變換void fft2(IplImage *src, IplImage *dst) / 實(shí)部、虛部IplImage *image_Re = 0, *im age_Im = 0, *Fourier = 0;/ int i, j;im age_Re = cvCreateIm age(cvGetSize(src), IPL_DEPTH_64F, 1);/Im aginary partim age_Im = cvCreateImage(cvGetSize(src), IPL_DEPTH_64F, 1);/2 channels

16、 (image_Re, image_Im)Fourier = cvCreateImage(cvGetSize(src), IPL_DEPTH_64F, 2);/ Real part conversion from u8 to 64f (double) cvConvertScale(src, im age_Re, 1, 0);/ Im aginary part (zeros) cvZero(im age_Im);/ Join real and imaginary parts and stock them in Fourier image cvMerge(im age_Re, image_Im,

17、0, 0, Fourier);/ Application of the forward Fourier transform cvDFT(Fourier, dst, CV_DXT_FORWARD); cvReleaseImage(&image_Re);cvReleaseImage(&image_Im); cvReleaseImage(&Fourier);/*/src IPL_DEPTH_64F/dst IPL_DEPTH_8U*/void fft2shift(IplImage *src, IplImage *dst)IplImage *image_Re = 0, *ima

18、ge_Im = 0;int nRow, nCol, i, j, cy, cx;double scale, shift, tmp13, tmp24;image_Re = cvCreateImage(cvGetSize(src), IPL_DEPTH_64F, 1); /Imaginary partimage_Im = cvCreateImage(cvGetSize(src), IPL_DEPTH_64F, 1); cvSplit( src, im age_Re, image_Im, 0, 0 );/ 具體原理見(jiàn)岡薩雷斯數(shù)字圖像處理 p123/ Com pute the magn itude of

19、 the spectrumMag = sqrt(ReA2 +/ 實(shí)部/ 虛部Im A2)/ 計(jì)算傅里葉譜cvPow( im age_Re, image_Re, 2.0);cvPow( im age_Im, image_Im, 2.0); cvAdd( im age_Re, image_Im, image_Re);cvPow( im age_Re, image_Re, 0.5 );/ 對(duì)數(shù)變換以增強(qiáng)灰度級(jí)細(xì)節(jié) ( 這種變換使以窄帶低灰度輸入圖像值映射/ 一寬帶輸出值，具體可見(jiàn)岡薩雷斯數(shù)字圖像處理 p62)/ Com pute log(1 + Mag);cvAddS( im age_Re, cvS

20、calar(1.0), image_Re ); / 1 +MagcvLog( im age_Re, image_Re ); / log(1 +Mag)/Rearrange the quadrants of Fourier image so that the origin is atthe im agecenternRow = src->height;nCol = src-> width;cy = nRow/2; / image centercx = nCol/2;/CV_IMAGE_ELEM 為 OpenCV 定義的宏，用來(lái)讀取圖像的像素值，這一部分就是進(jìn)行中心 變換for( j

21、= 0; j < cy; j+ )for( i = 0; i < cx; i+ )/ 中心化，將整體份成四塊進(jìn)行對(duì)角交換tmp13 =CV_IMAGE_ELEM( image_Re, double, j, i);CV_IMAGE_ELEM( image_Re, double, j, i) = CV_IMAGE_ELEM( im age_Re, double, j+ cy, i+ c x);CV_IMAGE_ELEM( image_Re, double, j+cy, i+ cx) = t mp13;tmp24 = CV_IMAGE_ELEM( image_Re, double, j,

22、 i+ cx);cx) =cy, i);cy, i) = t mp24;CV_IMAGE_ELEM( image_Re, double, j, i+CV_IMAGE_ELEM( image_Re, double, j+CV_IMAGE_ELEM( image_Re, double, j+ / 歸一化處理將矩陣的元素值歸一為 0,255 /(f(x,y)-minVal)/(maxVal- minVal)*255 double minVal = 0, m axVal = 0;/ Localize minimum and m aximum valuescvMinMaxLoc( image_Re, &

23、amp;minVal, &maxVal );/ Normalize image (0 - 255) to be observed as an u8 image scale = 255/(maxVal - minVal);shift = - minVal * scale; cvConvertScale(image_Re, dst, scale, shift); cvReleaseImage(&image_Re); cvReleaseImage(&image_Im);*/int m ain()IplImage *src; / 源圖像IplImage *Fourier; /

24、傅里葉系數(shù)IplImage *dst ;IplImage *ImageRe;IplImage *ImageIm;IplImage *Image;IplImage *ImageDst; double m ,M;double scale;double shift;src = cvLoadImage("D:main.jpg",0); / 加載源圖像，第二個(gè)參數(shù)表示將輸入的圖片轉(zhuǎn) 為單信道Fourier = cvCreateImage(cvGetSize(src),IPL_DEPT H_64F,2); dst = cvCreateIm age(cvGetSize(src),IPL_

25、DEPTH_64F,2);ImageRe = cvCreateImage(cvGetSize(src),IPL_DEPT H_64F,1);ImageIm = cvCreateIm age(cvGetSize(src),IPL_DEPTH_64F,1);Image = cvCreateImage(cvGetSize(src),src->depth,src-> nChannels); ImageDst = cvCreateIm age(cvGetSize(src),src -> depth,src-> nChannels); fft2(src,Fourier); / 傅里葉變換 fft2shift(Fourier, Image); / 中心化 cvDFT(Fourier,dst,CV_DXT_INV_SCALE);/ 實(shí)現(xiàn)傅里葉