多通道快速GPU 光线投射算法

林仁回，冯前进，罗博博，贠照强

（南方医科大学 生物医学工程学院，广东 广州510515）

0 引 言

直接体绘制 （direct volume rendering，DVR）［1］是高维图像分析中一种常用的可视化技术，光线投射法 （ray casting，RC）是直接体绘制的经典算法，重建图像质量高，但计算量很大，过程耗时，难以满足实时应用要求。因此，研究高效的RC方法是医学图像可视化的热点问题。

目前快速RC 方法研究主要集中在硬件加速技术与高效算法2方面。如VolumePro［2］系列专业图形卡，可以实现实时的DVR，但此类专业图形卡的价格过于昂贵，限制了它的广泛使用；随着通用图形处理器 （graphics proce－ssing unit，GPU）的发展，大量研究开始集中于采用GPU的通用计算管线完成快速RC［3］。在高效算法研究方面，主要包括降采样技术［4］、射线提前终止法［5］、差分采样［6］、虚拟采样［7］、空域跳跃技术［8］等。这些方法的核心都是利用沿光线方向上相邻数据间的相关性，降低采样点的个数，以减少计算量。实际上在三维数据场中与光线垂直方向上相邻的数据之间也具有很强的相关性，而当前的方法并未充分利用这一特性。为此，本文提出一种新的快速RC 算法，通过对成像平面降采样，减少光线数量，生成低分辨率重建图像，然后通过图像插值方法［9］得到高分辨率的重建图像。同时本文结合GPU 的通用计算管线，给出了算法的GPU 并行计算方案。

1 方 法

1.1 算法概述

现有的光线投射算法为了降低算法的运算量，多为沿光线方向进行降采样，通过减少采样点的数量降低计算量。具体做法如图1 （a）所示 （以平行投影为例），沿着光线前进的方向按照固定的步长 （通常大于三维数据体素大小）对体数据进行采样，把相邻2 个采样点 （图1 中的C1和C2）之间包含的体素属性信息视为与当前采样点的属性高度相关，通过插值方法对跳过的体素信息进行估计。本文方法根据奈奎斯特采样定理，利用光线垂直方向上相邻的数据之间的相关性，在沿光线的方向以体素大小为间隔进行采样 （如图1 （b）所示），而对光线进行降采样，通过减少入射光线生成低分辨率重建图像，再利用图像插值重建算法生成高分辨率图像。

图1 正常的光线投射算法与本文方法采样方式对比

本文方法分2 步完成 （如图2 所示）：①光线降采样RC；②高分辨率重建。通道一为一次GPU 运算，主要为渲染分辨率下的光线投射算法；通道二为一次GPU 运算，主要由渲染分辨率重建得到显示分辨率下的结果。其中显示分辨率是最终重建结果图像的分辨率，初始情况下显示分辨率与渲染分辨率相同，之后进行正常采样下的光线投射算法，得到渲染分辨率下的重建结果，并计算渲染时间，如果每秒帧率小于25帧，系统自动降低渲染分辨率，重新进行光线投射算法，当每秒帧率大于等于25帧时，记录下此时的渲染分辨率，之后均以此渲染分辨率进行重建。为了得到显示分辨率下的结果，将渲染分辨率下重建的结果重新作为管道线的输入，利用图像插值重建算法进行处理得到显示分辨率下的重建结果。

1.2 光线降采样RC

1.2.1 原理

图3描述了本文方法通道一光线降采样。具体可分为以下步骤：

（1）其中成像平面区域EFGH 表示显示分辨率大小，ABCD 表示渲染分辨率大小。初始条件下这2 个分辨率相同，为了达到实时交互的需要，渲染分辨率ABCD 要小于显示分辨率EFGH。

图2 多通道快速GPU 光线投射算法流程

图3 光线降采样RC

（2）渲染分辨率下，根据视点和渲染分辨率下成像平面ABCD 中任意一个像素点P生成一条投射光线，沿光线前进方向计算与三维数据场的交点，如图3中的点I；

（3）从入射点I开始，沿着光线的方向，在三维数据场内进行采样，利用传输函数将数据值转换为光学信息，进行光线合成，直到出射点J停止计算；

（4）针对成像平面ABCD 的每个像素点均生成一条光线穿过数据场，重复步骤 （3），到成像平面上的所有点计算完成结束；

（5）利用OpenGL 离屏绘制技术［10］（frame buffer object，FBO）将重建结果保存到纹理对象。

1.2.2 GPU 实现

光线降采样RC涉及大量的插值运算和颜色混合运算，计算量巨大，且每一条光线都是相对独立的对三维体数据进行采样，非常适合利用GPU 的流式并行架构进行加速［11］。因此本文编写了片元着色程序，针对每条光线使用相同的片元着色程序进行颜色值和不透明度值的积分运算，只是各自处理的数据有所不同。具体实现框架如图4所示。

图4 光线降采样RC实现框架

OpenGL提供的FBO 离屏绘制技术允许执行离屏渲染、动态地更新纹理图像。首先将体数据读出后以三维纹理形式一次性全部装入纹理内存，然后创建一个帧缓存对象和纹理对象，并对它们的状态进行初始化。帧缓存对象并不提供存储空间，所以它需要关联纹理对象用于存储渲染分辨率图像，纹理的大小就是渲染分辨率图像的大小。为了在重建图像质量和绘制效率之间取得一个好的平衡点，渲染分辨率图像的分辨率往往需要根据最终显示图像的分辨率进行动态调整，直到满足实时帧速为止。这些初始化工作都是在CPU 端进行的。由于OpenGL 是一种状态机，OpenGL的相关属性只需要初始化一次就可以保持当前的状态进行重复渲染。为了执行离屏渲染，每次执行光线降采样RC 之前，先把绘图缓冲区绑定到应用程序创建的帧缓存对象，然后触发GPU 执行片元着色程序，并利用提前光线终止技术加速。当光线降采样RC 执行完毕后渲染分辨率图像就会存储到帧缓存对象所关联的纹理。

1.3 高分辨率重建

1.3.1 原理

图像放大几乎都是采用内插值方法，即在原有图像像素的基础上在像素点之间采用合适的插值算法插入新的元素。传统的图像插值算法主要包括最近邻插值、双线性插值、双三次插值、B 样条插值与Catmull－Rom 样条插值等［12］。最近邻插值会产生明显的锯齿边缘和马赛克现象。双线性插值法具有平滑功能，能有效地克服最近邻插值法的不足，但会退化图像的高频部分，使图像细节变模糊。高阶插值，如双三次插值和三次样条插值等比低阶插值效果好。Catmull－Rom 样条是插值分段三次曲线，具有C1连续，插值计算仅依赖于邻近的4个节点，计算效率高，插值效果比较好，因此本文采用其进行高分辨率图像插值重建。Catmull－Rom 样条函数由式 （1）［13］给出

双三次插值利用了最近的16个像素进行插值产生新的像素，图5描述了利用最近的4＊4个像素 （从F（p－1，q－1）到F（p＋2，q＋2），其中F（p，q）是最邻近F（p′，q′）的像素点）进行插值产生新的像素F（p′，q′）。被插值像素F（p′，q′）由式 （2）［14］给出，其中F（p＋m，q＋n）表示点（p＋m，q ＋n）处的像素值，R（x）代表一个双三次插值函数。

图5 双三次插值图

为了获得更好的图像重建效果，本文采用式 （2）对渲染分辨率图像进行双三次插值重建，其中式 （2）中的R（x）由式 （1）替代。为了利用GPU 的可编程管道线，我们编写了基于Catmull－Rom 插值算法的片元着色程序，将光线降采样RC的重建结果重新输入到GPU 之中，利用片元着色程序对渲染分辨率图像进行插值处理，填补丢失的图像细节信息，得到显示分辨率下的重建结果。

1.3.2 GPU 实现

高分辨率重建是在二维图像空间进行插值重建，其重建时间远远低于在三维图像空间进行光线降采样RC 所耗费的时间。为了有效利用GPU 的纹理插值功能和流式并行计算能力提高图像的重建速度，本文同样实现了Catmull－Rom 插值算法片元着色程序。渲染分辨率图像有2个用途。其一是绑定为片元着色程序的2D 采样器，作为Catmull－Rom 插值算法的纹理提取目标；其二是将渲染分辨率图像绑定为贴图纹理，借助OpenGL 纹理映射功能实现显示。在渲染最终图像之前切换到窗口系统提供的帧缓冲区，并把视口大小调整到最终显示图像的分辨率，利用纹理贴图技术触发GPU 执行Catmull－Rom 插值算法片元着色程序，对渲染分辨率图像进行插值重建得到高分辨率图像，最后绘制到帧缓冲区。

2 实验结果分析与讨论

本文的实验测试平台如下：CPU AMD Phenom （tm）II X6 1055TProcessor 2.80GHz，RAM 16GB DDR3 1333 MHz，NVIDIA GeForce GTX 275 896 MB显存；编程环境为Windows 7操作系统下Visual Studio 2010，使用C＋＋和OpenGL编程实现，片元着色程序使用nVidia Cg［15］渲染语言；实验中使用的数据信息为：头部CT 数据，长宽高为512×512×460，数据类型是unsigned short，体素的间距为0.49×0.49×0.70；脚部CT 数据，长宽高为512×512×163，数据类型是unsigned short，体素的间距为0.51×0.51×1.25；为了测试在不同数据大小下本文算法的速度，对头部CT 数据进行了降采样处理，分别生成了256×256×230和128×128×115大小的测试数据。

图6为不同算法和参数下头部CT 数据的三维重建效果对比，图7为不同算法和参数下脚部CT 数据的三维重建效果对比。统一为800×800显示分辨率，图6和图7中第1行为渲染的整体效果，第2行为整体渲染效果的局部放大图，对应于第1 行中的红色方框内的区域。其中图6 （a）和图7 （a）是正常采样频率下的光线投射算法重建结果；图6和图7中图 （b）、图 （c）、图 （d）分别为本文算法重建效果，图6 （b）、图7 （b）采用的渲染分辨率为400×400，图6 （c）、图7 （c）采用的渲染分辨率为300×300，图6d﹑图7d 采用的渲染分辨率为200×200；图6 （e）、图7 （e）为降采样下的重建结果，二者的采样步长分别是图6 （a）、图7 （a）正常采样的6倍。从图6各图的局部放大效果图中可以看出，图6 （a）的重建效果最好，图6（e）的重建质量最差，由于采样频率的降低，造成重建结果中存在大量的伪影和失真；从图6 （b）到图6 （d），图像重建依次有所下降，主要表现在图像细节的模糊，但与图6 （e）相比，其图像质量有了极大的提高。从图7各图的局部放大效果图中可以看出，图7 （a）的重建效果最好，图7 （e）的重建质量最差，图像已经严重失真。从图7 （b）到图7 （d），由于渲染分辨率的顺序递减，图像重建质量也依次有所下降，主要表现在图像细节的模糊和边缘锯齿化，但与图7 （e）相比，其图像质量也有了很大的提高。

与图6、图7相对应的不同重建方法和参数下速度的对比见表1：其中表的第1列表示不同的人体部位CT 数据；其中表的第2列表示不同的重建方法和参数；第3列是重建时间，计算一帧需要的时间，单位为毫秒；第4 列为帧率，表示每秒渲染的帧数；第5列是帧率倍数，计算方法为用当前方法帧率除以对应的正常采样方式下的帧率，表示与正常采样方法相比，速度提升的倍数。正常采样步长下，即使采用GPU 进行加速，光线投射算法计算速度仍然非常慢；而降采样处理后得到的重建图像，由于采样步长增大到了原来的6倍，头部重建速度达到了正常采样的3.4倍，帧速达到了24 帧／秒，基本达到了实时要求，而脚部重建速度达到了正常采样的5.4倍，帧速达到了38帧／秒，已经满足实时要求，但与图6和图7其它方法相比，图像的细节信息丢失严重，且重建结果严重失真，已经无法识别人体的结构部位。表1中第5列的第2、3、4行和第7、8、9行分别为本文方法采用不同渲染分辨率下头部、脚部重建速度的测试结果，随着渲染分辨率的降低，重建速度也不断的增加，其中在200×200渲染分辨率下的重建速度分别是正常采样速度的3.0 倍、4.7 倍，基本与降采样速度相当。

图6 不同算法头部重建效果，其中 （e）的采样步长是其它算法的1／6

图7 不同算法脚部重建效果，其中 （e）的采样步长是其它算法的1／6

表1 不同人体部位CT 数据在不同采样方式下的重建时间比较

通过对重建图像质量和速度的对比分析，本文方法相对于正常采样和降采样方法，可以在图像质量和重建速度中取得折衷，虽然重建结果中图像的细节有所模糊，但与降采样方法结果相比，其结果更易于被临床所接受。

本文方法在800×800显示分辨率和300×300渲染分辨率情况下，对不同大小的头部CT 数据进行光线投射算法的重建速度结果见表2，其中第1列为数据的维数，表示数据的长宽和高；其它2 列和表1 相同。可以看到针对256＊256＊230的数据，采用本文的方法基本可以达到了58帧／秒，完全满足了用户实时交互的需求。

表2 300×300渲染分辨率下不同大小的头部CT 数据重建时间比较

3 结束语

本文提出了一种多通道快速GPU 光线投射算法，主要将超分辨率图像重建算法与光线投射算法相结合，以满足人机交互过程中对重建图像质量和速度的双重需求。实验结果表明，该方法与正常采样光线投射算法相比，在不同渲染分辨率下，重建速度有了很大的提高，同时重建结果中出现了一定的模糊现象，但与降采样方法重建结果失真的情况相比，重建质量有了极大的提高，可以在重建图像质量和速度上取得有效的折衷。后期工作我们会针对医学数据体绘制的特点，研究体数据矢量场的变化特征，控制采样光线在三维体数据上的分布，使得体数据梯度变化率大的区域入射光线密集，梯度变化率小的区域入射光线稀疏，从而实现更快速更高质量的的三维体绘制。

［1］Bruckner S，Grller ME.Instant volume visualization using maximum intensity difference accumulation ［J］.Computer Graphics Forum，2009，28 （3）：775－782.

［2］Karadayi K，Managuli R，Kim Y.Three－dimensional ultrasound：From acquisition to visualization and from algorithms to systems［J］.IEEE Reviews in Biomedical Engineering，2009，2：23－39.

［3］Jeong WK，Beyer J，Hadwiger M，et al.Ssecrett and neurotrace：Interactive visualization and analysis tools for large－scale neuroscience data sets ［J］.Computer Graphics and Applications，IEEE，2010，30 （3）：58－70.

［4］Dorgham OM，Laycock SD，Fisher MH.GPU accelerated generation of digitally reconstructed radiographs for 2－D／3－D image registration ［J］.IEEE Transactions on Biomedical Engineering，2012，59 （9）：2594－2603.

［5］Zhang Q，Eagleson R，Peters TM.Volume visualization：A technical overview with a focus on medical applications ［J］.Journal of Digital Imaging，2011，24 （4）：640－664.

［6］Knoll A，Hijazi Y，Westerteiger R，et al.Volume ray casting with peak finding and differential sampling ［J］.IEEE Transactions on Visua－lization and Computer Grap－hics，2009，15 （6）：1571－1578.

［7］Lee B，Yun J，Seo J，et al.Fast high－quality volume ray cas－ting with virtual samplings［J］.IEEE Transactions on Visualization and Computer Grap－hics，2010，16 （6）：1525－1532.

［8］Liu B，Clapworthy G，Dong F，et al.Octree rasterization：Accelerating high－quality out－of－core GPU volume rendering［J］.IEEE Transactions on Visualization and Computer Graphics，2013，19 （10）：1732－1745.

［9］Yang J，Wright J，Huang TS，et al.Image super－resolution via sparse representation ［J］.IEEE Transactions on Image Processing，2010，19 （11）：2861－2873.

［10］Dave Shreiner，The Khronos OpenGL ARB Working Group.OpenGL programming guide ［M ］. Boston： Addison－Wesley，2009.

［11］Smelyanskiy M，Holmes D，Chhugani J，et al.Mapping high－fidelity volume rendering for medical imaging to CPU，GPU and many－core architectures［J］.IEEE Transactions on Visualization and Computer Graphics，2009，15 （6）：1563－1570.

［12］CHEN Liping，LIU Wenhua.Catmull－Rom image expansion based on Canny edge detection ［J］.Application Research of Computers，2012，29 （10）：3993－3996 （in Chinese）. ［陈利平，刘文华.基于Canny边缘检测的Catmull－Rom 图像放大 ［J］.计算机应用研究，2012，29 （10）：3993－3996.］

［13］Burger W，Burge MJ.Principles of digital image processing：Core algorithms［M］.Germany：Springer，2009.

［14］Hosseini SM，Farsi H，Yazdi HS.Best clustering around the color images ［J］.Int J Comput Electr Eng，2009，1 （1）：20－24.

［15］ZHANG Shu，CHU Yanli.High performance computing using GPU CUDA ［M］.Beijing：China Water Power Press，2009：126－167 （in Chinese）. ［张舒，褚艳利.GPU 高性能运算之CUDA ［M］.北京：中国水利水电出版社，2009：126－167.］