GPU并行計(jì)算在FD-OCT成像中加速插值計(jì)算的應(yīng)用

易 俊，秦曉萌，岑穎珊，劉碧旺，韓定安，王茗祎，周月霞

（佛山科學(xué)技術(shù)學(xué)院 物理與光電工程學(xué)院，廣東 佛山 528200）

0 引言

圖1 FD-OCT系統(tǒng)組成Fig.1 FD-OCT system composition

頻域光學(xué)相干層析成像技術(shù)（Frequency Doman Optical Coherence Tomography, FD-OCT）是一種高軸向分辨率的新興無創(chuàng)三維光學(xué)成像技術(shù)[1]。自提出至今已經(jīng)得到了迅速發(fā)展，能實(shí)現(xiàn)對(duì)生物組織的功能成像，目前在醫(yī)學(xué)領(lǐng)域有著廣泛的應(yīng)用。近年來，隨著FD-OCT 技術(shù)的發(fā)展及超高速CMOS 線陣相機(jī)的更新，F(xiàn)D-OCT 系統(tǒng)的采集速度基本能滿足實(shí)時(shí)成像的要求。在FD-OCT 中，成像計(jì)算的關(guān)鍵步驟是插值計(jì)算。由于圖像重建過程中需要對(duì)每一線（即線陣相機(jī)的每一次曝光采集到的數(shù)據(jù)）數(shù)據(jù)都要進(jìn)行插值計(jì)算并進(jìn)行正向傅里葉變換。一幀完整的二維圖像至少由500 條線組成，數(shù)據(jù)吞吐量較大，從而加重了重建圖像中數(shù)據(jù)處理環(huán)節(jié)的計(jì)算負(fù)擔(dān)，在Matlab 語言下，基于CPU 算法的程序一幅由500 線組成的FD-OCT 圖像實(shí)現(xiàn)3 次樣品插值需要至少60ms 以上。其中插值計(jì)算占據(jù)了成像三之二的計(jì)算時(shí)間。因此，提高插值的計(jì)算速度，是提高FDOCT 成像速度的關(guān)鍵核心問題。

為了解決FD-OCT 系統(tǒng)在處理數(shù)據(jù)方面遇到的問題，相關(guān)研究人員提出了許多的解決方法：利用多核CPU 對(duì)FD-OCT 成像系統(tǒng)的數(shù)據(jù)進(jìn)行并行處理；也有相關(guān)研究人員在FD-OCT 系統(tǒng)中增加硬件模塊來加快數(shù)據(jù)處理速度，這樣使得系統(tǒng)的成本大大地增加了[2]。隨著GPU 技術(shù)的飛速發(fā)展，它不僅在傳統(tǒng)的圖像處理方面發(fā)揮作用，也早已被用于磁共振成像，系統(tǒng)用于數(shù)據(jù)處理的加速。Watanabe 和Itagaki 將GPU 應(yīng)用于FD-OCT 系統(tǒng)，利用其對(duì)FD-OCT 數(shù)據(jù)做了一維正向傅里葉變換，達(dá)到了每秒8 幀[3]。本文研究將經(jīng)濟(jì)實(shí)惠的商用GPU 應(yīng)用到FD-OCT 系統(tǒng)中，大幅度地加快了系統(tǒng)的數(shù)據(jù)處理速度。在每線像素點(diǎn)為1024，每幀1000 線的實(shí)驗(yàn)條件下，成像重建速度達(dá)到52 幀每秒，實(shí)現(xiàn)了FD-OCT 的實(shí)時(shí)成像。

1 系統(tǒng)組成和工作原理

1.1 系統(tǒng)的組成

本文用到的FD-OCT 系統(tǒng)組成如圖1 所示，主要包括超輻射發(fā)光激光二極管（λ0=1310nm，Δλ=400nm），樣品臂、參考臂、光譜儀和計(jì)算機(jī)（包含數(shù)據(jù)采集卡、圖像采集卡和GPU）5 個(gè)部分。其中數(shù)據(jù)采集卡為美國國家儀器公司的PCIe-6711，圖像采集卡采用的是來自美國國家儀器公司PCIe-1429；GPU 采用了英偉達(dá)公司的Ge Force GTX1060 顯卡（NVIDA CUDA 核心數(shù)1280 個(gè)，3GB 顯存）。

1.2 系統(tǒng)的工作原理

FD-OCT 系統(tǒng)的主要構(gòu)成是麥克爾遜干涉儀和光譜儀。寬帶光源的光被光纖耦合器一分為二，其中一束射向樣品臂，另一束射向參考臂；樣品臂上的后向散射光和參考臂的反射光在耦合器處相遇，產(chǎn)生相干光，這一部分即為麥克爾遜干涉儀。干涉光通過光譜儀分光，而后該干涉光被CCD 相機(jī)接收。根據(jù)單色光干涉理論可知，分光鏡處的光強(qiáng)可以表示為：

其中，IR和IS是參考光和樣品光的直流信號(hào)，AR為參考光的振幅，AS為樣品光的振幅，zj為等光程面的探測(cè)深度，Δφ 為相位差。當(dāng)干涉光經(jīng)光柵分光后，CCD 相機(jī)的每個(gè)線陣單元接收到的強(qiáng)度信號(hào)可以表示為：

由于光柵方程：

m 為衍射級(jí)次，λ 為衍射波長，φ 為入射角，θ 為衍射角以及d 為光柵常數(shù)。根據(jù)該光柵方程對(duì)入射光譜進(jìn)行分光[4]。干涉光通過光柵分為不同波長的光，之后衍射到線陣CCD 上。因此，相機(jī)記錄到的信號(hào)是一個(gè)以波長λ 為參變量的光強(qiáng)信號(hào)，根據(jù)干涉公式（2），實(shí)際需要獲得的信號(hào)是以波矢k 為參變量的信號(hào)。由于將波長信號(hào)插值到波矢信號(hào)是FD-OCT 成像的關(guān)鍵步驟，需要對(duì)相機(jī)采集到的光譜信號(hào)進(jìn)行插值，獲得真正的干涉信號(hào)。占據(jù)了成像的三分之二的計(jì)算時(shí)間，因而提高插值的計(jì)算速度是成像的關(guān)鍵核心問題。本文應(yīng)用了基于GPU 的CUDA 語言編寫了一套三次樣條插值算法，加速了整個(gè)插值計(jì)算過程。

1.3 三次樣條插值

三次樣條插值（Cubic Spline Interpolation），是一種應(yīng)用廣泛的插值函數(shù) 。三次樣條插值函數(shù)是由多段的低次多項(xiàng)式組成，并且它的曲線非常光滑；得到插值點(diǎn)處的函數(shù)值及提供（n-1）個(gè)邊界節(jié)點(diǎn)處的導(dǎo)數(shù)信息就的三次樣條插值函數(shù)[5]。

1.3.1 三次樣條插值函數(shù)的定義[6-8]

圖2 三次樣條插值求解流程圖Fig.2 Three-spline interpolation solution flowchart

假設(shè)有以下兩個(gè)節(jié)點(diǎn)x:a=x0＜x1＜...＜xn=b，y:y0y1...yn。樣條曲線是一個(gè)分段定義的函數(shù)。在n+1 個(gè)數(shù)據(jù)點(diǎn)中，共有n 個(gè)區(qū)間，三次樣條插值函數(shù)滿足以下條件：1）在每個(gè)區(qū)間[xi,xi+1]上，S(x)都是三次多項(xiàng)式；2）滿足S(xi)=yi；3）S(x)的導(dǎo)數(shù)S'(x)及二階導(dǎo)數(shù)S"(x)在區(qū)間[a,b]上都是連續(xù)的，即S(x)是光滑的曲線[9]。所以n 個(gè)三次多項(xiàng)式分段可以寫作：

其中ai，bi，ci，di代表曲線擬合的4 個(gè)系數(shù)。

1.3.2 三次樣條函數(shù)的求解

求解三次樣條插值函數(shù)是利用n+1 個(gè)數(shù)據(jù)節(jié)(x0,y0),(x1,y1),(x2,y2),…,(xn,yn)先計(jì)算步hi=xi+1-xi，令mi=Si"(xi)，由插值條件Si(xi+1)=yi+1以及Si(x)在樣點(diǎn)xi處二階導(dǎo)數(shù)的條件，增加自然邊界條件：

得到如下矩陣方程：

解矩陣方程（7），求出二次微分值mi[10]，從而求得系數(shù)ai，bi，ci，di并將系數(shù)代入獲得擬合曲線，獲得待求插值點(diǎn)的數(shù)值，最終完成三次樣條插值函數(shù)求解。

2 基于GPU并行計(jì)算的數(shù)據(jù)處理模塊的實(shí)現(xiàn)

2.1 CUDA：并行計(jì)算架構(gòu)

CUDA（Compute Unified Device Architecture）是一種并行計(jì)算架構(gòu)[11]。CUDA 將C 語言作為編程語言，并提供大量的計(jì)算指令，使得效率更高的密集數(shù)據(jù)計(jì)算解決方案能夠在GPU 突出計(jì)算能力的基礎(chǔ)上建立起來[12]。

CUDA 構(gòu)架分為兩部分：Host 和Device。通常而言，Host 指的是主機(jī)，Device 指的是GPU[13]。在CUDA 構(gòu)架中，CPU 作為主機(jī)，負(fù)責(zé)執(zhí)行程序的串行部分；GPU 作為協(xié)助處理器主要負(fù)責(zé)對(duì)密集且大量的數(shù)據(jù)進(jìn)行并行計(jì)算[14]。用CUDA 編寫好的程序叫做核（kernel）函數(shù)。CUDA 允許程序員定義稱為核的C 語言函數(shù)或使用CUDA 指令，在調(diào)用此類函數(shù)時(shí)，它將由N 個(gè)不同的CUDA 線程并行執(zhí)行N 次，執(zhí)行核的每個(gè)線程都會(huì)被分配一個(gè)獨(dú)特的線程ID，可通過內(nèi)置的threadIdx 變量在內(nèi)核中訪問此ID[15]。在 CUDA 程序中，主程序在調(diào)用GPU 內(nèi)核之前，務(wù)必對(duì)核執(zhí)行配置操作，以此確定線程塊數(shù)和每個(gè)線程塊中的線程數(shù)以及共享內(nèi)存大小。

2.2 數(shù)據(jù)處理流程

在FD-OCT 系統(tǒng)處理過程中，首先，CCD 相機(jī)采集到的原始數(shù)據(jù)經(jīng)過圖像采集卡傳輸?shù)紺PU 內(nèi)存上，然后再傳輸?shù)紾PU 顯存上面進(jìn)行數(shù)據(jù)處理；其中數(shù)據(jù)處理的主要過程包括對(duì)采集數(shù)據(jù)進(jìn)行去背景、色散補(bǔ)償、波長空間（λ空間）到波數(shù)空間（k 空間）的轉(zhuǎn)換、插值、FFT 變換、取模和取對(duì)數(shù)。最后，把處理好的數(shù)據(jù)再傳輸?shù)紺PU 內(nèi)存上并由計(jì)算機(jī)進(jìn)行顯示。由于插值過程計(jì)算復(fù)雜且系統(tǒng)采集到的圖像數(shù)據(jù)每列是相互獨(dú)立的，處理的時(shí)候也是分開一列一列數(shù)據(jù)進(jìn)行處理的，所以可以基于GPU 強(qiáng)大的并行計(jì)算能力對(duì)數(shù)據(jù)進(jìn)行插值計(jì)算，并把CPU 計(jì)算插值部分代碼改寫成在GPU 上執(zhí)行的kemel 函數(shù)，從而加快數(shù)據(jù)處理的速度，使成像時(shí)間減少。

2.3 程序設(shè)計(jì)

首先，定義整個(gè)數(shù)據(jù)處理過程中需要的變量，并給他們分配對(duì)應(yīng)大小的內(nèi)存空間；其中Host 端用malloc 函數(shù)開空間，Device 端用CUDA 內(nèi)部指令cudaMalloc 開空間；利用cudaMemcpy 將原始數(shù)據(jù)以及運(yùn)算所需中間變量從主機(jī)內(nèi)存復(fù)制到設(shè)備內(nèi)存。根據(jù)分配的線程空間借用對(duì)應(yīng)ID 進(jìn)行數(shù)據(jù)的調(diào)用和處理，其中包括調(diào)整數(shù)據(jù)類型以及去除背景噪聲。其次，由于CUDA 編程沒有自帶插值的庫函數(shù)，本文編寫了一套三次樣條插值算法并用CUDA 語言編寫成在GPU 上運(yùn)行的kernel 函數(shù)。

圖2 為三次樣條插值函數(shù)求解的流程圖，首先計(jì)算出步長，然后通過計(jì)算得到三對(duì)角矩陣，利用CUDA 自帶的庫函數(shù)求解三對(duì)角矩陣[16]；最后用求解矩陣得到的結(jié)果去計(jì)算多項(xiàng)式擬合系數(shù)，并將系數(shù)代入得到擬合曲線，獲得待求插值點(diǎn)的數(shù)值，完成三次樣條插值求解。

按照三次樣條插值的原理以及流程圖，以下為本文編

寫的基于GPU 三次樣條插值具體程序代碼：

//計(jì)算步長

__global__ void Kernel001(float *K, float *ks, int *index,float *h, float *dx, float *x31, float *xn)

{

h[idx] = K[idx] - ks[index[idx]]；

}

//得到三對(duì)角矩陣

__global__ void Kernel002(float *dx, float *x31, float *xn,float *A, float *B, float *C)

{

int idx = threadIdx.x + blockDim.x * blockIdx.x；

C[idx] = x31[0]； B[idx] = dx[1]； A[idx] = 0；

C[idx] = dx[idx - 1]； B[idx] = 2 * (dx[idx] + dx[idx - 1])；

A[idx] = dx[idx]；

C[idx] = 0； B[idx] = dx[idx - 2]；

A[idx] = xn[0]；

}

//三對(duì)角矩陣求解

cusparseCreate(&cusparseH)；

cusparseSgtsvInterleavedBatch(cusparseH,0, pixel, d_dl0,d_d0, d_du0, d_b, line*pic, d_work)；

以上兩條函數(shù)其專門用于求解三對(duì)角矩陣源于CUDA自帶的庫函數(shù)CUSPARSE LIBRARY。第一條用于創(chuàng)建句柄，第二條負(fù)責(zé)運(yùn)算求解矩陣。調(diào)用GPU 自帶函數(shù)求解矩陣，極大地提高了數(shù)據(jù)處理的速度以及精確度。

//用求解矩陣結(jié)果計(jì)算多項(xiàng)式擬合系數(shù)

_global__ void Kernel0061(float *divdif, float *s, float *dx,float *Y, float *dd, float *cc, float *bb, float *aa)

{

int idx = threadIdx.x + blockDim.x * blockIdx.x；

int idy = threadIdx.y + blockDim.y * blockIdx.y；

int idz = threadIdx.z + blockDim.z * blockIdx.z；

if (idy ＜ line && idx ＜ pixel - 1 && idz ＜ pic)

{ aa[idx + (pixel - 1)*idy + idz*s_frame] = Y[idx +pixel*idy + idz*frame]；

bb[idx + (pixel - 1)*idy + idz*s_frame] = s[idx + pixel*idy+ idz*frame]；

cc[idx + (pixel - 1)*idy + idz*s_frame] = 2 * (divdif[idx +(pixel - 1)*idy + idz*s_frame] - s[idx + pixel*idy + idz*frame]) /dx[idx] - (s[idx + 1 + pix-el*idy + idz*frame] - divdif[idx + (pixel- 1)*idy + idz*s_frame]) / dx[idx]；

dd[idx + (pixel - 1)*idy + idz*s_frame] = ((s[idx + 1 +pixel*idy + idz*frame] - divdif[idx + (pixel - 1)*idy + idz*s_frame]) / dx[idx] - (divdif[idx + (pixel - 1)*idy + idz*s_frame] -s[idx + pixel*idy + idz*frame]) / dx[idx]) / dx[idx]；

}

}

//系數(shù)代入得到擬合曲線，獲得待求插值點(diǎn)數(shù)值

__global__ void Kernel007(float *h, int *index, float *Y,float *dd, float *cc, float *bb, float *aa)

{ int idx = threadIdx.x + blockDim.x * block-Idx.x；

int idy = threadIdx.y + blockDim.y * blockIdx.y；

int idz = threadIdx.z + blockDim.z * blockIdx.z；

if (idx ＜ pixel && idy ＜ line && idz ＜ pic) {

Y[idx + pixel*idy + idz*frame] = h[idx] * (h[idx] * (h[idx] *dd[index[idx] + idy*(pixel - 1) + idz*s_frame] + cc[index[idx] +idy*(pixel - 1) + idz*s_frame]) + bb[index[idx] + idy*(pixel - 1)+ idz*s_frame]) + aa[index[idx] + idy*(pixel - 1) + idz*s_frame]；

}

}

然后，應(yīng)用CUDA 自帶的cuFFT 庫函數(shù)對(duì)插值后的數(shù)據(jù)做快速傅里葉變化，并改寫成kernel 函數(shù)；同時(shí)添加對(duì)快速傅里葉變換后數(shù)據(jù)進(jìn)行取模和取對(duì)數(shù)處理的核函數(shù)，實(shí)現(xiàn)FD-OCT 系統(tǒng)數(shù)據(jù)在GPU 完成所有數(shù)據(jù)處理環(huán)節(jié)。最后，將最終處理完的數(shù)據(jù)再次用cudaMemcpy 從顯存復(fù)制到內(nèi)存并在計(jì)算機(jī)上顯示和釋放所有顯存以及內(nèi)存空間。

3 程序驗(yàn)證以及實(shí)驗(yàn)結(jié)果分析

本系統(tǒng)采用Microsoft Visual Studio2015 中集成CUDA Toolkit 64 bit 9.2 和Nvidia Driver for Win-dows7 64bit 為開發(fā)環(huán)境，分別用基于CPU 模式的Matlab 語言和基于GPU 模式的CUDA 語言對(duì)采集到的人體表皮圖像數(shù)據(jù)進(jìn)行數(shù)據(jù)處理得到圖像。Matlab 語言是基于CPU 數(shù)據(jù)處理單元的串行數(shù)據(jù)處理方法，而CUDA 語言則用到了GPU 突出的并行計(jì)算能力來實(shí)現(xiàn)成像。

3.1 判斷依據(jù)

為了比較兩種模式成像速度的快慢，可以通過比較兩種模式下處理數(shù)據(jù)和顯示每幀圖總用時(shí)。因此，對(duì)于兩種模式，在設(shè)計(jì)程序時(shí)添加了幀率計(jì)算及顯示和記錄，以及通過CUDA 的計(jì)時(shí)函數(shù)和Matlab 的計(jì)時(shí)函數(shù)可得出數(shù)據(jù)處理每個(gè)環(huán)節(jié)所用時(shí)間以及總用時(shí)，這方便比較兩者的成像速度。

3.2 實(shí)驗(yàn)結(jié)果及分析

本文用B 掃描模式對(duì)人體表皮進(jìn)行了成像實(shí)驗(yàn)（采集到的每幀圖像數(shù)據(jù)大小為1000 線×1024 像素/線×12 字節(jié)/像素），實(shí)驗(yàn)成像效果如圖3 所示，圖(a)為基于CPU模式的Matlab 語言B 掃描成像截面圖，圖(b)為基于GPU模式的CUDA 語言B 掃描成像截面圖，圖(c)對(duì)應(yīng)于(a)圖中白色虛線的強(qiáng)度圖，圖(d)對(duì)應(yīng)于(b)圖中白色虛線的強(qiáng)度圖。由圖3 可知，基于CPU 模式的Matlab 語言和基于GPU 模式的CUDA 語言處理數(shù)據(jù)得到的圖像質(zhì)量無明顯差異，且經(jīng)過強(qiáng)度圖對(duì)比，數(shù)據(jù)相同。

圖3 人體表皮B掃描實(shí)驗(yàn)圖(1000Lines)Fig.3 Human epidermis B scan experiment (1000Lines)

圖4 不同大小圖片在基于CPU模式的Matlab語言和基于GPU模式的CUDA語言的成像時(shí)間對(duì)比Fig.4 Comparison of image times for pictures of different sizes in the CPU-based matlab language and the GPU-based CUDA language

除此之外，本文還改變了每幀圖像的線數(shù)，得到兩種模式的成像時(shí)間如圖4 所示。從圖中可以得出，在不同線數(shù)每幀圖的情況下，基于GPU 模式的CUDA 語言的成像時(shí)間都要比基于CPU 模式的Matlab 要短15 倍左右。特別是隨著數(shù)據(jù)量的變大，應(yīng)用GPU 進(jìn)行插值計(jì)算以及成像的速度更快，用時(shí)更短；而基于CPU 的Matlab 語言則隨著數(shù)據(jù)量的增加變得更慢，進(jìn)一步說明了把GPU 應(yīng)用到FD-OCT系統(tǒng)中的插值計(jì)算以及并行數(shù)據(jù)處理極大地提高了FDOCT 成像的速度。

4 結(jié)論

本文在沒有改變?nèi)魏斡布腇D-OCT 系統(tǒng)的基礎(chǔ)上，將基于GPU 的CUDA 語言應(yīng)用到系統(tǒng)的數(shù)據(jù)處理過程中，特別是針對(duì)三次樣條插值巨大且重復(fù)的數(shù)據(jù)計(jì)算，很好地利用了GPU 突出的并行計(jì)算能力，將插值計(jì)算以及需要并行計(jì)算的數(shù)據(jù)處理過程用CUDA 進(jìn)行改寫，使得整個(gè)的成像時(shí)間縮短了一個(gè)數(shù)量級(jí)。在國內(nèi)已有的相關(guān)研究中，劉巧艷，劉銳，朱珊珊等人也把GPU 應(yīng)用到OCT 成像系統(tǒng)的數(shù)據(jù)處理中，取得了不錯(cuò)的成績。本文著重是應(yīng)用GPU 對(duì)整個(gè)數(shù)據(jù)處理過程的插值計(jì)算部分進(jìn)行加速和優(yōu)化以達(dá)到整個(gè)數(shù)據(jù)處理過程加速，速度比單一的基于CPU 的Matlab語言處理速度加快了將近15 倍。所以說明GPU 加速模式很好地解決了傳統(tǒng)基于CPU 的Matlab 語言插值計(jì)算慢以及不能實(shí)時(shí)成像的問題，實(shí)現(xiàn)了FD-OCT 系統(tǒng)實(shí)時(shí)的2D 成像。但是，整個(gè)程序的圖像顯示環(huán)節(jié)沒有實(shí)現(xiàn)并行化處理，整個(gè)成像速度還有提升的空間，需要進(jìn)一步探索，更大程度地發(fā)揮GPU 的突出并行計(jì)算能力。本文研究為FD-OCT系統(tǒng)三維實(shí)時(shí)成像打下了堅(jiān)實(shí)的基礎(chǔ)。