一种基于欧拉流体模拟算法的心脏血液流动示意显示方法

摘要

本发明提供一种基于欧拉流体模拟算法的心脏血液流动示意显示方法，包括了三个步骤：心脏血液模型初始化阶段，预处理心血管的表面模型，剖分出欧拉网格，并初始化欧拉算法的时间边界；血液流体模拟阶段，实现实时的欧拉流体模拟算法，依照流体的物理模型仿真心血管内血液流动的情况；血液流场示意性显示阶段，在心血管空间内对流体进行分析，利用简洁明了的符号示意性地表示流场的特征。本发明完全基于GPU来模拟心血管内血液的流动，并可对血液流场进行分析与示意性显示，具有基于物理真实模型、实时性好的特点。

说明书

技术领域

本发明涉及一种基于欧拉流体模拟算法的心脏血液流动示意显示方法。

背景技术

随着计算机科学和虚拟现实技术的不断发展，如今计算机已经被用来模拟各种各样的场 景。这些场景的模拟有着巨大的应用前景，例如协助其他科研领域的实验、辅助工业设计与 测试、在电影和游戏中制作特效、或者在实践教学过程中提供生动的可交互环境等等。在现 代医学领域，计算机仿真技术的应用需求非常巨大，尤其是在人体器官的物理建模、远程会 诊、虚拟手术等方面。这既是计算机仿真技术与医学相结合共同发展的机遇与动力，也是对 计算机仿真技术的挑战。目前计算机仿真技术已经初步应用于虚拟手术训练、远程会诊、手 术规划及导航、远程协作手术等方向。在相关领域中，众多学者已经开展了大量研究。

在医学影像领域中，血液流动的模拟是十分重要的问题，但同时它又是比较复杂的问题。 在图形学方面，流体的仿真与绘制得到了广泛的研究。近些年来，基于物理的流体的仿真技 术得到了越来越多的关注。目前，基于物理的流体仿真方法大致分为三种：基于粒子系统的 拉格朗日（Lagrange）方法、基于网格的欧拉（Euler）方法和基于微观模型和力学方程的晶 格玻尔兹曼（Lattice Boltzmann）方法。

本发明通过对欧拉流体模拟方法进行研究和实现，在此基础上设计一种示意性的显示方 法表达流场的特征，最终完成了一个实时的心脏血液流动过程示意显示工具。这种对血液流 动的模拟显示方式像三维动态的教科书插图一样，可以清晰地表示出心血管区域内流场的物 理性质，方便医生之间更为直观地表达和交流，也给医生与病患者之间良好的沟通铺平了道 路，还可以用来进行医学专业和公共卫生方面的教学。

发明内容

本发明解决的技术问题是：克服了现有血液流动显示工具的不严谨，提供了一种基于流 体物理模型的心脏血液流动示意显示工具，并利用GPU的高度并行性能对物理模型的仿真 与分析进行加速，满足了实时显示的需求。

本发明采用的技术方案为：一种基于欧拉流体模拟算法的心脏血液流动示意显示方法， 包括了以下三个步骤：

步骤（1）心脏血液模型初始化阶段：预处理心血管的表面三角面片模型，剖分出欧拉 算法需要的三维交错网格，并将网格建立在索引数据结构之上，在保证存取性能的基础上优 化内存占用，并在网格中初始化欧拉算法的时间边界；

步骤（2）血液流体模拟阶段：在三维交错网格中利用欧拉流体模拟算法，依照流体的 物理模型迭代仿真心血管内血液的流动情况，并将仿真求解步骤中涉及的网格存储于显存 中，以GPU的高度并行性能达到实时的模拟；

步骤（3）血液流场示意性显示阶段：根据步骤（2）中计算得到的血液流动情况，利用 八叉树的数据结构模型在心血管空间内对流体进行分析，利用简洁的符号示意性地表示流场 的特征，并对显示结果进行绘制。

本发明的原理在于：

（1）通过欧拉流体模拟算法，将流体状态映射到流场，在GPU上实现速度对流、浓度 扩散、压强投影等步骤，可以得到基于物理的流体仿真结果。基于物理的方法有计算量大、 仿真速度慢的特点，本发明将该算法移植到GPU，利用其高度并行的性能，并进行数据存取 和计算的优化，达到实时的时间要求。

（2）为了达到示意显示的目的，本发明对于流体仿真的结果进行分析。利用八叉树数 据结构模型自适应的特性，可以在不同粒度下提取流场的特征，并且一定程度上在提高运算 速度的基础上保证细节、过滤噪声。

本发明与现有技术相比的优点在于：

1、本发明提出了基于欧拉流体模拟算法的心脏血液流动示意显示方法。一方面显示的 血流是通过流体物理模型仿真得到，可以保证其物理上的真实性；另一方面由于采用了GPU 来进行并行计算加速，算法的运算效率提高，可以进行实时的显示。

2、本发明提出的流场分析方法利用八叉树的模型，在表面附近增加细节，在其他区域 采用较低分辨率即可，此举能在分析结果正确的前提下减少运算、提升速度、获得更多细节。

附图说明

图1为基于欧拉流体模拟算法的心脏血液流动示意显示方法的处理流程图；

图2为欧拉流体模拟算法过程示意图；

图3为欧拉流体模拟算法的迭代处理流程；

图4为三维交错网格示意图；

图5为所采用的心血管模型；

图6为血液仿真矢量场（图中为速度场）分析结果的展示；

图7为血液仿真标量场（图中为浓度场）分析结果的展示；

图8为本发明的操作界面。

具体实施方式

图1给出了本发明的总体处理流程，图8展示了本发明的操作界面，下面结合其他附图 及具体实施方式进一步说明本发明。

本发明提供一种基于欧拉流体模拟算法的心脏血液流动示意显示方法，主要步骤介绍如 下：

1.心脏血液模型初始化

心脏血液模型初始化主要包括三部分：心血管模型的引入、网格的优化、数据的调度。

1.1、心血管模型的引入

心血管模型引入的要点在于模型向网格的转化。输入数据是心血管内表面的三角面片模 型（图5），而执行欧拉流体模拟算法需要交错网格（图4）的支持。本发明设计了一个剖分 算法，将心血管内表面三角面片模型转化为分辨率为128³的交错网格，使其适应欧拉流体模 拟算法。本发明设计的模型剖分的算法具体内容如下：

读入所有三角面片数据，找到三角面片模型的AABB包围盒（Axis-Aligned Bounding  Box，AABB）；

（1）将三角面片模型的AABB包围盒划分为128³的网格G₀，此时任一个三角面片中的 任一顶点必定会落在G₀的一个网格单元中；

（2）建立一个三角面片的栈，全部三角面片入栈；

（3）弹出栈顶三角面片tr，将tr的三个顶点落在G₀中的网格单元分别标记为边界；

（4）若tr的三个顶点没有落在G₀同一个网格单元，则将tr沿其最长边的中线划分成两 个三角面片tr₁和tr₂，将tr₁和tr₂压入栈；

（5）若栈非空，执行（4），否则执行（7）；

（6）将G₀中坐标为（0,0,0）的网格单元标记为空单元；

（7）对在G₀中标记为空单元的网格单元执行Flood Fill算法（即每一个与空单元相邻的 没有标记的单元都被标记为空单元，反复执行此步骤直到空单元的数量不再增加）；

（8）将此时没有标记的网格单元标记为内部单元。

1.2、网格的优化

欧拉法的缺点之一是数据量巨大。由于网格结构的原因，网格中每一个单元，不论其中 是否存在流体，都要占用一定的存储空间。这种情况导致了大量无用的空间被占据，也凭空 浪费了很多存取带宽和CUDA线程。实际上本发明所用的心血管模型只占用了1283分辨率 网格的一小部分，大多数网格单元都被标记为“空单元”，即网格是稀疏的。为了满足实时 的要求，需要对网格占用的存储空间进行优化。本发明用改进的数据结构重新组织网格，目 的是在不影响访问效率的情况下减少其对存储空间的需求。

首先在交错网格中统计各个属性需要分配的数量。此时只统计欧拉流体模拟算法需要更 新的属性数量，而忽略存在于“空单元”内的属性，因为欧拉流体模拟算法不会更新“空单 元”，为“空单元”内的属性分配空间是一种浪费。然后，按照各个属性所需的数量为其分 配连续的存储空间。最后，按照“非空单元”的坐标与其属性的一一对应关系建立一个索引 函数。此函数的输入是一空间坐标，输出是该坐标属性所在地址，该函数可以在正比于常数 的时间内给出结果。向函数中输入“空单元”的坐标会得到空地址。

用上述方法代替传统的三维数组来表示网格可以大幅度削减其占用的空间大小，而索引 结构能够在常数的时间内存取，达到了在不影响访问效率的前提下减少网格对存储空间需求 的要求。存储空间的减少意味着数据传输与缓存效率的提高，同时也会在CUDA启动Kernel 函数时节约大量的线程（此时只需为“非空网格”分配线程），加快程序的执行速度。

1.3、数据的调度

本发明实现的工具需要在CPU和GPU端同时运行。其中CPU逻辑能力运算强，负责 程序的组织、心血管模型的引入以及流场示意显示中速度场的分析工作；GPU的并行计算能 力强，负责欧拉流体模拟算法的执行、流场示意显示中浓度场的更新和最终绘制。

在CPU和GPU端同时运行的程序存在的问题是，CPU的寻址空间位于内存，GPU的 寻址空间位于显存，CPU端的程序不能直接访问到显存，GPU端的Kernel函数也不能访问 内存。内存与显存中的数据传输不方便，需要通过调用CUDA提供的接口实现，速度也相 对较慢。因此需要在各个已经实现的算法之间穿插数据调度步骤，实现CPU与GPU端的数 据同步。

本发明实现的欧拉流体模拟算法有一个特点，其运行所需的所有数据都保存在显存中， 在CPU不加干预的情况下，该算法可以持续执行而不需要内存与显存的数据交换。因此只 需在算法执行开始之前对显存中的数据进行初始化，便再也不需要因为该算法执行数据同 步。

在流场示意显示的部分中，对浓度场的更新与欧拉流体模拟算法类似，同时也保持了无 需数据同步的特性。对速度场的特征分析是在CPU中执行的，需要从显存中将八叉树叶结 点上最新的速度数据同步到内存。由于速度场显示的数据量不大，分析阶段完成后，在CPU 端调用OpenGL的接口进行绘制。

因此本发明的数据调度流程主要是：程序初始化阶段在内存中引入模型、分配空间、初 始化流场数据，然后一次性拷贝到显存中；在之后的每一个时间步中，将速度场的数据从显 存同步到内存，以便分析。

2.血液流体模拟

在欧拉法中，存储各种流体属性的空间网格之间的关联度非常低，每个网格单元基本只 与相邻的单元有关系，这种结构非常适合大规模的并行。下面介绍本发明根据此特性将欧拉 法在GPU上并行实现的细节。

本发明对空间采用了分辨率为128³的网格，即在空间中共有2097152个观察点，这对于 N-S方程的求解是一个非常巨大的计算量。在GPU上实现欧拉流体模拟算法的主要目的是 为了利用GPU的高度并行性，加快计算过程。通过NVIDIA提供的CUDA平台，我们可以 编写kernel函数供显卡并行执行。不幸的是，在GPU中执行的代码只能访问到显存空间， 不能访问到内存空间，而显存和内存的相互交流需要在CUDA的kernel函数外进行，并且 由于受制于带宽限制，大规模的数据传输会占用比较多的时间，产生延迟。因此本发明先将 所需的空间在显卡上分配出来，初始的情况在内存组织好后一次性拷进显存，随后算法运行 期间不再进行内存与显卡之间的数据通信。这样可以极大地提高算法的执行速度。

欧拉流体模拟算法每一个时间步中的步骤主要是计算外力项、计算速度对流项、计算粘 性扩散项、执行投影（如图2）。在GPU中实现时，需要加入一些对GPU的控制步骤。算法 实现的整体流程如图3所示。

2.1、速度对流

在GPU上实现速度对流的步骤为：

（1）在显存中建立一个临时的速度场；

（2）为每一个网格单元分配一个线程，该线程的任务是在原速度场中“回退”找到该 网格单元下一时刻速度值的来源，将来源处的速度值填入临时速度场对应位置（隐式积分）；

（3）在所有线程结束以后，用临时速度场更新原速度场。

2.2、粘度扩散

在GPU上实现粘性扩散的步骤为：

（1）在显存中建立一个临时的速度场；

（2）为每一个网格单元分配一个线程，该线程利用该网格单元周围一邻域内的速度值 来计算自己的新速度，并填入临时速度场对应位置；

（3）在所有线程结束以后，交换临时速度场和原速度场；

反复步骤（2）、（3）若干次，直到速度场收敛到一个较为稳定的情况。

2.3、压强投影

在GPU上实现压强投影的步骤为：

（1）利用当前的速度场求出速度散度场；

（2）在显存中建立一个原压强场和一个临时压强场，将原压强场中所有位置的压强置 为0；

（3）为每一个网格单元分配一个线程，该线程利用该网格单元周围一邻域内的压强值 来计算自己的新压强，并填入临时压强场对应位置；

（4）在所有线程结束以后，交换临时压强场和原压强场；

（5）反复步骤（3）、（4）若干次，直到压强场收敛到一个较为稳定的情况；

（6）为每一个网格单元分配一个线程，该线程从速度场中减去当前压强场的梯度，获 得压强投影后的速度场；

（7）等待所有线程结束。

3.血液流场示意性显示

3.1、矢量场的示意性显示

本发明在八叉树的数据结构上自行设计了矢量场的显示方法。矢量场的显示方法要求是 拓扑无关的、自适应的、实时的。八叉树方法是与拓扑结构信息无关的方法。这种方法将网 格单元按照位置临近、属性相似的标准进行划分，然后利用划分后的结果进行示意显示。经 过一系列的优化，该方法可以在效率上达到实时。下面介绍本发明提出的矢量场显示方法的 实现（以速度场为例）。

建立与网格相对应的八叉树，其根结点表示整个网格，根结点的八个子结点分别表示将 整个网格切分成的八个体积相等的立方体，依此递归下去。由于网格的分辨率是1283，因此 八叉树的最大深度为7。每一个结点都包含其对应的立方体中各个网格单元的速度信息。对 于每一个结点是否该继续向下生长，制定如下标准：

（1）如果当前结点深度已达到7，则停止继续向下生长，否则参考（2）；

（2）心血管模型并不是充满整个网格的，只是占据了网格的一部分。如果当前结点中 空网格单元所占的比例高于一定的阈值，则继续向下生长，否则参考（3）；

（3）制定一个衡量结点所对应的立方体中各个网格单元速度相似度的算法，并手动设 定一个与当前深度有关的阈值。如果当前结点所对应的立方体内速度大致相同，则停止继续 向下生长；反之，若当前结点所对应的立方体内速度差异大于阈值，则继续向下生长。

八叉树从根结点开始按照如上三条标准生长。最终其所有的叶子结点互不相交，完全覆 盖住心血管模型，并且其中的每一个都代表一个速度近似的立方体区域。若用箭头表示速度， 则用每一叶子结点中网格单元的位置插值出箭头的位置，用网格单元的速度插值出箭头的大 小和方向，在模型的相应位置画出箭头即可。效果如图6所示。

3.2、标量场的示意性显示

本发明的标量场（如浓度场）显示方法是将空间网格上的标量场数据映射为颜色，调用 OpenGL的接口绘制出来。效果如图7所示。

4.工具运行结果分析

4.1、运行效果

图5从两个不同的角度展示了本发明所使用的心血管模型。该模型来源于真实心血管的 CT扫描数据，输入数据形式是心血管内表面的三角面片。

图6在心血管模型中示意显示了欧拉流体模拟算法产生的速度场。其中黄色的箭头表示 速度场八叉树叶结点上的速度。从图中可以方便地观察心血管模型中流体的速度场。在特征 粒度比较粗的结果中，大量的网格单元聚成少量的箭头，在大体上描述了速度场，而忽略掉 了细节信息；在特征粒度比较细的结果中，八叉树向下细化层次较深，箭头数量较多，直观 地描述了速度场的细节信息。

图7展示了浓度场在模型中随时间的发展变化。图中用红色表示浓度的大小，灰色表示 浓度为0的区域。根据这幅图所示，浓度从入口进入血管，逐渐随着流体运动到整个模型中。

4.2、时间效率分析

表1列出了本发明所实现工具的CPU与GPU版本在仅执行流体模拟功能时的效率对比， 其中帧数的单位为帧/秒（fps）。从表中可以看出，相比CPU的串行计算，GPU强大的并行 计算能力能使流体模拟的速度至少提高一个数量级。尤其是在网格分辨率较高时（如128³）， CPU执行流体模拟的速度很难保证实时，而GPU执行流体模拟的速度比CPU版本快40倍 以上，仍然能维持在实时水平。下文讨论的工具执行速度均为GPU版本的速度。

表1CPU与GPU仅执行流体模拟功能效率对比

表2工具各种功能速度统计

表2列出了本发明所实现工具的各种功能执行速度统计。从表中可以看出，本发明所实 现的工具在各项功能中都达到了实时的速度要求（不低于24帧/秒）。

对于欧拉流体模拟算法，每个步骤完全在GPU上实现，充分利用了GPU强大的并行计 算的能力。在代码编写的过程中，我们十分关注算法速度的问题，尽可能地进行优化，使得 数据的访问次数和算法计算量最大程度地减小，并在一些不显著的地方进行了近似处理。实 验证明本发明实现的欧拉流体模拟算法速度非常快，帧数能够普遍达到90帧/秒以上，相当 于每一个时间步的更新只需要不到11毫秒的时间。

浓度场的显示速度在90帧/秒左右。除了欧拉流体模拟算法的求解以外，主要是浓度场 的更新过程占用了时间。本发明的浓度场更新过程与欧拉流体模拟算法十分类似，也是利用 GPU强大的并行能力求解网格上关于浓度离散形式的微分方程。另外，对于浓度场的绘制， 本发明是将浓度映射成模型表面的颜色，再利用CUDA与OpenGL互操作的措施批量进行 绘制。这个数据组织的过程也需要耗费一定的时间。

速度场的显示速度在27帧/秒左右。此功能相对较慢，原因是速度场的显示需要显存到 内存的数据同步，需要CPU根据速度场数据维护八叉树结构，还需要CPU遍历八叉树的叶 结点进行绘制。因此速度场的分析和显示过程相比其他基于GPU的功能来说有些耗时。

对浓度场和速度场的同时显示速度在27帧/秒左右，能够达到实时的速度要求。虽然速 度场的分析和浓度场的更新都是耗时的工作，但是由于这两项工作分别在CPU和GPU上同 时完成，彼此互不影响，因此并行之后并没有带来过多的时间增加。