馬 靜，邸江磊，肖 鋒

(1.西安工業(yè)大學(xué) 計算機科學(xué)與工程學(xué)院，陜西 西安 710032;2.西北工業(yè)大學(xué) 理學(xué)院，陜西 西安 710072)

0 引 言

數(shù)字全息術(shù)是顧德門[1]在1967年提出的一種全新的全息成像方法，以CCD等光電耦合器件取代傳統(tǒng)的干版記錄全息圖，并由計算機以數(shù)字的形式對全息圖進行再現(xiàn)。隨著全息技術(shù)的發(fā)展，數(shù)字全息技術(shù)在顯微三維物體表面形變測量、粒子場測試、流場測定[2-3]、生物樣品顯微測量等方面有了新的研究[4-5]。數(shù)字全息不需要成像透鏡，靈敏度高，是一種比較理想的微小物體數(shù)字三維測量方法[6-7]。

數(shù)字全息將全息圖記錄的物光波和參考光波相干涉疊加時產(chǎn)生的干涉條紋數(shù)字化，然后在計算機中通過再現(xiàn)算法進行處理，從而得到物體的三維強度像和相位像。早期的單次傅里葉變換算法應(yīng)用較為廣泛，該算法使再現(xiàn)像素大小與再現(xiàn)距離呈正比例變化。在光學(xué)的標量衍射理論和傅里葉變換方法的基礎(chǔ)上，又發(fā)展出了三種最常用的再現(xiàn)算法：菲涅耳變換法、卷積法和角譜法[8]。

文中基于數(shù)字全息技術(shù)，通過卷積算法，實現(xiàn)了數(shù)字全息三維重構(gòu)，針對重構(gòu)計算量大、重構(gòu)時間過長等問題，利用OpenMP并行程序?qū)?shù)字全息三維重構(gòu)過程進行了提速，并對其進行驗證。

1 基于卷積的數(shù)字全息重構(gòu)原理

數(shù)字全息重構(gòu)是利用電荷耦合器件CCD獲得干涉全息圖后，模擬生成參考光，然后利用計算機技術(shù)和相關(guān)算法對全息圖進行處理，恢復(fù)物體的振幅和相位信息，重構(gòu)物體的三維模型。通過數(shù)字全息重構(gòu)，可以直接得到被記錄物體再現(xiàn)像的復(fù)振幅分布，從而獲得被記錄物體的表面亮度和形貌分布。而全息中采用的CCD曝光時間短，分辨率高，記錄和再現(xiàn)過程都比經(jīng)典光全息方法更快捷。

卷積法是基于快速傅里葉變換的數(shù)值再現(xiàn)方法[9]。利用卷積法重構(gòu)圖像的尺寸不隨重建距離變化。衍射積分可以看作是物波函數(shù)與自由空間脈沖響應(yīng)函數(shù)的卷積。

g(x',y',ξ,η)=

(1)

利用解卷積的方法就可以重構(gòu)原物像，即

u'=F-1{F[h·r]·F[g]}

(2)

對于離軸全息而言，需要通過離軸物理光路，利用干涉原理將物體信息以干涉條紋的形式記錄下來，通過CCD相機以圖像的形式保存下來，形成全息圖[10]。

設(shè)照射到全息平面上的物光波和參考光波分別為O(x,y)和C(x,y)，則全息圖平面的干涉條紋強度為H(x,y)：

O(x,y)C*(x,y)+O*(x,y)C(x,y)

(3)

其中，C*(x,y)和O*(x,y)分別表示物光波和參考光波的共軛光波。

上述全息圖需要通過CCD記錄離散化的光強分布。設(shè)CCD的像素數(shù)為Nx×Ny，光敏面尺寸為Lx×Ly，則像元尺寸分別為Δx=Lx/Nx，Δy=Ly/Ny，CCD記錄的離散圖像為：

(4)

CCD記錄的離散化的干涉強度即為數(shù)字全息圖。

根據(jù)卷積理論，物光的復(fù)振幅U(x,y)可根據(jù)全息圖平面的干涉條紋強度、參考光波以及脈沖響應(yīng)函數(shù)表示為：

U(xi,yi)=F-1{F[H(x,y)C(x,y)]F[g(x,y,xi,yi)]}

(5)

其中，H(x,y)為全息圖平面的干涉條紋強度；C(x,y)為參考光波；g(x,y,xi,yi)為脈沖響應(yīng)函數(shù)。

脈沖響應(yīng)函數(shù)與具體的物理光路有關(guān)，計算方法如下：

g(x,y,xi,yi)=

g(xi-x,yi-y)

(6)

其中，λ為波長。

通過上述方法可以得到物光的復(fù)振幅[11-13]。物光的復(fù)振幅中攜帶了物體的三維信息，如果物體起伏超出了波長范圍，則需要進行解包裹，從而提取物體的三維分布信息。

2 數(shù)字全息三維重構(gòu)過程

數(shù)字全息三維重構(gòu)過程如圖1所示。數(shù)字全息處理的對象為重構(gòu)硬件系統(tǒng)采集的全息圖。根據(jù)重構(gòu)硬件系統(tǒng)組成的不同，對全息圖的處理也不同。

圖1 數(shù)字全息三維重構(gòu)過程

對全息圖的處理需要使用參考光，如果參考光是平面光，可直接用采集到的全息圖進行后續(xù)處理，如果是非平面光，則需用參考光乘以全息圖。處理后的圖像要進行傅里葉變換，從而得到其頻譜分布圖。

頻譜分布圖包含零級區(qū)域和正負一級區(qū)域。其中物體信息的頻譜分布在正負一級中心附近，而位于圖像中心最強的零級則沒有任何物體信息，因此需要進行選頻和移頻。選頻是選擇正一級或負一級頻譜區(qū)域，移頻是將正一級或負一級頻譜中心移到圖像中心。這一步需要用戶參與，為了便于用戶進行選頻操作，文中根據(jù)頻譜圖的特征實現(xiàn)了正一級或負一級頻譜中心及頻譜區(qū)域的自動生成。

正一級或負一級頻譜中心自動生成的基本思想就是令圖像的極大值為零級中心，次極大值和第三極大值為正一級或負一級頻譜中心。由于數(shù)字全息三維重構(gòu)的物理系統(tǒng)中都使用了激光，因此全息圖存在大量激光散斑，從而導(dǎo)致次極大值有可能出現(xiàn)在零級中心區(qū)域附近，而不是正一級或負一級中心。為此，首先采用均值濾波器對需要顯示的頻譜圖像進行濾波處理，然后在尋找次極大值和第三極大值，這樣可以有效消除激光散斑的影響，準確識別出正負一級中心。有了零級中心、正負一級中心，就可以生成正負一級區(qū)域。

文中生成的正負一級中心及其區(qū)域，只是為用戶提供初始依據(jù)，用戶還可以手動對其進行修正。另外需要說明的是，文中只是對顯示給用戶的頻譜圖進行了濾波處理，實際的頻譜圖并沒有改變。

根據(jù)脈沖響應(yīng)函數(shù)得到傳遞函數(shù)，將傳遞函數(shù)與移頻后的圖像相乘，對其乘積進行傅里葉逆變換。逆變換的圖像是一個復(fù)數(shù)圖像，其相位包含了物體信息，此時相位在0～2 pi之間，需要再進行解包裹處理。

解包裹采用最小Lp范數(shù)法，這種方法不依賴路徑的全局算法，計算量大，但對誤差點的控制較好。最小范數(shù)法是要得到解包裹相位Φ的相鄰像素相位差和包裹相位Ψ的相鄰像素相位差最小范數(shù)Lp。最小Lp范數(shù)解要使式(7)中的J取到極小值。

(7)

范數(shù)p的選取是關(guān)鍵。當p?2時，解包裹曲面與真實梯度面相差很大；當p<2時，結(jié)果與局部梯度較匹配，但是權(quán)重的作用變大；p=2是目前使用最多的最小范數(shù)解包裹方法，也就是最小二乘法，可將解包裹過程轉(zhuǎn)化為泊松方程，再用迭代方法求解[14-15]。

解包裹后的相位代表實際物體高程，需要根據(jù)標定值得到物體相位信息。

上述三維數(shù)字全息重構(gòu)算法是一個復(fù)雜耗時的過程。在重構(gòu)過程中，一般需要選擇合適的頻譜區(qū)域，并對重構(gòu)距離、波長等重構(gòu)參數(shù)進行設(shè)定。

3 三維重構(gòu)的并行處理

數(shù)字全息三維重構(gòu)過程要進行大量的圖像處理工作，運算量非常大，采用并行處理算法可加快運算速度，從而提高數(shù)字全息三維重構(gòu)的實時性。

OpenMP是一個在共享存儲的多處理機上編寫并行程序的應(yīng)用編程接口。對于同步共享變量、合理分配負載等任務(wù)都提供了有效的支持。

OpenMP的共享存儲模型最底層是處理器集合，這些處理器可以訪問內(nèi)存中的同一個共享位置，因而可以通過共享變量進行交互和同步。OpenMP可根據(jù)需要設(shè)置包含子句項,在沒有其他約束條件下，子句可以無序，也可以任意地選擇。#pragma omp parallel for[子句#]是最頻繁使用的編譯指導(dǎo)語句，可搭配firstprivate，if，lastprivate，private，reduction，schedule等子句使用。

在進行數(shù)字全息三維重構(gòu)時，一般有兩種應(yīng)用情況，一種是針對單幅圖像的重構(gòu)，另一種是針對視頻圖像的重構(gòu)。一種并行設(shè)計思想是針對不同任務(wù)進行并行處理，另一種思想是針對某一圖像重構(gòu)的不同任務(wù)進行并行處理。為了兼顧對圖像和視頻的重構(gòu)任務(wù)，采用第二種并行設(shè)計思想。就是對每一幅圖像數(shù)字全息三維重構(gòu)的各個環(huán)節(jié)采用并行算法進行處理。

根據(jù)數(shù)字全息三維重構(gòu)過程，采用并行處理算法的環(huán)節(jié)包括：全息圖與參考光的相乘；全息圖的傅里葉變換；頻譜圖移頻；傳遞函數(shù)的生成；傳遞函數(shù)與移頻后圖像的相乘；相位解包裹；解包裹后相位轉(zhuǎn)換為物高以及顯示圖像生成等。

上述任務(wù)的共同特點是都需要對圖像進行處理，需要采用二維循環(huán)來完成。因此文中采用針對循環(huán)運算的并行處理方法，在需要進行并行運算的for語句前增加以下OpenMP編譯指令：

#pragma omp parallel for schedule(kind,[size])

其中kind表示模式，共有三種模式，分別為static、dynamic和guided。static模式下為每個線程分配size次循環(huán)任務(wù)，按照線程數(shù)量一次性完成任務(wù)分配。如果size數(shù)值過大，則可能導(dǎo)致出現(xiàn)有的內(nèi)核任務(wù)很重，有的內(nèi)核沒有任務(wù)的情況。

dynamic模式也是為每個內(nèi)核分配size次循環(huán)任務(wù)，但是分配過程是動態(tài)的，哪個內(nèi)核處于空閑狀態(tài)就給它分配size次循環(huán)任務(wù)。

guided模式按照從大到小的方式給每個內(nèi)核分配任務(wù)，先分配較多任務(wù)，然后逐漸減小任務(wù)量，最小任務(wù)量為size次循環(huán)。

為了檢測并行算法以及各種并行運算模式的運行速度，在配有AMD Phenom II X6 1055T六盒處理器的臺式機上，對軟件運行時間進行了測試。所采用全息圖的分辨率為1 024×768，測試時對全息圖進行400次相同的重構(gòu)，計算重構(gòu)時間的均值及其方差，上述過程多次測量選擇中值，結(jié)果如表1所示。

表1 測量結(jié)果

從表1可以看出，未采用并行處理算法時，每次重構(gòu)的平均時間為0.369 s，在三種不同的并行方法處理下，重構(gòu)時間縮減到0.23 s左右，單幅重構(gòu)圖像時間提高約1.56倍。不同參數(shù)下，三種并行方法重構(gòu)的時間相差0.002～0.004 s，采用guided方式，方差相比而言略小。

4 處理結(jié)果

根據(jù)上述重構(gòu)算法，用C++語言開發(fā)了數(shù)字全息三維重構(gòu)軟件。選取一幅通過CCD采集的分辨率板全息圖，如圖2所示，該幅分辨率板全息圖像的主要參數(shù)如下：波長632.8 nm；像素尺寸4.65 μm；重構(gòu)距離118.025 mm。

圖2 分辨率板全息圖

根據(jù)上述圖像參數(shù)對分辨率板全息圖進行重構(gòu)，重構(gòu)結(jié)果如圖3所示。由于分辨率板為相位型物體，因此相位圖就已經(jīng)反映了分辨率板的形狀。圖中的三維重構(gòu)圖是沒有進行解包裹的結(jié)果。

圖3 分辨率板重構(gòu)圖像

從頻譜圖中可以看出，在零級周圍有一個正一級中心，還有一個負一級中心，無論是正一級還是負一級都可以重構(gòu)出相位信息，其對應(yīng)的重構(gòu)距離一個為正，另一個為負。

5 結(jié)束語

利用OpenMP技術(shù)，根據(jù)卷積的離軸數(shù)字全息重構(gòu)的基本原理，設(shè)計并開發(fā)了數(shù)字全息三維重構(gòu)軟件。采用并行運算提高了三維重構(gòu)的運行速度，并對OpenMP實現(xiàn)并行運算的方式進行了對比分析，結(jié)果表明各種并行模式的平均重構(gòu)時間基本相同，但是guided模式的方差相對較小。

[1] GOODMAN J W,LWRENCE R W.Digital image formulation from electronically detected holograms[J].Applied Physics Letters,1967,11(3):77-79.

[2] FLEMING C P，ECKERT J，HALPERN E F，et al.Depth resolved detection of lipid using spectroscopic optical coherence tomography[J].Biomedical Optics Express,2013,4(8):1269-1284.

[3] LUCENTE M.Interactive three-dimensional holographic displays:seeing the future in depth[J].ACM SIGGRAPH Computer Graphics,1997,31(2):63-67.

[4] DI CAPRIO G,DARDANO P,COPPOLA G,et al.Digital holographic microscopy characterization of superdirectivebeam by metamaterial[J].Optics Letters,2012,37(7):1142-1144.

[5] ALIEVA T,BASTIAANS M J.On fractional Fourier transform moments[J].IEEE Signal Processing Letters,2000,7(11):320-323.

[6] LIEBLING M,BLU T,UNSER M.Fresnelets:new multriesolution wavelet bases for digital holography[J].IEEE Transactions on Image Processing,2003,12(1):29-43.

[7] 趙 潔,王大勇，李 艷，等.數(shù)字全息顯微術(shù)應(yīng)用于生物樣品相襯成像的實驗研究[J].中國激光,2010,37(11):2906-2911.

[8] 張志會，王華英，劉佐強，等.基于快速傅里葉變換的相位解包裹算法[J].激光與光電子學(xué)進展,2012,49(12):59-65.

[9] 王華英，劉佐強，廖 薇，等.基于最小范數(shù)的四種相位解包裹算法比較[J].中國激光,2014,41(2):122-127.

[10] 郭恒光，瞿 軍．基于小波域雙譜分析的磨粒圖像多尺度形狀特征提取[J]．計算機應(yīng)用與軟件，2016，33(9)：224-226．

[11] 賴建新,胡長軍,趙宇迪,等.OpenMP任務(wù)調(diào)度開銷及負載均衡分析[J].計算機工程，2006，32(18)：58-60.

[12] 馬旭龍，林 峰.基于OpenMP的快速并行分層算法[J].計算機輔助設(shè)計與圖形學(xué)學(xué)報，2015，27(4):747-753.

[13] 鄒賢才,李建成,汪海洪,等.OpenMP并行計算在衛(wèi)星重力數(shù)據(jù)處理中的應(yīng)用[J].測繪學(xué)報,2010,39(6):636-641.

[14] 奎 因.并行程序設(shè)計C、MPI與OpenMP[M].北京:清華大學(xué)出版社,2005.

[15] 潘哲郎，李仕萍，鐘金鋼.用數(shù)字全息層析成像技術(shù)測量毛細管的內(nèi)徑及壁厚[J].光學(xué)精密工程，2013，21(7):1643-1650.