OpenMP并行程序在數字全息三維重構中的應用

馬 靜，邸江磊，肖 鋒

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

0 引 言

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

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

文中基于數字全息技術，通過卷積算法，實現了數字全息三維重構，針對重構計算量大、重構時間過長等問題，利用OpenMP并行程序對數字全息三維重構過程進行了提速，并對其進行驗證。

1 基于卷積的數字全息重構原理

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

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

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

(1)

利用解卷積的方法就可以重構原物像，即

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

(2)

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

設照射到全息平面上的物光波和參考光波分別為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記錄離散化的光強分布。設CCD的像素數為Nx×Ny，光敏面尺寸為Lx×Ly，則像元尺寸分別為Δx=Lx/Nx，Δy=Ly/Ny，CCD記錄的離散圖像為：

(4)

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

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

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)為脈沖響應函數。

脈沖響應函數與具體的物理光路有關，計算方法如下：

g(x,y,xi,yi)=

g(xi-x,yi-y)

(6)

其中，λ為波長。

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

2 數字全息三維重構過程

數字全息三維重構過程如圖1所示。數字全息處理的對象為重構硬件系統采集的全息圖。根據重構硬件系統組成的不同，對全息圖的處理也不同。

圖1 數字全息三維重構過程

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

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

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

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

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

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

(7)

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

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

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

3 三維重構的并行處理

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

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

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

在進行數字全息三維重構時，一般有兩種應用情況，一種是針對單幅圖像的重構，另一種是針對視頻圖像的重構。一種并行設計思想是針對不同任務進行并行處理，另一種思想是針對某一圖像重構的不同任務進行并行處理。為了兼顧對圖像和視頻的重構任務，采用第二種并行設計思想。就是對每一幅圖像數字全息三維重構的各個環節采用并行算法進行處理。

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

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

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

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

dynamic模式也是為每個內核分配size次循環任務，但是分配過程是動態的，哪個內核處于空閑狀態就給它分配size次循環任務。

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

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

表1 測量結果

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

4 處理結果

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

圖2 分辨率板全息圖

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

圖3 分辨率板重構圖像

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

5 結束語

利用OpenMP技術，根據卷積的離軸數字全息重構的基本原理，設計并開發了數字全息三維重構軟件。采用并行運算提高了三維重構的運行速度，并對OpenMP實現并行運算的方式進行了對比分析，結果表明各種并行模式的平均重構時間基本相同，但是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] 趙 潔,王大勇，李 艷，等.數字全息顯微術應用于生物樣品相襯成像的實驗研究[J].中國激光,2010,37(11):2906-2911.

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

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

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

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

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

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

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

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