随着医疗卫生领域中虚拟现实技术的不断发展, 虚拟手术系统降低手术成本, 制定手术方案, 训练外科医生等方面发挥了重大的作用.

近年来, 国内外学者对医学图像的三维可视化展开了研究并建立了特定的虚拟手术系统. Lee 等^[1]通过CT扫描的方法构建了下颌和上颌的三维数字模并开发了一个骨骼复杂模型引导的虚拟正颌外科系统. Maini等^[2]利用虚拟现实技术在髋臼骨折手术前为患者特制预轮廓钢板模板, 有助于简化手术程序, 提高手术的精确性. Marinković等^[3]将简化的几何非线性共旋转有限元公式用于模拟虚拟现实手术中的高度非线性的软组织形变以提高数值运算的精度. Rahman等^[4]讨论了利用三维虚拟手术与自编语言矫治器综合治疗特发性髁突吸收及继发颌骨畸形的方法. 汪军等^[5]搭建了一种具有力反馈手感的胆囊切除虚拟手术仿真训练平台. 肖亮等^[6]将计算机三维手术计划系统用于肝脏外科教学中, 取得了良好的效果. 莫建清等^[7]针对膝关节镜手术术前规划的需求, 设计并完成膝关节镜手术规划及训练系统的开发工作. 高洪林等^[8]采用三维数字化方法建立了踝关节骨折数字化虚拟手术设计系统. 黄铭明等^[9]使用Unity3D和HTC VIVE开发了一种下颌骨骨折虚拟手术培训系统.

真实感绘制和实时性研究是虚拟手术领域广泛讨论的热点问题. Dong等人^[10]提出了一种基于二维纹理贴图的体纹理合成方法. 该方法在一些纹理细节较少且周期性特征较强的体纹理合成中具有良好的合成效果及合成速度, 但在一些纹理细节丰富且无明显结构性特征的纹理合成中效果不佳. 潘翔^[11]提出了一种基于复用计算的体纹理合成方法, 采用复用计算的肝脏体纹理合成速度要明显高于非复用计算. 但随着肝脏体模型表面三角面片数目的增大, 将肝脏体纹理从纹理空间映射到肝脏体模型的过程的耗时仍是尚未解决的一大难题. 综合国内外学者在纹理合成领域的相关研究成果来看, 但多数体纹理是通过单张二维样图纹理来合成的, 在物体内部无法表现出丰富的纹理细节, 因此很难绘制出具有较高真实感的体纹理. 另外, 体纹理的绘制需要对庞大的纹理数据量进行计算和处理, 因此很难在保证绘制效果的同时兼顾绘制的速度. 尤其在虚拟手术领域, 需要在建立精细的手术模型及场景的同时还能够满足虚拟手术的实时性和交互性要求, 这对现有的硬件性能是一个极大的挑战.

本文的研究基于计算机统一设备架构(Compute Unified Device Architecture, CUDA), 以具有高度并行计算能力的GPU为硬件基础, 对体纹理合成与映射及肝脏模型切割的方法进行改进, 使其与GPU的硬件特性相结合, 以此来提高虚拟肝脏手术系统中肝脏体纹理合成与映射以及模型切割的速度, 从而兼顾其对真实感与实时性的要求.

1 需求分析 1.1 功能需求

在虚拟肝脏手术的研究中, 采用体纹理的绘制方法, 不仅可以提高虚拟肝脏模型的真实感表现, 还能够更好的展现切割操作时模型内部切面的纹理细节. 本系统以美国可视化人体数据集(Visible Human Dataset, VHP)为数据源, 提取连续的肝脏横断面合成样本体纹理, 对肝脏模型进行三维重建, 对建立的肝脏体模型进行不同程度的网格消减处理, 得到实验参考模型; 同时对连续的肝脏横断面的轮廓数据进行最优感兴趣区域(Region Of Interest, ROI)的提取, 再进行纵向叠加排列, 最后生成样本体纹理. CUDA加速的肝脏体纹理空间合成需要根据肝脏体模型的大小计算出能够包含该模型的体纹理空间的大小, 从样本体纹理中的随机位置取出特定大小的体纹理对纹理空间进行填充. CUDA加速的肝脏体纹理映射需要将肝脏体纹理从纹理空间映射到模型的表面形成精细的纹理以提高真实感. CUDA加速的肝脏体切割需要采用鼠标模拟手术器械来对虚拟肝脏模型进行切割的操作.

1.2 系统运行环境

系统开发环境是系统集成开发的基础, 主要包括了系统开发的硬件平台、软件平台, 其中软件平台包括了操作系统、集成开发环境以及一些开放的工具库.

虚拟手术中肝脏体实时切割系统基于CUDA架构, 以VS2012为平台, 使用OpenGL编写了虚拟肝脏的3D图形视觉渲染程序. 实验用的PC机采用Windows操作系统并配置了Intel® Core^TM i7-2670QM 2.2 GHz 的CPU, 8 GB内存, NVIDIA GeForce GT 540 MB 的显卡及2 GB显存.

2 系统总体设计

模块化设计是从应用的需求出发, 将系统拆分出不同的功能块, 然后以功能块为单位进行程序设计. 模块化设计的目的是为了降低程序复杂度, 使程序的设计、调试和维护等操作简单化, 同时增加程序的复用性.

虚拟手术中肝脏体实时切割系统在开发上遵循模块化的设计思想. 系统从整体上可分为4大模块, 分别是数据处理模块、交互模块、渲染模块、切割模块, 各大模块又分别由不同的子模块构成. 不同模块负责处理不同的任务, 具有不同的功能. 系统总体设计图如图1所示.

数据处理模块是构建虚拟手术中肝脏体实时切割系统的基础. 该模块负责提取VHP数据集中连续横断面图像, 采用区域的图像分割算法从连续的横断面图像中提取出肝脏轮廓的区域, 然后对肝脏轮廓进行层间插值, 利用VTK (Visualization ToolKit)工具包结合面绘制方法构建出肝脏的三维模型. 交互设计是系统设计中必不可少的环节, 交互模块主要负责处理人机交互, 提高系统的易用性和可控性. 渲染模块是利用虚拟现实和计算机三维成像技术为外科医生呈现出的具有真实感的虚拟肝脏模型. 切割操作是虚拟肝脏手术中最基本的操作之一, 医务工作者需要通过对虚拟肝脏模型的切割来模拟虚拟肝脏手术的进行. 当切割操作完成后, 医务工作者需要实时地观察到切割面的肝脏纹理.

2.1 数据处理模块

数据处理模块包括肝脏轮廓提取子模块和肝脏模型三维重建子模块.

肝脏轮廓提取子模块: 肝脏轮廓提取是指从连续的横断面图像中提取出肝脏轮廓的区域. 三维重建是指利用单视图或多视图建立三维物体在计算机中的数学模型, 它是虚拟现实中进行分析、处理和操作的基础. 在虚拟肝脏手术的研究中, 肝脏模型的三维重建是后续研究工作的基础.

肝脏模型三维重建子模块: 肝脏模型三维重建需要对肝脏轮廓进行层间插值, 利用VTK工具包结合面绘制方法构建出肝脏的三维模型. 在基于样图的二维纹理合成中, 首先需要获取二维的样图作为合成的参考. 在三维空间中, 合成体纹理也需要获取能够反映全局特征的纹理信息的三维样本块, 即样本体纹理. 在本文的研究中, 需要以肝脏的轮廓中提取出的ROI区域为二维纹理样本, 来合成包含丰富纹理信息的三维纹理块.

2.2 交互模块

交互设计是系统设计中必不可少的环节, 交互模块主要负责处理人机交互, 提高系统的易用性和可控性. 交互模块包括鼠标交互子模块、键盘交互模块和GUI参数设置模块.

鼠标交互子模块: 负责接收和处理鼠标事件, 将事件转换为不同的指令传递给不同的模块以执行相应的任务: 传递给渲染模块, 用于改变摄像头的位置以控制观察的角度及距离; 传递给形变计算模块, 用于选中并拉动模型顶点以设置外力; 传递给四面体切割算法模块, 用于控制手术刀模型, 形成分割面.

键盘交互子模块: 负责接收和处理键盘事件, 将事件转换为指令专递给相应的模块. 键盘事件主要针对渲染模块, 不同的按键对应不同的指令, 用于控制仿真场景中的物体, 如地板、场景盒等, 是否进行渲染, 或者控制肝脏模型的渲染属性, 如是否显示包围盒, 网格线等.

GUI参数设置子模块: 主要针对形变仿真的参数, 负责提供设置接口, 使用户能够对这些参数进行交互式地调整.

2.3 渲染模块

渲染模块包括样本体纹理生成子模块、肝脏体纹理合成与映射子模块.

样本体纹理生成子模块: 在基于样图的二维纹理合成中, 首先需要获取二维的样图作为合成的参考. 在三维空间中, 合成体纹理也需要获取能够反映全局特征的纹理信息的三维样本块, 即样本体纹理. 在本文的研究中, 需要以肝脏的轮廓中提取出的ROI区域为二维纹理样本, 来合成包含丰富纹理信息的三维纹理块.

肝脏体纹理合成与映射子模块: 为了将体纹理样本映射到肝脏体模型上形成纹理细节丰富且真实感强的肝脏体纹理, 需要建立精细的肝脏体纹理空间. 本文通过合成一个能够包含肝脏模型大小的体纹理空间, 将体纹理从纹理空间映射到模型空间中, 来生成具有真实感的肝脏模型.

2.4 切割模块

切割模块的主要功能是完成切割操作. 切割操作是虚拟肝脏手术中最基本的交互方式之一, 通过切割不仅需要将模型分割成两个部分, 还要将纹理从肝脏体纹理空间映射到切割面上, 形成真实的切割面纹理. 本系统采用鼠标来模拟虚拟手术器械进行切割操作,切割模块包括碰撞检测子模块和模型网络重构模块.

碰撞检测子模块在切割仿真的不同阶段起着不同的作用: 在切割未执行阶段, 碰撞检测主要用于确定切割操作的起始时间和起始位置; 而在切割操作完成之后, 碰撞检测则用于检测仿真物体与其他场景物体的接触情况, 进而做出相应的处理, 产生必要的碰撞反应.

模型网络重构子模块将用户通过鼠标在屏幕上划动形成连续的切割线量化为若干切割线段, 然后将每一条线段投影到肝脏模型上与模型的表面产生交点, 连接交点形成一个有界平面. 接着通过径向基函数来拟合线段集所产生的有界平面来生成一个新的平面, 即切割面.

3 关键技术 3.1 贝塞尔曲线

贝塞尔曲线^[12]由起始点、锚点(即终止点)和控制点组成, 通过调整控制点, 可以使其形状发生变化. 贝塞尔曲线有以下几种形式, 用于绘制各类二维平面上的曲线形状.

(1)一阶贝塞尔曲线:

$B({{t}}) = (1 - {{t}}){{{p}}_0} + {{t}}{{{p}}_1},{{t}} \in [0,1]$

(1)

其中, B(t)为t时刻下的点的坐标值, p₀为起始点, p₁为终止点.

(2)二阶贝塞尔曲线:

$B\left( {{t}} \right) = {\left( {1 - {{t}}} \right)^2}{{{p}}_0} + 2{{t}}\left( {1 - {{t}}} \right){{{p}}_1} + {{{t}}^2}{{{p}}_2},{{t}} \in \left[ {0,1} \right]$

(2)

二阶贝塞尔曲线由p₀、p₁和p₂三个定点确定, 在二维平面上表现为一条平滑的抛物线.

(3) n阶贝塞尔曲线:

通过低阶贝塞尔曲线的参数形式, 可以推导出基于给定点p₀、p₁、p₂直至p_n 的n阶贝塞尔曲线通式为:

$B({{t}}) = \sum\limits_{{{i}} = 0}^{{n}} {\left( {{}_{{i}}^{{n}}} \right){{{p}}_{{i}}}{{(1 - {{t}})}^{{{n}} - {{i}}}}{{{t}}^{{i}}}} ,{{t}} \in [0,1]$

(3)

本文采用人工确定关键点结合贝塞尔曲线拟合的方式来提取肝脏的轮廓, 为了获得较好的实验效果, 需要找到合适的控制点. 首先沿着肝脏的轮廓间隔地指定拟合的关键点, 关键点的位置尽量贴合肝脏轮廓的边缘, 依次连接轮廓上各关键点组成一个多边形. 连接该多边形各边上的中点, 得到中点的连线段. 在改中点的连线段上找出一关键点, 使得该点分线段的比与相邻两边的长度比相等, 最后平移中点的连线段, 使该线段上的关键点与对应的顶点重合, 从而完成曲线拟合. 曲线拟合示意图如图2所示.

3.2 区域生长法

区域生长法的基本思想是根据像素的相似性生长的原则来划分集合构成分割的区域. 首先在预分割的区域内确定一个种子点, 该点必须是待分割区域内具有代表性的点, 种子点可以是单个像素点, 也可以是相似的像素点区域. 在种子点的邻域内寻找与该点相似的像素, 相似性的判定的原则是计算周边像素点与种子点彩色、梯度、灰度等量间的距离是否在差值范围内, 将符合条件的像素点包含到种子点的区域内, 然后以新区域内的像素点作为新的种子点继续生长, 直到找不到满足条件的像素点为止^[13].

3.3 层间插值法

假设相邻的断层图像为第k–1层和第k+1层, 需要通过插值法计算出第k层的图像轮廓数据. 但由于肝脏的轮廓包含了大量的点, 直接对其进行形状插值将引入大量的计算, 为了简化计算, 通常采用的方法是在原始的轮廓中寻找一个点, 使剩下的点构成的多边形的面积与原来最接近. 以此类推, 分别找到面积差值最小的m个点. 本文采用了等分圆周的方法来简化轮廓点数, 从而获得最逼近原始轮廓的多边形. 首先假设原始轮廓上的点数为N, 经过简化后轮廓上的最终点数为M(M<N), 计算相邻轮廓 ${C_{k - 1}}$ 和 ${C_{k + 1}}$ 的质心 ${\overline m _{k - 1}}$ 、 ${\overline m _{k + 1}}$ :

$ \overline {{m_{k - 1}}} = \left\{ {\begin{aligned} & {{{\overline x }_{k - 1}} = \frac{1}{{{N_{k - 1}}}}\sum\limits_{i = 1}^{{N_{k - 1}}} {{x_i}} }\\ & {{{\overline y }_{k - 1}} = \frac{1}{{{N_{k - 1}}}}\sum\limits_{i = 1}^{{N_{k - 1}}} {{y_i}} } \end{aligned}} \right. $

(4)

$ \overline {{m_{k + 1}}} = \left\{ {\begin{aligned} & {{{\overline x }_{k + 1}} = \frac{1}{{{N_{k + 1}}}}\sum\limits_{i = 1}^{{N_{k + 1}}} {{x_i}} }\\ & {{{\overline y }_{k + 1}} = \frac{1}{{{N_{k + 1}}}}\sum\limits_{i = 1}^{{N_{k + 1}}} {{y_i}} } \end{aligned}} \right. $

(5)

过 ${C_{k + 1}}$ 的质点 ${\overline m _{k + 1}}$ 作一条平行于Z轴的直线, 平移 ${C_{k - 1}}$ 使该直线经过其质点 ${\overline m _{k - 1}}$ , 然后将 ${C_{k - 1}}$ 投影到 ${C_{k + 1}}$ 所在的平面上, 记其投影轮廓为 $C_{k - 1}'$ . 接下来以质点为起点, 任意角度为起始角度, 间隔角度为 ${360^ \circ }/{{M}}$ 向外作射线, 射线与 $C_{k - 1}'$ 和 ${C_{k + 1}}$ 的交点即为新的轮廓点, 只要M值足够大, 就能够使新的轮廓点连接成的多边形尽可能地逼近原始轮廓线, 分别记简化后的轮廓为 $C_{k - 1}^{''}$ 和 $C_{k + 1}'$ , 假设第k–1, k, k+1层Z方向上的坐标为 ${z_{k - 1}}$ , ${z_k}$ 和 ${z_{k + 1}}$ , 则计算 ${C_k}$ 的插值公式为:

${C_k} = \dfrac{{{z_{k + 1}} - {z_k}}}{{{z_{k + 1}} - {z_{k - 1}}}} \cdot C_{k - 1}^{''} + \dfrac{{{z_k} - {z_{k - 1}}}}{{{z_{k + 1}} - {z_{k - 1}}}} \cdot C_{k + 1}'$

(6)

3.4 基于CUDA的并行计算

基于CUDA的并行计算主要分为以下6个步骤, 计算流程如图3所示.

步骤1. 在利用CUDA进行并行计算之前, 需要在GPU上为需要进行运算的数据及运算结果申请相应的显存空间.

步骤2. 在GPU中分配好相应的显存空间后, 将内存地址中需要在GPU上进行运算的数据拷贝至显存空间中.

步骤3. 在利用CUDA进行GPU计算时, 需要启动合适的线程数来覆盖并行计算的规模. 例如进行网格节点着色计算时, 定义Kernel函数名为facePointKernel, 其线程表达式为facePointKernel<<<(faceNum+1023) /1024,1024>>>(…), 其中faceNum为表面网格节点数, <<<(faceNum+1023)/1024,1024>>>表示启动的线程块数为(faceNum+1023)/1024, 每个线程块启动1024条线程, (…)中为需要转移到GPU上计算的数据.

步骤4. 设置Kernel函数即CUDA的核函数. 所有的并行计算都在Kernel中进行, 为了保证计算的并行性, 需要建立线程索引, 使索引的每一个线程都运行一次Kernel函数中的计算.

步骤5. 在Kernel中完成计算后, 需要将计算的结果传送回主机, 将数据从GPU的显存拷贝到CPU的内存中.

步骤6. 在CUDA上完成计算后, 也需要释放动态分配的显存空间, 以免造成泄漏.

4 实现

在虚拟肝脏手术中, 需要利用人体连续横断面中肝脏轮廓的数据来生成肝脏模型以及肝脏的样本体纹理块. 本文采用了美国可视化人体项目的VHP图片集来作为研究的数据来源, 因此对VHP数据集进行肝脏轮廓的提取就成为了研究工作的基础.

对VHP数据集进行肝脏轮廓提取是对一系列图像进行处理的过程, 因此在进行轮廓提取的操作之前, 需要先将数据集导入到系统中. 为了确保数据集能够顺序地在系统中进行处理, 在导入前对数据集采用依顺序编号命名的原则, 使系统在导入数据集的过程中能够按照各图片的文件名进行先后次序的排列. 系统利用OpenCV的cvLoadImage()函数来实现图像数据的读取, 由于VHP数据集的图像数据较大, 包含了许多无用的信息, 为了避免进行过多的处理, 需要先对读取的图像进行裁剪处理,根据肝脏所在区域大小将图像裁剪成分辨率为600×630的图像集. 在如图4所示的系统界面中, 能够查看导入到系统后经过裁剪的每一张肝脏横断面图像.

肝脏轮廓提取的精度直接影响到三维重建后肝脏模型的还原程度, 为了提高图像分割的精度, 本系统采用人工参与肝脏轮廓关键点的指定, 再用贝塞尔曲线拟合曲线的方法来标记出肝脏的轮廓, 由于导入到系统中的图像为连续的肝脏横断面图像, 因此相邻两层之间的轮廓差异性较小, 可以采用当前图像指定的关键点作为下一张图像的起始点, 再进行细微的调整即可. 随着肝脏轮廓的变化, 只需根据需要对关键点的数量进行增加或删除, 便可以方便地标记出VHP数据集中各横断面上的肝脏轮廓. 贝塞尔曲线拟合的结果如图5所示.

通过人工拟定关键点及贝塞尔曲线的拟合,在连续的横断面图像中标记出肝脏区域的轮廓, 然后对肝脏的区域进行分割提取. 图像分割结果如图6所示.

利用分割出的肝脏轮廓数据, 采用层间插值及面绘制的方法对肝脏模型进行三维重建. 层间插值的目的是为了建立更精细的肝脏模型. 根据VHP数据集提供的数据, 其相邻的横断面之间的距离为1 mm, 为了改善断层的过渡情况, 提高模型的精度, 在相邻两层切面之间插入额外的两层切面, 插值面的轮廓形状由相邻的两切面计算得到. 本系统采用基于多边形逼近的形状插值法来进行层间的插值, 首先利用等圆周的方法来简化原轮廓上点的数量, 以获得最大程度逼近原始轮廓面积的多边形, 再利用3.3中的插值公式进行插值轮廓的计算. 在进行层间插值后, 构建模型的切面的总层数为188+187×2=562层.

本系统采用VTK工具包结合Marching Cubes面绘制的方法, 利用562层切面轮廓进行肝脏模型的三维重建. 在尚未进行网格数简化的情况下, 即网格简化率为1时, 未进行网格削减的肝脏表面网格模型如图7所示.

未经过网格简化处理的模型包含了过多的顶点, 这将给后续的映射计算及切割带来大量的数据计算, 造成系统实时性的降低. 为了避免这种情况的发生, 需要对初始模型进行一定程度的网格削减工作. 本系统采用了VTK封装的VTKDecimatePro网格削减类, 通过设置不同的简化率来减少初始模型表面的网格数目, 网格简化率为0.3时所生成的肝脏体模型如图8所示.

5 切割效果与分析

本文实现了肝脏体模型的切割和切面体纹理的快速绘制. 本文分别在GT 540M和GTX 650的GPU以及Intel® Core^TM i7-2670QM 2.2 GHz的CPU上进行了实验, 并选用了顶点数14 663, 三角面片数18 676的肝脏模型作为切割模型, 针对3组不同形状的切割进行了实验数据的记录, 并与CPU下的切面纹理的映射速度进行了对比, 切割效果如图9所示.

由图9可以看出, 利用鼠标模拟虚拟手术刀具在肝脏模型上画出的曲线, 可以将肝脏模型切割成任意形状, 且能够看到切割产生的新平面上的纹理. 表1所示为图9中3组不同切割形状的模型在GPU和CPU下的切割时间对比.

6 结论

本文的主要工作包括肝脏模型的肝脏轮廓提取和三维重建. 在肝脏轮廓提取的工作中, 首先分析了各图像分割算法的利弊, 提出了一种人工干预的曲线拟合方法来标记肝脏区域的轮廓, 接着通过区域生长法提取出该轮廓区域. 然后采用面绘制的方法, 结合层间插值技术, 利用VTK工具包进行肝脏模型的三维重建工作, 建立了肝脏网格表面模型,并引入CUDA技术, 提高运算效率, 为医学诊断和手术规划提供了辅助工具.