CN103247073A - 基于树状结构的三维脑血管模型构造方法 - Google Patents

Abstract

本发明提供一种基于树状结构的三维脑血管模型构造方法，首先，从CT或者MRA设备获取三维脑血管体数据场；采用分割算法分离脑血管和背景噪声；其次，计算脑血管骨架线；根据骨架线构建树状脑血管拓扑结构；再次，根据骨架线并采用弹性球算法计算每个控制点脑血管半径；最后，将构建的树状结构脑血管进行三维显示。本发明提供的构造方法符合脑血管的空间拓扑结构信息，具有血管显示精度高，结果误差小的优点。能够对脑血管病变区域进行检测，并能针对不同大小的显示窗口采用多尺度方式进行绘制。

Description

基于树状结构的三维脑血管模型构造方法

技术领域

本发明属于医学领域，具体涉及一种基于树状结构的三维血管模型构造方法。

背景技术

现有的血管造型技术大致可以分为两类：自由模型（model-free）的方法和基于模型（model-based）的方法。

自由模型的方法中最典型最常见的曲面重建方法是MC（Marching Cubes）算法。该方法通过选取一个合适的阈值并用线性插值的办法来计算出一个等值面从而将空间分为两个部分达到重建的目的。然而线性插值的使用和将空间按阈值分割成两部分的做法也过于简单，因此重建的效果并不理想。在重建之后需要通过平滑来消除曲面的锯齿效应，简单的拉普拉斯平滑方法会破坏细小的分枝。Taubin提出了低通滤波的方法，Vollmer对拉普拉斯平滑进行了改进，都取得了较好的效果。约束弹性曲面网（CESN）通过将平滑的初始曲面的顶点约束在它们所属的单元内较好地平衡了精确性和光滑性取得了理想的效果。

基于模型的方法假设血管是管状的结构，利用各种几何形状和构造方法来逼近血管的横截面来达到重构血管的目的，常用的是采用管状或者球状表示血管模型，该方法不仅定义了三维实体模型内的所有点，并精确定义了其中心线（骨架），十分有利于对模型进行实时操控、变形、演化，是一个十分适合用于表达构造血管这样的管状物体的构造方法。

发明内容

针对现存的方法并没有考虑血管拓扑结构、在处理过程中破坏细小分支的缺陷，本发明提出了一种基于树状结构的三维脑血管模型的构建方法，该模型采用树状结构符合脑血管的拓扑特性，采用球结合中心线方式表示单根血管，该描述方法能够对脑血管病变区域进行检测，并能针对不同大小的显示窗口采用多尺度方式进行绘制。

为实现上述目的，本发明采用了以下的技术方案：

基于树状结构的三维血管模型构造方法，包括以下步骤：

（1）从CT或者MRA设备获取三维脑血管体数据场；

（2）采用分割算法分离脑血管和背景噪声；

（3）计算脑血管骨架线；

（4）根据骨架线构建树状脑血管拓扑结构；

（5）根据骨架线采用弹性球算法计算每个节点脑血管半径；

（6）将构建的树状结构脑血管进行三维显示。

优选地，所述步骤（2）包括：

（2.1）采用高斯滤波对三维脑血管体数据场进行平滑处理，经过MIP投影得到MIP图像，采用二维OTSU算法借助MIP图像获得三维血管种子点；

（2.2）定义全局与局部信息相结合的区域增长规则，通过区域增长算法对血管进行粗分割，得到区域增长血管轮廓；

（2.3）采用Catt扩散模型对三维脑血管体数据场进行各项异性滤波，采用局部自适应C-V模型，将初步分割结果作为自适应活动轮廓模型的初始轮廓线进行二次分割。

优选地，所述步骤（3）中是采用Hessian矩阵的方法计算脑血管走势。

优选地，所述步骤（5）是通过弹性球运动力方程的构造，寻找弹性球中心线。

优选地，所述步骤（2.3）通过以下公式实现：

$E^{\mathrm{LCV}}(c_1,c_2,c_1^{\′},c_2^{\′},\mathrm{\φ})=\mathrm{\μ}{\∫}_{\mathrm{\Ω}}\mathrm{\δ}\left(\mathrm{\φ}(x,y)\right)|\▿\mathrm{\φ}(x,y)|\mathrm{dxdy}+V{\∫}_{\mathrm{\Ω}}H\left(\mathrm{\φ}(x,y)\right)\mathrm{dxdy}+$

${\mathrm{\λ}}_1{\∫}_{\mathrm{\Ω}}|u_0(x,y)-c_1{|}^{2}H\left(\mathrm{\φ}(x,y)\right)\mathrm{dxdy}+{\mathrm{\λ}}_2{\∫}_{\mathrm{\Ω}}{|u_0(x,y)-c_2|}^2(1-H\left(\mathrm{\φ}(x,y)\right))\mathrm{dxdy}+$

${\mathrm{\λ}}_1^{\′}{\∫}_{\mathrm{\Ω}}{|u_0(x,y)-c_1^{\′}|}^{2}H\left(\mathrm{\φ}(x,y)\right)\mathrm{dxdy}+{\mathrm{\λ}}_2^{\′}{\∫}_{\mathrm{\Ω}}{|u_0(x,y)-c_2^{\′}|}^2(1-H\left(\mathrm{\φ}(x,y)\right))\mathrm{dxdy}---\left(1\right)$

式（1）、式（2）中，

$c_1=\frac{{\∫}_{{\mathrm{\Ω}}_1}{u}_0(x,y)H\left(\mathrm{\φ}(x,y)\right)\mathrm{dxdy}}{{\∫}_{{\mathrm{\Ω}}_1}H\left(\mathrm{\φ}(x,y)\right)\mathrm{dxdy}}$ $c_2=\frac{{\∫}_{{\mathrm{\Ω}}_2}{u}_0(x,y)(1-H\left(\mathrm{\φ}(x,y)\right))\mathrm{dxdy}}{{\∫}_{{\mathrm{\Ω}}_2}(1-H\left(\mathrm{\φ}(x,y)\right))\mathrm{dxdy}}$

$c_1^{\′}=\frac{{\∫}_{{\mathrm{\Ω}}_1^{\′}}{u}_0(x,y)H\left(\mathrm{\φ}(x,y)\right)\mathrm{dxdy}}{{\∫}_{{\mathrm{\Ω}}_1^{\′}}H\left(\mathrm{\φ}(x,y)\right)\mathrm{dxdy}}$ $c_2^{\′}=\frac{{\∫}_{{\mathrm{\Ω}}_2^{\′}}{u}_0(x,y)(1-H\left(\mathrm{\φ}(x,y)\right))\mathrm{dxdy}}{{\∫}_{{\mathrm{\Ω}}_2^{\′}}(1-H\left(\mathrm{\φ}(x,y)\right))\mathrm{dxdy}}.$

优选地，所述步骤（3）所述获得脑血管中心线和半径信息的具体求取过程如下：

设置弹性力方程，决定弹性球中心走向：

$F(i,j,k)=\underset{m=i-r}{\overset{i+r}{\mathrm{\Σ}}}\underset{n=j-r}{\overset{j+r}{\mathrm{\Σ}}}\underset{p=k-r}{\overset{k+r}{\mathrm{\Σ}}}\mathrm{Seg}(m,n,k)---\left(3\right)$

式（3）中：

式（4）中，V(m,n,k)表示分割结果，Position(i,j,k)表示(i,j,k)位置像素的矢量位置，r表示当前计算的血管半径，其计算过程如下：

Step1：对于当前中心线点位端点，设置r=1；

Step2：计算弹性力；

Step3：如果弹性力为零，设置r=r+1，重复Step2，否则进行下一步；

Step4：沿力的矢量方向移动中心点，重新计算弹性力；

Step5：如果弹性力为零且出现双相切，搜索结束，否则下一步；

Step6：记录新修订的中心线以及半径信息，返回Step3。

优选地，所述脑血管树状拓扑结构构建方法如下：

采用节点方式表示每一个分支血管结构，数据结构定义如下：

Struct Node

{

父亲节点，

{中线点位置，半径，圆面法向量}，

左孩子节点，

右孩子节点

}

骨架上节点分类策略如下：

端点：骨架线开始或者终止的位置，端点的邻域关系为：只有一个邻接点；

普通点：构成骨架线的基本点，普通点周围存在两个临接点；

分叉点：骨架线的分支点，分叉点是二分结构，分叉点周围存在三个邻接点；

树状拓扑结构构建：

在上述节点分类的基础上，采用深度优先的侧率构建血管树状拓扑结构，算法如下：

Step1:将当前观察节点为起始点，直接建立一新树结构，并把当前点作为根节点存入；

Step2：进入下一个待观察节点；

Step3：当前观察点为普通点，加入到正在处理的节点结构中，返回Step2；否则转向Step4；

Step4：若观察点为分叉点，建立两个新节点分别作为该节点的左孩子节点和右边孩子节点，针对左孩子节点和右边孩子节点，分别递归进入Step2；

Step5：当处理完所有的节点后，血管树的拓扑结构构建完毕。

本发明所述的基于树状结构的三维脑血管模型构造方法，不仅能够有效分割脑血管粗大分支，而且还能精确提取脑血管细小结构。本发明构建的三维脑血管模型采用树状结构符合脑血管的拓扑特性，采用球结合中心线方式表示单根血管，该描述方法能够对脑血管病变区域进行检测，并能针对不同大小的显示窗口采用多尺度方式进行绘制。

附图说明

图1为本发明所述的基于树状结构的三维脑血管模型构造方法的流程图；

图2为分割算法分离脑血管和背景噪声的流程图；

图3所提取的脑血管骨架线的结果；

图4弹性球算法获取中心线上点对应的半径；

图5树状脑血管中每个节点数据的结构设计；

图6三维脑血管模型构造结果显示。

具体实施方式

以下结合附图对本发明的优选实施例进行详细描述，但本发明并不仅仅限于这些实施例。本发明涵盖任何在本发明的精髓和范围上做的替代、修改、等效方法以及方案。为了使公众对本发明有彻底的了解，在以下本发明优选实施例中详细说明了具体的细节，而对本领域技术人员来说没有这些细节的描述也可以完全理解本发明。

如图1所示，本发明提供的构造方法，包括如下的步骤：

步骤S601：获取医学图像的原始数据，从CT或者MRA设备获取三维脑血管体数据场；

步骤S602：采用分割算法将血管数据和背景分离，将血管数据与背景分离。步骤S602的具体实施方法如图2所示，首先采用高斯滤波对三维脑血管体数据场进行平滑处理，经过MIP投影得到MIP图像，采用二维OTSU算法借助MIP图像获得三维血管种子点；然后，采用自定义全局与局部信息相结合的区域增长规则，通过区域增长算法对血管进行粗分割，得到区域增长血管轮廓；最后，采用Catt扩散模型对三维脑血管体数据场进行各项异性滤波，采用局部自适应C-V模型，将初步分割结果作为自适应活动轮廓模型的初始轮廓线进行二次分割，获得脑血管骨架线如图3所示。所述局部自适应C-V模型参数设置如式（1）、式（2）所示，

$E^{\mathrm{LCV}}(c_1,c_2,c_1^{\′},c_2^{\′},\mathrm{\φ})=\mathrm{\μ}{\∫}_{\mathrm{\Ω}}\mathrm{\δ}\left(\mathrm{\φ}(x,y)\right)|\▿\mathrm{\φ}(x,y)|\mathrm{dxdy}+V{\∫}_{\mathrm{\Ω}}H\left(\mathrm{\φ}(x,y)\right)\mathrm{dxdy}+$

${\mathrm{\λ}}_1{\∫}_{\mathrm{\Ω}}|u_0(x,y)-c_1{|}^{2}H\left(\mathrm{\φ}(x,y)\right)\mathrm{dxdy}+{\mathrm{\λ}}_2{\∫}_{\mathrm{\Ω}}{|u_0(x,y)-c_2|}^2(1-H\left(\mathrm{\φ}(x,y)\right))\mathrm{dxdy}+$

${\mathrm{\λ}}_1^{\′}{\∫}_{\mathrm{\Ω}}{|u_0(x,y)-c_1^{\′}|}^{2}H\left(\mathrm{\φ}(x,y)\right)\mathrm{dxdy}+{\mathrm{\λ}}_2^{\′}{\∫}_{\mathrm{\Ω}}{|u_0(x,y)-c_2^{\′}|}^2(1-H\left(\mathrm{\φ}(x,y)\right))\mathrm{dxdy}---\left(1\right)$

其中，

$c_1=\frac{{\∫}_{{\mathrm{\Ω}}_1}{u}_0(x,y)H\left(\mathrm{\φ}(x,y)\right)\mathrm{dxdy}}{{\∫}_{{\mathrm{\Ω}}_1}H\left(\mathrm{\φ}(x,y)\right)\mathrm{dxdy}}$ $c_2=\frac{{\∫}_{{\mathrm{\Ω}}_2}{u}_0(x,y)(1-H\left(\mathrm{\φ}(x,y)\right))\mathrm{dxdy}}{{\∫}_{{\mathrm{\Ω}}_2}(1-H\left(\mathrm{\φ}(x,y)\right))\mathrm{dxdy}}$

$c_1^{\′}=\frac{{\∫}_{{\mathrm{\Ω}}_1^{\′}}{u}_0(x,y)H\left(\mathrm{\φ}(x,y)\right)\mathrm{dxdy}}{{\∫}_{{\mathrm{\Ω}}_1^{\′}}H\left(\mathrm{\φ}(x,y)\right)\mathrm{dxdy}}$ $c_2^{\′}=\frac{{\∫}_{{\mathrm{\Ω}}_2^{\′}}{u}_0(x,y)(1-H\left(\mathrm{\φ}(x,y)\right))\mathrm{dxdy}}{{\∫}_{{\mathrm{\Ω}}_2^{\′}}(1-H\left(\mathrm{\φ}(x,y)\right))\mathrm{dxdy}}$

轮廓曲线C将图像划分为内部区域Ω_in和外部区域Ω_out；以当前演变观察的像素位置为圆点，半径为r的球将Ω_in分为Ω₁,Ω'₁，将Ω_out分为Ω₂,Ω'₂。C₁、C₁'、C₂、C'₂分别表示区域Ω₁、Ω'₁、Ω₂、Ω'₂的均值。U₀(x,y)表示点（x,y）均值。

H(φ)的一个变体函数定义：

$\mathrm{H\ϵ}=\frac{1}{2}(1+\frac{2}{\mathrm{\π}}\arctan \left(\frac{\mathrm{\φ}}{\mathrm{\ϵ}}\right))---(3-16)$

实验结果表明本发明所述的算法不仅能够有效分割脑血管粗大分支，而且还能精确提取脑血管细小结构。

步骤S603：计算脑血管骨架线，获得脑血管中心线和半径信息，获取中心线上点对应的半径如图4所示，所使用的是弹性球算法，其具体求取过程如下：

设置弹性球弹性力方程，决定弹性球中心走向：

$F(i,j,k)=\underset{m=i-r}{\overset{i+r}{\mathrm{\Σ}}}\underset{n=j-r}{\overset{j+r}{\mathrm{\Σ}}}\underset{p=k-r}{\overset{k+r}{\mathrm{\Σ}}}\mathrm{Seg}(m,n,k)---\left(3\right)$

式（3）中，

r表示当前计算的血管半径，

式（4）中，V(m,n,k)表示分割结果，Position(i,j,k)表示(i,j,k)位置像素的矢量位置；其计算过程如下：

Step1：对于当前中心线点位端点，设置r=1；

Step2：计算弹性力；

Step3：如果弹性力为零，设置r=r+1，重复Step2，否则进行下一步；

Step4：沿力的矢量方向移动中心点，重新计算弹性力；

Step5：如果弹性力为零且出现双相切，搜索结束，否则下一步；

Step6：记录新修订的中心线以及半径信息，返回Step3。

步骤S604：根据骨架线构建树状脑血管拓扑结构。具体的，采用节点方式表示每一个分支血管的结构，如图5所示，其数据结构定义如下：

Struct Node

{

父亲节点，

{中线点位置，半径，圆面法向量}，

左孩子节点，

右孩子节点

}

将血管的节点分为三类：端点、普通点和分叉点，定义如下：

端点：位于骨架线开始或者终止的位置，端点的邻域关系为：只有一个邻接点。

普通点：构成骨架线的基本点，普通点周围存在两个临接点。

分叉点：骨架线的分支点，因为脑血管中很少存在多分叉结构，故一般情况下，分叉点是二分结构，分叉点有三个邻接点。

根据上述节点分类的基础上，采用深度优先的侧率构建血管树状拓扑结构，算法如下：

Step1:将当前观察节点为起始点，直接建立一新树结构，并把当前点作为根节点存入；

Step2：进入下一个待观察节点；

Step3：当前观察点为普通点，加入到正在处理的节点结构中，返回Step2；否则转向Step4；

Step4：若观察点为分叉点，建立两个新节点分别作为该节点的左孩子节点和右边孩子节点，针对左孩子节点和右边孩子节点，分别递归进入Step2；

Step5：当处理完所有的节点后，血管树的拓扑结构构建完毕。

步骤S605：根据骨架线并采用弹性球算法计算每个节点脑血管半径，参考具体步骤（5）的详细描述。

步骤S606：在计算机中将构建的血管模型可视化，显示结果如图6所示。

本发明优选实施例只是用于阐述本发明，且并没有详尽叙述所有的实施方式，本发明不限制该发明仅为所述的具体实施方式。显然，根据本说明书的内容，能够作很多的修改和变化。本说明书选取并具体描述这些实施例，是为了更好地解释本发明的原理和实际应用，从而使所属技术领域技术人员能很好地利用本发明。本发明仅受权利要求书及其全部范围和等效物的限制。

Claims (7)

1.基于树状结构的三维脑血管模型构造方法，其特征在于，包括以下步骤：

（1）从CT 或者MRA设备获取三维脑血管体数据场；

（2）采用分割算法分离脑血管和背景噪声；

（3）计算脑血管骨架线；

（4）根据骨架线构建树状脑血管拓扑结构；

（5）根据骨架线采用弹性球算法计算每个节点脑血管半径；

（6）将构建的树状结构脑血管进行三维显示。

2.如权利要求1所述的构造方法，其特征在于，所述步骤（2）包括：

（2.1）采用高斯滤波对三维脑血管体数据场进行平滑处理，经过MIP投影得到MIP图像，采用二维OTSU算法借助MIP图像获得三维血管种子点；

（2.2）定义全局与局部信息相结合的区域增长规则，通过区域增长算法对血管进行粗分割，得到区域增长血管轮廓；

（2.3）采用Catt扩散模型对三维脑血管体数据场进行各项异性滤波，采用局部自适应C-V模型，将初步分割结果作为自适应活动轮廓模型的初始轮廓线进行二次分割。

3.如权利要求1所述的构造方法，其特征在于，所述步骤（3）中是采用Hessian矩阵的方法计算脑血管走势。

4.如权利要求1所述的构造方法，其特征在于，所述步骤（5）是通过弹性球运动力方程的构造，寻找弹性球中心线。

5.如权利要求2所述的构造方法，其特征在于，所述步骤（2.3）通过以下公式实现：

$E^{\mathrm{LCV}}(c_1,c_2,c_1^{\′},c_2^{\′},\mathrm{\φ})=\mathrm{\μ}{\∫}_{\mathrm{\Ω}}\mathrm{\δ}\left(\mathrm{\φ}(x,y)\right)|\▿\mathrm{\φ}(x,y)|\mathrm{dxdy}+V{\∫}_{\mathrm{\Ω}}H\left(\mathrm{\φ}(x,y)\right)\mathrm{dxdy}+$

${\mathrm{\λ}}_1{\∫}_{\mathrm{\Ω}}|u_0(x,y)-c_1{|}^{2}H\left(\mathrm{\φ}(x,y)\right)\mathrm{dxdy}+{\mathrm{\λ}}_2{\∫}_{\mathrm{\Ω}}{|u_0(x,y)-c_2|}^2(1-H\left(\mathrm{\φ}(x,y)\right))\mathrm{dxdy}+$

${\mathrm{\λ}}_1^{\′}{\∫}_{\mathrm{\Ω}}{|u_0(x,y)-c_1^{\′}|}^{2}H\left(\mathrm{\φ}(x,y)\right)\mathrm{dxdy}+{\mathrm{\λ}}_2^{\′}{\∫}_{\mathrm{\Ω}}{|u_0(x,y)-c_2^{\′}|}^2(1-H\left(\mathrm{\φ}(x,y)\right))\mathrm{dxdy}---\left(1\right)$

式（1）、式（2）中：

$c_1=\frac{{\∫}_{{\mathrm{\Ω}}_1}{u}_0(x,y)H\left(\mathrm{\φ}(x,y)\right)\mathrm{dxdy}}{{\∫}_{{\mathrm{\Ω}}_1}H\left(\mathrm{\φ}(x,y)\right)\mathrm{dxdy}}$ $c_2=\frac{{\∫}_{{\mathrm{\Ω}}_2}{u}_0(x,y)(1-H\left(\mathrm{\φ}(x,y)\right))\mathrm{dxdy}}{{\∫}_{{\mathrm{\Ω}}_2}(1-H\left(\mathrm{\φ}(x,y)\right))\mathrm{dxdy}}$

$c_1^{\′}=\frac{{\∫}_{{\mathrm{\Ω}}_1^{\′}}{u}_0(x,y)H\left(\mathrm{\φ}(x,y)\right)\mathrm{dxdy}}{{\∫}_{{\mathrm{\Ω}}_1^{\′}}H\left(\mathrm{\φ}(x,y)\right)\mathrm{dxdy}}$ $c_2^{\′}=\frac{{\∫}_{{\mathrm{\Ω}}_2^{\′}}{u}_0(x,y)(1-H\left(\mathrm{\φ}(x,y)\right))\mathrm{dxdy}}{{\∫}_{{\mathrm{\Ω}}_2^{\′}}(1-H\left(\mathrm{\φ}(x,y)\right))\mathrm{dxdy}}.$

6.如权利要求3所述的构造方法，其特征在于，所述步骤（3）所述获得脑血管中心线和半径信息的具体求取过程如下：

设置弹性力矢量方程，决定弹性球中心走向：

$F(i,j,k)=\underset{m=i-r}{\overset{i+r}{\mathrm{\Σ}}}\underset{n=j-r}{\overset{j+r}{\mathrm{\Σ}}}\underset{p=k-r}{\overset{k+r}{\mathrm{\Σ}}}\mathrm{Seg}(m,n,k)$

式（3）中：

r表示当前计算的血管半径，

式（4）中，V(m,n,k)表示分割结果，Position(i,j,k)表示(i,j,k)位置像素的矢量位置，其计算过程如下：

Step1：对于当前中心线点位端点，设置r=1；

Step2：计算弹性力；

Step3：如果弹性力为零，设置r=r+1，重复Step2，否则进行下一步；

Step4：沿弹性力的矢量方向移动中心点，重新计算弹性力；

Step5：如果弹性力为零且出现双相切，搜索结束，否则下一步；

Step6：记录新修订的中心线以及半径信息，返回Step3。

7.如权利要求2所述的构造方法，其特征在于，所述树状脑血管拓扑结构的构建方法如下：

采用节点方式表示每一个分支血管结构，数据结构组织如下：

Struct Node

{

父亲节点，

{中线点位置，半径，圆面法向量}，

左孩子节点，

右孩子节点

}

骨架上节点分类如下：

端点：骨架线开始或者终止的位置，端点的邻域关系为：只有一个邻接点；

普通点：构成骨架线的基本点，普通点周围存在两个临接点；

分叉点：骨架线的分支点，分叉点是二分结构，分叉点周围存在三个邻接点；

树状拓扑结构构建：

在上述节点分类的基础上，采用深度优先的侧率构建血管树状拓扑结构，算法如下：

Step1:将当前观察节点为起始点，直接建立一新树结构，并把当前点作为根节点存入；

Step2：进入下一个待观察节点；

Step3：当前观察点为普通点，加入到正在处理的节点结构中，返回Step2；否则转向Step4；

Step4：若观察点为分叉点，建立两个新节点分别作为该节点的左孩子节点和右边孩子节点，针对左孩子节点和右边孩子节点，分别递归进入Step2；

Step5：当处理完所有的节点后，血管树的拓扑结构构建完毕。

Images