肖 波,鄭華東,2，劉柯健，李 飛，高智方

(1.上海大學(xué) 精密機(jī)械工程系,上海 200444；2.教育部先進(jìn)顯示與系統(tǒng)應(yīng)用集成重點(diǎn)實(shí)驗(yàn)室，上海 200444)

引言

計(jì)算機(jī)生成全息技術(shù)(CGH)是一種利用具有光學(xué)全息數(shù)學(xué)模型通過計(jì)算機(jī)輔助生成數(shù)字全息圖的技術(shù)。盡管CGH可通過計(jì)算機(jī)技術(shù)快速實(shí)現(xiàn)，但仍需要大量耗時(shí)的計(jì)算復(fù)核，這對計(jì)算機(jī)生成全息圖的速率產(chǎn)生了很大的瓶頸。為了解決這些問題，已經(jīng)提出了各種方法。賈甲提出通過使用壓縮表法降低計(jì)算復(fù)雜度和減少內(nèi)存存儲(chǔ)空間以提高速度[1]，Yoshikawa提出使用泰勒展開近似方法降低計(jì)算復(fù)雜度，以提高計(jì)算速度[2]。H. Shimobaba利用現(xiàn)場可編程門陣列(FPGA)的硬件和圖形處理單元(GPU)硬件進(jìn)行加速處理[3-4]，蔣曉瑜教授等人利用查表法對電源法數(shù)據(jù)壓縮加速[5]。2018年，國防科技大學(xué)辛誠提出利用CUDA(compute unified device architecture)實(shí)現(xiàn)基于GPU環(huán)境的matlab計(jì)算全息圖的加速計(jì)算像息圖快速計(jì)算[6]。寧苒提出采取GPU的matlab計(jì)算全息圖[7]。

雖然Shimobaba提出的硬件加速算法可以大幅提高運(yùn)算速度，但是該方法具有以下限制：開發(fā)FPGA板的成本很高，開發(fā)時(shí)間長且需要專業(yè)的FPGA的處理知識(shí)。而通過matlab調(diào)用庫函數(shù)計(jì)算無法充分發(fā)揮GPU底層計(jì)算資源，還有進(jìn)一步的優(yōu)化空間。同時(shí)，由于已經(jīng)有很多方法對計(jì)算全息的點(diǎn)源法和三角面片法進(jìn)行了加速處理，本文使用層析法計(jì)算三維物體全息圖的方式以CUDA C進(jìn)行程序編寫，通過對線程和存儲(chǔ)架構(gòu)優(yōu)化，并通過流處理方式掩蓋數(shù)據(jù)傳輸?shù)难訒r(shí)時(shí)間，對菲涅爾衍射算法進(jìn)行最大幅度的優(yōu)化以獲得最快的計(jì)算速率，對計(jì)算全息進(jìn)行了數(shù)值模擬和光學(xué)再現(xiàn)以獲得三維物體計(jì)算全息圖的再現(xiàn)像。對GPU加速計(jì)算分析表明，使用GPU并行計(jì)算比傳統(tǒng)CPU串行運(yùn)算速度有大幅提高，并且隨著處理的切片數(shù)據(jù)量的增加優(yōu)勢更加明顯。

1 計(jì)算全息計(jì)算原理

1.1 全息的基本原理

計(jì)算全息術(shù)是一種利用光的干涉和衍射特性顯示的技術(shù)。光全息顯示可以用圖1來說明。利用從物體反射或者透射的攜帶了物體信息的物光波傳播到記錄平面(xh,yh)，其在記錄平面上與參考光波干涉疊加的復(fù)振幅分布稱為物光O(xh,yh)，照射到記錄平面的另一束光稱為參考光R(xh,yh)，全息圖就是物光和參考光干涉的條紋分布；全息再現(xiàn)過程也利用了光的衍射原理，當(dāng)光波照射到全息記錄介質(zhì)時(shí)，將再現(xiàn)和還原之前保存的原物的記錄信息。

圖1 光全息記錄光路原理圖

1.2 基于層析法的計(jì)算全息原理

首先通過計(jì)算機(jī)軟件模擬三維物體數(shù)據(jù)模型可以得到光波的數(shù)學(xué)描述。根據(jù)處理過程的不同，分為點(diǎn)源法[8]、三角面片法[9]和層析法[10]。然后通過對物體或者波面進(jìn)行離散抽樣出樣點(diǎn)計(jì)算出物光波在全息平面上的復(fù)振幅分布。最后可以把全息平面上的光波復(fù)振幅分布編碼成全息圖透過率變化，在空間光調(diào)制器上進(jìn)行再現(xiàn)[11-13]。

本文使用的是層析法進(jìn)行全息圖計(jì)算，該方法認(rèn)為一個(gè)三維物體可以離散為一系列垂直于Z軸的多層平行的截面，如圖2所示。每一個(gè)平面圖形是一個(gè)二維數(shù)組，可以方便地利用快速二維傅里葉變換計(jì)算其衍射光波的分布，然后計(jì)算每個(gè)平面的衍射光波，最后將所有平面的物光波疊加并抽取其相位信息，最終平面上疊加的相位信息就是我們需要得到的全息圖數(shù)據(jù)結(jié)果。

圖2 層析法計(jì)算三維物體全息圖的原理圖

將物體沿法線方向分為N層，第n層平面圖與全息圖平面距離為zn，振幅透射率為fn(x,y)，則其在全息記錄平面上的衍射分布菲涅爾近似為[14]

(1)

其中：k為波數(shù)，k=2π/λ；xh、yh為全息圖平面坐標(biāo)，On(xh,yh)表示物面第n層的復(fù)振幅。

三維物體各個(gè)層面衍射到全息面的復(fù)振幅之和為

(2)

考慮到人眼結(jié)構(gòu)的分辨率[15]，為了得到連續(xù)的再生像，相鄰之間的距離Δz必須小于人眼的縱向分辨距離。人眼的最小分辨距離公式為

(3)

人眼的分辨角ε=2.9×10-6z2rad，瞳距為Le=65 mm，則可以計(jì)算出相應(yīng)三維物體的最小分割層數(shù)。

2 CUDA并行計(jì)算

2.1 GPU計(jì)算原理

CUDA(compute unified device architecture)是NVIDIA推出的一個(gè)基于GPU通用計(jì)算的并行計(jì)算處理平臺(tái)和編程模型[16]。利用圖形處理器(GPU)多處理單元的編程特性，通過大量的線程并行實(shí)現(xiàn)計(jì)算性能的顯著提高。一個(gè)GPU具有多個(gè)流多處理器簇(streaming multiprocessors, SM)，如圖3，每一個(gè)SM都分別設(shè)置有獨(dú)立訪問紋理內(nèi)存(texture memory)、常量內(nèi)存(constant memory)和全局內(nèi)存(global memory)及它們的總線，并通過里面的多個(gè)流處理器(streaming processors, SP)共享控制邏輯和指令緩存[17]。全局內(nèi)存具有千兆字節(jié)(GB)的數(shù)據(jù)傳輸帶寬，但是具有較長的讀取時(shí)間延時(shí)，SM能夠在4個(gè)時(shí)鐘周期內(nèi)訪問共享緩存，訪問全局緩存則需要400～600個(gè)時(shí)鐘周期，為了提高處理和計(jì)算速度和處理瓶頸，通常需要將GPU和CPU進(jìn)行流并行處理異步計(jì)算[18]。

2.2 基于層析法快速計(jì)算

為了實(shí)際使用全息圖，快速計(jì)算CGH非常重要，通常CGH需要消耗大量的計(jì)算時(shí)間[19]。對于1 024×1 024分辨率的10 000個(gè)點(diǎn)光源的3D物體，全息圖需要大約十億次計(jì)算，為了減少計(jì)算時(shí)間，我們使用了GPU并行處理。圖4顯示了中央處理單元(CPU)和GPU的結(jié)構(gòu)差異。

圖3 SM內(nèi)部組成結(jié)構(gòu)圖

圖4 CPU架構(gòu)和GPU架構(gòu)示意圖

如上所述，GPU的并行計(jì)算能力遠(yuǎn)超過CPU，因?yàn)镚PU擁有比CPU多幾個(gè)數(shù)量級的運(yùn)算核心。因此許多研究已經(jīng)使用了GPU快速計(jì)算復(fù)雜算法，在本文中，我們使用了CUDA是因?yàn)镃UDA比OpenCL更容易實(shí)現(xiàn)CGH運(yùn)算算法。如圖5所示，先通過對三維物體計(jì)算分層，提取出每一層次的截面信息，然后分配線程塊和內(nèi)存，對各個(gè)截面層次的圖像信息進(jìn)行迭代并行計(jì)算，并同步進(jìn)行各個(gè)分層的全息圖的合成，最終得出三維物體的衍射全息信息。

圖5 GPU并行處理步驟

在一定條件下，可知菲涅爾衍射近似公式(4)可以看作為一個(gè)二維傅里葉變換的過程：

(4)

其中：F{·}表示傅里葉變換；x0、y0為衍射面的坐標(biāo)；x、y為觀察面的坐標(biāo)。如果對公式做變形推導(dǎo)，可以得出只要求一次的傅里葉變換的離散求解，稱為S-FFT算法(single fast Fourier transform algorithm)，如(5)式所示。其中：Δx、Δy是衍射面x方向和y方向上的抽樣間隔；Δx、Δy是觀察面x方向和y方向上的抽樣間隔；m0、n0、m、n分別為衍射面和抽樣面的抽樣點(diǎn)數(shù)，且m0=1,2,…,M0，n0=1,2,…,N0。

(5)

通過利用GPU流多處理器優(yōu)化算法執(zhí)行和數(shù)據(jù)傳輸過程(其示意圖如圖6所示)，使得全息圖運(yùn)算時(shí)間和物點(diǎn)信息的傳輸時(shí)間達(dá)到穩(wěn)定的平衡點(diǎn)。同時(shí)還可以使用CUDA線程優(yōu)化的方法來優(yōu)化硬件資源的分配，實(shí)現(xiàn)的關(guān)鍵代碼如圖7所示，通過CUDA流處理掩蓋傳輸延遲，提高運(yùn)算速率[20-21]，還可以優(yōu)化算法的執(zhí)行順序進(jìn)行異構(gòu)計(jì)算等方式進(jìn)一步提高其計(jì)算速度。由于GPU計(jì)算瓶頸主要有帶寬、內(nèi)存和浮點(diǎn)數(shù)運(yùn)算能力，而全息圖計(jì)算時(shí)，主要性能瓶頸存在于帶寬。因此將菲涅爾變換算法分為2個(gè)核函數(shù)分開進(jìn)行流處理，以覆蓋存儲(chǔ)傳輸時(shí)的延遲。將菲涅爾算法分為傳輸運(yùn)算二維傅里葉變換和結(jié)果累加計(jì)算兩個(gè)部分。

圖6 GPU流處理計(jì)算示意圖

GPU并行處理層析平面數(shù)據(jù)示意圖如圖8所示，根據(jù)二維圖像像素點(diǎn)劃分為一個(gè)個(gè)線程矩陣進(jìn)行計(jì)算。

圖7 GPU端部分加速代碼

3 實(shí)驗(yàn)結(jié)果與討論

本次實(shí)驗(yàn)采用的硬件平臺(tái)為Intel Core i5-3470，運(yùn)行內(nèi)存為8 G，實(shí)驗(yàn)顯卡為筆記本顯卡，參數(shù)為NVIDIA GeForce GT1050Ti，顯卡的計(jì)算能力為6.1，顯存4 G。如表1所示。

實(shí)驗(yàn)顯卡為筆記本獨(dú)立顯卡有768個(gè)核心，最大計(jì)算頻率為1.62 GHz，內(nèi)存處理頻率為3 504 MHz。

圖8 基于層析法的CUDA并行化示意圖

表1 硬件環(huán)境

采用層析法計(jì)算三維物體全息圖的算法結(jié)果如圖9和圖10所示，可以知道硬件并行加速計(jì)算三維模型的數(shù)據(jù)的正確性，同時(shí)通過光電全息平臺(tái)進(jìn)一步驗(yàn)證了并行計(jì)算程序的可行性。在此基礎(chǔ)上，分別通過GPU和CPU運(yùn)算采集其計(jì)算速度，結(jié)果見表2。

圖9 300 mm距離下512×512分辨尺寸再現(xiàn)效果

圖10 層析法計(jì)算三維物體計(jì)算全息圖結(jié)果數(shù)值仿真及光電再現(xiàn)驗(yàn)證

表2顯示了不同分辨率尺寸下的三維模型在一定分層數(shù)量時(shí)所花費(fèi)的時(shí)間和加速比。可以看出，經(jīng)過優(yōu)化后，隨著圖片切片的增加，在GPU端流處理并行運(yùn)行節(jié)省的時(shí)間也大幅增加，增加了1個(gè)數(shù)量級。

表2 CPU和GPU計(jì)算方式的計(jì)算時(shí)間和計(jì)算效率

4 結(jié)論

本文通過對層析法進(jìn)行GPU的加速研究，在GPU并行計(jì)算方法下，對三維模型分層切片同時(shí)優(yōu)化硬件資源和算法，再使用流處理的傳輸策略有效加快了計(jì)算全息圖生成速度。分析了不同分辨率和不同分層情況下的加速計(jì)算效率，基于GPU的全息圖并行加速計(jì)算方法比基于CPU串行計(jì)算方法的計(jì)算效率提高66倍左右。該方法為三維顯示計(jì)算全息圖的快速計(jì)算提供了有益的參考。