傅里葉變換輪廓術(shù)的Matlab仿真實現(xiàn)

吳應(yīng)山，張啟燦

(四川大學 電子信息學院，四川 成都 610064)

?

傅里葉變換輪廓術(shù)的Matlab仿真實現(xiàn)

吳應(yīng)山，張啟燦

(四川大學 電子信息學院，四川 成都 610064)

在傅里葉變換輪廓術(shù)測量方法中，測量系統(tǒng)從一個角度投影結(jié)構(gòu)光場，系統(tǒng)中的成像裝置從另一個角度獲取由物體高度調(diào)制后的變形光場，并通過傅里葉變換、頻域濾波和傅里葉逆變換恢復(fù)出物體的高度信息。FTP測量方法中，關(guān)于傅里葉變換，相位展開的相關(guān)知識，涉及大量復(fù)雜的數(shù)學運算使其抽象難以理解。針對這一問題，文中設(shè)計了基于Matlab的仿真實驗。根據(jù)相應(yīng)原理編寫了仿真代碼，運用Matlab的數(shù)學運算和可視化功能，模擬整個測量流程，完成了仿真實驗。這有助于學習和理解FTP測量方法原理及相關(guān)知識。

傅里葉變換輪廓術(shù)；Matlab；相位展開

隨著計算機技術(shù)的發(fā)展，三維數(shù)字化技術(shù)逐步成熟，并廣泛用于各個領(lǐng)域[1-3]，傅里葉變換輪廓術(shù)是1983年由M.Takeda和K.Mutoh將傅里葉變換用于三維物體的測量中而提出的三維面形測量技術(shù)[4]，只需采集一或兩幅變形條紋圖，就可以進行三維重構(gòu)，實現(xiàn)對物體輪廓的測量。所以測量速度快、易于實現(xiàn)，適合測量動態(tài)物體的三維面形[5-7]。在實際應(yīng)用中，可使用MFC與Opengl開發(fā)三維測量軟件[8]進行測量。Matlab是Mathworks公司推出的一款數(shù)學軟件，它集數(shù)值分析、矩陣計算、信號處理和圖形顯示于一身，是一種簡潔、高效的編程工具[9]。運用Matlab Gui設(shè)計數(shù)據(jù)分析界面[10]，可提高數(shù)據(jù)分析的效率。

1 傅里葉變換輪廓術(shù)

FTP(Fourier Transform Profilometry)測量光路如圖 1所示，投影儀和CCD的光軸相交于O點，O為參考平面的原點，建立如圖1所示的坐標軸。A和B分別是投影儀的出瞳和CCD的入瞳，d為A與B間的距離，L是B到參考平面R的距離，h為物面上H點到R的距離，光心連線AB與參考平面R平行，在FTP中，當投影正弦光柵產(chǎn)生的結(jié)構(gòu)光場經(jīng)過物體表面的反射，由成像系統(tǒng)得到的變形條紋可表示為

I(x,y)=a(x,y)+b(x,y)cos(2πf0x+φ(x,y))

(1)

其中，a(x,y)是背景光強；b(x,y)是條紋對比度；f0是投影光柵的頻率；j(x,y)是由物體高度分布h(x,y)引起的相位調(diào)制。為得到相對高度分布和消除系統(tǒng)測量誤差，對參考平面進行測量，對于物體高度h(x,y)=0的參考變形條紋可表示為

圖1 FTP光路圖

I0(x,y)=α(x,y)+b(x,y)cos(2πf0x+φ0(x,y))

(2)

在圖1中，AC投影到參考平面的位置為C，經(jīng)光線CB在CCD陣列上成像，光強為I0，調(diào)制相位為φ0(x,y)，當投影在被測物體上時，對同一條光線AC而言，由于光線CB和物面產(chǎn)生交點H，通過CCD觀察到的H點的光強為I，相當于參考面上D點在CCD中的光強，由物體高度h(x,y)產(chǎn)生的調(diào)制相位為φ(x,y)，則二者的相位差為[1]

(3)

由ΔAHB與ΔDHC的相似關(guān)系與式(3)得到

(4)

為便于進行傅里葉變換，將式(1)和式(2)改寫為指數(shù)形式

I(x,y)=a(x,y)+q(x,y)ej2πf0x+q*(x,y)e-2πf0x

(5)

I0(x,y)=a(x,y)+q0(x,y)ej2πf0x+q0*(x,y)e-2πf0x

(6)

上式中,*表示取共軛。對式(5)進行二維傅里葉變換可得

F(u,v)=A(u,v)+Q(u-f0,v)+Q*(u+f0,v)

(7)

其中，A(u,v)對是a(x,y)的二維傅里葉變換，Q(u-f0,v)，Q*(u+f0,v)分別是q(x,y)和q*(x,y)的二維傅里葉變換，其頻譜示意圖如圖2所示。

圖2 傅里葉頻譜分布示意圖

圖2中虛線包圍的頻譜為包含物體高度變化的有用信息，選用適當?shù)膸V波窗口將基頻分量Q(u-f0,v)濾出，再計算其傅里葉逆變換。即

(8)

其中，F(xiàn)-1表示二維傅里葉逆變換，同樣對參考條紋也進行二維傅里葉變換并濾波，傅里葉逆變換處理后可得

(9)

通過式(8)與式(9)得到s(x,y)與s0(x,y)進行下列運算，獲取相位差

(10)

式中，atan{·}表示反正切函數(shù)；imag[·]表示取復(fù)數(shù)虛部；real[·][11]表示取復(fù)數(shù)實部，得到的相位差通過相位展開技術(shù)，獲取連續(xù)相位，根據(jù)式(4)恢復(fù)出物體的高度信息h(x,y)。

2 相位展開方法

假定滿足抽樣定理的要求，即抽樣頻率大于最高空間頻率的2倍。則任何兩個相鄰抽樣點之間的連續(xù)相位差變化在(-π,π)內(nèi)，用Matlab中的atan2函數(shù)將式(10)改寫為

(11)

經(jīng)轉(zhuǎn)換后得到的二維相位呈現(xiàn)不連續(xù)的截斷分布，須用相位展開技術(shù)，消除相位分布中的截斷跳變展開成連續(xù)的相位，才能恢復(fù)高度信息。

1982年Itoh分析了一維相位展開問題，建立了嚴謹?shù)臄?shù)學模型，并且得到包裹相位圖的相位差分的再包裹就是實際相位圖的相位差分，相位展開實際上是對包裹相位差分再包裹值的積分過程[12-14]，此思想就是路徑跟蹤算法的基礎(chǔ)。在無噪聲或者其他干擾的理想情況下，只需沿著相位圖的行和列，逐點將相位展開即可。對于二維相位圖x方向展開后的相位值用up表示[15]

up(xk,y)=wp(xk,y)-2nkπ，nk∈Z

(12)

其中，wp(xk,y)是展開前的截斷相位；(xk,y)為采樣點的像素坐標，根據(jù)任何兩個相鄰抽樣點之間的連續(xù)相位變化<π，即

-π

(13)

根據(jù)式(12)與式(13)，得到式(14)

-π<[wp(xk,y)-wp(xk-1,y)]-2nkπ+2nk-1π<π

(14)

將式(14)加π,再除以2π得到

0<[wp(xk,y)-wp(xk-1,y)]/2π+0.5-nk+nk-1<1

(15)

對不等式各項同時向下取整得

nk=floor{[wp(xk,y)-wp(xk-1,y)]/2π+0.5}+nk-1

(16)

floor是Matlab中的向下取整函數(shù)，根據(jù)此分析，假定第一點截斷相位為已知的相位，則n0=0，編寫uphase1d函數(shù)，輸入不連續(xù)的截斷相位wp，輸出連續(xù)相位up。

function up = uphase1d(wp)

n= zeros(size(wp));

for i=2:length(wp)

n(i) = floor((wp(i)-wp(i-1))/(2*pi)+0.5)+n(i-1);

end

up = -2*pi*n + wp;

end

同理，上述方法對也對y方向的相位展開適用。假定二維截斷相位圖(1,1)點的相位為已知相位，由此點完成第一行的相位展開，然后再以展開的第一行各點為基準點，逐列展開，最后得到連續(xù)分布的二維相位。

function up = uphase2d(wp)

[r,c]= size(wp);

wp(1,:) = uphase1d(wp(1,:));%第一行的相位展開

for n = 1:c

wp(:,n) =uphase1d(wp(:,n));%各列的相位展開

end

up = wp;

end

3 計算機仿真實驗

3.1 FTP仿真實驗流程圖

根據(jù)傅里葉變換輪廓術(shù)的基本原理，設(shè)計仿真實驗流程圖如圖3所示。

3.2 在Matlab中仿真FTP測量方法

根據(jù)圖1所示，設(shè)置系統(tǒng)參數(shù)L=800 mm，d=350 mm，模擬的被測物體大小為512 pixel × 512 pixel的peaks函數(shù)，如圖4所示。

圖4 模擬物體

投影周期為p=8 pixel正弦參考條紋，當參考條紋投影到被測物體時，由于物體高度的變化，使參考條紋的相位發(fā)生變化將發(fā)生形變。分別對參考條紋與變形條紋做二維傅里葉變換，將得到它們的頻譜，用鼠標選取濾波范圍，將自動生成矩形濾波窗。根據(jù)式(10)獲取截斷相位差，如圖5所示。

圖5 截斷相位

圖6 展開后的連續(xù)相位

調(diào)用uphase2d函數(shù)將二維截斷相位展開成二維連續(xù)相位，如圖6所示。根據(jù)式(4)進行恢復(fù)物體高度信息，還原后的物體如圖7所示。

圖7 恢復(fù)后的物體

圖8 誤差分布圖

用恢復(fù)后的物體減去模擬的物體，得到其誤差的分布圖如圖8所示。其誤差分布較大，是由于用矩形濾波窗進行頻譜濾波時，會產(chǎn)生頻譜泄漏造成的。恢復(fù)后最大誤差為2.73 mm，平均誤差為0.154 mm。

FTP仿真實驗的Matalb代碼如下

N = 512;L = 800;d = 350;p = 8; f0=1/p; a = 127.5; b = 127.5;

x = 0:1:N-1;[X,Y] = meshgrid(x);object =5.6*peaks(length(x));

I0 = a+b*cos(2*pi*f0*X);I = a+b*cos(2*pi*f0*(X +(d/L*object)));%產(chǎn)生變形條紋

F0 = fftshift(fft2(I0)); F = fftshift(fft2(I)); %二維傅里葉變換

rect_window_filter= zeros(size(F));

figure,contour(abs(F),80);r = round(getrect);%用鼠標選取濾波范圍

rect_window_filter(r(2):r(4)+r(2),r(1):r(3)+r(1)) = 1; %生成矩形濾波窗

Q = rect_window_filter.*F; Q0 = F0.*rect_window_filter;%濾波

s0 = ifft2(ifftshift(Q0)); s = ifft2(ifftshift(Q));%逆傅里葉變換

Deltaphi= atan2(imag(s.*conj(s0)),real(s.*conj(s0))); %獲取相位差

up = uphase2d(Deltaphi); %截斷相位展開成連續(xù)相位

recovery_object = L*up./(2*pi*f0*d-up); %恢復(fù)物體高度

4 結(jié)束語

本文對FTP的基本原理進行了介紹、并對理想條件下相位展開方法進行了推導(dǎo)，文中的仿真實驗采用的是遠心光路，平行光投影方式，用鼠標選取濾波范圍生成矩形濾波窗進行濾波。在相位展開方法上，假定無噪聲無干擾、滿足抽樣定理的理想條件下，根據(jù)相鄰像素相位差為依據(jù)進行相位展開。在Matlab中恢復(fù)了模擬物體的三維形貌，完成了整個仿真實驗。

[1] 蘇顯渝,李繼陶.信息光學[M].北京:科學出版社,2013.

[2] 蘇顯渝,張冠申,陳澤先,等.鞋楦三維面形光電自動測量系統(tǒng)[J].光學工程,1989(6):1-5.

[3] 金國藩,李景鎮(zhèn).激光測量學[M].北京:科學出版社,1998.

[4] Takeda M,Mutoh K.Fourier transform profilometry for the automatic measurement of 3-D object shapes[J]. Applied Optics,1983,22(24):3977-3982.

[5] 張啟燦,蘇顯渝.動態(tài)三維面形測量的研究進展[J].激光與光電子學進展,2013,50(1):1-14.

[6] 張啟燦.動態(tài)過程三維面形測量技術(shù)研究[D].成都:四川大學,2005.

[7] 李劍,蘇顯渝,陳峰,等.基于傅里葉變換輪廓術(shù)的動態(tài)爆轟過程研究[J].激光雜志,2005,26(2):47-48.

[8] 廖宏剛,劉榮.基于OpenGL的三維真實感地形的實現(xiàn)[J].電子科技,2013,26(9):164-165.

[9] 胡新艷,霍文曉,李愛濤.Matlab在《信號與系統(tǒng)》課程教學中的應(yīng)用[J].科技信息,2009(27):136-137.

[10] 李婭麗,忻尚芝,鄭春雷.基于Matlab/GUI的數(shù)據(jù)分析界面設(shè)計[J].電子科技,2016,29(4):88-91.

[11] 劉衛(wèi)國.Matlab程序設(shè)計與應(yīng)用[M].北京:高等教育出版社,2006.

[12] Itoh K.Analysis of the phase unwrapping algorithm[J].Applied Optics,1982,21(14):2470-2470.

[13] 吳明云.二維相位展開算法的研究[D].天津:天津大學,2012.

[14] 肖楓. InSAR相位解纏算法的研究[D].上海:同濟大學,2008.

[15] Macy W W. Two-dimensional fringe-pattern analysis[J].Applied Optics,1983,22(23):3898-3901.

Fourier Transform Profilometry Simulation in Matlab

WU Yingshan,ZHANG Qican

(School of Electronic Information, Sichuan University, Chengdu 610064, China)

In the method of Fourier transform profilometry (FTP), the measurement system obtains the deformation of the object height modulated light field by the imaging device, and restores the height information of the object through the Fourier transform, frequency domain filtering and Fourier inverse transformation. The Fourier transform and phase unwrapping knowledge in the FTP measurement method involves a large number of complex mathematical calculations, making it abstract and difficult to understand. In view of this situation, we design simulation experiments based on Matlab and prepare the corresponding simulation codes. The entire measurement process is simulated by using Matlab numerical calculation and visualization function, which helps to understand the principles of the FTP measurement.

Fourier transform profilometry; Matlab; phase unwrapping

2016- 08- 31

國家重大儀器專項基金(2013YQ490879)

吳應(yīng)山(1985-)，男，碩士。研究方向：動態(tài)三維測量。張啟燦(1974-)，男，博士，教授，博士生導(dǎo)師。研究方向：三維傳感等。

10.16180/j.cnki.issn1007-7820.2017.06.003

TN29

A

1007-7820(2017)06-009-04