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

Abstract

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

Description

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

技术领域

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

背景技术

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

在医学影像领域中，血液流动的模拟是十分重要的问题，但同时它又是比较复杂的问题。在图形学方面，流体的仿真与绘制得到了广泛的研究。近些年来，基于物理的流体的仿真技术得到了越来越多的关注。目前，基于物理的流体仿真方法大致分为三种：基于粒子系统的拉格朗日（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 BoundingBox，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上同时完成，彼此互不影响，因此并行之后并没有带来过多的时间增加。

Claims (4)

1.一种基于欧拉流体模拟算法的心脏血液流动示意显示方法，包括以下三个步骤：

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

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

步骤（3）、血液流场示意性显示阶段：根据步骤（2）中计算得到的血液流动情况，利用八叉树的数据结构模型在心血管空间内对流体进行分析，利用符号示意性地表示流场的主要特征，并对显示结果进行绘制,所述符号包括箭头或颜色指示,所述主要特征包括速度场、浓度场。

2.根据权利要求1所述的基于欧拉流体模拟算法的心脏血液流动示意显示方法，其特征在于：步骤（1）中所述的初始化阶段，将依照模型创建的三维交错网格中有效数据存储在索引表中，而对无效数据不进行存储；在后续的数据处理与计算中，只需对每一索引项分配显存、线程资源，省去了无效数据的资源占用；在初始化阶段将模型数据一次性调入显存，后续的模拟过程中不再进行内存和显存之间的数据调度，只在需要分析流体时进行数据调度，这保证了高效的运算；

依照模型创建三维交错网格的过程为：

读入所有三角面片数据，找到三角面片模型的AABB包围盒（Axis-Aligned BoundingBox，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）将此时没有标记的网格单元标记为内部单元；

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

3.根据权利要求1所述的基于欧拉流体模拟算法的心脏血液流动示意显示方法，其特征在于：步骤（2）中所述的血液流体模拟阶段，是在三维交错网格中实现的，速度对流步骤模拟了流体按照其速度场的运动，粘性扩散步骤模拟了流体的粘性所导致的运动，压强投影步骤模拟了流体的不可压缩性；利用双缓存技术，每一个步骤利用一块缓存作为数据源，另一块缓存作为输出，步骤结束以后交换两块缓存，从而简化数据调度流程，整个算法完全基于GPU实现，并达到性能上的实时；

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

4.根据权利要求1所述的基于欧拉流体模拟算法的心脏血液流动示意显示方法，其特征在于：步骤（3）中所述的血液流场示意性显示阶段基于八叉树的数据结构，利用其自适应的特点可以寻找到流场的特征，八叉树根据流场中的近似度自适应地生长，在近似度高的位置深度较小，描述低频的趋势，在近似度低的位置深度较大，描述高频的细节；通过八叉树自适应的生长过程，其叶节点将流场划分成若干互不相交的子区域，对于每一个子区域里的数据进行显示，所述显示方式包括箭头描述或颜色映射；

八叉树自适应生长的的标准为：

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

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

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

Images