基于改进Marching Cubes算法的血流实时仿真研究

王　娜，陈国栋，陈　怡

(1.福建师范大学福清分校 电子与信息工程学院，福建 福清 350300；2.福州大学 物理与信息工程学院，福建 福州 350108)



基于改进Marching Cubes算法的血流实时仿真研究

王娜1，陈国栋2*，陈怡2

(1.福建师范大学福清分校 电子与信息工程学院，福建 福清 350300；2.福州大学 物理与信息工程学院，福建 福州 350108)

针对传统的Marching Cubes算法空体元检测时间过多影响执行效率的问题，设计了一种针对流体表面绘制的Marching Cubes改进算法。在算法中，首先检测了规则点阵的密度，然后通过设定阈值将粒子密度低于阈值的区域与密度高于阈值的区域分离，仅将密度较高的区域使用简化版Marching Cubes算法绘制。仿真实验证明，与球形渲染算法和Marching Cubes算法相比，本文提出的算法减少了对空体元的访问，提高了显示的质量，从而使得整体绘制算法符合实时渲染的要求。

Marching Cubes(MC)算法; 体元; 血流; 等值面; 阈值

血流仿真作为虚拟手术系统中的一个重要组成部分，随着虚拟手术的进行而实时产生，仿真的效果将直接决定虚拟手术系统的整体真实感。作为虚拟手术系统中的关键技术，血液流动仿真也得到了快速的发展。总体而言，血液流动仿真作为流体仿真领域中的一员，经历了早期的基于构造的仿真和后期的基于物理的仿真。基于构造的血流仿真计算量小，适合早期计算机硬件水平比较落后的情况，但其突出的问题在于真实感不强。因此，对血流仿真领域的研究已经渐渐地转向使用基于物理的方法。该方法相比基于构造的仿真方法而言，真实感比较强。

近年来，国内外众多学者一直在努力对于基于物理的血流仿真领域的研究。文[1]结合立方插值拟质点的方法来模拟血液流动，该方法可以使用较为稀疏的网格对流体进行精确的模拟，然而当网格出现大的变形时会造成网格扭曲而导致精度降低。文[2]将弹簧模型与SPH(光滑粒子流体动力学)相结合实现了血液与血管的相互作用。文[3]利用流体动力学工具结合3D技术对动脉血管中的血流进行了仿真。文[4]、[5]采用稳定的半拉格朗日方法对血液流动进行模拟。文[6]提出了一种结合力反馈机制的渗血模拟仿真模型，该模型对传统的拉格朗日法进行了简化。

综上可知，血液流动现象仿真的主要问题集中于实时性与真实感的矛盾上。因此，权衡真实感与实时性这对矛盾将是血流模拟中必须考虑的问题。

1　粒子位置信息的预处理

为了提高计算的实时性，需要将粒子的位置信息进行预处理，具体做法为构建一个可以包围所有粒子的长方体，并将长方体均匀分割成大量的立方体，根据粒子的位置对其近邻的立方体顶点进行标记，最终构建出一个规则的三维点阵。

由于采用直接计算粒子与立方体顶点之间的距离寻找与粒子最近的立方体顶点的方法计算量大，本文通过粒子的三个维度的坐标与立方体边长a之间进行比较，找到最近邻的立方体顶点，具体算法如下：

由于长方体分割为立方体时是均匀分割，因此，可以计算出与xi最近的立方体顶点的x轴坐标xa与xb，其中xa

(1)

(2)

(3)

(4)

图1　寻找最近邻的立方体顶点

2　Marching Cubes算法

在三维空间数据场中构造等值面有很多种方法，其中最有代表性的就是Marching Cubes算法。Marching Cubes算法简称为MC算法，是W.E.Lorenson与H.E.Cline在1987年的时候提出的[7]。由于该算法具有容易实现、通用性强的优点，因此得到较为广泛的应用。

2.1Marching Cubes算法的主要思想

Marching Cubes算法的主要思想是把离散的三维规则数据场构建成一系列相邻的立方体，采用线性插值计算出等值面与立方体棱边的交点，将交点按一定方式连接成几个相邻的三角面片，作为等值面在局部体元内的一个逼近表示。

Marching Cubes算法的流程如下：

(1) 将三维规则数据场读入内存中；

(2) 根据三维规则数据场，从前往后，从左往右构建体元。每个体元的8个角点分成两层，分别取自上下两层的规则数据场；

(3) 根据等值面所采用的阈值，将体元的每个角点分成两类，其一为当体元角点的函数值大于阈值时，将该角点标记为黑色，其二为当体元角点的函数值小于阈值时，对角点不进行标记。对一个体元的所有角点进行遍历，即可得到该体元的状态表；

(4) 根据状态表，即可得到与等值面相交的体元边界；

(5) 利用线性插值的方法计算出体元与等值面的交点位置；

(6) 利用中心差分算法求得体元各角点位置的法向，再通过线性插值计算三角形各个顶点位置的法向量；

(7) 遍历所有层的数据；

(8) 根据三角面片的顶点坐标和该位置上的法向量绘制出等值面图象。

2.2Marching Cubes算法存在的缺陷

传统的Marching Cubes算法需要通过在相邻切片间构建体元，每个体元上共有8个角点，因此每个体元的角点状态共有28=256种，其相对应的等值面分布也有256种。但是考虑到旋转对称和反转对称的因素，因此可以将等值面的种类缩小为15种[8]。如果对每个体元进行处理，根据顶点位置与等值面阈值比对构造出关系表，再由线性插值得到多个三角形，将得到的多个三角形网格形成等值面。一个医学序列图像所包含的三角形切片有可能数以百万，进行等值面绘制的操作对计算机硬件的要求比较高，由于数据量太大，耗费时间长，处理代价高，实时性差。研究表明，在一般情况下，真正与等值面相交的体元只占总数量的很小部分(至多10%)，算法执行过程中30%～70%的时间是花费在空体元的检测上。

3　本文所采用的方法

本文所使用的方法是一种基于Marching Cubes和球形绘制的混合算法，使Marching Cubes算法中检测的空体元数量减少，从而提高算法的效率。算法主要思想为，将预处理得到的立方体顶点序列进行分层检测，从x，y，z三个维度进行检测，根据设定的阈值，分离出标记为“1”的顶点密度大与密度小的区域，并针对这两种区域进行独立处理，密度小的区域中标记为“1”的顶点直接使用圆球进行绘制，对于密度大的区域使用Marching Cubes算法进行绘制，在绘制过程中，通过简化等值面与立方体边界角点的计算公式以达到降低计算量的目的。

(1) 区域分离

先从x方向进行区域分离。首先设置阈值Nmin，在后面的检测中，将以该阈值作为x方向分割的依据。x方向的区域分离分成两个步骤，即x轴坐标从最小值增加的方向和x轴坐标从最大值减少的方向。

首先从x增加方向进行检测与分离。让x的值从最小值xmin开始增加，每增加一层，x的值增加Δx，计算每层中被标注为1的立方体顶点数n。如果n小于被设定的阈值Nmin，那么将该层分割出来，作为球形绘制的区域。当x增加到某个值xmin+k·Δx时该层标注为1的立方体顶点数n的值不小于阈值Nmin，这时停止该方向的检测。

然后，从x减少的方向进行检测与分离。让x的值从最大值xmax开始减少，每减少一层，x的值减少Δx，计算每层中被标注为1的立方体顶点数n。如果n小于被设定的阈值Nmin，那么将该层分割出来，作为球形绘制的区域。当x减少到某个值xmax-l·Δx时该层标注为1的立方体顶点数n的值不小于阈值Nmin，这时停止该方向的检测。

对每次得到的长方体经过三个维度上的分离，每次都要设置新的阈值，按上述步骤进行处理，得到的分离结果如图2所示，将原长方体划分成以下七个长方体。其中前六个长方体中的点将使用圆球直接对其进行绘制。长方体九中的点将使用Marching Cubes算法进行绘制。

图2　分离结果

(2)使用球形进行绘制

用OpenGL绘制球形只要直接调用GLUT工具包中的glutSolidSphere()函数即可。该函数的原型为void glutSolidSphere(GLdouble radius , GLint slices , GLint stacks)。其中，参数radius代表球形的半径，slices代表以竖直轴上线段为直径分布的圆周线的条数，stacks代表围绕在竖直轴周围的线的条数。slices与stacks这两个参数决定了圆球绘制的精度，这两个值越大那么绘制出来的圆球越精细，但同时带来的是计算量的上升。使用glutSolidSphere( )函数绘制圆球时，需要将其球心的位置平移至要绘制的位置，使用glTranslatef( )函数进行绘制位置的平移。

(3)构建体元

在(1)中，将长方体进行了一系列的分离，最终确定了使用Marching Cubes算法进行绘制的区域，即(1)中所得到的长方体九，将长方体中的立方体当成体元即可。

(4)角点的标记

将粒子预处理后转化为立方体顶点点阵，在立方体顶点中添加一位标志位。立方体顶点的标志位为“1”时，该角点标记为黑色；立方体顶点的标志位为“0”时，该角点不标记。

(5)确定与等值面相交的体元边界

由2.2节可知，等值面的种类可缩小为15种[8]。根据这十五种等值面分布类型可以构建出一个查找表，查找表的长度为256，记录着所有情况下等值面的连接方式，因此根据体元的角点状态就可以从查找表中寻找到与等值面相交的体元边界以及等值面的连接方式。

(6)计算体元和等值面交点位置

Marching Cubes算法中，计算体元与等值面的交点位置是其计算量的重要组成部分之一，因此这一步的计算量将很大程度上影响整体的绘制速度。本文选用计算量较小的中点插值计算公式：

(5)

其中，P代表等值点坐标，P1、P2代表两个端点的坐标。

(7)计算体元和等值面交点处的法向量

与上一步的计算公式相对应，本文计算等值点的法向量也是采用中点插值的计算公式：

(6)

其中，N代表等值点处的法向量，N1、N2代表两个端点处的法向量。

(8)遍历所有层的数据

(9)根据三角面片的顶点坐标和该位置上的法向量绘制出等值面图象

4　实验结果及分析

上述算法已经在移动工作站(2.8 GHz Intel处理器，4.0G内存，NVIDIA QuadroFX 770M显卡，显存512M)上利用标准OpenGL图形库实现。我们将粒子分布设置为长方体的形状，长方体分为两个部分，其中一个部分的粒子密度比较大，一个部分的粒子密度比较小。使用三种方法进行绘制，这三种方法分别是球形绘制、Marching Cubes算法、本文所用的方法。每种方法绘制的结果分成两种形式，一种为俯视图，一种为整体图。球形绘制的结果如图3(b)与4(a)所示，从其俯视图可以发现，使用球形绘制时，密度比较小的区域绘制结果比较好，但是密度较大的区域会出现明显的凹凸，表面不够平整，效果欠佳。而从其整体图也可以发现这个问题，即表面的光滑度不好。使用传统的Marching Cubes算法进行绘制的结果如图3(a)和4(b)所示，从图中可以发现，Marching Cubes算法绘制出来的结果中密度较大的区域绘制结果比较光滑，但是边缘与真实情况差异较大。本文所使用的方法如图3(c)和4(c)所示，该方法对粒子密度较大的区域采用Marching Cubes算法进行绘制，而对边缘的粒子采用球形绘制方法进行绘制。

图3　绘制结果的俯视图

图4　绘制结果的整体效果图

最后，我们将该方法应用于血流绘制，如图5所示。其中，图5(a)是利用本算法绘制出来的全局视角的血流图。我们可以看到在血流最靠近边缘的位置上有部分离散点，这些点就是使用球形绘制的；中间的部分比较平整，是使用Marching Cubes算法进行绘制的。图5(b)是同一次实验中将视角改变后所看到的结果。由于在改变视角的同时，粒子还处于运动情况，因此与图5(a)并不是同一帧。将视角转变为俯视的情况，从图中依然可以看出边缘是球形绘制，中间区域是Marching Cubes算法绘制。图6分别表示血流在第40帧、第50帧、第60帧时的绘制结果。从图中可以看到，本文的绘制算法结合了球形绘制和Marching Cubes算法，绘制出来的结果比较真实可靠。

(a)正常视图　　　　　　(b)俯视图(不同帧)图5　不同视角的血液流动现象绘制

(a)第40帧　　　　(b)第50帧　　(c)第60帧图6　不同帧的绘制结果

5　结论

本首先将粒子的位置进行预处理，使得分布杂乱无章的粒子转变为立方体顶点点阵。得到立方体顶点点阵之后，根据实时性的要求，将传统的面

绘制算法Marching Cubes算法与传统的球形绘制方法相结合，形成一种新的绘制方法。该方法首先将绘制区域划分成七个区域，其中六个区域采用传统的球形绘制算法进行绘制，另外一个区域采用精简版的Marching Cubes算法进行绘制。实验表明，该方法可以实时、真实地绘制出分布杂乱无章的粒子群。

[1] Rianto S, Li L. Fluid dynamic visualisations of cuttings-bleeding for virtual reality heart beating surgery simulation[C]// Australasian Computer Science Conference, Brisbane, Australia, January. Australian Computer Society, Inc. 2010:53-60.

[2]Qin J, Pang W M, Nguyen B P, et al. Particle-based simulation of blood flow and vessel wall interactions in virtual surgery[C]// Symposium on Information and Communication Technology, Soict 2010, Hanoi, Viet Nam, August. 2010:128-133.

[3]Reorowicz P, Obidowski D, Klosinski P, et al. Numerical simulations of the blood flow in the patient-specific arterial cerebral circle region.[J]. Journal of Biomechanics, 2014, 47(7):1642-51.

[4] 黄雷,肖双九,顾力栩,等.虚拟手术训练系统的血流模拟[J].计算机应用与软件,2011,28(1):65-68.

[5] 郑广超. 虚拟手术系统中模拟手术场景的渲染和平台的构建[D]. 上海:上海交通大学, 2008.

[6] 施鹏, 熊岳山, 徐凯, 等. 虚拟肝脏手术中实时动态渗血效果模拟[J]. 计算机应用, 2013, 33(10): 2911-2913.

[7] Lorensen W E, Cline H E. Marching cubes: A high resolution 3D surface construction algorithm[C]//ACM Siggraph Computer Graphics. ACM, 1987, 21(4): 163-169.

[8] Mor A B. Progressive cutting with minimal new element creation of soft tissue models for interactive surgical simulation[D]. Pittsburgh：Carnegie Mellon University, 2001.

(责任编辑：曾晶)

Research on Real-Time Blood Flow Simulation Based on the Improvement Marching Cubes Algorithm

WANG Na1， CHEN Guodong2*，CHEN Yi2

(1. Department of Electronics and Information Engineering, Fuqing Campus of Fujian Normal University, Fuqing 350300, China；2. Institute of Physics and Information Engineering, Fuzhou University,Fuzhou 350108, China)

The conventional Marching Cubes algorithm takes too much time in detecting null voxels dataset, which leads to the low efficiency. To resolve this issue and to conduct visual simulation of fluid surface images, an improved algorithm was proposed. The density of the rule lattices was detected and a threshold value was set. The particle density below the threshold region was isolated from the one above the threshold region. Only the regions with high density were simplified version of Marching Cubes algorithm. The simulation experiment results show that compared with Spherical Rendering algorithm and Marching Cubes, the Marching Cubes algorithm proposed in this paper decimates the access to null voxels with good rendering performances and it enables the real-time rendering of the model.

Marching Cubes (MC) algorithm; voxels; blood flow; isosurface; threshold

1000-5269(2016)02-0084-04

10.15958/j.cnki.gdxbzrb.2016.02.19

2015-10-07

国家自然科学基金(61471124)；福建省自然科学基金项目(2016J01293)；福建省科技计划重点项目(2011H0027)；福建省中青年教师教育科研项目(JA15574)

王娜(1978-)，女，副教授，硕士，研究方向：计算机图形学和虚拟现实，Email：studyres@126.com .

陈国栋，Email：fzucgd@126.com.

TP391.41

A