應用CUDA的數字全息三維粒子位移場并行重建及實現

朱 崗,楊 巖

(重慶理工大學 機械工程學院， 重慶 400054)

傳統并行計算需要使用者有較高的專業素質，多數是命令行的操作，不適宜于普通研究人員的開發。2006年，英偉達公司推出了統一計算設備架構CUDA(compute unified device architecture)，它是以GPU(graphics processing unit)作為數據并行計算設備的軟硬件體系。與傳統的網格和集群等并行架構相比，GPU的并行不僅表現出節能、體積小、性價比高等優勢，而且加速比更優。目前，CUDA已廣泛應用于科研和工程應用等領域，包括視頻與音頻處理、醫學成像、物理效果模擬、流體力學模擬、光線追蹤、CT圖像重組、石油天然氣勘探、產品設計、計算生物學和地震分析等[1～7]。

在數字全息測量中，對特定平面波前的重建是通過計算機對光波的模擬衍射來實現的,特別是在如流場分析等多焦平面對象的重建過程中，必須對每個不同記錄距離的平面進行計算。重建平面越多，其軸向分辨率越高，計算結果越準確，但是整個計算過程也更加耗時。為了提高計算速度，使重建過程能達到實時重建的目標，國內外研究人員已經采用了多種手段來提高計算速度。在單平面重建方面，Shimobaba等[8]利用GPU進行單平面的實時重建，獲得一幅512像素×512像素全息圖的重建像的時間為41.6 ms。朱竹青等[9]構建了實時數字全息再現系統，以旋轉的骰子為對象，實時拍攝全息圖，并進行單平面重建，重建幀速為20幀/s，重建圖大小為512像素×512像素。

在多平面快速重建方面，Abe等[10]采用16個FPGA實現全息圖的快速重建，其系統能在3.3 s時間內對1 024像素×1 024像素全息圖進行100個平面的重建。Cheong等[11]通過硬件加速與軟件優化相結合的方法來優化系統，重建一副200像素×200像素的圖像僅需要3 ms，若采用3個獨立的計算線程則處理速度可以達到8 幀/s。

本研究的目標是利用GPU的并行計算能力開展數字全息中多焦平面目標的快速重建，并將其應用于流體流場的實際測試中。

1 基本原理

1.1 數字全息基本原理

數字全息包括光學記錄和數字再現2個階段。數字全息光學記錄是利用CCD記錄參考光波與通過物體的反射或透射物光波干涉形成的全息圖。在光學重建中，若采用參考光的共軛光波照射全息圖，則實像出現在原物體的位置，而虛像出現在CCD相反方向的相同距離處。菲涅爾-基爾霍夫衍射實像平面的光場分布為[12]

O(ξ,η)=

(1)

(2)

重建平面與全息平面位置關系如圖1所示。

在數字重建中有菲涅爾衍射重建算法、卷積重建算法、角譜重建算法和小波重建算法等多種方法。卷積重建算法由于分辨率與全息圖相同，無衍射距離限制，計算速度快，因此受到廣泛應用。

令：

g(ξ,η,x,y)=

(3)

那么，式(1)則變為

O(ξ,η)=

(4)

在數字重建中，數字脈沖響應函數為

g(k,l)=

(5)

這樣，利用卷積的性質有：

(6)

圖1 重建平面與全息平面位置關系

1.2 并行計算

CUDA 是一種基于新的并行編程模型和指令集架構的通用計算架構，其結構上包括主機端與設備端。CPU被稱為主機端，GPU被稱為設備端，GPU作為CPU的協處理器。在計算任務分配中，在CPU上執行邏輯性強的串行計算任務，在GPU上執行大量的并行計算任務。

設備代碼在設備端GPU上執行，被稱為內核(kernel)，內核不是一個完整的程序，而是任務中能分解為并行執行的步驟的集合。每個內核函數包括2個層次的并行，即網格中線程塊間的并行和線程塊內線程間的并行。當運行至設備代碼時，調用內核函數，啟動多個線程。每個線程塊都包括多個線程，對不同的數據并行執行同一指令(SIMD)，實現線程塊內的并行。同時，在每個線程網格中，不同線程塊間也并行執行。每個線程具有私有的寄存器和板載局部存儲器，同一線程塊內線程具有片內共享儲存器。所有線程均能訪問全局、常量和紋理存儲器。GPU包括多個流多處理器(SM)，每個流多處理器又包含大量的流處理單元(SP)。線程以塊為單位分配到流多處理器上，隨后線程塊內的線程又被發送到流處理單元。線程調度器可以自動根據線程塊的數量和每個線程塊的線程數將1個或多個線程塊分配到1個流多處理器(SM)中。流多處理器(SM)中以32個線程組成一個warp作為最小的調度和執行單位。CUDA編程模型如圖2所示。

圖2 CUDA編程模型

2 實現方法

2.1 三維粒子位移場串行計算

利用CPU進行流體位移場計算的基本步驟如下：① 在獲得了t1時刻的全息圖之后，需要在光軸方向進行多個平面的數值重建，對每個重建像均需要進行圖像濾波、二值化、粒子三維空間位置提取。同時，在完成多個平面的粒子提取之后，還要進行重復粒子刪除處理，以獲得t1時刻粒子三維空間位置。② 對t2時刻的全息圖進行相同的處理，獲得t2時刻粒子三維空間位置。③ 對相鄰2個時刻的2幅全息圖所提取的粒子進行粒子配對，獲得t1和t2時刻粒子的位移矢量場。通過除以拍攝間隔時間，還可以獲得三維速度矢量場。具體流程如圖3所示。從該流程中可以發現整個處理過程是利用CPU串行執行的，這樣將會花費大量的計算時間，無法達到程序的實時處理要求。

圖3 位移矢量場獲取流程

2.2 利用GPU實現程序并行化

在流體的位移矢量場計算過程中，在獲得某時刻的全息圖后，需要利用數值計算程序對不同重建距離的圖像進行數字聚焦。不同距離的多個重建圖像的數值計算過程相互獨立，為并行計算的開展提供了可行性。為了實現流體中位移場的快速并行計算，除了需要采用GPU硬件之外，還需要對原有Matlab程序進行并行化處理，以便發揮GPU的并行計算能力。整個并行計算過程主要包括執行粗粒度的并行與細粒度的并行2個層次。

2.2.1 粗粒度的并行化

粗粒度的并行主要包括1張全息圖多個不同位置的圖像并行處理與2個連續全息圖序列執行的并行。如果計算卡具有多個GPU,則可以利用多線程技術控制多個GPU，每個GPU又分別設置多個流。即利用狀態機實現多個流同時進行不同距離的多個平面的數據處理，實現三維粒子速度場并行重建。具體執行過程如下：

1) 利用第1個線程的hallo流在第1個GPU上對t1時刻所獲得的全息圖進行二維離散傅里葉變換等預處理。

2) 將結果以異步方式發送給其他線程控制的GPU。

3) 多個流在多個GPU上進行多個不同位置的圖像并行處理(包括圖像重建、濾波、二值化、連通域處理等)。

4) 各個流等待同步。

5) 利用第1個GPU的hallo流對t2時刻所獲得的全息圖進行二維離散傅里葉變換等預處理。同時，利用第2個GPU的“0”號流完成t1時刻全息圖的重復非聚焦粒子刪除、粒子配對等后處理。

6) 第1個線程的hallo流，將t2時刻全息圖的預處理結果以異步方式發送給其他線程控制的GPU。

7) 各個GPU的多個流(除第2個GPU上的“0”號流外)進行多個不同位置的圖像并行處理。

8) 第2個GPU上的“0”號流完成t1時刻的后處理后，首先查詢t2時刻的多個不同位置的圖像并行處理是否完成。如果完成，則所有流同步；如果未完成，則轉入查詢t2時刻全息圖的預處理結果是否收到。如果t2時刻的預處理尚未收到，則等待；如果預處理已收到，但是多個不同位置的圖像的并行處理尚未完成，則其加入圖像的并行處理。

9) 各個流等待同步。

10) 利用第1個GPU的hallo流對t3時刻所獲得的全息圖進行二維離散傅里葉變換等預處理。同時，利用第2個GPU的“0”號流完成t2時刻全息圖的重復非聚焦粒子刪除、粒子配對等后處理。

11) 重復執行上述步驟，完成全息圖序列的處理。

具體流程如圖4所示。

圖4 利用狀態機實現位移場重建程序設計

2.2.2 細粒度的并行化

細粒度的并行主要是指將核函數分為多個Block處理，利用大量的核并行計算。細粒度的并行化需要將濾波、二值化、連通域提取、刪除非聚焦重復粒子等多個原有串行執行的程序進行并行處理，以便能在GPU上運行。

這里，以刪除非聚焦重復粒子為例，說明如何對串行程序進行并行化處理。利用CPU從tn時刻全息圖提取到粒子位置及大小信息，這些粒子既包括了聚焦粒子，也包括了離焦粒子。為了將離焦粒子刪除掉，按照粒子重建圖像的規律，以二值化后獲得的粒子大小信息為判斷準則。如果2顆粒子的三維空間位置的歐式距離小于3倍粒子直徑，則認為它們是同一顆粒子的聚焦像或者離焦像，并認為粒子面積較大的粒子為聚焦粒子，而面積較小的粒子為非聚焦粒子，此時保留聚焦粒子的信息，刪除非聚焦粒子的信息。如果采用CPU對重復粒子進行刪除，那么串行算法是從第1顆粒子開始，往后進行逐一比較，直到與最后一顆粒子進行比較，完成第1顆粒子的比較循環，然后繼續從第2顆粒子開始依次跟后面的粒子進行比較。如果獲得了n顆粒子，若利用CPU串行計算，其循環次數為M，其計算效率非常低。

M=(n-1)·(n-2)·(n-3)…3·2·1=

(n-1)!

(7)

如果采用GPU進行并行計算，可以將n顆粒子分為J組，每組有N顆粒子。這樣每一組的粒子序號分別為1～N,N+1～2N,…(J-1)N～JN。將粒子序號按組從左向右排列，同時將粒子序號從上往下排列，形成一個二維矩陣，如圖5所示。這樣橫坐標中任意一組(第x組)便可以跟縱坐標中的任意一組(第y組)同時進行比較。具體運行中將其化為一個Block塊，其總的Block塊數為

B=J+(J-1)+(J-2)…+3+2+1=

(8)

在一個Block內，橫坐標為第x組，其粒子序號為(x-1)N～xN，縱坐標為第y組,其粒子序號為(y-1)N～yN，能分別進行同時比較。具體在GPU運行中，可以利用一個SM中的多個ALU進行并行計算。這樣可以大大提高計算效率。

圖5 刪除重復粒子并行計算模型

3 實驗

為了分析重建平面數量、分辨率對重建速度的影響，分別配置了基于CPU和GPU的3種硬件系統，以比較其計算速度。具體實驗裝置配置如下：

1) Window7系統，Matlab r2014a軟件(用于CPU測試),CPU型號為Inter(R) Core(TM) i7 CPU 930(多核并行)，48G內存；

2) 1片Nvidia Tesla 40C；

3) 1片Nvidia Tesla K80。

K40C具有1個開普勒架構核心(GK110)，開啟了15組流多處理器，具有2 880個流處理器，專用存儲器總容量為12 GB，顯存帶寬為288 GB/s，其單精度計算能力可達4.29萬億次/s，而雙精度計算能力可達1.43萬億次/s。K80采用了2個開普勒架構核心(GK210)，每一個核都只開啟了15組流多處理器陣列中的13組，即K80具有4 992個流處理器，顯存是2組384位的GDDR5，容量為24 GB，其單精度計算能力可達8.74萬億次/s，而雙精度計算能力可達2.91萬億次/s。

實驗主要從重建平面數的變化及重建像分辨率的變化2方面進行。

3.1 重建平面數對重建速度的影響

固定分辨率為512像素×512像素，重建平面數從500幅開始，每次增加100幅，直到1 000幅。通過CPU與GPU分別進行多焦平面重建，分別獲取其運行時間，得出GPU的加速效果，即加速比。運行時間如表1所示。運行時間曲線如圖6所示。

表1 不同重建平面數量3種系統運行時間對比 s

由表1可知：在不同的重建平面數時，采用GPU并行計算的效率均遠高于CPU，特別是在重建平面數為1 000時，K80、K40C和CPU930的重建時間分別0.866、1.8和180.733 s。隨著重建幅數的增加，K40C加速比從84倍增加到100倍，而K80含2片GPU，對K40C的加速比維持在2倍左右。每增加100幅，CPU運算時間延長20～30 s，而K40C在大約200 ms以內，K80在大約100 ms以內，差異浮動比較小。

圖6 3種系統運行時間與重建平面數的關系

3.2 重建像分辨率對重建速度的影響

重建平面數取為512幅，實驗采用的全息圖分辨率為1 296像素×966像素，故實驗選擇從512像素×512像素分辨率開始，每次增加128像素×128像素，直到896像素×896像素。不同分辨率對重建時間的影響如表2所示。運行時間曲線如圖7所示。

表2 不同分辨率對重建時間的影響 s

由表2可知：隨著全息圖的分辨率大小的增加，CPU運算時間急劇增加；GPU相對變化差異小，尤其是K80。

綜上所述，GPU相對于CPU來說，不僅運算速度快，并且運算性能更平滑，時間范圍更具可控性。

3.3 維粒子位移場并行計算測試

實驗光路如下：氦氖氣體連續激光器發出波長為632.8 nm的激光(索雷博HNL210L-EC)，通過針孔濾波器進行濾波，再通過準直鏡(西格瑪 DLB-30-300PM)形成平行光，通過裝有純凈水的樣品池，其光程為10 mm，在其中散布有直徑為51.7 μm的三聚氰胺甲醛樹脂微粒(武漢華科微科)。其同軸全息圖被CCD(BOBCAT ICL-B1310)接收，CCD的分辨率為1 296像素×966像素，像元3.75 μm，其記錄距離為21.5～31.5 mm。圖8為拍攝的三聚氰胺甲醛樹脂微粒的同軸全息光路示意圖，圖9為某時刻拍攝到的三聚氰胺微粒的數字全息圖。

圖7 3種系統運行時間與重建分辨率的關系

圖8 同軸粒子場光路示意圖

圖9 三聚氰胺全息圖

將并行程序用于實際三維流體的全息圖的快速重建實驗中。利用個人筆記本電腦進行計算，其處理器為Core i5-4210U@1.7 GHz，內存為8 GB，操作系統為Windows 10，軟件為Matlab(r2014a)。截取512像素×512像素的全息圖，完成512個重建平面的速度場計算需要10.9 s。為了實現快速計算，利用計算統一設備架構，將工作站作為主機,選用K80作為設備端，其處理器為Xeon E5-2630@2.3 GHz，內存為64 GB，操作系統為Windows 7，軟件為Visual studio 2014。采用單線程單GPU計算時間需要0.46 s，而采用雙線程雙GPU則其計算時間為0.268 s。運行時間如表3所示，獲得的位移場分布如圖10所示，高度方向已折算為10 mm粒子深度范圍所對應的像素數量。

表3 個人電腦與帶GPU的工作站計算時間比較 s

圖10 微粒的位移矢量圖

通過以上結果可知：利用并行重建方法能正確獲得流體中示蹤微粒位移矢量場；在圖像分辨率為512像素×512像素、重建512個平面時，利用K80采用雙線程雙GPU的速度比采用CPU的速度快40倍，其全息圖的處理速度達到3.7幀/s。

4 結束語

利用基于CUDA的架構將原有串行的重建步驟進行二級并行化處理。同時結合多線程技術實現數字全息圖粒子3維位移場的快速重建。并將其應用到實際獲得的全息圖序列中，其處理的單幅全息圖為512像素×512像素，重建512個平面的時間為268 ms，即3.7幀/s，運算速度為采用CPU進行串行計算的40.6倍。

通過測試發現，每個GPU設置不同的流的數量對實驗結果也有影響。分別將流的數目設置為4、8和16，發現當流的數目設置為8時，其運算速度最快。另外，如果減少重建平面的數量，增加GPU的數量或者采用更高性能的計算卡，那么其處理全息圖的幀率可以得到進一步提升，完全可能達到對全息圖進行實時處理的要求。

參考文獻：

[1] 盧文龍，王建軍，劉曉軍.基于CUDA的高速并行高斯濾波算法[J].華中科技大學學報(自然科學版),2011，39(5)：10-13.

[2] 朱奭，常晉義.基于CUDA的數字化放射圖像重建算法[J].計算機應用研究，2014，31(5)：1577-1580.

[3] 江順亮，黃強強，董添問，等.基于CUDA的離散粒子系統模擬仿真及其實現[J].南昌大學學報(工科版)，2011，33(3)：290-294.

[4] 周煜坤，陳清華，余瀟.基于CUDA 的大規模流體實時模擬[J].計算機應用與軟件，2015，32(1)：143-148.

[5] 岳田爽，趙懷慈，花海洋.基于CUDA的光線追蹤優化算法研究與實現[J].計算機應用與軟件，2015，32(1)：161-162.

[6] 王玨，曹思遠，鄒永寧.利用CUDA技術實現錐束CT圖像快速重建[J].核電子學與探測技術，2010，30(3)：315-320.

[7] 劉金碩.基于CUDA的并行程序設計[M].北京：科學出版社，2016.

[8] SHIMOBABA T,SATO Y,MIURA J,et al.Real-time digital holographic microscopy using the graphic processing unit[J].Opt Express,2008,16(16):97035.

[9] 朱竹青，孫敏，聶守平,等.基于GPU 的數字全息實時再現系統設計及實驗研究[J].光電子·激光，2012，23(3)：325-329.

[10] ABE Y,MASUDA N,WAKABAYASHI H,et al.Special purpose computer system for flow visualization using holography technology[J].Opt Express,2008,16(11):94254.

[11] CHEONG F C,SUN B,DREYFUS R,et al.Flow visualization and flow cytometry with holographic videomicroscopy[J].Opt Express,2009,17(15):109416.

[12] ULF S,WERNER J.Digital Holography[M].Heidelberg：Springer，2005.