【作 者】劉巧艷，李躍杰，徐秋晶，趙金城，王立偉，高用賀

1 北京協(xié)和醫(yī)學(xué)院，北京市，100730

2 中國(guó)醫(yī)學(xué)科學(xué)院生物醫(yī)學(xué)工程研究所，天津市，300192

3 天津市眼科醫(yī)學(xué)設(shè)備技術(shù)工程中心，天津市，300384

0 引言

光學(xué)相干層析成像技術(shù)（Optical Coherence Tomography），自上世紀(jì)90年代被成功用于眼科疾病診斷[1]后，得到了迅速發(fā)展：由時(shí)域OCT階段發(fā)展到頻域OCT階段[2]；由組織結(jié)構(gòu)成像向組織功能成像發(fā)展[3-4]；由眼科診斷拓展到心血管、皮膚、口腔、組織工程等領(lǐng)域中的應(yīng)用[5-6]。

隨著超高速CMOS線陣掃描相機(jī)的發(fā)展，頻域OCT光譜譜線轉(zhuǎn)換及線采樣率已經(jīng)可以達(dá)到300 k線/s[7]，為臨床OCT系統(tǒng)實(shí)時(shí)成像提供了前提。目前影響眼科OCT系統(tǒng)實(shí)時(shí)成像的技術(shù)瓶頸是需要先將采樣數(shù)據(jù)進(jìn)行頻譜域空間(λ空間)到波數(shù)空間(K空間)變換，進(jìn)行插值變換和FFT變換，然后再將變換后的數(shù)據(jù)進(jìn)行2D或3D成像。由于成像的數(shù)據(jù)量很大，特別是進(jìn)行C模式掃描成像(如眼底視網(wǎng)膜en-face成像模式)時(shí)，需要先將獲得的3D圖像數(shù)據(jù)進(jìn)行處理后，再將得到的數(shù)據(jù)成像。因此，如何提高數(shù)據(jù)處理的速度，成為眼科OCT系統(tǒng)實(shí)現(xiàn)實(shí)時(shí)成像的關(guān)鍵。

為了克服眼科OCT系統(tǒng)實(shí)時(shí)成像的技術(shù)瓶頸，近年來(lái)研究人員進(jìn)行了大量研究，并提出了一些相關(guān)的解決方法。有研究人員將多核CPU引入到OCT系統(tǒng)進(jìn)行采樣數(shù)據(jù)的并行處理，對(duì)于1024個(gè)采樣點(diǎn)模式OCT系統(tǒng)，非均勻K空間數(shù)據(jù)處理的速度可達(dá)到80 k線/s[8]，均勻K空間數(shù)據(jù)的處理速度可達(dá)到207 k線/s[9]。也有研究人員通過(guò)在原有系統(tǒng)中增加DSPs[10]或者FPGAs[11]硬件模塊，加快實(shí)現(xiàn)數(shù)據(jù)處理過(guò)程。隨著圖形處理單元（GPU）技術(shù)的發(fā)展，其運(yùn)算能力和可編程性得到大幅度提高，它除了在傳統(tǒng)的圖像處理領(lǐng)域應(yīng)用繼續(xù)保持優(yōu)勢(shì)外，作為通用的并行計(jì)算處理器已被廣泛用于科學(xué)計(jì)算、工程、金融以及其他領(lǐng)域中。目前，有研究人員將GPU引入到OCT系統(tǒng)中，利用它強(qiáng)大的并行計(jì)算能力來(lái)加速采樣數(shù)據(jù)處理過(guò)程，解決OCT系統(tǒng)實(shí)時(shí)成像的技術(shù)瓶頸。Kang Zhang和Jin U.Kang將GPU引入到OCT系統(tǒng)中，在專業(yè)圖像工作站平臺(tái)上采用線性插值算法，在1024個(gè)采樣點(diǎn)模式下數(shù)據(jù)處理速度達(dá)到680 k線/s，2048個(gè)采樣點(diǎn)模式下數(shù)據(jù)處理速度可達(dá)到320 k線/s，達(dá)到了3D實(shí)時(shí)成像對(duì)于處理速度的要求[12]?，F(xiàn)階段將GPU應(yīng)用于OCT系統(tǒng)還處于實(shí)驗(yàn)室研究階段，為獲得更高處理速度，往往需要配置高性能的專用圖像工作站（配置多核處理器）和GPU，硬件成本較高。

本課題在不改變實(shí)驗(yàn)室現(xiàn)有硬件平臺(tái)的條件下，把低成本GPU引入到我們正在開發(fā)應(yīng)用的眼科OCT儀器中，實(shí)現(xiàn)儀器性能的大幅度提高，解決眼科OCT系統(tǒng)實(shí)時(shí)成像的問(wèn)題。

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

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

1.1 系統(tǒng)組成

本系統(tǒng)組成如圖1所示，主要包括眼科OCT信號(hào)采集系統(tǒng)和信號(hào)處理系統(tǒng)兩大部分。其中眼科OCT信號(hào)采集系統(tǒng)采用寬帶超亮發(fā)光二極管（SLD）(λ0=840 nm，△λ= 50 nm)作為系統(tǒng)光源，以2048像素CMOS線陣相機(jī)（采樣速率為70 k線/s）作為光譜儀的檢測(cè)器。信號(hào)處理系統(tǒng)由計(jì)算機(jī)和GPU組成。計(jì)算機(jī)CPU為 Intel Celeron Dual-Core E3400 @ 2.60 GHz，2 G內(nèi)存；GPU采用 NVIDIA 公司的GeForce GTX 460顯卡（336個(gè)流處理器，1 GB顯存）。

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

光源發(fā)出的光經(jīng)過(guò)50:50的光纖分束器后，被均勻分成兩束光，分別進(jìn)入OCT系統(tǒng)的參考臂和樣品臂。從樣品臂反射回來(lái)的信號(hào)光和從參考臂返回的參考光再次經(jīng)過(guò)光纖分束器匯合后發(fā)生干涉。包含樣品不同深度信息的干涉信號(hào)光譜，經(jīng)光譜儀的CMOS線陣掃描相機(jī)采集，并通過(guò)相機(jī)數(shù)據(jù)線傳輸?shù)接?jì)算機(jī)，由計(jì)算機(jī)內(nèi)的圖像采集卡對(duì)干涉信號(hào)光譜進(jìn)行A/D轉(zhuǎn)換，并將轉(zhuǎn)換結(jié)果作為采樣數(shù)據(jù)存入到計(jì)算機(jī)內(nèi)存中。將采樣數(shù)據(jù)通過(guò)PCIE×16總線傳輸?shù)紾PU顯存，借助GPU強(qiáng)大的并行數(shù)據(jù)處理能力進(jìn)行數(shù)據(jù)處理，并將處理好的結(jié)果數(shù)據(jù)送回計(jì)算機(jī)進(jìn)行圖像顯示。顯示的圖像包含了檢測(cè)樣品不同深度的結(jié)構(gòu)信息。

2 基于GPU加速技術(shù)的數(shù)據(jù)處理流程

2.1 CUDA架構(gòu)

CUDA(Compute Unified Device Architecture)是一種由NVIDIA公司推出的通用并行計(jì)算架構(gòu)，該架構(gòu)使GPU能夠解決復(fù)雜的計(jì)算問(wèn)題。在CUDA架構(gòu)下，開發(fā)人員可以通過(guò)CUDA C語(yǔ)言（CUDA C語(yǔ)言是對(duì)標(biāo)準(zhǔn)C語(yǔ)言的一種簡(jiǎn)單擴(kuò)展）對(duì)GPU編程[13]。

CUDA架構(gòu)中，CPU作為主機(jī)（Host），GPU作為協(xié)處理器或者設(shè)備（Device）。在一個(gè)系統(tǒng)中可以存在一個(gè)主機(jī)和多個(gè)設(shè)備。 CPU 主要負(fù)責(zé)進(jìn)行邏輯性強(qiáng)的事物處理和串行計(jì)算，GPU 則專注于執(zhí)行高度線程化的并行處理任務(wù)。CPU 、GPU 各自擁有相互獨(dú)立的存儲(chǔ)器地址空間：主機(jī)端的內(nèi)存和設(shè)備端的顯存；在CUDA程序中，將運(yùn)行在GPU上一個(gè)可以被并行執(zhí)行的步驟稱為kernel（內(nèi)核函數(shù)）。

圖2 CPU-GPU系統(tǒng)數(shù)據(jù)處理流程圖Fig.2 Signal processing flow chart of CPU-GPU hybrid system architecture

2.2 數(shù)據(jù)處理過(guò)程

本系統(tǒng)的數(shù)據(jù)處理過(guò)程如圖2所示，其中實(shí)箭頭所指方向代表數(shù)據(jù)在不同設(shè)備間的流動(dòng)方向，空箭頭所指方向代表數(shù)據(jù)在GPU內(nèi)部的流動(dòng)方向。系統(tǒng)數(shù)據(jù)處理過(guò)程包括對(duì)采樣數(shù)據(jù)進(jìn)行預(yù)處理、頻譜域空間（λ空間）到波數(shù)空間（K空間）變換、插值變換、FFT變換和POST-FFT變換。

在頻域OCT系統(tǒng)中，采樣數(shù)據(jù)是通過(guò)對(duì)OCT的光路系統(tǒng)掃描由相機(jī)采集到的，掃描一次得到一列數(shù)據(jù)（一個(gè)A-SCAN）。處理時(shí)是一列一列數(shù)據(jù)進(jìn)行處理的。針對(duì)每列數(shù)據(jù)彼此相互獨(dú)立、可以并行處理的特點(diǎn)，利用CUDA架構(gòu)將OCT系統(tǒng)整個(gè)數(shù)據(jù)處理過(guò)程改寫成適合在GPU上執(zhí)行的kernel函數(shù)，大大提高了數(shù)據(jù)處理速度，從而達(dá)到系統(tǒng)實(shí)時(shí)成像的要求。

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

首先確定系統(tǒng)數(shù)據(jù)處理過(guò)程中的串行部分和并行部分，選擇合適的算法，并按照算法確定數(shù)據(jù)和任務(wù)的劃分方式，將每個(gè)需要并行實(shí)現(xiàn)的步驟映射為一個(gè)滿足CUDA兩層并行模型的內(nèi)核函數(shù)。其中頻譜域空間（λ空間）到波數(shù)空間（K空間）變換模塊在主程序中執(zhí)行；數(shù)據(jù)預(yù)處理、插值運(yùn)算、POST-FFT運(yùn)算需要重新改寫成適合在GPU上執(zhí)行的kernel函數(shù)；FFT運(yùn)算部分采用CUDA自帶的CUFFT庫(kù)函數(shù)來(lái)進(jìn)行計(jì)算，CUDA4.1版本目前已經(jīng)支持FFT雙精度運(yùn)算。整個(gè)程序設(shè)計(jì)如下：

(1)給采樣數(shù)據(jù)和計(jì)算過(guò)程中的各個(gè)變量分配空間，包括內(nèi)存空間和顯存空間，在Host端CPU程序中，準(zhǔn)備好計(jì)算要用到的數(shù)據(jù)，將數(shù)據(jù)從主機(jī)內(nèi)存?zhèn)鬏數(shù)斤@存中。

(2)數(shù)據(jù)預(yù)處理包括數(shù)據(jù)類型轉(zhuǎn)換和去噪運(yùn)算。去噪運(yùn)算是將每一列數(shù)據(jù)組成的數(shù)組（一個(gè)A-SCAN）都減去一列同大小的噪聲數(shù)組。

(3)插值運(yùn)算實(shí)現(xiàn)了在進(jìn)行FFT變換之前，采樣數(shù)據(jù)K空間分布的均勻化。系統(tǒng)最常用的插值算法包括最鄰近插值算法、線性插值算法和三次樣條插值算法。

(4)利用CUDA自帶的CUFFT庫(kù)函數(shù)對(duì)插值運(yùn)算結(jié)果進(jìn)行FFT變換。

(5)POST-FFT運(yùn)算包括對(duì)FFT運(yùn)算結(jié)果取模、取對(duì)數(shù)并進(jìn)行歸一化。POST-FFT運(yùn)算結(jié)果存儲(chǔ)在體積數(shù)組中。

(6)根據(jù)不同成像平面的需要，抽取體積數(shù)組數(shù)據(jù)作為GPU結(jié)果數(shù)據(jù)。

(7)將結(jié)果數(shù)據(jù)從顯存?zhèn)鬏數(shù)街鳈C(jī)內(nèi)存中。

2.4 優(yōu)化分析

在用CUDA對(duì)GPU上的kernel函數(shù)進(jìn)行改寫時(shí)，采用了一些優(yōu)化方法：

(1)為了使內(nèi)存和顯存之間的數(shù)據(jù)傳輸速度更快，在分配主機(jī)端內(nèi)存的時(shí)候采用了pinned memory。

(2)使用流運(yùn)算來(lái)隱藏主機(jī)端和設(shè)備端數(shù)據(jù)通信的時(shí)間。

(3)合理利用讀取速度更快的shared memory來(lái)存放計(jì)算過(guò)程中的中間變量。

(4)將已經(jīng)計(jì)算好的在插值運(yùn)算過(guò)程中值保持不變的采樣數(shù)據(jù)的等間隔和非等間隔橫坐標(biāo)值存放在常數(shù)寄存器里面。

(5)根據(jù)系統(tǒng)資源進(jìn)行g(shù)rid和block維度設(shè)計(jì)。

(6)使用CUDA profiler對(duì)CUDA程序進(jìn)行性能測(cè)試，對(duì)耗時(shí)較長(zhǎng)的模塊進(jìn)行優(yōu)化。

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

本系統(tǒng)采用Microsoft Visual Studio2010中集成CUDA Toolkit 32 bit 4.1、CUDA SDK 32 bit 4.1和Nvidia Driver for Windows732 bit為開發(fā)環(huán)境，對(duì)采樣數(shù)據(jù)進(jìn)行B掃描模式和C掃描模式成像。

3.1 B掃描模式成像

B掃描模式成像圖像能提供視網(wǎng)膜斷層結(jié)構(gòu)，能清晰地顯示視網(wǎng)膜各層細(xì)微結(jié)構(gòu)及病理改變，并作出定性或定量分析，目前已成為視網(wǎng)膜疾病和青光眼強(qiáng)有力的診斷工具[14]。

我們采用100幀共計(jì)195 MB數(shù)據(jù)（每幀數(shù)據(jù)大小為500線×2048 像素/線×2 字節(jié)/像素）進(jìn)行B掃描模式成像。分別采用線性插值算法和三次樣條插值算法，利用CUDA提供的計(jì)時(shí)函數(shù)分別對(duì)CPU模式和CPU-GPU模式下系統(tǒng)單幀B掃描模式圖像成像時(shí)間進(jìn)行計(jì)時(shí)(計(jì)算100幀圖像成像時(shí)間取平均)，實(shí)驗(yàn)結(jié)果如表1所示。 從表1可知，采用GPU+CPU模式執(zhí)行成像數(shù)據(jù)處理的速度較CPU模式執(zhí)行同樣數(shù)據(jù)處理的速度提高超過(guò)一個(gè)數(shù)量級(jí)，其中采用線性插值算法速度提高了60倍，采用三次樣條插值算法速度提高了35倍。實(shí)驗(yàn)成像效果圖如圖3所示，其中(a)為采用線性插值算法在CPU模式下系統(tǒng)成像圖像；(b)為采用線性插值算法在CPU+GPU模式下系統(tǒng)成像圖像；(c)為采用三次樣條插值算法在CPU模式下系統(tǒng)成像圖像；(d)為采用三次樣條插值算法在CPU+GPU模式下系統(tǒng)成像圖像。由圖3可知，采用相同的插值算法，在CPU模式和CPU-GPU模式下系統(tǒng)成像圖像質(zhì)量無(wú)差異，采用三次樣條插值算法成像圖像質(zhì)量比采用線性插值算法成像圖像質(zhì)量要好。

圖3 視網(wǎng)膜B掃描成像圖像Fig.3 B-scan imaging images of the retina

插值運(yùn)算是系統(tǒng)數(shù)據(jù)處理過(guò)程中計(jì)算量最大、耗時(shí)最長(zhǎng)的一個(gè)環(huán)節(jié)。插值算法的選取不僅影響系統(tǒng)加速比，也影響系統(tǒng)成像速度和成像圖像質(zhì)量。采用線性插值算法，成像速度快，成像圖像質(zhì)量差；采用三次樣條算法，成像速度慢、成像圖像質(zhì)量好。

表1 CPU和GPU 成像速度對(duì)比Tab.1 Comparison of CPU and GPU imaging speed

3.2 C掃描模式成像

C掃描模式成像圖像能直觀地顯示眼底視網(wǎng)膜血管和黃斑等眼底組織的結(jié)構(gòu)信息，臨床上可用于診斷眼底病變，如黃斑病變[15]、眼底神經(jīng)組織結(jié)構(gòu)變化等，可用于對(duì)視網(wǎng)膜下新生血管的深度和層次進(jìn)行準(zhǔn)確定位[16]。

進(jìn)行C掃描模式成像時(shí)，我們先將采樣數(shù)據(jù)進(jìn)行B掃描模式處理，并將處理后的結(jié)果數(shù)據(jù)組織成3D紋理數(shù)組存放在顯存之中。根據(jù)不同的需要，對(duì)3D數(shù)組進(jìn)行切片提取顯示en-face單層圖像，或者通過(guò)多切片疊加取平均顯示眼底組織結(jié)構(gòu)信息。

我們采用的3D數(shù)據(jù)塊為480幀共計(jì)204 MB數(shù)據(jù)（每幀數(shù)據(jù)大小為249線×896 像素/線×2 字節(jié)/像素）進(jìn)行C掃描模式成像。3D數(shù)據(jù)塊平均成像速度為1.8 s，能快速實(shí)現(xiàn)視網(wǎng)膜en-face單層切片成像或en-face多切片疊加平均成像。成像效果如圖4所示，其中圖(a)～(c)分別顯示了視網(wǎng)膜en-face單層切片圖像；沿縱軸方向看，(b)切片比圖4(a)切片深約30 μm，(c)切片比(b)切片深約30 μm。圖4(a)～(c)清晰地顯示了視網(wǎng)膜血管、微血管、黃斑等眼底組織細(xì)微結(jié)構(gòu)。視網(wǎng)膜en-face單層切片成像的軸向分辨率可達(dá)到5 μm。與眼底相機(jī)只能對(duì)眼底復(fù)合結(jié)構(gòu)成像相比，視網(wǎng)膜en-face單層切片成像可以對(duì)眼底視網(wǎng)膜下的細(xì)微組織的深度和層次進(jìn)行準(zhǔn)確定位。圖4(d)～(f)分別顯示了視網(wǎng)膜11個(gè)en-face 切片疊加后再取平均的成像圖像。沿縱軸方向看，圖4(d)～(f)、分別以圖4(a)～(c)的切片為中心，前后各取5切片疊加后取平均所成圖像。從圖4(d)～(f)能清晰地觀察到視網(wǎng)膜血管、微血管、黃斑等眼底組織細(xì)微結(jié)構(gòu)，與圖4(a)～(c)單切片成像圖相比，圖4(d)～(f)包含了更多的復(fù)合結(jié)構(gòu)信息。

圖4 視網(wǎng)膜en-face成像圖像Fig.4 En-face imaging images of the retina

實(shí)驗(yàn)結(jié)果表明：在成像圖像質(zhì)量不變的前提下，在眼科OCT系統(tǒng)中，采用GPU+CPU模式實(shí)現(xiàn)成像數(shù)據(jù)處理的速度，較傳統(tǒng)的基于CPU平臺(tái)的串行計(jì)算和成像模式執(zhí)行同樣數(shù)據(jù)處理的速度提高超過(guò)一個(gè)數(shù)量級(jí)。

4 結(jié)語(yǔ)

本文利用實(shí)驗(yàn)室現(xiàn)有標(biāo)準(zhǔn)OCT系統(tǒng)平臺(tái)，在未增加硬件的基礎(chǔ)上，利用計(jì)算機(jī)通用顯卡GPU，并將基于GPU的統(tǒng)一計(jì)算設(shè)備架構(gòu)（CUDA）引入到眼科OCT系統(tǒng)成像中的數(shù)據(jù)處理過(guò)程，借助GPU強(qiáng)大的并行數(shù)據(jù)處理能力和浮點(diǎn)計(jì)算能力，用CUDA對(duì)OCT系統(tǒng)數(shù)據(jù)處理過(guò)程進(jìn)行改寫，使得眼科OCT系統(tǒng)的成像速度較基于CPU平臺(tái)處理成像速度提高了數(shù)十倍，達(dá)到了臨床2D實(shí)時(shí)成像的要求，為眼科3D實(shí)時(shí)成像打下了基礎(chǔ)。

[1]D.Huang,E.A.Swanson,C.P.Lin,et al.Optical coherence tomography[J].Science,1991,254(5035): 1178-1181.

[2]R.Leitgeb, C.K.Hitzenberger, A.F.Fercher.Performance of fourier domain vs.time domain optical coherence tomography[J].Opt Express,2003,11(8): 889–894.

[3]Ruikang，K.Wang.In vivo structural and flow imaging: US,US20100027857A1[P],2010-2-4.

[4]B.Cense,T.C.Chen,B.H.Park,et al.Thickness and birefringence of healthy retinal nerve fiber layer tissue measure measured with polarization-sensitive optical coherence tomography[J].Invest Ophth Vis Sci,2004,45(12): 2606-2612.

[5]N.D.Gladkova,G.A.Petrova,N.K.Nikulin,et al.In vivo optical coherence tomography imaging of human skin: norm and pathology[J].Skin Res Technol,2000,6(1): 6-16.

[6]L.Vabre,A.Dubois,A.C.Boccara.Thermal-light full-field optical coherence tomography[J].Opt Lett,2002,27(7): 530-532.

[7]B.Potsaid,I.Gorczynska,V.J.Srinivasan,et al.Ultrahigh speed spectral/Fourier domain OCT ophthalmic imaging at 70,000 to 312,500 axial scans per second[J].Opt Lett,2008,16(19): 15149–15169.

[8]G.Liu,J.Zhang,L.Yu,et al.Real-time polarization-sensitive optical coherence tomography data processing with parallel computing[J]. Appl Opt,2009,48(32): 6365–6370.

[9]J.Probst,P.Koch,G.Huttmann.Real-time 3D rendering of optical coherence tomography volumetric data[C].Proc SPIE,2009,7372,73720Q.

[10]A.W.Schaefer,J.J.Reynolds,D.L.Marks,et al.Real-time digitalsignal processing-based optical coherence tomography and Doppler optical coherence tomography[J].IEEE Trans Biomed Eng,2004,54: 186-190.

[11]T.E.Ustun,N.V.Iftimia,R.D.Ferguson,et al.Real-time processing for Fourier domain optical coherence tomography using a field programmable gate array[J].Rev Sci Instrum,2008,79:114301-114310.

[12]Kang Zhang,Jin U.Kang,Real-time 4D signal processing and visualization using graphics processing unit on a regular nonlinear-k Fourier-domain OCT system[J], Opt Express,2010,18(11): 11772-11784.

[13]Jason Sanders,Edward Kandrot.GPU高性能編程CUDA實(shí)戰(zhàn)[M].北京: 機(jī)械工業(yè)出版社,2011

[14]J.S.Schuman,M.R.Hee,C.A.Puliafito.et al.Quantification of nerve fiber layer thickness in normal and glaucomatous eyes using optical coherence tomography[J].Arch Ophth,1995,113(5): 586-596.

[15]M.Altaweel,M.Ip.Macular hole: improved understanding of pathogenesis,staging,and management based on optical coherence tomography[J].Semin Ophth,2003,18(2): 58-66.

[16]紀(jì)淑興,張軍軍,唐健,等.中心性滲出性脈絡(luò)膜視網(wǎng)膜病變的光學(xué)相干斷層掃描圖像[J].中華眼底病雜志,2002,18(2): 121-124.