APP下载

GPU在眼科FD-OCT系统实时图像显示及数据处理中的应用

2013-12-05刘巧艳李跃杰徐秋晶赵金城王立伟高用贺

中国医疗器械杂志 2013年1期
关键词:扫描模式插值数据处理

【作 者】刘巧艳,李跃杰,徐秋晶,赵金城,王立伟,高用贺

1 北京协和医学院,北京市,100730

2 中国医学科学院生物医学工程研究所,天津市,300192

3 天津市眼科医学设备技术工程中心,天津市,300384

0 引言

光学相干层析成像技术(Optical Coherence Tomography),自上世纪90年代被成功用于眼科疾病诊断[1]后,得到了迅速发展:由时域OCT阶段发展到频域OCT阶段[2];由组织结构成像向组织功能成像发展[3-4];由眼科诊断拓展到心血管、皮肤、口腔、组织工程等领域中的应用[5-6]。

随着超高速CMOS线阵扫描相机的发展,频域OCT光谱谱线转换及线采样率已经可以达到300 k线/s[7],为临床OCT系统实时成像提供了前提。目前影响眼科OCT系统实时成像的技术瓶颈是需要先将采样数据进行频谱域空间(λ空间)到波数空间(K空间)变换,进行插值变换和FFT变换,然后再将变换后的数据进行2D或3D成像。由于成像的数据量很大,特别是进行C模式扫描成像(如眼底视网膜en-face成像模式)时,需要先将获得的3D图像数据进行处理后,再将得到的数据成像。因此,如何提高数据处理的速度,成为眼科OCT系统实现实时成像的关键。

为了克服眼科OCT系统实时成像的技术瓶颈,近年来研究人员进行了大量研究,并提出了一些相关的解决方法。有研究人员将多核CPU引入到OCT系统进行采样数据的并行处理,对于1024个采样点模式OCT系统,非均匀K空间数据处理的速度可达到80 k线/s[8],均匀K空间数据的处理速度可达到207 k线/s[9]。也有研究人员通过在原有系统中增加DSPs[10]或者FPGAs[11]硬件模块,加快实现数据处理过程。随着图形处理单元(GPU)技术的发展,其运算能力和可编程性得到大幅度提高,它除了在传统的图像处理领域应用继续保持优势外,作为通用的并行计算处理器已被广泛用于科学计算、工程、金融以及其他领域中。目前,有研究人员将GPU引入到OCT系统中,利用它强大的并行计算能力来加速采样数据处理过程,解决OCT系统实时成像的技术瓶颈。Kang Zhang和Jin U.Kang将GPU引入到OCT系统中,在专业图像工作站平台上采用线性插值算法,在1024个采样点模式下数据处理速度达到680 k线/s,2048个采样点模式下数据处理速度可达到320 k线/s,达到了3D实时成像对于处理速度的要求[12]。现阶段将GPU应用于OCT系统还处于实验室研究阶段,为获得更高处理速度,往往需要配置高性能的专用图像工作站(配置多核处理器)和GPU,硬件成本较高。

本课题在不改变实验室现有硬件平台的条件下,把低成本GPU引入到我们正在开发应用的眼科OCT仪器中,实现仪器性能的大幅度提高,解决眼科OCT系统实时成像的问题。

1 系统组成和工作原理

图1 FD-OCT成像系统组成Fig.1 FD-OCT imaging system configuration

1.1 系统组成

本系统组成如图1所示,主要包括眼科OCT信号采集系统和信号处理系统两大部分。其中眼科OCT信号采集系统采用宽带超亮发光二极管(SLD)(λ0=840 nm,△λ= 50 nm)作为系统光源,以2048像素CMOS线阵相机(采样速率为70 k线/s)作为光谱仪的检测器。信号处理系统由计算机和GPU组成。计算机CPU为 Intel Celeron Dual-Core E3400 @ 2.60 GHz,2 G内存;GPU采用 NVIDIA 公司的GeForce GTX 460显卡(336个流处理器,1 GB显存)。

1.2 系统工作原理

光源发出的光经过50:50的光纤分束器后,被均匀分成两束光,分别进入OCT系统的参考臂和样品臂。从样品臂反射回来的信号光和从参考臂返回的参考光再次经过光纤分束器汇合后发生干涉。包含样品不同深度信息的干涉信号光谱,经光谱仪的CMOS线阵扫描相机采集,并通过相机数据线传输到计算机,由计算机内的图像采集卡对干涉信号光谱进行A/D转换,并将转换结果作为采样数据存入到计算机内存中。将采样数据通过PCIE×16总线传输到GPU显存,借助GPU强大的并行数据处理能力进行数据处理,并将处理好的结果数据送回计算机进行图像显示。显示的图像包含了检测样品不同深度的结构信息。

2 基于GPU加速技术的数据处理流程

2.1 CUDA架构

CUDA(Compute Unified Device Architecture)是一种由NVIDIA公司推出的通用并行计算架构,该架构使GPU能够解决复杂的计算问题。在CUDA架构下,开发人员可以通过CUDA C语言(CUDA C语言是对标准C语言的一种简单扩展)对GPU编程[13]。

CUDA架构中,CPU作为主机(Host),GPU作为协处理器或者设备(Device)。在一个系统中可以存在一个主机和多个设备。 CPU 主要负责进行逻辑性强的事物处理和串行计算,GPU 则专注于执行高度线程化的并行处理任务。CPU 、GPU 各自拥有相互独立的存储器地址空间:主机端的内存和设备端的显存;在CUDA程序中,将运行在GPU上一个可以被并行执行的步骤称为kernel(内核函数)。

图2 CPU-GPU系统数据处理流程图Fig.2 Signal processing flow chart of CPU-GPU hybrid system architecture

2.2 数据处理过程

本系统的数据处理过程如图2所示,其中实箭头所指方向代表数据在不同设备间的流动方向,空箭头所指方向代表数据在GPU内部的流动方向。系统数据处理过程包括对采样数据进行预处理、频谱域空间(λ空间)到波数空间(K空间)变换、插值变换、FFT变换和POST-FFT变换。

在频域OCT系统中,采样数据是通过对OCT的光路系统扫描由相机采集到的,扫描一次得到一列数据(一个A-SCAN)。处理时是一列一列数据进行处理的。针对每列数据彼此相互独立、可以并行处理的特点,利用CUDA架构将OCT系统整个数据处理过程改写成适合在GPU上执行的kernel函数,大大提高了数据处理速度,从而达到系统实时成像的要求。

2.3 程序设计

首先确定系统数据处理过程中的串行部分和并行部分,选择合适的算法,并按照算法确定数据和任务的划分方式,将每个需要并行实现的步骤映射为一个满足CUDA两层并行模型的内核函数。其中频谱域空间(λ空间)到波数空间(K空间)变换模块在主程序中执行;数据预处理、插值运算、POST-FFT运算需要重新改写成适合在GPU上执行的kernel函数;FFT运算部分采用CUDA自带的CUFFT库函数来进行计算,CUDA4.1版本目前已经支持FFT双精度运算。整个程序设计如下:

(1)给采样数据和计算过程中的各个变量分配空间,包括内存空间和显存空间,在Host端CPU程序中,准备好计算要用到的数据,将数据从主机内存传输到显存中。

(2)数据预处理包括数据类型转换和去噪运算。去噪运算是将每一列数据组成的数组(一个A-SCAN)都减去一列同大小的噪声数组。

(3)插值运算实现了在进行FFT变换之前,采样数据K空间分布的均匀化。系统最常用的插值算法包括最邻近插值算法、线性插值算法和三次样条插值算法。

(4)利用CUDA自带的CUFFT库函数对插值运算结果进行FFT变换。

(5)POST-FFT运算包括对FFT运算结果取模、取对数并进行归一化。POST-FFT运算结果存储在体积数组中。

(6)根据不同成像平面的需要,抽取体积数组数据作为GPU结果数据。

(7)将结果数据从显存传输到主机内存中。

2.4 优化分析

在用CUDA对GPU上的kernel函数进行改写时,采用了一些优化方法:

(1)为了使内存和显存之间的数据传输速度更快,在分配主机端内存的时候采用了pinned memory。

(2)使用流运算来隐藏主机端和设备端数据通信的时间。

(3)合理利用读取速度更快的shared memory来存放计算过程中的中间变量。

(4)将已经计算好的在插值运算过程中值保持不变的采样数据的等间隔和非等间隔横坐标值存放在常数寄存器里面。

(5)根据系统资源进行grid和block维度设计。

(6)使用CUDA profiler对CUDA程序进行性能测试,对耗时较长的模块进行优化。

3 实验结果

本系统采用Microsoft Visual Studio2010中集成CUDA Toolkit 32 bit 4.1、CUDA SDK 32 bit 4.1和Nvidia Driver for Windows732 bit为开发环境,对采样数据进行B扫描模式和C扫描模式成像。

3.1 B扫描模式成像

B扫描模式成像图像能提供视网膜断层结构,能清晰地显示视网膜各层细微结构及病理改变,并作出定性或定量分析,目前已成为视网膜疾病和青光眼强有力的诊断工具[14]。

我们采用100帧共计195 MB数据(每帧数据大小为500线×2048 像素/线×2 字节/像素)进行B扫描模式成像。分别采用线性插值算法和三次样条插值算法,利用CUDA提供的计时函数分别对CPU模式和CPU-GPU模式下系统单帧B扫描模式图像成像时间进行计时(计算100帧图像成像时间取平均),实验结果如表1所示。 从表1可知,采用GPU+CPU模式执行成像数据处理的速度较CPU模式执行同样数据处理的速度提高超过一个数量级,其中采用线性插值算法速度提高了60倍,采用三次样条插值算法速度提高了35倍。实验成像效果图如图3所示,其中(a)为采用线性插值算法在CPU模式下系统成像图像;(b)为采用线性插值算法在CPU+GPU模式下系统成像图像;(c)为采用三次样条插值算法在CPU模式下系统成像图像;(d)为采用三次样条插值算法在CPU+GPU模式下系统成像图像。由图3可知,采用相同的插值算法,在CPU模式和CPU-GPU模式下系统成像图像质量无差异,采用三次样条插值算法成像图像质量比采用线性插值算法成像图像质量要好。

图3 视网膜B扫描成像图像Fig.3 B-scan imaging images of the retina

插值运算是系统数据处理过程中计算量最大、耗时最长的一个环节。插值算法的选取不仅影响系统加速比,也影响系统成像速度和成像图像质量。采用线性插值算法,成像速度快,成像图像质量差;采用三次样条算法,成像速度慢、成像图像质量好。

表1 CPU和GPU 成像速度对比Tab.1 Comparison of CPU and GPU imaging speed

3.2 C扫描模式成像

C扫描模式成像图像能直观地显示眼底视网膜血管和黄斑等眼底组织的结构信息,临床上可用于诊断眼底病变,如黄斑病变[15]、眼底神经组织结构变化等,可用于对视网膜下新生血管的深度和层次进行准确定位[16]。

进行C扫描模式成像时,我们先将采样数据进行B扫描模式处理,并将处理后的结果数据组织成3D纹理数组存放在显存之中。根据不同的需要,对3D数组进行切片提取显示en-face单层图像,或者通过多切片叠加取平均显示眼底组织结构信息。

我们采用的3D数据块为480帧共计204 MB数据(每帧数据大小为249线×896 像素/线×2 字节/像素)进行C扫描模式成像。3D数据块平均成像速度为1.8 s,能快速实现视网膜en-face单层切片成像或en-face多切片叠加平均成像。成像效果如图4所示,其中图(a)~(c)分别显示了视网膜en-face单层切片图像;沿纵轴方向看,(b)切片比图4(a)切片深约30 μm,(c)切片比(b)切片深约30 μm。图4(a)~(c)清晰地显示了视网膜血管、微血管、黄斑等眼底组织细微结构。视网膜en-face单层切片成像的轴向分辨率可达到5 μm。与眼底相机只能对眼底复合结构成像相比,视网膜en-face单层切片成像可以对眼底视网膜下的细微组织的深度和层次进行准确定位。图4(d)~(f)分别显示了视网膜11个en-face 切片叠加后再取平均的成像图像。沿纵轴方向看,图4(d)~(f)、分别以图4(a)~(c)的切片为中心,前后各取5切片叠加后取平均所成图像。从图4(d)~(f)能清晰地观察到视网膜血管、微血管、黄斑等眼底组织细微结构,与图4(a)~(c)单切片成像图相比,图4(d)~(f)包含了更多的复合结构信息。

图4 视网膜en-face成像图像Fig.4 En-face imaging images of the retina

实验结果表明:在成像图像质量不变的前提下,在眼科OCT系统中,采用GPU+CPU模式实现成像数据处理的速度,较传统的基于CPU平台的串行计算和成像模式执行同样数据处理的速度提高超过一个数量级。

4 结语

本文利用实验室现有标准OCT系统平台,在未增加硬件的基础上,利用计算机通用显卡GPU,并将基于GPU的统一计算设备架构(CUDA)引入到眼科OCT系统成像中的数据处理过程,借助GPU强大的并行数据处理能力和浮点计算能力,用CUDA对OCT系统数据处理过程进行改写,使得眼科OCT系统的成像速度较基于CPU平台处理成像速度提高了数十倍,达到了临床2D实时成像的要求,为眼科3D实时成像打下了基础。

[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实战[M].北京: 机械工业出版社,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].中华眼底病杂志,2002,18(2): 121-124.

猜你喜欢

扫描模式插值数据处理
滑动式Lagrange与Chebyshev插值方法对BDS精密星历内插及其精度分析
认知诊断缺失数据处理方法的比较:零替换、多重插补与极大似然估计法*
基于低频功率数据处理的负荷分解方法
ILWT-EEMD数据处理的ELM滚动轴承故障诊断
不同扫描方案联合不同锐利卷积核对冠状动脉支架显示的影响
基于pade逼近的重心有理混合插值新方法
混合重叠网格插值方法的改进及应用
双光能X射线骨密度仪测量腰椎骨密度不同扫描模式的对比研究*
更 正
基于希尔伯特- 黄变换的去噪法在外测数据处理中的应用