GPU并行計算在FD-OCT成像中加速插值計算的應用

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

（佛山科學技術學院 物理與光電工程學院，廣東 佛山 528200）

0 引言

圖1 FD-OCT系統組成Fig.1 FD-OCT system composition

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

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

1 系統組成和工作原理

1.1 系統的組成

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

1.2 系統的工作原理

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

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

由于光柵方程：

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

1.3 三次樣條插值

三次樣條插值（Cubic Spline Interpolation），是一種應用廣泛的插值函數 。三次樣條插值函數是由多段的低次多項式組成，并且它的曲線非常光滑；得到插值點處的函數值及提供（n-1）個邊界節點處的導數信息就的三次樣條插值函數[5]。

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

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

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

其中ai，bi，ci，di代表曲線擬合的4 個系數。

1.3.2 三次樣條函數的求解

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

得到如下矩陣方程：

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

2 基于GPU并行計算的數據處理模塊的實現

2.1 CUDA：并行計算架構

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

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

2.2 數據處理流程

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

2.3 程序設計

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

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

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

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

//計算步長

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

{

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

}

//得到三對角矩陣

__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]；

}

//三對角矩陣求解

cusparseCreate(&cusparseH)；

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

以上兩條函數其專門用于求解三對角矩陣源于CUDA自帶的庫函數CUSPARSE LIBRARY。第一條用于創建句柄，第二條負責運算求解矩陣。調用GPU 自帶函數求解矩陣，極大地提高了數據處理的速度以及精確度。

//用求解矩陣結果計算多項式擬合系數

_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]；

}

}

//系數代入得到擬合曲線，獲得待求插值點數值

__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]；

}

}

然后，應用CUDA 自帶的cuFFT 庫函數對插值后的數據做快速傅里葉變化，并改寫成kernel 函數；同時添加對快速傅里葉變換后數據進行取模和取對數處理的核函數，實現FD-OCT 系統數據在GPU 完成所有數據處理環節。最后，將最終處理完的數據再次用cudaMemcpy 從顯存復制到內存并在計算機上顯示和釋放所有顯存以及內存空間。

3 程序驗證以及實驗結果分析

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

3.1 判斷依據

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

3.2 實驗結果及分析

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

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

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

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

4 結論

本文在沒有改變任何硬件的FD-OCT 系統的基礎上，將基于GPU 的CUDA 語言應用到系統的數據處理過程中，特別是針對三次樣條插值巨大且重復的數據計算，很好地利用了GPU 突出的并行計算能力，將插值計算以及需要并行計算的數據處理過程用CUDA 進行改寫，使得整個的成像時間縮短了一個數量級。在國內已有的相關研究中，劉巧艷，劉銳，朱珊珊等人也把GPU 應用到OCT 成像系統的數據處理中，取得了不錯的成績。本文著重是應用GPU 對整個數據處理過程的插值計算部分進行加速和優化以達到整個數據處理過程加速，速度比單一的基于CPU 的Matlab語言處理速度加快了將近15 倍。所以說明GPU 加速模式很好地解決了傳統基于CPU 的Matlab 語言插值計算慢以及不能實時成像的問題，實現了FD-OCT 系統實時的2D 成像。但是，整個程序的圖像顯示環節沒有實現并行化處理，整個成像速度還有提升的空間，需要進一步探索，更大程度地發揮GPU 的突出并行計算能力。本文研究為FD-OCT系統三維實時成像打下了堅實的基礎。