CN103699716A - 一种个性化三维医学图像驱动的器官虚拟显示方法 - Google Patents

Abstract

本发明涉及一种个性化三维医学图像驱动的器官虚拟显示方法，包括以下步骤:（1）将用于大尺度形变仿真模拟的拉格朗日显式计算模型进行了扩展；（2）通过用户透明地构建材质敏感的代理物理仿真结构来驱动个性化三维医学图像的非线性形变和切割仿真；（3）对手术仿真框架中的主要处理阶段进行了基于CUDA的全程并行加速，仿真过程可以生理准确的展现器官任意切面的解剖结构。

Description

一种个性化三维医学图像驱动的器官虚拟显示方法

技术领域

在虚拟手术仿真方面，本文以基于非线性有限元的形变和切割仿真、个性化医学数据的直接使用、器官内部解剖结构的真实呈现作为基本目标展开研究：手术仿真环境必须能够提供与真实手术相似的虚拟环境，其中人体软组织及器官在不同体位下的自变形和在手术器械接触、拉伸、切割等条件下的大尺度变形，以及切割过程的物理仿真，是虚拟手术和术中导航研究的重点和难点之一。人体软器官的行为和材质结构很复杂，一般具有不均匀性、各向异性、准不可压缩性、非线性、塑性和粘弹性等性质。目前，对于软组织的变形主要采用基于连续力学的有限元分析法，以保证变形的精确性与真实性。因此，如何基于器官的物理仿真模型，实时精确的仿真器官的大尺度变形及其对切割等交互操作是虚拟手术仿真必须研究的内容。现有的手术仿真方法尚不能直观、有效的传达个性化人体器内部的组织结构、材质组成及其视觉特性等信息。因此，在手术仿真技术换代发展需求的牵引下，如何从个性化的CT、MRI等三维医学影像数据出发，并在有限元仿真模型的驱动下，进行个性化数器官体纹理动态合成与绘制是个性化手术仿真必须解决的问题。此外，虚拟手术的有限元仿真以及仿真结果的逼真呈现需要大量的运算，为了保证仿真的可交互性，如何结合CUDA进行全程并行加速也是本文需要着重研究的一个重要问题。

背景技术

个性化三维医学图像驱动的器官形变和切割仿真相关方面为三维图像的形变操控、基于有限元的形变和切割仿真以及基于CUDA的通用计算三个方面。

三维图像的形变操控

直接以三维图像为对象进行形变操控，可最大限度的描述实体外观及其内部结构的逼真性。许多早期的工作一般都是通过基于图像的形变技术来解决这个问题。在对三维图像任意剖面的内部结构进行探究的需求的驱动下，后来一些研究者相继提出了一些更为复杂的几何形变操控策略。如：Masutani等人通过构建分片的多边形平面几何机构，并以此作为三维图像分片纹理映射的载体，借助多边形顶点的位移控制实现了对三维图像形变的近似模拟；Correa等人则通过在三维图像剖面上散布代理的几何控制点，通过操控这些控制点相对高效的实现了三维图像的形变控制。但是这些方法仍属于基于几何的自由形变的范畴，因而大都缺乏物理真实感。最近，Nakao等人通过将代理点算法进行扩展，通过自动生成四面体代理结构实现了较为精确的三维图像形变的物理仿真。后来，又将优化算法融入其中改进了提升了四面体代理结构的构建精度，进一步提升了仿真的精确性。但是，由于这种四面体代理结构的生成非常依赖于对三维图像内部曲率的检测，并需要辅以一定的后处理，因此，代理结构的构建只能在预处理阶段完成，这导致当三维图像发生大尺度形变时，很难实时的对代理几何结构进行同步优化更新来保证网格结构的质量。此外，该方法也忽略了对切割导致的拓扑结构变更及其自碰撞的处理，因此，它在虚拟手术中的实际应用仍受到很大的限制。

基于有限元的形变和切割仿真

目前有限元对物体形变进行物理建模和力学仿真的有效性已经的得到了广泛的认可。起初，有限元方法大都只能模拟线性形变，只有一小部分基于共旋模型的有限元方法可以处理几何非线性形变的情况。最近，为了对物体的大尺度形变进行仿真，Christian将基于共旋模型的有限元方法进行了扩展，并在采用规则六面体有限元网格的基础上，借助基于多重网格的偏微分方程组加速求解算法实现了线弹性物体大尺度形变的实时仿真，但是从物体材质的本构关系方面来说，该方法仍属于线性有限元的范畴。Miller等人提出的TLED模型对材质非线性进行了严谨的数学建模，并可通过显式的时间积分来对有限元的控制方程进行快速求解，并且基于GPU的并行实现对具有10000个自由度的物体模型已经可以实现实时的形变仿真，同时，它作为一种可变形模型也已经被应用于三维图像非刚性配准领域。后来，通过在TLED中融入对粘弹性的支持，实现了对物体各向异性粘弹性形变模拟的仿真。但是，目前TLED方法在切割仿真方面的扩展还比较欠缺，尚需进行进一步的研究。由于在因切割引起拓扑改变时，需要对有限元网格及其附属物理属性进行动态更新，目前基于有限元方法的物体精确切割仿真还很难满足交互性的要求。其中，当切割发生时，最简单的处理方法是将虚拟工具所接触到单元网格直接删除掉，但这一般会导致仿真效果发生较为严重的失真。基于有限元细分的方法虽然可以较为有效的来应对因切割导致的拓扑改变，但这种方法往往会产生一些病态的仿真单元，从而影响到仿真的稳定性。后来，尽管Wicke等人通过用不规则仿真单元来代替原来的规则仿真单元，对该方法进行了改进，但这种改进需要付出极高的计算代价。与上述方法思路不同的是，Molino等人提出了一种虚拟节点算法来处理应对拓扑更新问题，在切割发生时，他们主要通过对被切割单元进行复制，并同时将被切割单元所剩余的材质成分赋给副本单元，避免了对被切割单元的直接分割处理。然而，由于切割后所得到的每一个片段都要求至少包含原始有限元网格中的一个有效单元，因此这种方法所能允许的切割次数严重受限。最近，基部自适应六面体剖分的有限元方法在线弹性物体切割仿真方面取得了成功应用。此外，通过从三维图像中提取的多分辨率几何模型，Jerabkova等人通过在最精细层次上去除体素单元并对相应几何模型进行同步更新，对物体切割取得了交互式的仿真效率。因此，这些工作为本发明将用户透明地规则单元网格自适应构建与个性化医学图像相结合，进行非线性有限元驱动的三维图像形变和切割仿真带来了很大的启发。

基于CUDA的通用计算

现代GPU具有高度并行化的体系结构设计，可以以比CPU高多个数量级的效率来执行浮点运算。特别地，近几年随着GPU统一计算框架CUDA的出现，一个CUDA核可以同时执行成千上万的轻量级线程，并且在片上缓存、共享缓存以及全局内存等多种GPU存储形式的辅助下，这些线程具有极为灵活的数据访问和同步协调处理能力。因而，各种计算复杂度较高的通用计算任务都开始尝试着借助CUDA来提升计算效率。对CUDA通用计算方法及其应用进行完备综述已经超越了本文的范畴。

发明内容

一种个性化三维医学图像驱动的器官虚拟显示方法，其特征在于包括以下步骤:

（1）将用于大尺度形变仿真模拟的拉格朗日显式计算模型（Total Lagrangian ExplicitDynamic algorithm，TLED）进行了扩展；

（2）通过用户透明地构建材质敏感的代理物理仿真结构来驱动个性化三维医学图像的非线性形变和切割仿真；

（3）对手术仿真框架中的处理阶段进行了基于CUDA的全程并行加速。

附图说明

图1三维医学图像非线性有限元驱动下的形变和切割仿真流程图。

图2基于高斯混合模型的三维图像分割结果及其与三维水平集方法的对比。

图3为用户透明的代理结构自适应生成过程示意图。

图4为规则背景栅格辅助的代理仿真结果拓扑更新和碰撞检测处理。

图5为基于MLS拟合的位移场上采样计算方法示意图。

图6为基于CUDA的动态三维图像体绘制算法的数据流。

图7为局部固定的三维图像在重力作用下的形变仿真结果。

图8为三维图像在重力和地面限制条件共同作用下的形变仿真结果。

图9为局部固定的女性三维图像在重力作用下的多次切割仿真结果。

图10为局部固定的男性三维图像在重力作用下的多次切割仿真结果。

具体实施方式

下面结合附图以及具体实施例进一步说明本发明。

一种个性化三维医学图像驱动的器官形变和切割仿真，具体实施如下：

1.物理仿真代理结构的创建

1.1基于高斯混合模型的三维图像自动分割

针对特定器官来说，为了分析其物理拓扑结构，本发明需要首先对三维医学图像进行分割。众所周知，由于图像颜色的相似性，用户交互较少的自动或半自动分割方法很难实现医学图像的精确分割，鉴于本发明研究内容的着力点在于统一的个性化手术仿真框架，因此，给定个性化的三维医学图像，本发明采用分割效果相对较好的高斯混合模型（GMM）来快速地对三维图像进行自动的多标识分割。

图2给出了多组图像的GMM算法分割结果及其与三维水平集算法的对比，其中图中的第一列给出了原始三维医学图像，第二列给出了GMM算法的多标识分割结果。在此基础上，本发明根据分割结果，通过指定目标器官所属的标识类型，实现了对目标器官的拾取，拾取结果被罗列在图中的第三列中，而图的第四列则是三维水平集方法的分割结果。通过对比本发明可以发现，GMM算法要比三维水平集方法相对精确一些，特别对于具有较弱颜色对比度的图像（如男性心脏三维图像和女性心脏三维图像）来说，这种效果就更为明显。而在分割效率方面，根据本发明的实验统计，对2562*70大小的三维图像来说，GMM一般只需几十秒钟便可完成多标识分割，这完全可以满足个性化数据即时应用的需要。

1.2用户透明代理仿真单元的自适应生成

基于多标识分割结果，本发明通过用户透明地构建一种代理网格结构来来刻画器官的物理属性分布。如图3(A)所示，首先基于自适应八叉树将三维图像所在的工作空间进行层次化的剖分，对当前层次的某个子空间来说，其是否需要继续进行细分主要取决于该子空间所包含的器官材质的分布情况。而图3(A)中用不同颜色标记的每个子空间的材质类型则主要通过该子空间所包含体素的分割标识来进行确定。在材质敏感的空间剖分的基础上，本发明继续将每个处于叶节点的子空间用四面体网格进行进一步的剖分。为了避免相邻四面体单元的边出现交叉现象，如图3(B)所示，本发明定义了两种不同的四面体剖分模板。通过分析八叉树叶节点的的空间邻接关系，本发明可以很容易地确定每个处于叶节点的子空间具体应该采用哪种模板来进行四面体剖分。图3(C)以分割得到的肝脏为例，给出了材质敏感的四面体网格代理结构生成结果的一个实例，其中不同颜色的四面体代表不同的材质类型。

1.3基于TLED的物理仿真模型构建

给定特定器官的密度ρ，杨氏模量E，泊松比ν等材质系数，本发明首先根据四面体单元的体积来计算每个仿真单元的质量，然后，将仿真单元的质量非均匀的分配到其所属的四面体顶点上。同时，本发明可以计算出用于显式积分求解的时间步长△t，计算公式为：

$\left\{\begin{array}{c}\mathrm{\Δt}=\frac{L_e}{c}\\ c=\sqrt{\frac{E(1-v)}{\mathrm{\ρ}(1+v)(1-2v)}}\end{array}\right.,$

其中L_e代表所有四面体的最短边长，而c则被定义为材质的波动传播速度。在TLED模型中，t时刻的形变主要通过相对于0时刻的变形梯度^t ₀X来进行度量，^t ₀X能够很好地反映由非线性形变引起的材质拉伸和旋转程度，它可被定义为：

$X_{\mathrm{ij}}^0{}_t=\frac{{\∂}^t{x}_{i|}}{{\∂}^0{x}_j}.$

基^t ₀X，本发明进而可以分别计算得到右柯西变形张量（Right Cauchy-Green deformationtensor）^t ₀X和格林-拉格朗日应变张量（Green-Lagrange strain tensor）^t ₀E如下：

$C_0^t=X_T^0{}_{t}X_0^t,$

$E_0^t=\frac{1}{2}(C_0^t-I),$

这也是TLED模型所采用的应变度量方式，而第二类皮奥拉-基尔霍夫应力模型（thesecond Piola-Kirchhoff stress）^t ₀S则被用来度量应力的大小：

$S_0^t=\frac{\mathrm{\ρ}_0}{\mathrm{\ρ}_t}\left(X_t^0\right)\left(T_t\right)\left(X_T^t{}_0\right),$

这里^tT代表形变物体单位面积上所受的柯西应力，并且⁰ _tX=(⁰ _tX)^-1，而物体质量密度的变化率⁰ρ/^tρ则受控于：

$\frac{\mathrm{\ρ}_0}{\mathrm{\ρ}_t}=J_t=\det \left(X_0^t\right).$

基于上述定义，TLED将连续形式的平衡公式定义为：

${\∫}_{V_0}S_{\mathrm{ij}}^0{}_t{\mathrm{\δ}}_0^t{E}_{\mathrm{ij}}{d}^{0}V={\∫}_{V_0}f_i^B{}_0{}^t{\mathrm{\δ}}_{u_i}{d}^{0}V+{\∫}_{S_f^0}f_i^S{}_0{}^t{\mathrm{\δ}}_{u_i}^S{d}^{0}S,$

其中^t ₀f^B _i代表体积力，^t ₀f^S _i则是对应于位移δ^S _ui作用在曲面⁰S_f上的表面张力。对于离散形式的四面体仿真单元来说，本发明可以采用下列形函数来在四面体内部对相关物理量进行插值计算。

$\left\{\begin{array}{c}{h}_1=(1-r)(1-s)(1-t)/8\\ h_2=(1+r)(1-s)(1-t)/8\\ h_3=(1+s)(1-t)/4\\ h_4=(1+t)/2\end{array}\right.$

其中h1、h2、h3和h4分别是四面体四个顶点的形函数。这样本发明可以得到标准的有限元平衡公式如下：

$M\stackrel{\·\·}{U}+\mathrm{DU}+K\left(U\right)U=R.$

由于TLED模型采用显式的时间积分来对公式进行求解，因此，有限元刚度矩阵无需再进行动态装配，可直接通过下列公式在仿真单元层面进行计算。

$K\left(U\right)U=\underset{e}{\mathrm{\Σ}}F_t\left(e\right),$

$F_t={\∫}_{V_0}\left(B_L^T{}_0{}^t\right)\left(\hat{S}_0^t\right)d^{0}V,$

这里的t₀^S是第二类皮奥拉-基尔霍夫应力的向量形式表示，^t ₀B_L是应变-位移矩阵，它的第i个子矩阵可基于下列公式通过用变形梯度^t ₀X对₀B_L0进行变换得到。

而₀B⁽ⁱ⁾ _L0则可通过公式所定义的形函数的空间导数计算得到，其计算公式为：

假设本发明已经计算得到了所有节点所受的应力^tF和外力^tR，那么每个仿真单元顶点在当前时刻的位移量可通过下列公式计算得到。

$U_i^{t+\mathrm{\Δt}}=\frac{R_i^t-F_i^t+\frac{2M_{\mathrm{ii}}}{\mathrm{\Δ}{t}^2}\left(U_i^t\right)+(\frac{D_{\mathrm{ii}}}{2\mathrm{\Δt}}-\frac{M_{\mathrm{ii}}}{\mathrm{\Δ}{t}^2})\left(U_i^{t-\mathrm{\Δt}}\right)}{D_{\mathrm{ii}}/\left(2\mathrm{\Δt}\right)+M_{\mathrm{ii}}/\mathrm{\Δ}{t}^2},$

并且，其中的一些项可以被预计算，具体如下：

$A_i=\frac{1}{D_{\mathrm{ii}}/\left(2\mathrm{\Δt}\right)+M_{\mathrm{ii}}/\mathrm{\Δ}{t}^2},$

$B_i=\frac{\frac{2M_{\mathrm{ii}}}{\mathrm{\Δ}{t}^2}}{D_{\mathrm{ii}}/\left(2\mathrm{\Δt}\right)+M_{\mathrm{ii}}/\mathrm{\Δ}{t}^2}=\frac{2M_{\mathrm{ii}}}{\mathrm{\Δ}{t}^2}{A}_i,$

$C_i=\frac{\frac{D_{\mathrm{ii}}}{2\mathrm{\Δt}}-\frac{M_{\mathrm{ii}}}{\mathrm{\Δ}{t}^2}}{D_{\mathrm{ii}}/\left(2\mathrm{\Δt}\right)+M_{\mathrm{ii}}/\mathrm{\Δ}{t}^2}=\frac{2D_{\mathrm{ii}}}{2\mathrm{\Δt}}{A}_i-\frac{B_i}{2}.$

因此，基于上述推导过程，本发明可以为生成的代理仿真结构进行物理属性的绑定，并对其中的一些中间物理量进行预计算。

2.代理结构的物理行为仿真

2.1基于CUDA的形变仿真

给定生成的代理仿真结构，预计算结果Ai、Bi、Ci及积分的时间步长△t，本发明可以通过对顶点位移进行迭代更新来实现变形仿真，因此，公式可被重写为：

$U_i^{t+\mathrm{\Δt}}=A_i(R_i^t-F_i^t)+B_{i}U_i^t+C_{i}U_i^{t-\mathrm{\Δt}}.$

鉴于外力^tR_i的计算比较简单，在每个时间步长中本发明只需根据公式来重新计算^tF_i。从预计算得到的₀B_L0中本发明可以发现，_tF_i实质上间接依赖于t₀^S和^t ₀X。而对于^t ₀^S，本发明可以借助下面的超弹性模型来对其进行计算。

$\hat{S}=\mathrm{\μ}({\mathrm{\δ}}_{\mathrm{ij}}-C_{\mathrm{ij}}^{-1}{}_0{}^t)_0^t+\mathrm{\λ}\left(J_t\right)(J_t-1)\left(C_{\mathrm{ij}}^{-1}{}_0{}^t\right),$

其中μ和λ是Lamé常量，δij是Kronecker’s delta，而^t ₀C和^tJ则分别来自于公式的定义。对于^t ₀X来说，本发明无法直接基于公式对其进行计算，只能通过形函数插值的方式得到：

$X_0^t=u_e^T{}_t\∂H+I,$

其中^tu_e是由四面体单元的顶点位移构成的4×3矩阵，而

是由相应形函数导数构成的4×3矩阵，I是一个3×3的单位矩阵。对体积为⁰V的单个四面体单元来说，公式中的积分项^tF可通过下列公式以离散的形式进行计算。

$F_t=V_0\left(B_L^T{}_0{}^t\right)\left(\hat{S}_0^t\right).$

因此，本发明最终可以得到一个12×1向量来表示^tF，它包含当前四面体单元四个顶点的受力大小。所以，对于被多个四面体共享的顶点i来说，通过对相应的受力情况进行累加，本发明可以得出当前时刻其所受应力^tF_i的大小。根据公式，这种逐顶点的受力计算方式非常适合用CUDA进行并行加速。

2.2基于CUDA的切割和碰撞处理仿真

如图4(A)所示，在切割仿真过程中，由于时间间隔△t_f一般都比较短，因此在这个过程中手术刀的运动方向可认为是不变的，而手术刀扫过的曲面也可以近似地看作是平面，可用一个平行四边形

来进行描述，其中p1p2和p3p4分别代表手术刀在当前时刻和上一时刻所处的位置。给定一个四面体单元，假设pmpl是它的六条棱边中的一条，并假设它与平行四边形

相交于p_i点，本发明可以得出下列几何关系：

$\left\{\begin{array}{c}{p}_i=p_1+r(p_3-p_1)+s(p_2-p_1)\\ p_i=p_m+t(p_l-p_m)\end{array}\right.,$

其中r，s，t代表pi的参数坐标，本发明进而可以得到：

$r(p_3-p_1)+s(p_2-p_1)+t(p_m-p_l)=p_m-p_1.$

假设pa=p3-p1，pb=p2-p1，pc=pm-pl，pd=pm-p1并将公式以矩阵的形式进行重写，交点pi的参数坐标r，s，t可通过下列公式计算得到。

$\left[\begin{array}{c}r\\ s\\ t\end{array}\right]={\left[\begin{array}{ccc}{p}_{a1}& p_{b1}& p_{c1}\\ p_{a2}& p_{b2}& p_{c2}\\ p_{a3}& p_{b3}& p_{c3}\end{array}\right]}^{-1}\left[\begin{array}{c}{p}_{d1}\\ p_{d2}\\ p_{d3}\end{array}\right],$

并且，当且仅当r∈(0，1)，s∈(0，1)，t∈(0，1)同时成立时，边pmpl才被确定为切割边。若将四面体单元的六条边都进行这样的判别计算，本发明可以得到当前四面体所包含的切割边的数量。然后，本发明通过将Forest等人提出的拓扑更新算法基于CUDA进行并行优化，在切割发生时，便可对被切割单元实现动态的细分更新。在此基础上，本发明需要对相关仿真单元的体积0V，质量分布，A_i，B_i，C_i，以及相应的形函数导数进行计算更新。同时，显式积分的时间步长△t也需要根据拓扑改变后的代理结构进行重新计算。此外，由于切割的仿真计算过程可以以逐个单元的形式进行处理，因此，本发明可以使用CUDA来对其进行并行加速，相关算法的伪代码请参见算法2。值得注意的是，为了便于CUDA实现，在实际实现过程中，被切割的原始单元只是被标识为无效单元，而并不需要真正的将其从相应的数据列表中删除。对于自碰撞，主要借助规则的背景栅格来辅助CUDA进行并行处理。在初始化阶段本发明会将工作区域剖分为较为稀疏的空间栅格单元。为了尽早的剔除大量的不可能发生碰撞的四面体，首先根据它们的位置分布将其分配到相应的栅格单元中，并将Ganovelli算法进行基于CUDA的优化，来进行四面相交判断。如图30(C)所示，该算法的核心思想是：假设e是面f0，f1的共享边，而p0，p1代表由这些面扩展得到的半平面，并且W表示由此得到的能够包含当前四面体的楔形体，那么如果其他的四面体单元与W相交的话，就不可能存在一个能够包含e的分割平面。所以，将四面体单元按其所属的背景栅格单元进行分组，本发明便可以栅格单元为单位对其所包含的四面进行并行碰撞检测。如图30(B)所示，发生碰撞的四面体单元已用红色进行了标记。若某个四面体单元被检测到会与其他单元发生碰撞，可利用公式沿其运动方向的逆方向对四面体四个顶点的位移进行后退调整，这相当于在检测到碰撞时，为相关单元施加了一种惩罚力。

$U_i^t=U_i^t+(U_i^{t-\mathrm{\Δt}}-U_i^t)/2.$

3.仿真结果的逼真绘制

3.1基于CUDA的动态三维纹理坐标计算方法

到目前为止，在仿真循环的每个时间周期内，本发明都可以得到当前时刻四面体网格顶点的位移及其形变以后的位置坐标。然而，如图中的细点所示，由四面体网格顶点构成的标量场太过稀疏，根本无法精确的描述形变后的三维图像。因此，为了得到图5中粗点所示的颜色信息可随形变发生变化的稠密标量场，本发明需要执行上采样操作。很明显，若直接对由四面体网格顶点颜色构成的稀疏标量场进行上采样，仍然不可能得到精确的效果。对于稠密标量场中的每一个采样点p，本发明需要首先对四面体网格顶点位移所构成的向量场进行上采样，来得到p的位移向量，然后在此基础上，相对于原始三维图像，计算p点的三维纹理坐标，并据此从原始三维图像中拾取相应的颜色值作为该点的颜色。本发明采用基于移动最小二乘拟合的方法来对稀疏的位移场进行上采样。类似于对自碰撞的处理，本发明也采用同样的规则背景栅格来辅助基于CUDA的局部拟合计算，这样，在局部拟合过程中可极大的提高对邻域体顶点的搜索效率。首先计算出p点所处的栅格单元的索引号，然后以该栅格单元为中心，在它的3×3×3邻域栅格单元所构成的区域内搜索出这些栅格单元所包含的四面体网格顶点的集合Sv。因此，这也保证了相邻的局部拟合会具有一定的重叠性。正如本发明在5.4.2节中所描述的，当切割发生时，原来的一些四面体单元会被标识为无效单元，因为它们可能已经被一些更小、更精准的四面体所替代了，因此在顶点集合Sv的生成过程中，需要把这些无效四面体所包含的顶点剔除掉。

3.2基于CUDA的动态三维图像体绘制方法

个性化虚拟手术的形变和切割仿真效果的逼真展现很大程度上主要依赖于能否对器官切面及其内部解剖结构进行生理准确的绘制，这等价于要求对时变的三维图像进行动态的体绘制。尽管本发明已经可以基于CUDA完全在显存中对三维图像的拓扑形变进行物理仿真，但是传统的体绘制方法仍然需要在每一个仿真周期内重新将时变的三维图像拷贝到内存，然后使用图形API将其作为三维纹理再次加载到GPU来进行体绘制，其中GPU和CPU间的数据交换是影响绘制效率的最大瓶颈。如图32所示，本发明采用像素缓存对象(PBO)作为媒介载体，可直接将形变后的三维图像数据从CUDA缓存映射到纹理缓存，并且这个过程可以直接通过DMA通道来实现，而无需任何的CPU干预。对最终显示结果中的每一个像素，由于以GPU光线投射算法为代表的完全基于图像的体绘制算法，需要进行大量的三维纹理拾取并对拾取的颜色进行混合，对已经被有限元计算占据了大量计算资源的虚拟手术仿真来说，这种方法对仿真结果的绘制很难满足可交互性的需求。因此，需要动态的创建图元来辅助体绘制。Christof等人通过将三维图像所占据的空间分解为一组可分别进行纹理映射的多边形图元序列，实现了三维图像的高效体绘制。如图32所示，首先视点相关的为分解后的每一个切片创建一个六边形图元，其中图32)中以虚线表示的顶点是冗余顶点，这样做的目的是为了使每一个切片图元都可以统一的用六边形来进行表示，从而便于CUDA的处理。然后，为了使这些图元在GPU上的物理存储地址可被直接映射为CUDA缓存，并尽可能的减少图元绘制命令的调用次数，本发明将生成的六边形图元以顶点缓存对象(VBO)的形式进行组织。在这种完全基于GPU的数据传递和绘制机制的保证下，本发明可实现时变三维图像的可交互绘制。

Claims (5)

1.一种个性化三维医学图像驱动的器官虚拟显示方法，其特征在于包括以下步骤:

（1）将用于大尺度形变仿真模拟的拉格朗日显式计算模型（Total Lagrangian ExplicitDynamic algorithm，TLED）进行了扩展；

（2）通过用户透明地构建材质敏感的代理物理仿真结构来驱动个性化三维医学图像的非线性形变和切割仿真；

（3）对手术仿真框架中的处理阶段进行了基于CUDA的全程并行加速。

2.根据权利要求1所述的个性化三维医学图像驱动的器官虚拟显示方法，其特征是：所述虚拟手术仿真计算框架，将虚拟手术仿真各个阶段的处理流程进行了无缝集成，并快速地直接利用个性化的三维医学图像进行手术仿真，而无需任何繁琐的人工几何重建处理。

3.根据权利要求1所述的个性化三维医学图像驱动的器官虚拟显示方法，其特征是：将用于大尺度形变仿真模拟的拉格朗日显式计算模型（Total Lagrangian Explicit Dynamicalgorithm，TLED）进行了扩展，使其能够同步处理切割仿真；具体的扩展过程为：给定特定器官的密度ρ，杨氏模量E，泊松比v以及材质系数，首先根据四面体单元的体积来计算每个仿真单元的质量，然后，将仿真单元的质量非均匀的分配到其所属的四面体顶点上；同时，计算出用于显式积分求解的时间步长△t，计算公式为：

$\left\{\begin{array}{c}\mathrm{\Δt}=\frac{L_e}{c}\\ c=\sqrt{\frac{E(1-v)}{\mathrm{\ρ}(1+v)(1-2v)}}\end{array}\right.,$

其中L_e代表所有四面体的最短边长，而c则被定义为材质的波动传播速度；

在TLED模型中，t时刻的形变通过相对于0时刻的变形梯度^t ₀X来进行度量，^t ₀X能够很好地反映由非线性形变引起的材质拉伸和旋转程度，它被定义为：

$X_{\mathrm{ij}}^0{}_t=\frac{{\∂}^t{x}_i}{{\∂}^0{x}_j}.$

$E_0^t=\frac{1}{2}(C_0^t-I),$

这也是TLED模型所采用的应变度量方式；同时，将自碰撞检测及其响应处理；也纳入到了该扩展模型，并采用并行计算方法来进行加速。

4.根据权利要求1所述个性化三维医学图像驱动的器官虚拟显示方法，其特征是：通过用户透明地构建材质敏感的代理物理仿真结构来驱动个性化三维医学图像的非线性形变和切割仿真，并在此基础上，采用了基于代理结构稀疏位移场上采样的动态三维纹理坐标计算方法，具体过程如下：在仿真循环的每个时间周期内，得到当前时刻四面体网格顶点的位移及其形变以后的位置坐标；对于稠密标量场中的每一个采样点p，需要首先对四面体网格顶点位移所构成的向量场进行上采样，来得到p的位移向量，然后在此基础上，相对于原始三维图像，计算p点的三维纹理坐标，并据此从原始三维图像中拾取相应的颜色值作为该点的颜色；

采用基于移动最小二乘拟合的方法来对稀疏的位移场进行上采样；类似于对自碰撞的处理，采用同样的规则背景栅格来辅助基于CUDA的局部拟合计算，这样，在局部拟合过程中可极大的提高对邻域体顶点的搜索效率；首先计算出p点所处的栅格单元的索引号，然后以该栅格单元为中心，在它的3×3×3邻域栅格单元所构成的区域内搜索出这些栅格单元所包含的四面体网格顶点的集合Sv；因此，这也保证了相邻的局部拟合会具有一定的重叠性；当切割发生时，原来的一些四面体单元会被标识为无效单元，因为它们可能已经被一些更小、更精准的四面体所替代了，因此在顶点集合Sv的生成过程中，需要把这些无效四面体所包含的顶点剔除掉；并且，整个动态三维纹理坐标计算过程可借助CUDA进行并行计算。实现了器官切面及其内部解剖结构的高效逼真绘制，进一步提高了虚拟手术仿真的真实感。

5.根据权利要求1所述的个性化三维医学图像驱动的器官虚拟显示方法，其特征是：对手术仿真框架中的主要处理阶段进行了基于CUDA的全程并行加速，并从GPU数据结构组织，算法流程控制，CPU/GPU同步控制方面进行了重新设计或优化，包括切割仿真、有限元结构的拓扑更新、自碰撞处理、时变三维图像的实时体绘制。

Images