基于二級并行架構的顆粒全息圖快速重建方法

浦世亮，金其文，陳曉鋒，毛 慧，吳學成*

(1.杭州海康威視數字技術股份有限公司，杭州 310051；2.浙江大學 能源清潔利用國家重點實驗室，杭州 310027)

引 言

全息術作為一種3維照像技術，自發明以來經過記錄介質的更新換代，由光學全息發展為數字全息，并廣泛地應用在流場和顆粒場測量領域[1-2]。數字全息技術測量顆粒場的優勢在于可同時測量3維空間中的顆粒位置、形貌、粒徑參量，結合粒子圖像/軌跡測速(particle image/trackingvelocimetry,PIV/PTV)還可進一步測量速度、旋轉等參量。然而數字全息技術的重建過程需要大量的計算時間，全息圖數量越多、分辨率越高、重建截面越多，計算耗時也迅速增加。這限制了全息技術應用于工業過程中顆粒場實時在線監測[3]。

為了提高全息重建的速度，國內外的研究人員已經開展了相關的研究工作。加快重建速度一般可以從兩方面著手：一方面是采用性能更強的計算硬件，比如采用圖形處理單元(graphics processing unit,GPU)[4-7]和大規模可編程門陣列(field-programmable gate array,FPGA)[8]；另一方面是優化重建算法和并行架構，比如采用多線程編譯框架(open multi-processing,OpenMP)。其中，GPU技術由于具有開發過程更為簡單和可移植性更強的優點[9]，已經被廣泛使用，一般需要搭配NVIDIA顯卡的統一計算設備架構(compute unified device architecture,CUDA)[10]。

快速傅里葉變換(fast Fourier transform,FFT)是全息重建中的主要耗時算法，MASUDA等人[8]設計的全息PTV系統，FFT 由FPGA芯片完成，對于256×256的全息圖，重建100個截面的時間為266ms。ZHANG等人[11]比較了CPU和GPU進行全息重建的計算性能，使用OpenMP多線程技術進行算法優化，重建速度顯著增加。YAMAMOTO等人[12]發展了一種用于高速全息3維成像的專用計算機，包括4個計算模塊，重建速度相比于普通計算機提高了68倍。ZHU等人[7]基于GPU并行技術，開發了數字全息實時再現系統，每秒平均可以處理20張分辨率為512×512的全息圖。MA等人[13]將OpenMP應用于數字全息3維重構，比較了不同并行方法的加速效果，證明了OpenMP可有效提高全息重建速度。

為了滿足數字全息技術在工業應用過程中的實時性需求，快速獲取3維顆粒場信息，本文中提出了一種基于OpenMP圖片級并行和基于CUDA像素級并行的二級并行架構。以小波重建方法為例，分析了小波重建的計算流程，在此基礎上提出了計算流程改進、并行計算思路和二級并行實現方式。以煤粉顆粒全息圖為實驗對象，測試了二級并行架構的重建準確性和加速性能。

1 理論基礎與并行實現

數字同軸全息3維顆粒場測量主要包括顆粒全息圖的記錄和數值重建兩個過程。在全息記錄過程中，一束平行激光照射顆粒場，一部分光由顆粒散射作為物光，另一部分未經散射的光作為參考光，物光與參考光干涉后，由數字相機記錄形成顆粒全息圖。數字全息重建是指用計算機模擬光的傳播，采用全息重建算法獲得全息圖中的3維顆粒場信息。菲涅耳積分重建、角譜重建、分數傅里葉變換重建和小波重建是常見的全息重建方法。小波重建具有全息重建圖像信噪比高和圖像背景均勻的優點，本文中利用小波重建方法進行全息重建提取顆粒信息[14]，并對基于小波重建的并行設計進行分析介紹。根據小波母函數Ψ0(x,y)=sin(x2+y2)，構造小波函數：

(1)

式中,a為尺度參量，與重建距離有關，其表達式為：

(2)

式中，λ是激光波長。

對構造的小波函數進行適當調整，引入高斯窗函數exp[-(x2+y2)/(a2σ2)]使小波函數在偏離中心點處的值迅速趨于零，同時引入一個調零參量K,使其均值為零，K(σ)的表達式為：

(3)

式中,σ為窗函數帶寬因子，代表窗函數的寬度，依賴于全息圖采樣特性，其表達式為：

(4)

式中,N和δ分別為相機分辨率和像素尺寸，ε為一個數值極小的常數。最終得到校正的小波函數可表示為：

Ψ(x,y)=

(5)

全息圖Iholo(x,y)重建圖像的光強I(x,y,z)可以表示為：

I(x,y,z)=

(6)

Fig.1 Flow chart of reconstruction with wavelet method

全息重建是在深度方向進行切分，得到一系列不同深度z處的重建圖像。計算過程中有些變量與z有關，如a和σ，有些變量與z無關，如(x2+y2)，對于與z無關的變量可作為局部變量保存，避免多次計算。全息圖的重建過程具有二級特征，第1級為重建截面之間的獨立性，第2級為像素點之間的獨立性，因為重建截面之間和像素點之間沒有數據的交換，得到結果像素的計算過程對于每個像素是獨立的。這使得并行計算非常適用于全息圖的快速重建，而且是“圖片級+像素級”的二級并行，二級并行架構的具體實現如圖 2所示。利用OpenMP實現圖片級并行，利用CUDA實現像素級并行。除了用于單張全息圖的并行重建，當面對海量的全息圖數據而這些全息圖重建參量一致時，可以將每個重建截面的F[Ψ(x,y)]作為局部變量保存，同樣可以避免重復計算，加快批處理速度。

Fig.2 Framework of two-level parallel reconstruction

OpenMP是一個應用程序接口(application programming interface,API)，用于在多個處理器上編寫并行程序實現多線程之間的并行計算、內存共享和負載分配。OpenMP使用最頻繁的編譯指導語句為#pragma omp parallel for[子語句]，用OpenMP實現圖片級并行的主要代碼如圖3所示。主程序中firstprivate(z)表示重建深度z為私有變量，變量image為全息圖的傅里葉變換F[Iholo(x,y)]，為公有變量。CUDA并行架構程序包括主機端和設備端兩部分，分別在GPU和CPU上執行，主要流程是顯存申請(cudaMalloc)、數據拷貝(cudaMemcpy)、數據運算、數據拷貝。數據運算由核函數在GPU上運行，采用__global__ static void function進行函數定義，再用function函數可設置線程塊和線程數量。此外CUDA中有針對傅里葉變換的專用函數CUFFT，可直接調用。本文中針對小波重建編寫的一些子函數及其功能如圖3所示。

Fig.3 Some codes for OpenMP and CUDA

2 實驗與結果

2.1 實驗裝置

本文中采用煤粉顆粒全息圖作為實驗對象，測試二級并行架構的準確性和加速效果。實驗裝置如圖4所示。采用同軸數字全息系統記錄煤粉顆粒流的全息圖。激光器發出的激光束經過光強衰減、濾波、擴束、準直后照射煤粉顆粒場，CCD相機記錄煤粉顆粒全息圖，其中激光器波長為532nm，相機像素大小為5.5μm、分辨率可調。高性能工作站用于重建煤粉顆粒全息圖，工作站配備了2塊驍龍處理器(Intel Xeon E5-2680 v4)和4塊英偉達顯卡(NVIDIA Tesla P100)。實驗中獲得了一系列不同分辨率的煤粉顆粒全息圖。

Fig.4 Experimental setup

2.2 重建準確性驗證

在二級并行架構重建準確性的驗證中，比較了并行算法與單線程算法重建不同分辨率的煤粉顆粒全息圖的差異。圖5a為分辨率為1000×1000的煤粉顆粒全息圖;圖5b為單線程重建結果;圖5c為并行重建結果。重建深度z范圍為9cm至10cm，重建截面間隔為1mm。

Fig.5 Reconstructed images of coal powder hologram with different calculation method

計算圖5b和圖5c對應位置的像素值的方差，作為重建準確性的判別依據。計算得到圖5中二級并行與單線程重建結果圖的方差S=0.078，由于像素值的大小在0～255之間，所以可以基本認為圖5b和圖5c是一致的。分別處理多張不同分辨率和重建參量的全息圖，計算兩種重建算法的結果方差，結果如表1所示。表中，zmin為重建深度的最小值，zmax為重建深度的最大值，interval為重建間隔，variance為方差計算結果，隨著全息圖分辨率的增加，方差逐漸增大，但都小于0.1，表明不同方法重建圖差別極小。結果說明基于二級并行架構的全息重建程序與單線程全息重建程序在重建結果上是基本一致的，證明了二級并行架構的重建準確性。

Table 1 Comparisonofreconstructionresultsbetweensinglethreadandtwo-level parallel

2.3 重建加速效果對比

在二級并行架構重建加速效果的驗證中，首先同時使用單線程重建算法、二級并行架構重建算法以及單線程調用CUDA的并行算法重建了從100×100～5000×5000的6種不同分辨率的全息圖，重建截面數都為40。3種方法計算耗時以及加速比(speedup ratio)如圖6所示。speedup ratio 1為單線程耗時與單線程調用CUDA耗時的比值，speedup ratio 2為單線程耗時與二級并行耗時的比值。如圖6所示，重建低分辨率的全息圖時，二級并行架構的加速效果不明顯，主要是因為計算量小，單線程重建程序也有足夠的計算能力在短時間內處理完成，而并行計算能力沒有充分體現出來。當全息圖分辨率增大時，重建耗時明顯減少。對于100×100,200×200,500×500,1000×1000,2000×2000和5000×5000的全息圖，二級并行架構的加速比分別達到了22.8倍、30.7倍、36.1倍、42.7倍和48.3倍，充分顯示了其并行加速能力。全息圖的分辨率增大到1000×1000時，加速比的增加開始變得緩慢。這主要是由于一個重建截面的內存動態申請過程不能通過并行處理來加速，因此動態申請內存所消耗的時間占了很大的比例，加速比的增加也就逐漸平緩。單線程調用CUDA的重建耗時結果也證明了全息圖分辨率增大時動態申請內存的耗時會大大增加，分辨率大于1000×1000后，單線程調用CUDA的加速比迅速下降。在實際工程應用中，例如工業粉體粒度的在線監測，往往要求數字全息技術的結果反饋時間在分鐘級，假設全息圖分辨率為1000×1000，每張全息圖捕捉的顆粒數目為1000。60000顆顆粒能夠具有充分的代表性，那么至少需要處理60張全息圖。如果要在2min內反饋結果，則一張全息圖的處理時間要少于2s，而本文中采用的并行架構處理時間為1.41s，能夠滿足要求。

Fig.6 Time consumption of the three methods for holograms with different resolution

其次，選擇分辨率為1000×1000的全息圖，測試3種算法在重建截面數目不同時的計算耗時和加速比，結果如圖7所示。隨著重建截面的增多，單線程重建耗時急劇增加。二級并行算法的加速效果受限于線程數量，重建截面較少時，計算耗時主要來源于線程通信，重建截面數目影響較小，隨著重建截面數目的增加，線程趨于滿負荷狀態，此時二級并行算法的計算耗時隨之也會有較大的增加。而對于單線程調用CUDA，加速比都在20左右，重建截面數大于200后略有下降，這說明CUDA的并行效果穩定，單線程調用CUDA的重建耗時與單純單線程的重建耗時為近似線性的關系。以上不同分辨率和不同截面數下的實驗結果表明，采用二級并行架構重建全息圖，當分辨率過大或重建截面過多時，加速效果會受到內存申請和線程通信的限制，在實際應用過程中應加以考慮。

Fig.7 Time consumption of the three methods for reconstruction with different plane numbers

3 結 論

提出了一種基于二級并行架構的顆粒全息圖快速重建方法，該架構結合了OpenMP多線程技術和CUDA技術，利用OpenMP實現圖片級并行，利用CUDA實現像素級并行。以煤粉顆粒全息圖為實驗對象，驗證了二級并行架構全息重建的準確性和加速效果。結果表明二級并行架構在保證了重建準確性的同時，可極大提高重建速度。對于分辨率為5000×5000的全息圖，可實現48.3倍的加速比。但隨著全息圖分辨率增大，加速比的增加會趨于平緩；隨著重建截面的增多，加速比則趨于穩定。