CN103236058A - 获取四维心脏图像感兴趣体积的方法 - Google Patents

Abstract

一种图像处理技术领域的获取四维心脏图像感兴趣体积的方法，通过读取组织结构的体数据并进行体素化处理，经过三角面片集合向交互式多边形选择区域投影，得到基于体数据的心肌四维可视化数据；然后通过对网格化投影面中选择区域的精确提取、网格扩充处理，得到包含在感兴趣区域内的所有面元，最后通过八叉树编码，一致性提取得到任意时序图像的体素集合，从而实现基于体数据的心血管系统四维可视化数据。本发明可以自动的获得提取心脏感兴趣区域的轮廓所对应的局部组织体数据集，最后实现心脏四维图像的局部提取和可视化，由于提取的局部体数据集所占空间明显小于整体体数据集所占空间，所以提取出的图像堆栈所占有的空间与提取前的图像堆栈相比会大大减少，并由此显著降低了四维可视化所消耗的时间。

Description

获取四维心脏图像感兴趣体积的方法

技术领域

本发明涉及的是一种医学图像处理技术领域的方法，具体是一种在心脏体数据表面画出多边形感兴趣区域来对心脏的感兴趣体积（VOI，volume of interest）进行提取的方法。

背景技术

世界心脏联合会宣布，目前全世界每年因心脏疾病死亡的人数占死亡总人数的三分之一。因此，对心脏疾病的提前预防和治疗就变得十分重要。截至目前为止已经出现了很多针对人类心脏中一些解剖组织的分割算法，最常见的就是针对冠状动脉、心房和心室等组织结构的分割。由于这些分割算法的分割对象具有针对性，所以不可以用一种统一的方法进行分割。因为心脏的组织结构十分复杂，包括心房、心室、冠状动脉、血管、心包以及瓣膜等组织结构，心脏的连续搏动使得心脏影像具有连续改变的位移特性，很难进行四维图像的提取和定位，这就要求寻找到一种算法既能够完成局部任意感兴趣体积的三维提取，也能够方便选定任意三维组织进行四维可视化。这样不但满足了计算机辅助治疗和临床需求，在科研和临床也有很重要的意义。

伴随三维医学成像技术的快速发展，使得观察组织结构内部或局部的立体结构成为了可能，从而可以观察到人体局部器官详细的结构和功能。虽然医学图像三维重建技术已经很成熟，可是感兴趣体积分割、四维可视化和量化分析等技术仍然不成熟。在对感兴趣体积的提取之前必须要保证被识别和分割区域的可靠性，使得体积测量和心脏局部体积三维模型堆栈的四维可视化可以顺利进行。对心脏三维医学图像的任意区域实施提取是非常重要的研究工作，因为它是实现心脏四维可视化和手术导航的基础。介于心脏解剖结构及其运动生理特点的复杂性，使得心脏边界的选择和局部体积的分割成为了一个难题。

中国专利文献号为：CN200580040748.3，公开日：2007.10.31，发明名称为：“感兴趣体积的选择”中公开了一种在对象数据集中创建感兴趣体积（VOI）提取的方法，该对象数据集被格式化成对象数据切片，其中至少两个对象数据切片存在与该对象数据切片内的感兴趣体积相关的部分，这些切片上的相关点构成感兴趣区域的轮廓线。定义与至少两个感兴趣区域相切的表面，定义上述表面与轮廓线相交的曲线方程，通过该曲线方程，这些曲线可以限定该表面内感兴趣体积部分的边界，通过优选自动产生切片内感兴趣体积的轮廓线。但该技术只针对有一定封闭性的组织器官，同时曲线方程只可以被限定在边缘线附近，对组织内部和周围分散的血管组织来说用该方法就不再适合，不适合任意组织的任意区域的分割。对于切面和曲线的获取，该算法的复杂度也十分大。

美国专利文献号为：US2006/0170679A1，公开时间为2006.08.03，发明名称为“Representing a volume of interest as boolean combinations of multiple simple contour sets”中公开了一种用软件和硬件相结合来完成用多个简单轮廓集合以布尔表达式的形式代表感兴趣体积（VOI）的方法，在感兴趣体积（VOI）的提取过程中始终围绕单轮廓集合的选取上，该方法将体数据分成四层，每一层都会提取出相应的轮廓信息，该技术中还就轮廓可能产生的不同形状进行了分析和讨论，每一层的信息还要同时以布尔表达式的形式被硬件设备处理，该技术不能实现任意区域的四维提取和可视化，其实现过程是软硬件相结合来实现感兴趣体积（VOI）的提取，使得实现成本和复杂度无法满足现有操作需要。

Park，Sang Cheol；Kim，Won Pil；Zheng，Bin；等人2009.02.10，在“Medical Imaging2009-Image Processing”会议上发表的文章“Pulmonary airways tree segmentation from CTexaminations using adaptive volume of interest”中写到了一种结合自适应确定阈值的区域增长算法来提取肺部感兴趣体积（VOI）的分割，该文献主要用圆筒状感兴趣模型实现了肺支气管的三维分割和显示，这种区域增长算法为了降低在给定感兴趣体积内选择上的遗漏，采用迭代的方式来降低这种可能性的出现，这样做大大降低了算法的性能，使得分割一次感兴趣体积需要耗时大约一个小时十分钟。这种算法消耗时间太长，不适合临床应用。

Kim，Jinman；Cai，Weidong；Feng，Dagan等人在“Institute of Electrical and ElectronicsEngineers”期刊上发表的文章“A new way for multidimensional medical data management：Volumeof interest（VOI）-based retrieval of medical images with visual and functional features”中写到了一种应用模糊C-均值聚类分析算法实现感兴趣体积（VOI）的分割，然后将分隔出来的VOI存入数据库，最后实现VOI的图像检索，这里面应用的模糊C-均值聚类分析算法本身易陷入局部极小值和对初始值敏感的缺点，尤其在聚类数越大的情况下这种缺点越明显。

Guy Shechter，Cengizhan Ozturk，和Elliot R.McVeigh在2000.02.13，的“Image Display andVisualization”会议上发表的文章“Interactive Four Dimensional Segmentation of Multiple ImageSets”，详细记载了一种四维心脏数据的分割和可视化工具，该文献详细介绍了该交互式四维分割和可视化工具的功能和实现过程，该可视化软件在实现感兴趣区域的四维分割过程是建立在多边形网格的基础上得到包围感兴趣区域的曲面然后进行分割的，与本人的发明在网格的选取上和算法设计上有很大的区别，该文献在分割过程中还运用了B样条曲面插值算法对给定的轮廓进行细化。存在的不足是B样条曲面插值算法没有充分考虑待插值点周围领域范围的点对待插值点的影响，在处理图像分辨率上表现欠佳。除此之外，在分割心脏数据的过程中，只对有明显轮廓的序列图像才有效，有一定局限性。

中国专利文献号为：CN200510099178.7，公开日：2008.03.15，发明名称为：“从四维图像数据组中分割出解剖结构的方法”，中公开了一种用于以时间序列的心脏三维图像中分割出血管解剖结构的方法，其基本思想是从时间序列三维图像数据组中分割出冠状血管树，可以缩短分割同等数量的序列的三维数据所消耗的时间，这种方法可以高效的对序列三维图像进行一致性提取。这种方法主要依靠连续的三维图像之间不可能任意移动，这个特点不仅可以限制立体数据组中代表冠状血管树的一部分元素，也可以明显的提高分割速度。通过时间步骤中的计算结果可以显著的改善分割的可靠性。该技术应用于分割冠状血管而非心脏的任意区域，同时也不能完成对血管的任意局部进行再度分割的功能，对于实现诊断血管及心房等内部结构来说都不够灵活。

空间分割技术研究的核心是将一个给定三维图形按照特定的区域划分法则实施快速高效的分割。在三维显示技术中，所提供的模型通常都是空间立体模型，为了绘制的更加快捷，立体模型要预先离散为单个的面元模形，空间分割可以概括为对空间大量面元按照特定规则进行分区的问题，对立方体分区，平移和缩放变换不能改变面元与分区之间的空间位置关系，空间分割技术研究的核心是将一个给定三维图形按照特定的区域划分法则实施快速高效的分割。为对空间大量面元进行体素化，并按照特定规则对体素进行分区。将物体的几何形式表示转换成最接近该物体的体素表示形式，产生体数据集，其不仅包含模型的表面信息，而且能描述模型的内部属性。

空间分割技术不仅在数据较大的计算机图形显示中发挥着主导作用，近年来已经延伸到其他三维显示技术相关的领域。该方法将大的三维图形按空间位置关系拆解为较小的子模块，利用子模块之间的空间位置关系的固定性，能够快速的对较大的三维立体数据进行操作、处理。例如二叉空间分割BSP（Binary Space Partitioning）以其高效和递归性成为空间分割技术的主导，可是BSP的缺点是：当想选择一个将空间中的多边形分割成一对相等的子集合的面是很复杂的，因为有无数个可以挑选的面。所以在集合中挑选最适合的面实现起来是很困难的。

发明内容

本发明针对现有技术存在的上述不足，提出一种获取四维心脏图像感兴趣体积的方法，可以自动的获得提取心脏感兴趣区域的轮廓所对应的局部组织体数据集，最后实现心脏四维图像的局部提取和可视化，由于通过适当的转变后，提取出的图像堆栈所占有的空间与提取前的图像堆栈相比会大大减少，并由此显著降低了四维可视化所消耗的时间。

本发明是通过以下技术方案实现的，本发明通过读取组织结构的体数据并进行体素化处理，经过三角面片集合向交互式多边形选择区域投影，得到基于体数据的心肌四维可视化数据；然后通过对网格化投影面中选择区域的精确提取、网格扩充处理，得到包含在感兴趣区域内的所有面元，最后通过八叉树编码，一致性提取得到任意时序图像的体素集合，从而实现基于体数据的心血管系统四维可视化数据。

所述的组织结构包括：冠状动脉、心室、心房、瓣膜、心包以及血管等组织结构。

所述的体数据是指：组织结构的序列CT图像，采用最小二乘法，得到各采样点对应的近似切平面中心，同时计算其组成的协方差矩阵，计算其特征值，把最小的特征值所对应的单位特征向量，当作该近似切平面的法向量。对切平面进行法向量调整使其指向曲面的同侧，对其进行局部线性逼近，构造一有符号距离函数。最后，对距离函数进行插值，采用Marching Cubes算法来进行等值面的抽取，最终输出重建结果，得到由三角形面片组成的表面模型。

所述的交互式多边形选择区域投影是指：将三角面片集合向投影面做投影，得到投影区域，然后通过交互方式确定一个任意多边形后，判断多边形的每一条边与投影面中正方形网格的位置关系：

当该直线经过A(x_a，y_a)和B(x_b，y_b)，那么利用两点式直线的方程可表示为：由公式（5）通过连线的方程将x递增并向下取整，很容易得出x=1，2，3，4，…时y的值：由公式（6）再对y_n值向下取整，所得(x₁，y₁)，(x₂，y₂)，(x₃，y₃)，…一定是这条直线经过的正方体网格。但上述方法会遗漏某些特殊情况，即左上方的正方形网格也可能是该直线经过的正方形网格，这是由于只在x轴上递增的向下取整，才导致这些点没有被考虑在内。

当该直线的斜率和A点与第一个网格右上方顶点M(x_C，y_C)连线的斜率k做比较来得出结论。用公式（7）利用这种方法能够快速的计算出直线经过的所有正方形网格。推广到任意多边形所经过的网格也可以全部求出来。

所述的体素化（Voxelization）处理是指：将组织结构由三角形面片组成的表面模型转换成内部离散化的体素（Voxel）并产生体数据集（Volume Datasets）；所述的体素（Voxel）是指：二维像素在三维空间的转换，它们是一组分布在正交网格中心的立方体单元，不仅可以包含模型的表面信息，还能描述模型的内部属性。

所述的体素化处理包括：三角面片顶点的体素化、三角面片边的体素化、三角面片的体素化以及三角面片包围模型对模型内部的体素化。

所述的精确提取包括：获取多边形经过正方形网格、获取面元边缘点的坐标集合以及获取多边形的最小外接矩形。

所述的网格扩充是指：对最小外接矩形进行网格扩充循环递归，然后用射线法得到包含在感兴趣区域内的所有面元；

所述的八叉树编码是指：八叉树结构表达三维图像的方法，其定位码规则包括：

1）给定分辨率n，也即确定了坐标系统的大小，每个坐标轴的取值范围从0到2ⁿ-1，n的取值为3。

2）八叉树每个点的编码位置q_i是（0，1，2，…，7）八个数之一，q_i的个数取决于分辨率n。

3）整个编码按Z字形的方式进行，其方向取决于如何选择第一、二、三和第四维坐标，基于第一个规则，坐标（X，Y，Z）已知的八叉树点的编码可用下式表示：

$\left\{\begin{array}{c}X=c_{n-1}{2}^{n-1}+c_{n-2}{2}^{n-2}+...+c_0\\ Y=d_{n-1}{2}^{n-1}+d_{n-2}{2}^{n-2}+...+d_0\\ Z=e_{n-1}{2}^{n-1}+e_{n-2}{2}^{n-2}+...+e_0\end{array}\right.---\left(1\right)$

其中：系数c_i、d_i和e_i（i=0，…，n-2，n-l）的取值为0或1，可通过比较方程式两边的值来定。根据步骤2、3，q_i值可由下式确定：

$q_i=e_i{2}^2+d_i{2}^1+c_i(i=n-1,n-2,\·\·\·,0)---\left(2\right)$

所述的八叉树编码优选为线性八叉树，如线性八叉树的地址码有八进制和十进制数两种，其中：十进制地址码亦称Morton码。由于Morton码是自然数码，所以可以将二维数组转化成以Morton码为下标的一维数组。

所述的线性八叉树的自然数编码的表达式为：

$N=(c_0{8}^0+c_1{8}^1+...+c_{n-1}{8}^{n-1})+2^1(d_0{8}^0+d_1{8}^1+...+d_{n-1}{8}^{n-1})+2^2(e_0{8}^0+e_1{8}^1+...+e_{n-1}+8^{n-1})$ （3）

其中：c₀，d₁，c_n-1，d₀，d_n-1和e₀，e₁，e_n-1分别是行号和列号二进制化后由低位到高位的权；当图形、图像恢复时，可由N码经逆变换为下列十进制的行、列号：

$I_i=\underset{k=n-1}{\overset{k=0}{\mathrm{\Σ}}}\mathrm{MOD}(T_K,2)\·2^k---\left(4\right);$

其中：I_i表示维欧氏空间中第i个坐标轴的十进制坐标，T_K表示N码与8^k比值的取整结果，k表示迭代项，其中k为非负整数。

技术效果

与现有技术相比，本发明采用的四维图像空间分割利用心脏序列图像的顺序，一次性将心脏周期的所有三维图像分割成有一定顺序的体元集合，不但可以高效的对四维心脏序列图像标出感兴趣种子点，提取出关联的四维感兴趣序列体数据，还能够观察心脏任意多边形部分覆盖的体元集合（这些集合可以是：冠状动脉、心室、心房、瓣膜、心包以及血管等组织结构）的运动状态，这里的任意部分是指四维心脏图像上所有医生感兴趣的体积，是用户通过人机交互手段在四维图像上留下的种子点连接而成的多边形。这种提取方法对实现血管及心脏局部区域的四维可视化是十分有价值的。

附图说明

图1为本发明流程示意图。

图2为三维物体三角网格投影到正交网格面示意图。

图3为当y_n值经过向下取整示意图；

图中：直线在投影面上经过的正方体网格。

图4为本发明中y_n递增向下取整后，容易忽略网格的情况。

图5为投影面示意图；

图中：（a）为通过人工交互方式画出的凸多边形，（b）为包含在多边形内部的正方形网格，这些方格不与多边形相交，（c）为多边形包含和经过的所有正方形网格。

图6为正方形网格集合示意图；

图中：M_n为最小外接矩形，R₀为最小外接矩形包括的正方形网格所在的集合，R₁为交互多边形内的正方形网格所在的集合，R为交互多边形经过的所有正方形网格所在的集合。

图7为本发明人工交互方式示意图；

图中：对最小外接矩形内的每一个边缘网格进行N次网格划分，得到的任意多边形面元示意图。

图8为使用射线法判断点在多边形内部的四种特殊情况示意图。

图9为实施例效果示意图；

图中：（a）为本发明在VOI提取过程中三组横断面，矢状面，横断面终止面和VR显示效果图，其中：1、5、9为横断面、2、6、10为矢状面，3、7、11为横断面终止面，4、8、12为VR显示效果图；（b）中：1为通过任意多边形面元提取出的血管簇的三维模型，2为通过人工交互操作在四维图形中画出的任意多边形面元；（c）中：1为通过任意多边形面元提取出的整个心脏的三维模型，2为通过人工交互操作在四维图形中画出的任意多边形面元。

图10为本发明利用感兴趣多边形得到的最小立方体盒，对序列体数据进行一次性提取的原理图。

具体实施方式

下面对本发明的实施例作详细说明，本实施例在以本发明技术方案为前提下进行实施，给出了详细的实施方式和具体的操作过程，但本发明的保护范围不限于下述的实施例。

实施例1

如图1所示，本实施例包括以下步骤：

预处理：

第一步：读取序列二维切片（CT图像来自于回顾性心电门控技术采集的不同节段的切片集）对心脏组织结构的序列CT图像采用最小二乘法，得到各采样点对应的近似切平面中心，同时计算其组成的协方差矩阵，计算其特征值，把最小的特征值所对应的单位特征向量，作为该近似切平面的法向量。对切平面进行法向量调整使其指向曲面的同侧，对其进行局部线性逼近，构造一有符号距离函数。最后，对距离函数进行插值，采用Marching Cubes算法来进行等值面的抽取，最终输出重建结果，得到由三角形面片组成的表面模型（即读取体数据后将其转化成三角面片集合）

第二步：模型的体素化

1.1.0三维实体的表面模型的基本元素由三角形面片构成。为方便起见，初始时将模型的AABB包围盒（Axis-aligned bounding box）变换到三维空间的第1象限中，也即所有坐标值均为正。

所述的AABB包围盒是指：利用体积较大且特性简单的几何体，来近似的代替复杂的几何对象。立方体的包围盒用（x，y，z，l）表示，其中，x，y，z是该立方体的顶点坐标，l表示立方体的边长，此处所述的AABB包围盒指代的是包围所有三角形面片的最小立方体模型。

1.1.1三角面片顶点的体素化：通过求出D₀点所在体素的局部网格坐标，由公式（3）可以得出顶点所在体素，其中：D₀指代三角面片的顶点。

1.1.2三角面片边的体素化：处理三角面片的边，相当于从起点发出一条射线，然后跟踪该射线，找到该射线穿过的所有体素。一条射线从起点所在的体素发出，指向所处理边的终点，必与起点体素相交，交点为P，相交的可能有3种情况：与面相交、与棱边相交和与顶点相交。按上面的方法递归地找到与该边相交的体素，直到找到边终点所在的体素为止。

1.1.3三角面片的体素化：将三角面D₀D₁D₂投影到xoz、xoy和yoz平面时投影三角形为D_d0D_d1D_d2。前面已经找到了属于顶点和边的体素，扫描网格坐标的两个分量iu（k）和iv（k），当有网格坐标分量i0满足iu（k）<i0<iv（k）。就找到了属于三角面内部的体素网格。

1.1.4三角面片包围模型对模型内部的体素化：将模型表面体素化的操作进行完之后即可得到对模型体素表示的一个“外壳”，接下来要做的操作就是进行模型的内部体素化操作。

这里采用一种简单的方法：对于模型AABB中的所有空体素，以其中心位置以轴对齐方向来发射两条射线，这两条射线的方向相反，但基本方向都是轴对齐的。对于这两条射线利用空间模型与3D模型的相交位置，并得到相交点的法向量及到相交点的距离，然后根据这两点法向量之间的关系来判断得到当前体素是在3D模型的内部或是在3D模型的外部。将这样的操作施加于每一个空的体素之后就可以完成对3D模型内部的体素化操作。

但是将这样的操作施加于每一个空的体素速度比较慢，故而此处可以采用扫描的方法来进行加速处理；当判断得到某个体素的位置为模型内部后，就可以根据射线的方向及这两条射线与模型的交点处的距离来对当前体素相邻的体素进行扫描，这样不需要再做判断就可以标记出相邻体素的状态，这样就加速了整个模型内部的体素化操作。

第三步：三角面片集合向交互式多边形选择区域投影

在这个离散化的立方体体元集合中，为了使感兴趣体积能够通过手工画出的选择区域被提取出来，就要将三角面片集合向投影面做投影，得到投影区域，然后通过鼠标画出一个任意多边形后，判断多边形的每一条边与投影面中正方形网格的位置关系，如图2所示：

当该直线经过A(x_a，y_a)和B(x_b，y_b)，那么利用公式（5）通过连线的方程将x递增并向下取整，很容易得出x=1，2，3，4，…时y的值：由公式（6）再对y_n值向下取整，所得(x₁，y₁)，(x₂，y₂)，(x₃，y₃)，…一定是这条直线经过的正方体网格。但上述方法会遗漏某些特殊情况，即左上方的正方形网格也可能是该直线经过的正方形网格，这是由于只在x轴上递增的向下取整，才导致这些点没有被考虑在内。

当该直线的斜率和A点与第一个网格右上方顶点M(x_C，y_C)连线的斜率k做比较来得出结论。用公式（7）利用这种方法能够快速的计算出直线经过的所有正方形网格。推广到任意多边形所经过的网格也可以全部求出来。

第四步：网格化投影面中选择区域的精确提取

3.1.1获取多边形经过正方形网格

设直线j上有两个坐标分量都是非整数的点A（x_a，y_a）和B（x_b，y_b）。除那些包括端点A，B的网格外，直线经过的矩形网格都能够和直线有2个交点，而相连的网格公用1个交点，那么直线AB经过的矩形网格数目为M=N+1。但是仅仅知道经过的正方体网格的数目还远远不够，还要求出直线具体经过了哪些正方形网格。由于该直线经过A（x_a，y_a）和B（x_b，y_b），那么利用直线方程可表示为：

$(x_b-x_a)(y-y_a)-(x-x_a)(y_b-y_a)=0---\left(5\right)$

通过连线的方程（1）将x递增并向下取整，如图6所示，很容易得出x=1，2，3，4，…时y的值：

$y_n=(x_b-x_a)(y-y_a)-(x-x_a)(y_b-y_a),x_n=\mathrm{1,2,3}\&CenterDot;\&CenterDot;\&CenterDot;---\left(6\right)$

再对y_n值向下取整，所得（x₁，y₁），（x₂，y₂），（x₃，y₃），…是这条直线经过的第一次划分的正方体网格。如图3所示。

但上述方法会遗漏某些特殊情况，即图4a和图4b中的A网格可能是该直线经过的正方形网格，这是由于只在x轴上递增的向下取整，才导致这些点没有被考虑在内，如图4所示。

解决这个问题可以利用该直线的斜率和A点与图中红色中心点C（x_C， _C）连线的斜率k做比较来得出结论，即：

$k\max =\left\{\begin{array}{c}\mathrm{kAB}=\tan \mathrm{\&alpha;}=(\mathrm{yb}-\mathrm{ya})/(\mathrm{xb}-\mathrm{xa})\\ \mathrm{kAC}=\tan \mathrm{\&beta;}=(\mathrm{yc}-\mathrm{ya})/(\mathrm{xc}-\mathrm{xa})\end{array}\right.---\left(7\right)$

其中：k_max是指：直线AB与直线AC斜率比较后，其中斜率值较大的那条直线的斜率值，当k_AB＞0，那么判断当k_AB＞k_AC那么直线经过网格A；当k_AB＜0，那么当k_AB＜k_AC那么直线经过网格A；利用这种方法能够快速的计算出直线经过的所有正方形网格。

推广到任意多边形所经过的网格也可以全部求出来。这里所有经过的网格集合如图5所示，为： $R=\left\{R\right|\mathrm{Meshid}(x_n,y_n),n=\mathrm{1,2,3}\&CenterDot;\&CenterDot;\&CenterDot;\}---\left(8\right);$

其中：Meshid（）是指：直线经过的网格，R是指直线经过的所有网格所构成集合：

3.1.2获取面元边缘点的坐标集合

从集合R中任意取出一个多边形经过的正方形网格，它的四个端点坐标分别为G₁（e，f），G₂（e+1，f），G₃（e，f+1），G₄（e+1，f+1）。这样就能够得到所有多边形经过的正方形网格顶点（即边缘点）的集合为

$U=\left\{U\right|\mathrm{Gn}1(\mathrm{en},\mathrm{fn}),\mathrm{Gn}2(\mathrm{en},1,\mathrm{fn}),\mathrm{Gn}3(\mathrm{en},\mathrm{fn}+1),\mathrm{Gn}4(\mathrm{en}+1,\mathrm{fn}+1)\}---\left(9\right)$

其中：U对应于面元的边缘点的坐标集合，en表示多边形经过正方形网格的四个顶点的横坐标，fn表示多边形经过正方形网格的四个顶点的纵坐标，Gni表示多边形经过的所有正方形网格的四个顶点在投影面上的坐标（i=1、2、3、4…）。

3.1.3获取多边形的最小外接矩形

由于每个平面都是由两组平行线正交构成，那么可以把每个面上的所有直线正交后所得的交点看作一个以网格端点的集合。因为在包含所有网格点端点内的所有矩形中，最小外接矩形一定是一个以四个极点即左极点S₁（x₁，y₁），右极点S₂（x₂，y₂），上极点S₃（x₃，y₃），下极点S₄（x₄，y₄）的坐标为边缘的最小外接矩形，如图6，那么该最小外接矩形包括的正方形网格为：

$Z=\left\{Z\right|\mathrm{Meshid}(x_n,y_n),x_1\&le;x_n\&le;x_2,y_4\&le;y_n\&le;y_3\}---\left(10\right)$

第五步：对最小外接矩形进行网格扩充

对最小外接矩形进行网格扩充的循环递归，得到优化外接矩形，通过判断得知这里所有经过的正方形网格集合为：

$R^{\&prime;}=\{{\mathrm{Meshid}}^{*}(x^{*},y^{*}),n=\mathrm{1,2,3}\&CenterDot;\&CenterDot;\&CenterDot;\}---\left(11\right)$

所述的网格扩充循环递归的具体步骤为：当得到最小外接矩形后，通过改变网格间距可以将最小外接矩形划分成许多更小的正方形网格，取新的网格间距是初始网格间距的二分之一，通过坐标转换可以得到多边形的顶点在新坐标系下的坐标，由于网格细分，所以多边形端点更靠近细分后的网格边线；当完成一遍10次递归直至多边形的端点和优化外接矩形边缘重合，而多边形的边线在边缘网格内也将与边缘网格的边线重合时，得到优化外接矩形中多边形的轮廓线，其中：n=为递归次数，一般设置为10次，如图7所示。

第六步：用射线法得到包含在感兴趣区域内的所有面元

在优化外接矩形中的网格R₀的网格中心点为起点用射线法得到包含在感兴趣区域内的所有体元，在最小外接矩形M_n中找到包含在优化外接矩形内部的所有网格R₀和在优化外接矩形外的所有网格R*。

所述的射线法是指：首先遍历优化外接矩形中的每一行的网格找到每一行的网格中除去多边形的边经过的所有的网格，将每一行的空白网格全部挑出来并把这些网格顶点的坐标存在一个数组中，通过计算可以得到这些网格的中心点P的坐标；由于当从P作水平向右射线的话，如P在多边形内部，那么这条射线与多边形的交点个数必为奇数，当P在多边形外部，则交点个数必为偶数（0也在内）。所以，顺序考虑多边形的每条边，求出交点的总个数。如图8所示，还有一些特殊情况要考虑：

1）当射线正好穿过P₁或者P₂，那么这个交点会被算作2次，处理办法是当P的纵坐标与P₁，P₂中较小的纵坐标相同，则直接忽略这个顶点。

2）当边水平，则射线要么与其无交点，要么有无数个，这种情况也直接忽略这条边与射线的交点。

3）当边竖直，而P₀的横坐标小于P₁，P₂的横坐标，则必然相交。

4）再判断相交之前，先判断P是否在边P₁P₂的上面，当在边P₁P₂上则直接得出结论：P在多边形内部。

所述的射线法针对任意形状的多边形，包括：凹多边形、凸多边形、内部交叉多边形等任意形状的多边形。

第七步：用线性八叉树完成感兴趣体素集合的提取和显示

通过对步骤3.1.1的平面几何运算可以得到感兴趣多边形所经过和包含的所有正方形网格，通过对网格的坐标转换可以得到包围盒中对应的三角面片，通过三角面片和体素之间的关系可知：

Ⅰ、三角面片点的编码是通过求出D₀点所在体素的局部网格坐标，由公式（3）可以得出顶点所在体素网格的自然数编码N，由公式（4）经逆变换为十进制的行、列号，求出体素编码，通过索引编码就可以获得D₀所在的体素位置；

Ⅱ、处理三角面片的边与面相交、与棱边相交和与顶点相交不论是哪种相交情况，根据立方体体素自身的编码特点，均返回一个索引值，由此索引值，按照编码规则可以计算出与此索引值相邻的体素的编码，按上面的方法递归地找到与该边相交的体素，直到找到边终点所在的体素为止。

Ⅲ、前面已经找到了属于顶点和边的体素，扫描网格坐标的两个分量iu（k）和iv（k），当有网格坐标分量i0满足iu（k）<i0<iv（k）。就找到了属于三角面内部的体素网格。

由公式（3）可以得出三角面片所在的体素网格编码的自然数编码N，由公式（4）经逆变换为十进制的行、列号，求出体素编码，通过索引编码就可以获得D₀D₁D₂所在的体素位置。Ⅳ、模型内部的体素化：从中每一个空体素心位置以轴对齐方向来发射两条射线，得到其与3D模型的相交位置，并得到相交点的法向量及到相交点的距离，然后根据这两点法向量之间的关系来判断得到当前体素是在3D模型的内部或是在3D模型的外部。将这样的操作施加于每一个空的体素之后就可以完成对3D模型内部的体素化操作。。借助体素化转换可以将感兴趣区域转化成对应的体素集合，最后在最小立方体盒中利用线性八叉树结构的编码特性（Z字编码），根据26-邻接体素的拓扑关系，寻找感兴趣体积的体素集合，通过八叉树索引编码就可以获得该体素所在的体素位置，最后借助STL(标准模版库)的vector容器存储三维数据集索引，通过句柄传递实现时序三维体数据序列的查询和提取，根据26-邻接体素的拓扑关系，寻找感兴趣体积的体素集合。最后经过传输函数加以显示，最终实现了医学图像任意感兴趣体积的提取和四维可视化，如图9所示。

第八步：实现四维数据的一致性提取

由于四维图像是有时间顺序的三维数据序列，为了提高提取算法的高效性，本实施例中对序列体数据采用一致性提取方法，通过计算已经得到感兴趣体积的像素集合，得到包含在感兴趣多边形内的所有体素集合。由于任意相邻的两张四维数据之间都存在体素连续性，那么可以利用血管的这种生理特征来对感兴趣体积进行一致性提取，这样，其他的三维序列进行相同操作，如图13所示，得到其余图像相同选择区域所在的感兴趣体积，然后在相同位置提取出总的感兴趣体积集合为：Voi={V_n，1≤n≤N}。提取的顺序和序列体数据的原有顺序保持一致。最后将提取得到的血管序列三维图像进行四维可视化。如图10所示。

实施验证

1）表1中的Cardiac模型和表1中Head模型采用128的分辨率和256的分辨率时，三角面数、表面体素数和离散时间之间的关系，该算法在处理时间上有很大的提高。

表1模型参数及体素模型计算时间

2）表2中的不同物理尺寸的心脏局部模型，取不同三角面数的模型，在512分辨率下的处理时间相对较短，而且物理尺寸的增大对该算法的处理时间影响不大，但三角面数选取的越多，图像效果也比较清晰，表现出了该算法的在处理数据方面的优势。可见本算法采用较低的分辨率就可以比较好的逼近模型，并且具有较快的速度，如表2所示。

表2四个局部模型参数及26-邻接体素模型计算时间

Claims (10)

1.一种获取四维心脏图像感兴趣体积的方法，其特征在于，通过读取组织结构的体数据并进行体素化处理，经过三角面片集合向交互式多边形选择区域投影，得到基于体数据的心肌四维可视化数据；然后通过对网格化投影面中选择区域的精确提取、网格扩充处理，得到包含在感兴趣区域内的所有面元，最后通过八叉树编码，一致性提取得到任意时序图像的体素集合，从而实现基于体数据的心血管系统四维可视化数据。

2.根据权利要求1所述的方法，其特征是，所述的交互式多边形选择区域投影是指：将三角面片集合向投影面做投影，得到投影区域，然后通过交互方式确定一个任意多边形后，判断多边形的每一条边与投影面中正方形网格的位置关系：

当该直线经过A(x_a，y_a)和B(x_b，y_b)，那么利用两点式直线的方程可表示为：由（x_b–x_a）（y–y_a）–（x–x_a）（y_b–y_a）=0通过连线的方程将x递增并向下取整，很容易得出x=1，2，3，4，…时y的值：由y_n=（x_b–x_a）（y–y_a）–（x–x_a）（y_b–y_a），x_n=1，2，3…再对y_n值向下取整，所得(x₁，y₁)，(x₂，y₂)，(x₃，y₃)，…则为这条直线经过的正方体网格；

当左上方的正方形网格为该直线经过的正方形网格时，通过将直线的斜率和A点与第一个正方形网格右上方顶点M(x_C，y_C)连线的斜率k做比较来得出结论，即：

$k\max =\left\{\begin{array}{c}\mathrm{kAB}=\tan \mathrm{\&alpha;}=(\mathrm{yb}-\mathrm{ya})/(\mathrm{xb}-\mathrm{xa})\\ \mathrm{kAC}=\tan \mathrm{\&beta;}=(\mathrm{yc}-\mathrm{ya})/(\mathrm{xc}-\mathrm{xa})\end{array}\right.$ ，其中：k_max是指：直线AB与直线AC斜率比较后，其中斜率值较大的那条直线的斜率值，当k_AB＞0，那么判断当k_AB＞k_AC那么直线经过网格A；当k_AB＜0，那么当k_AB＜k_AC那么直线经过网格A；利用这种方法能够快速的计算出直线经过的所有正方形网格。

3.根据权利要求1所述的方法，其特征是，所述的体素化处理是指：将组织结构的由三角形面片组成的表面模型转换成内部离散化的体素并产生体数据集；所述的体素是指：二维像素在三维空间的转换，它们是一组分布在正交网格中心的立方体单元；

所述的体素化处理包括三角面片顶点的体素化、三角面片边的体素化、三角面片的体素化以及三角面片包围模型对模型内部的体素化。

4.根据权利要求1或3所述的方法，其特征是，所述的体素化处理具体包括以下步骤：

1.1.0）初始时将模型的AABB包围盒变换到三维空间的第1象限中，也即所有坐标值均为正；

1.1.1）三角面片顶点的体素化：通过求出D0点所在体素的局部网格坐标，由线性八叉树的自然数编码的表达式：

N＝(c₀8⁰+c₁8¹+...+c_n-18^n-1)+2¹(d₀8⁰+d₁8¹+...+d_n-18^n-1)+2²(e₀8⁰+e₁8¹+...+e_n-18^n-1)

得出顶点所在体素，其中：D₀指代三角面片的顶点；

1.1.2）三角面片边的体素化：根据任一三角面片边从起点所在的体素发出，指向所处理边的终点，与起点体素相交的交点为P、相交的情况包括：与面相交、与棱边相交和与顶点相交；递归地找到与三角形面片边相交的体素，直到找到边终点所在的体素为止；

1.1.3）三角面片的体素化：将三角面D₀D₁D₂投影到xoz、xoy和yoz平面时投影三角形为D_d0D_d1D_d2；对应已经获得的属于顶点和边的体素，扫描网格坐标的两个分量iu（k）和iv（k），当有网格坐标分量i0满足iu（k）<i0<iv（k），则获得属于三角面内部的体素网格；

1.1.4）三角面片包围模型对模型内部的体素化：

1.1.4.1）对于模型AABB中的所有空体素，以其中心位置以轴对齐方向来发射两条射线，这两条射线的方向相反，但轴对齐；

1.1.4.2）对于这两条射线利用空间模型与3D模型的相交位置，并得到相交点的法向量及到相交点的距离，然后根据这两点法向量之间的关系来判断得到当前体素是在3D模型的内部或是在3D模型的外部；

1.1.4.3)将这样的操作施加于每一个空的体素之后就可以完成对3D模型内部的体素化操作。

5.根据权利要求4所述的方法，其特征是，对于步骤1.1.4）中，当判断得到某个体素的位置为模型内部后，则根据射线的方向及这两条射线与模型的交点处的距离来对当前体素相邻的体素进行扫描，以实现整个模型内部的体素化操作的加速。

6.根据权利要求1所述的方法，其特征是，所述的精确提取包括：获取多边形经过正方形网格、获取面元边缘点的坐标集合以及获取多边形的最小外接矩形。

7.根据权利要求1或6所述的方法，其特征是，所述的精确提取具体包括：

3.1.1获取多边形经过正方形网格：设直线j上有两个坐标分量都是非整数的点A（x_a，y_a）和B（x_b，y_b）；除那些包括端点A，B的网格外，直线经过的矩形网格都能够和直线有2个交点，而相连的网格公用1个交点，那么直线AB经过的矩形网格数目为M=N+1；但是仅仅知道经过的正方体网格的数目还远远不够，还要求出直线具体经过了哪些正方形网格；由于该直线经过A（x_a，y_a）和B（x_b，y_b），那么利用直线方程可表示为：

（x_b–x_a）（y–y_a）–（x–x_a）（y_b–y_a）=0，

通过连线的方程 $\left\{\begin{array}{c}X=c_{n-1}{2}^{n-1}+c_{n-2}{2}^{n-2}+\&CenterDot;\&CenterDot;\&CenterDot;+c_0\\ Y=d_{n-1}{2}^{n-1}+d_{n-2}{2}^{n-2}+\&CenterDot;\&CenterDot;\&CenterDot;{+d}_0\\ Z=e_{n-1}{2}^{n-1}+e_{n-2}{2}^{n-2}+\&CenterDot;\&CenterDot;\&CenterDot;+e_0\end{array}\right.$ 将x递增并向下取整，得出x=1，2，3，4，…时y的值：y_n=（x_b–x_a）（y–y_a）–（x–x_a）（y_b–y_a），x_n=1，2，3…，再对y_n值向下取整，所得（x₁，y₁），（x₂，y₂），（x₃，y₃），…是这条直线经过的第一次划分的正方体网格；

当网格A为该直线经过的正方形网格时，则利用该直线的斜率和A点与图中红色中心点C（x_C，y_C）连线的斜率k做比较来得出结论，即： $k\max =\left\{\begin{array}{c}\mathrm{kAB}=\tan \mathrm{\&alpha;}=(\mathrm{yb}-\mathrm{ya})/(\mathrm{xb}-\mathrm{xa})\\ \mathrm{kAC}=\tan \mathrm{\&beta;}=(\mathrm{yc}-\mathrm{ya})/(\mathrm{xc}-\mathrm{xa})\end{array}\right.,$ 其中：k_max是指：直线AB与直线AC斜率比较后，其中斜率值较大的那条直线的斜率值，当k_AB＞0，那么判断当k_AB＞k_AC那么直线经过网格A；当k_AB＜0，那么当k_AB＜k_AC那么直线经过网格A；利用这种方法能够快速的计算出直线经过的所有正方形网格；

推广到任意多边形所经过的网格也可以全部求出，这里所有经过的网格集合为：

其中：Meshid（）是指：直线经过的网格，R是指直线经过的所有网格所构成集合；

3.1.2获取面元边缘点的坐标集合：从集合R中任意取出一个多边形经过的正方形网格，它的四个端点坐标分别为G₁（e，f），G₂（e+1，f），G₃（e，f+1），G₄（e+1，f+1），即能够得到所有多边形经过的正方形网格顶点，即边缘点的集合为：

$U=\left\{U\right|\mathrm{Gn}1(\mathrm{en},fn),\mathrm{Gn}2(\mathrm{en}+1,\mathrm{fn}),\mathrm{Gn}3(\mathrm{en},\mathrm{fn}+1),\mathrm{Gn}4(\mathrm{en}+1,\mathrm{fn}+1)\}$ ，其中：U对应于面元的边缘点的坐标集合，en表示多边形经过正方形网格的四个顶点的横坐标，fn表示多边形经过正方形网格的四个顶点的纵坐标，Gni表示多边形经过的所有正方形网格的四个顶点在投影面上的坐标（i=1、2、3、4…）；

3.1.3获取多边形的最小外接矩形：把每个面上的所有直线正交后所得的交点看作一个以网格端点的集合，对应最小外接矩形为一个以四个极点即左极点S₁（x₁，y₁），右极点S₂（x₂，y₂），上极点S₃（x₃，y₃），下极点S₄（x₄，y₄）的坐标为边缘的最小外接矩形，那么该最小外接矩形包括的正方形网格为： $Z=\left\{Z\right|\mathrm{Meshid}(x_n,y_n),x_1\&le;x_n\&le;x_2,y_4\&le;y_n\&le;y_3\}.$

8.根据权利要求1所述的方法，其特征是，所述的网格扩充是指：对最小外接矩形进行网格扩充循环递归，然后用射线法得到包含在感兴趣区域内的所有面元。

9.根据权利要求8所述的方法，其特征是，所述的网格扩充循环递归的具体步骤为：当得到最小外接矩形后，通过改变网格间距可以将最小外接矩形划分成许多更小的正方形网格，取新的网格间距是初始网格间距的二分之一，通过坐标转换可以得到多边形的顶点在新坐标系下的坐标，由于网格细分，所以多边形端点更靠近细分后的网格边线；当完成一遍递归直至多边形的端点和优化外接矩形边缘重合，而多边形的边线在边缘网格内也将与边缘网格的边线重合时，得到优化外接矩形中多边形的轮廓线，其中：n=为递归次数。

10.根据权利要求1所述的方法，其特征是，所述的一致性提取是指：通过计算已经得到感兴趣体积的像素集合，得到包含在感兴趣多边形内的所有体素集合；利用任意相邻的两张四维数据之间都存在体素连续性，对感兴趣体积进行提取；其他的三维序列进行相同操作，得到其余图像相同选择区域所在的感兴趣体积，然后在相同位置提取出总的感兴趣体积集合为：Voi={V_n，1≤n≤N}，最后将提取得到的血管序列三维图像进行四维可视化；

所述的一致性提取的顺序和序列体数据的原有顺序保持一致。

Images