GPU-Accelerated有限元方法建模光传输扩散光学断层扫描

文摘

我们引入一个GPU-accelerated有限元求解器的计算散射光传输媒体。提出模型的迭代方法计算最昂贵的组件在扩散光学层析图像重建,因此性能优化的解算器是提高效率的关键解决方案的逆问题。GPU提出解决使用CUDA实现评估图形硬件的稀疏线性系统产生的扩散方程的有限元公式。时域和频域的问题提出解决方案。比较CPU-based实现显示了显著的性能提升的图形加速解决方案,改进的双精度计算大约10倍,和因素超出20单精确计算。收益也证明是依赖于网格的复杂性,在最大的收益实现网格分辨率高。

1。介绍

扩散光学层析成像(DOT)是医学应用的功能成像形态有潜力提供散射和吸收参数分布的三维图像在活的有机体内,临床相关的生理参数,如组织和血液氧合状态和状态变化可以派生。应用包括大脑活动可视化(1,2),大脑氧化监测在婴儿3),和乳腺肿瘤检测(4]。

数据采集系统包括一个红外线交付系统,阐明了组织表面在不同的位置,和探测器测量透射光的表面位置。测量可以在连续波(CW)模式中,在输入脉冲超快时间分辨模式和时间分辨的探测器,或在频域模式下,使用调制光源和测量相移调制振幅在探测器的位置。

由于高水平的散射在大多数生物组织图像重建的点是一个不适定的非线性问题的解决方案通常需要制定提出模型的非均匀散射光传播的组织。经常利用光传输模型包括随机模型如蒙特卡罗模拟(5),或确定性模型如辐射传输方程(RTE) [6)或扩散方程(DE) (7]。数值解的方法包括有限差分、有限元、有限体积,或边界元方法。光传输模型是本文的有限元方法(FEM)扩散方程的解决方案。重建问题可以表示为一个非线性优化问题,在一个目标函数,定义为一个标准的测量数据和模型数据的区别对于一个给定的光学参数,是最小化,regularisation功能。重建方法包括方法,只需要提出模型的可用性,如马尔可夫链蒙特卡罗方法,它的一阶导数,如非线性共轭梯度方法,及其二阶导数,如牛顿型方法。

提出模型的迭代解决需要多个评估计算目标函数及其梯度。提出模型本身涉及到解决方案的大型线性系统与多个右手边。问题涉及高维参数空间导致耗时的正向模型的评估,这限制了在临床实践中重建方法的适用性。显著的性能改进应点医学成像的一个可行的工具。最近的事态发展在计算机硬件提供了利用并行计算的可能性。传统上,解决方案包括中央处理单元(CPU)基于适度并行系统共享内存访问(多处理器和多核实现)和大规模分布式并行系统限制了节点之间的数据传输(集群CPU实现)。最近,图形处理单元(GPU)的并行架构利用了通用的加速度计算,包括GPU密度方法的解决方案(8,9)或稀疏(10- - - - - -14线性系统。后者在有限元法的实现遇到了。

上下文中的扩散光学断层扫描和光学成像的相关领域,GPU-accelerated计算已经成功地用于实现蒙特卡罗光线传输模型(15- - - - - -17),计算独立光子轨迹和非常适合parallelisation由于缺乏进程间通信。加速度率超过300是可能的。Zhang et al。18)GPU加速应用于有限元计算相比,生物发光断层扫描和单线程和多线程CPU性能。他们报道GPU版本的显著的性能优势,但仅限于低网复杂性由于内存限制。在光学投影层析成像,基于gpu的重建方法一直受雇于Vinegoni et al。19]。渡边和Itagaki20.)用GPU实现实时可视化傅氏域光学相干断层扫描。

在本文中,我们正在调查潜在的GPU实现正向模型的点。我们现在统一计算设备架构(CUDA)版本之前给出的向前有限元解算器(21,22),使用尖端库(23]稀疏线性系统计算的图形处理器。CUDA是NVidia图形处理器的计算体系结构,可以通过一个应用程序接口(API)来解决。当前GPU硬件单精度算法的性能优化。我们调查的影响单精确计算的准确性提出模型的不同组合的光学参数。我们比较GPU的性能与同等解决CPU实现。我们可以实现显示显著的性能改进。提出的评价模型是最耗时的迭代元素逆动力学,因此解决前进中的任何效率提升直接转化为减少整体图像重建应用程序和运行时是一个重要的一步做点一个可行的成像在临床实践中的应用。

2。方法

2.1。有限元解算器

我们考虑扩散近似辐射传输方程(24,25)在稳态、时间或频率域的光传输模型组织向前发展。对于稳态问题,介质内的固定实光子密度计算连续波来源而引起的频域问题,源振幅调制,光子密度波的复值解分布。在时域问题,源被认为是delta-pulse时间,测量由时间色散的传输信号。给定一个紧凑的域有界的,扩散方程26在时间和频率域是由 分别在哪里是角调制频率源,和是空间不同的扩散和吸收系数,分别在哪里与散射系数。光速在中,是真正的和复数的光子密度字段。以下为简单起见,我们使用表示真正的或复数的性质。

Robin-type边界条件(27)适用于, 在哪里是一个真正的或复数源分布,折射率是一个边界反射术语合并tissue-air接口,表面点的表面正常吗。边界算子定义出射度通过是由Dirichlet-to-Neumann地图吗 测量数据的集合从源分布通过整合在测量资料从表面上看 的时域问题,时间分散配置文件的接收信号强度,频率域的问题,是由复杂的出射度值,通常表达的对数振幅吗和相移(28), 给定的一组数据所有测量数据的源分布,(1)(4)定义模型参数分布的地图吗测量值对于一个给定的域几何、调制频率、源分布,测量资料。

提出模型用有限元数值求解方法。域的一个部门到四面体元素定义的顶点节点提供了一种分段多项式的基础参数和光子密度。近似场在任何时候是由节点的插值系数使用分段多项式函数形状 分段多项式近似,连续参数,定义的节点系数,构建以同样的方式。应用伽辽金方法转换的连续问题(1)到一个维离散问题的节点字段值在所有节点,考虑到节点参数集。频域的问题,得到的线性系统是由 在哪里 ,是对称的稀疏矩阵由(7] 和右边是由 与的节点系数的扩张的基础。时域的解决问题的时间导数(1)时间由有限差分近似 的时间是近似在一组离散步骤吗有限差分方法和评估,给出的迭代 在哪里,时间的步骤,是一个控制参数,可用于选择隐式(),显式(),或者中间方案。这一步的长度由有限差分格式的稳定性问题。无条件稳定的隐格式,可以调整步长时间的曲率,允许增加指数衰减的尾巴的步长。

因此由有限元的解决问题(我)建筑系统的矩阵(9),(2)复数的线性问题的解决方案(7)或实值序列的线性问题(12),(iii)映射到测量(3)和(4)。主要成本计算是线性系统的解决方案,特别是在时域问题,而成本矩阵组装时间通常只有1 - 10%的单一线性的解决方案。线性系统与直接法可以得到解决,如柯列斯基分解的实值时域问题或LU分解的复数频域问题,或与迭代方法,如共轭梯度的实值问题和双共轭梯度的复数的问题。直接方法成为大规模的不切实际的问题,由于内存存储需求分解和增加计算时间。3 d节点密度高的问题,解决迭代通常使用。

2.2。GPU实现

重建的瓶颈问题是解决线性系统(7)或(12)。加速线性规划求解因此提高逆解算器性能的有效方法。我们已经向前有限元解算器的嵌入式图形processor-accelerated版本软件包到现有的烤面包(29日)光传输和重建提出之前(7]。GPU-accelerated代码使用NVidia的CUDA编程接口图形处理器硬件。实现利用尖端图书馆为稀疏线性代数和提供一个模板化的框架提供了共轭梯度(CG)和双共轭gradient-stabilised (BiCGSTAB)迭代解决稀疏线性系统。图书馆支持单和双精度计算如果支持的硬件。

我们使用压缩稀疏矩阵的行(CSR)格式存储。有替代存储格式如坐标、ELLPACK、或混合格式(11)可以提供更好的并行性能取决于矩阵填充结构,通常在更紧凑的存储的成本。然而,CSR格式构成良好的性能和通用性和之间的妥协是适合矩阵填充分布引起的非结构化有限元网格。

当线性问题的解决方案(7),我们扩展的复杂系统进真实系统的形式 尖端CG和BiCGSTAB解决者必须修改的提前终止迭代循环由于奇点的中间结果。提前终止条件在实践中偶尔发生在本文中考虑的问题,特别是由于单精度舍入错误。

主机和图形设备内存之间的数据流为一个求解步骤如图1。系统矩阵是聚集在主机内存对于一个给定的一组参数,结合源向量,复制到GPU设备内存。GPU解决所有眼睛,然后调用之后,预计的解决方案复制回主机内存。

有限差分时域问题的解决方案,整个迭代(12)可以评估在GPU最小主机和图形系统之间的通信,组成的初始复制系统矩阵和GPU,返回计算时间配置文件回主机。时域问题的数据流程图如图2。

2.3。单精度运算

GPU硬件通常是单精度浮点操作的优化。虽然GPU硬件与双精度能力是新兴,通常只有一小部分芯片基础设施致力于双精度运算,从而导致一个重大的性能代价。对于优化,因此有利于使用单精度运算,足够了。我们已经实现了有限元求解器在两个单引号和双精度对GPU和CPU平台。

当系统矩阵反映在单精度,必须注意在组装。系统矩阵是由单个元素的贡献(9)。全球系统中的顶点贡献矩阵局部元素的总和为所有元素相邻顶点值到顶点。在求和,精度损失可能发生如果加式之间的大小差异很大尾数相比精度的浮点表示。对于单精度运算,这可能是一个问题在特定的顶点附近有大量的元素,尤其是在3 d网格与四面体元素。期间损失精度矩阵组装可以减少如果贡献排序从最小到最高程度。然而,这会产生一个簿记开销,可以对性能的影响。相反,我们选择在双精度组装系统矩阵和的值映射到单一精度后组装。装配步骤是在主机上执行方面,性能的影响可以忽略不计,因为装配时间通常是小相比,解决时间。

比较的结果矩阵组装在单和双精度,我们执行一个有限元提出解决方案从系统单精度矩阵组装单和双精度。域是一个均匀的圆柱体的半径25毫米和高度50毫米,光学参数和毫米。一个点光源调制频率MHz是放置在缸套。网由83142年和444278年的四面体元素。

图3显示了单和双精度之间的差异提出解决方案在对数振幅(图3(一个)(图)和阶段3 (b))的复杂的光子密度场沿着一条线从源位置在圆柱的体积。实线代表系统的单精度误差矩阵已经聚集在双精度在映射到单一的精度,而虚线是误差引起的系统矩阵在单一的组装精度。可以看出,系统矩阵组装双精度可以显著降低解决方案错误,特别是从源的距离。

光学参数的影响在前进的单精度误差数据如图4。提出的解决方案是计算出三种不同组合的吸收和散射系数(我),(二)毫米⁻¹,毫米⁻¹,(iii)毫米⁻¹,毫米⁻¹。可以看出差异变得更加严重的光学参数值更高。结果是特别敏感的增加散射参数。由于衰减,对象内部的光子密度领域迅速衰减,导致大动态范围的数据。增加吸收和散射参数加剧这种效果,削弱单精确解的准确性,特别是在地区远离源。然而,值得注意的,温和的光学参数在一个典型的光学断层扫描范围,单精度的解决方案是准确的,最大相对误差来分别在对数振幅和相位。

3所示。结果

图形加速推进解决问题上执行一个NVidia GPU GTX 285。设备的技术规格表中列出1。设备都支持双以及单精度运算,所以收集的结果。性能比较,相同的模型计算也具有CPU-based串行实现的执行在2.0 GHz Intel Xeon处理器时钟与4 MB缓存和12 GB的内存。均匀柱面网格的有限元模型包括25毫米半径和高度50 mm在不同元素的决议。80 80源和探测器表面上。一起一个网格,生成的系统矩阵稀疏结构,如图5。

3.1。频域解算器

GPU和CPU运行时性能比较实现在各种情况下,我们评估了频域有限元模型使用不同网格的复杂性。提出解决方案计算的复值使用调制频率的问题兆赫和实值稳态问题。复数的问题,BiCGSTAB线性解算器被用来计算的线性系统(7),用CG解决实问题。我们测试了GPU的性能解决方案作为CG和BiCGSTAB函数收敛公差,没有预调节器或对角预调节器。结果所示的GPU性能的因素,给出的CPU和GPU运行时间的比率。图6显示了单精度(图的性能因素6(一))和双精度(图6 (b))计算。可以看出,GPU实现8和19个单精度之间的性能因素的计算,根据问题类型,在实值BiCGSTAB显示解决方案没有预处理改善14 - 19在最高,而复杂BiCGSTAB展品最小的解决方案没有预处理改善8 - 11.5。一般来说,较低的公差范围的性能系数下降。双精度的解决方案的性能因素是显著降低,介于3.7和4.7之间。这是由于这样的事实,虽然GPU性能显著降低双精度计算,解决CPU性能一般不受影响,事实上CPU性能是在双精度略高,因为它避免了铸造浮点缓冲区单和双精度之间的关系。性能下降因素低公差范围不存在双精度的结果。

下一个测试比较了CPU和GPU性能作为网格节点的密度函数和线性系统的大小。的性能因素提出解决应用于圆柱形网格的网格分辨率节点数的函数图所示7。在每个网格分辨率,我们解决了这两个实值预先处理这些CG稳态问题解决者,和一个复数的频域问题事先已经BiCGSTAB解算器,单精度(图7(一))和双精度(图7 (b))决议。所有求解结果是真实的或复杂的光子密度计算字段80源,为解决者宽容固定。可以看出,在所有情况下GPU性能改善和增加的大小线性系统。单精度解算器、性能因素的范围1 - 26为网格节点数在9000和分别为稳态问题,2到30之间的网格节点数9000分别为频域问题。注意,对于频域问题,无法计算性能因素的两个最大网格由于过度的CPU计算时间的解决方案。绝对的线性规划求解时间选择的情况下表所示2。可以看出,最大的网格分辨率,提出解决次CPU可以超过一个小时。这可以禁止迭代重建的临床应用,每一步的重建可能需要多个评估提出的问题来计算当前估计的目标函数及其梯度或执行一行搜索当前的搜索方向。相比之下,GPU时间这些问题通常需要2到10分钟,这对重建问题是可行的。

提供一个比较CPU-based并行解算器,我们也显示性能因素的一个共享内存线程版本的有限元解算器使用8线程,相比单线程串行实现。线程实现使用粗粒度的parallelisation策略,将线性问题的解决方案不同的右手边在工作线程可用。这种方法提供了更好的处理器利用率和减少通信开销的问题被认为是比精密战略parallelising迭代解算器本身。因为CPU实现没有显示出显著的性能差异单和双精度的解决方案,我们现在只有双精度的结果。图8显示了性能因素2,4,8线程使用CG解决实问题,和使用BiCGSTAB复数的问题解决者。CG解算器到达因素之间1.5(2线)和2.8(8线程),而BiCGSTAB解算器达到因素之间1.7(2线程)和4(8线程)。依赖网的复杂性没有标记为GPU的能手。

3.2。时域解算器

我们计算的有限差分实现时域问题(12)超过100 50秒的时间步长圆筒网格不同的复杂性。对于这些模拟,Crank-Nicholson计划()使用。信号强度时间资料计算80探测器位置为每个80源的位置。性能结果如图9。可以看出,GPU的性能改进实现再次强烈依赖于网格分辨率,从3至13倍双精度运算计算,和从6到17单精确计算。分辨率最高的网,总提出解决CPU运行时间约为8小时实现单和双精度而GPU运行时都是大约29和36分钟的单引号和双精度的解决方案,分别。

4所示。结论

我们已经开发出一种向前的GPU实现有限元模型漫射光传输,可以用作一个组件在一个迭代的非线性扩散光学层析成像重建方法。转发的效率解决重建性能有重大的影响,和减少重建时间是至关重要的在光学断层扫描一个可行的成像方式在临床诊断。

本文提供的模型支持真实和复数的问题,可以应用于稳态,时间,或者频域成像系统。线性系统的有限元discretisation解决与共轭梯度或双共轭梯度迭代求解稳定在GPU设备。我们已经表明,GPU的解算器可以实现显著的性能改进系列CPU实现的因素5至30日,根据网格的复杂性,公差极限,解算器类型。基于gpu的发展解决者提供更高的性能收益比CPU线程并行实现,用于比较。预计未来发展GPU硬件进一步提高性能。

我们已经表明,提出的问题一个精密线性解算器可以应用在临床应用典型的光学参数范围的光学参数。单精度计算收益率更高的性能尤其是图形计算平台。然而,在非常高的吸收和散射参数值,线性系统可能变得越来越坏脾气的,不再与单精度算法收敛。在这些情况下,需要双精度计算。

确认

这项工作是由EPSRC格兰特EP / E034950/1和欧盟第七框架计划(fp7/2007 - 2013)根据授权协议。201076年。

引用

1. d·a·博厄斯·d·h·布鲁克斯·e·l·米勒et al .,“与扩散光学断层扫描成像的身体,”IEEE信号处理杂志,18卷,不。6日,57 - 75,2001页。视图:出版商的网站|谷歌学术搜索
2. b·w·波格k·d·鲍尔森,c .银白杨和h·考夫曼,“绝对校准的近红外频域组织光谱吸收系数定量在新生儿头鬼,”《生物医学光学,5卷,不。2、185 - 193年,2000页。视图:出版商的网站|谷歌学术搜索
3. m .应付和d . t . deply”系统长期测量大脑血液和组织氧合的新生儿通过近红外透照,“医学和生物工程和计算机,26卷,不。3、289 - 294年,1988页。视图:谷歌学术搜索
4. d . j . Hawrysz和e . m . Sevick-Muraca”发展对诊断乳腺癌使用neer-infrared光学测量和荧光成像造影剂,”瘤形成,卷2,不。5,388 - 417年,2000页。视图:谷歌学术搜索
5. d·a·博厄斯j.p.斑鸠,j·j·斯托特和a·k·邓恩,“三维蒙特卡罗光子迁移代码通过复杂的异构媒体包括成年人的头,“光学表达,10卷,不。3、159 - 170年,2002页。视图:谷歌学术搜索
6. g . s . Abdoulaev和a . h . Hielscher”三维光学断层扫描的辐射传输方程,”电子杂志的成像,12卷,不。4、594 - 601年,2003页。视图:出版商的网站|谷歌学术搜索
7. m . Schweiger s . r . Arridge, Nissila,”高斯牛顿法在扩散光学层析图像重建,”物理学在医学和生物学,50卷,不。10日,2365 - 2386年,2005页。视图:出版商的网站|谷歌学术搜索
8. a . Moravansky和n Ag)、“稠密矩阵代数在GPU,”Direct3D ShaderX2w·f·恩格尔,Ed。p。2, Wordware出版,普莱诺,特克斯,美国,2003年。视图:谷歌学术搜索
9. n . Galoppo n . k . Govindaraju m·亨森和d . Manocha”LU-GPU:高效算法求解密度线性系统在图形硬件,”诉讼的ACM和IEEE 2005 SC会议,p。3, IEEE计算机协会,华盛顿特区,2005年12月美国。视图:出版商的网站|谷歌学术搜索
10. j·d·欧文斯d Luebke: Govindaraju et al .,“通用计算图形硬件上的调查,”《计算机图形学论坛34卷,第113 - 80页,布莱克威尔出版社,2007年3月。视图:谷歌学术搜索
11. n·贝尔和m .花环,“有效的稀疏矩阵向量乘法CUDA上,”技术。众议员极- 2008 - 004,NVIDIA公司,2008年版。视图:谷歌学术搜索
12. l . Buatois g . Caumon b·利维,”克朗彻并发数:一个有效的稀疏线性规划求解GPU,”第三届国际研讨会论文集高性能计算和通信,(HPCC ' 07)卷,4782在计算机科学的课堂讲稿施普林格,页358 - 371年,休斯顿,德克萨斯州,美国,2007年9月。视图:谷歌学术搜索
13. j . Krager r·韦斯,“线性代数运算符为GPU实现数值算法,”ACM交易图片22卷,第908 - 916页。视图:谷歌学术搜索
14. j . Bolz即农民、大肠Grinspun·施罗德“解决稀疏矩阵在GPU上:共轭梯度和多栅的,”ACM交易图片22卷,第917 - 924页。视图:谷歌学术搜索
15. 大肠Alerstam、t·斯文森和s . Andersson-Engels”并行计算与图形处理单元高速光子迁移的蒙特卡罗模拟,”《生物医学光学,13卷,不。6、文章ID 060504, 2008。视图:出版商的网站|谷歌学术搜索
16. 问:方舟子和d·a·博厄斯“蒙特卡罗模拟3 d光子迁移的浑浊的媒体加速图形处理单元,“光学表达,17卷,不。22日,第20190 - 20178页,2009年。视图:出版商的网站|谷歌学术搜索
17. n . Ren j .梁,x, j . Li b . Lu和j .田,“基于gpu的蒙特卡罗模拟光传播的复杂的多相组织,”光学表达,18卷,不。7,6811 - 6823年,2010页。视图:出版商的网站|谷歌学术搜索
18. b .张x杨,杨f . et al .,“CUBLAS和CULA基于GPU加速生物发光断层的自适应有限元框架”光学表达,18卷,不。19日,20201 - 20213年,2010页。视图:出版商的网站|谷歌学术搜索
19. c . Vinegoni l . Fexon p . f . Feruglio et al .,“高通量传输光学投影层析用低成本的图形处理单元,“光学表达,17卷,不。25日,第22332 - 22320页,2009年。视图:出版商的网站|谷歌学术搜索
20. 渡边y和t . Itagaki傅立叶域光学相干断层扫描系统上实时显示使用图形处理单元,“《生物医学光学,14卷,不。6、文章ID 060506, 2009。视图:谷歌学术搜索
21. s . r . Arridge m . Schweiger m . Hiraoka和d . t . deply”有限元方法建模光子在组织运输,”医学物理学,20卷,不。2、299 - 309年,1993页。视图:出版商的网站|谷歌学术搜索
22. m . Schweiger和s . r . Arridge”有限元方法在散射光的传播媒体:频域情况下,“医学物理学,24卷,不。6,895 - 902年,1997页。视图:出版商的网站|谷歌学术搜索
23. 尖端库,http://code.google.com/p/cusp-library/。
24. s . r . Arridge”在医学成像光学断层扫描逆问题,15卷,不。2,R41-R93, 1999页。视图:谷歌学术搜索
25. 公元前m . s . Patterson, b .机会和威尔逊,“时间分辨反射率和透射率的非侵入式测量组织光学特性,”应用光学,28卷,不。12日,第2336 - 2331页,1989年。视图:谷歌学术搜索
26. 石丸,波传播和散射随机介质,1卷,学术出版社,纽约,纽约,美国,1978年。
27. r . c . Haskell l . o . Svaasand t t -蔡,t·c·冯·m·s·麦克亚当斯和b . j . Tromberg”辐射传输的扩散方程,边界条件”美国光学学会杂志》上,11卷,不。10日,2727 - 2741年,1994页。视图:谷歌学术搜索
28. Nissila, k . Kotilahti k Fallstro, t . Katila”仪器精确测量的相位和振幅在光学断层扫描中,“审查的科学仪器,卷73,不。9日,第3312 - 3306页,2002年。视图:出版商的网站|谷歌学术搜索
29. m . Schweiger和s . r . Arridge“吐司重建计划,”http://toastplusplus.org。视图:谷歌学术搜索

版权

版权©2011年马丁Schweiger。这是一个开放的分布式下文章知识共享归属许可,它允许无限制的使用、分配和复制在任何媒介,提供最初的工作是正确引用。