基于有限元形变理论的血管壁弹性分析方法

摘要

一种基于有限元形变理论的血管壁弹性分析方法，属于医学图像分析领域，通过计算IVUS(血管内超声)图像中的血管壁在一定压力下的应变特征值来分析其弹性。采用超声探头悬停在血管腔内某一位置采集的若干个心动周期的图像，并将图像输入到计算机中进行处理。包括：采用模板匹配法，在模板图像中搜索出原始图像中标记点的匹配点；根据血管壁的组织特征，将血管壁划分为若干层；提出一种基于有限元法分割计算血管壁弹性的方法；根据三角形微元应变的特征，将血管壁划分为若干相连的三角形；考虑到三角形在形变时的相互影响，将影响因子也加入到应变值的计算中。改变了现有血管壁弹性分析方法所存在的局限性，可以准确的刻画血管壁病变的特征，为病变分析创造了条件。

说明书

技术领域

本发明涉及计算机医学图像分析领域，特别涉及一种基于有限元形变理论的血管壁弹性分析方法，通过对血管内超声(IVUS)图像中血管壁部分的处理分析，判别出血管壁的相对弹性，用于血管壁病变的分析和早期预测。

背景技术

血管内超声(Intravascular Ultrasound，IVUS)成像是近年来临床应用于诊断血管病变的一种新的诊断方式，它可以显示血管内腔、管壁和粥样斑块的组织形态学特征以及斑块的病理组成，如钙化、纤维组织、脂质核心及斑块破裂、离体等如图1。它比只能提供动脉长轴影像的冠状动脉造影成像更容易发现血管的早期病变。目前IVUS正越来越广泛地被用于评价冠状动脉粥样硬化特点及程度，被认为是诊断冠心病新的“金标准”，并用于指导选择介入性治疗方法、评估介入治疗效果等。在IVUS临床应用中，可以获取一系列的超声图像供大夫判断辨别血管的病变。

动脉血管弹性功能是反映人体血管系统好坏的重要指标，分析血管壁弹性对早期识别易损斑块、防治心血管疾病具有重要意义。而弹性又以一定力作用下血管壁产生的应变来衡量，研究管壁应变的精确计算方法成为判定血管弹性的关键。现有的计算应变的方法叫做IVUSE(血管弹性图)，该技术是基于组织受压产生变形的弹性原理，用血管内超声导管收集不同压力作用下冠状动脉血管壁和斑块的射频回波信号，经局部置换建立反映组织受牵拉情况的横截面弹性图。然而，该方法求出的是径向应变值，而平面内血管壁在管内的压力下还产生切向的应力，并且某处的应变结果还受相邻组织的影响。因此，该计算方法具有一定的局限性。

有限元是集合在一起能够表示实际连续域的离散单元。有限元的概念早在几个世纪前就已产生并得到了应用，例如用多边形(有限个直线单元)逼近圆来求得圆的周长。有限元法最初被称为矩阵近似方法，应用于航空器的结构强度计算，并由于其方便性、实用性和有效性而引起从事力学研究的科学家的浓厚兴趣。经过短短数十年的努力，随着计算机技术的快速发展和普及，有限元方法迅速从结构工程强度分析计算扩展到几乎所有的科学技术领域，成为一种丰富多彩、应用广泛并且实用高效的数值分析方法。

有限元方法与其他求解边值问题近似方法的根本区别在于它的近似性仅限于相对小的子域中。20世纪60年代初首次提出结构力学计算有限元概念的克拉夫(Clough)教授形象地将其描绘为：“有限元法＝Rayleigh Ritz法+分片函数”，即有限元法是Rayleigh Ritz法的一种局部化情况。不同于求解满足整个定义域边界条件的允许函数的Rayleigh Ritz法，有限元法将函数定义在简单几何形状(如二维问题中的三角形或任意四边形)的单元域上(分片函数)，且不考虑整个定义域的复杂边界条件，这是有限元法优于其他近似方法的原因之一。

基于有限元形变血管壁弹性分析是以线段微元为基础，利用三角形形变的特征，求出血管壁在应变时出现的所有形变情况。同时，巧妙的结合微元线应变与微元面应变，以及微元面间应变的影响。在不增加附加设备的情况下，充分利用现有IVUS设备提供的超声图像信息，以血管内超声图像为依据，把管壁区域进行三角形单元划分。只要划分的单元足够小，就可以把三角形区域近似看作一个应变均匀的整体，采用经典有限变形力学理论，结合血管壁受力情况，建立管壁应变的数学模型。计算出各个三角形的应变数值，即可得到血管壁各处的应变情况。再通过颜色编码把三角形应变数据转换为256色应变图像，使得管壁的应变显示形象直观。

模板匹配是图像识别中的一种常见的方法。它的基本思想是：要判读一幅图像中是否存在某种已知的图案，则以一幅与该图案有着相同尺寸和方向的小图像为模板，通过一定的方法可以在图中找到目标，确定其坐标位置。

在原始图像中标记出血管壁内外膜点，标记出来的点都是一些独立的点坐标。独立的像素点只有灰度的特征，在模板图像上搜索匹配点时很难定义搜索条件。因此，利用模板匹配法将像素点放在一个适当的背景下，利用背景的特征在模板图像中搜索匹配背景的位置，从而找到匹配点的位置。

发明内容

本发明的目的在于，通过提供基于有限元形变理论的血管壁弹性分析方法，以实现对血管壁的弹性分析，从而更加精确且简便的求出血管壁的应变分量。

本发明是采用以下技术手段实现的：

一种基于有限元形变理论的血管壁弹性分析方法，通过计算血管内超声图像中的血管壁在一定压力下的应变特征值来分析其弹性，采用超声探头悬停在血管腔内某一位置采集的若干个心动周期的图像，并将图像输入到计算机中进行处理等操作；采用模板匹配法，在模板图像中搜索出与原始图像中标记的点相匹配的点的位置；根据血管壁的组织特征，将血管壁划分为若干层；提出一种基于有限元法分割计算血管壁弹性的方法；根据三角形微元应变的特征，将血管壁划分为若干相连的三角形；将影响因子加入到应变值的计算中；

上述计算机处理通过超声机接口软件对血管内超声图像进行读/写处理，进入处理图像应变计算模块计算，并经显示器显示结果，包括下述步骤：

步骤1、计算机从超声机接口读入血管内超声图像信号，并保存在内存中；

步骤2、在保存在内存中的图像中，选择同一心动周期内不同血压下的两幅图像，一幅定义为原始图像，另一幅为模板图像；

步骤3、统计血管内超声图像中血管壁内外膜的特征，在原始图像中分别标记出内外膜点，要求内膜点的个数与外膜点的个数相同，且疏密程度相同；

步骤4、在模板图像中，使用模板匹配法搜索与原始图像中标记的点相匹配的点，模板匹配法的公式为：

$C(i,j)=\frac{\underset{m=1}{\overset{M}{\Sigma }}\underset{n=1}{\overset{N}{\Sigma }}{S}^{\mathrm{ij}}(m,n)\times T(m,n)}{\sqrt{\left\{\underset{m=1}{\overset{M}{\Sigma }}\underset{n=1}{\overset{N}{\Sigma }}\left[S^{\mathrm{ij}}{(m,n)}^2\right]\right\}}\times \sqrt{\left\{\underset{m=1}{\overset{M}{\Sigma }}\underset{n=1}{\overset{N}{\Sigma }}\left[T{(m,n)}^2\right]\right\}}}$

式中T(m，n)表示模板中的某个像素值；S(m，n)表示模板覆盖下的那块区域中的某个像素值：C(i，j)表示模板与模板覆盖区域的相似度。

前述的进入血管壁应变计算模块，包括以下步骤：

步骤21、根据步骤3标记的点，将血管壁标记点中临近的内膜与外膜点划分为一组。按照插值公式将血管壁划分问为若干层，插值公式为：

$\mathrm{arrsum}(x,y)=\mathrm{arrInner}(x,y)+\frac{\mathrm{num}-1}{\mathrm{nLayer}}[\mathrm{arrOuter}(x,y)-\mathrm{arrInner}(x,y)]$

式中nLayer代表血管壁划分的层数；arrsum表示包括内外膜在内的所有点序列；arrInner表示内膜点序列；arrOuter表示外膜点序列；num表示点的序列号；

步骤22、根据内外膜点的组合，将血管壁划分为相连的三角形，每个三角形与两个三角形相连；

步骤23、根据微元线段的形变公式：

ds^′2-ds²＝2ε_xxdx²+4ε_xydxdy+2ε_yydy²

式中ds、ds′表示微元线段变化前后的平方值。分量ε_xx、ε_yy、ε_xy通称为Green应变张量ε_ij。其中，ε_xx和ε_yy为正应变，分别表示变形前后沿x轴和y轴的微元线段的形态变化，ε_xy为角应变，表示微元线段的方向变化，这三个量代表了血管壁在压力作用下不同方向上产生的应变。设原始三角形ABC，形变后的匹配三角形为A′B′C′，它们对应点的坐标分别为(x_a，y_a)，(x_b，y_b)，(x_c，y_c)和(x_a′，y_a′)，(x_b′，y_b′)，(x_c′，y_c′)；计算出整个三角形的应变值，联立单个三角形三条边的应变公式如下：

$\left(\begin{array}{c}\frac{1}{2}({\overline{A\prime B\prime }}^2-{\overline{\mathrm{AB}}}^2)={ϵ}_{\mathrm{xx}}{(x_b-x_a)}^2+2{\gamma }_{\mathrm{xy}}(x_b-x_a)(y_b-y_a)+{ϵ}_{\mathrm{yy}}{(y_b-y_a)}^2\\ \frac{1}{2}({\overline{B\prime C\prime }}^2-{\overline{\mathrm{BC}}}^2)={ϵ}_{\mathrm{xx}}{(x_c-x_b)}^2+2{\gamma }_{\mathrm{xy}}(x_c-x_b)(y_c-y_b)+{ϵ}_{\mathrm{yy}}{(y_c-y_b)}^2\\ \frac{1}{2}({\overline{A\prime C\prime }}^2-{\overline{\mathrm{AC}}}^2)={ϵ}_{\mathrm{xx}}{(x_c-x_a)}^2+2{\gamma }_{\mathrm{xy}}(x_c-x_a)(y_c-y_a)+{ϵ}_{\mathrm{xy}}{(y_c-y_a)}^2\end{array}\right)$

按照所述的三角形的应变值的大小，用颜色填充到每个单元内。

前述的颜色填充由红到黄到绿的顺序排列成应变图，由静态填充和/或动态填充组成。

前述的静态填充为在程序中设定256个颜色等级，无应变的三角形单元的应变值为0，颜色为黄色；大于0逐渐变红；小于0逐渐变绿。

前述的动态填充根据最大与最小值的差值将颜色分为256个等级，从大到小的颜色变化是由红到黄、由黄再到绿，没有固定的0值颜色。

本发明一种基于有限元形变理论的血管壁弹性分析方法，与现有技术相比，具有以下明显的优势和有益效果：

本发明一种基于有限元形变理论的血管壁弹性分析方法，改变了现有血管壁弹性分析方法所存在的局限性，可以准确的刻画血管壁病变的特征，为病变分析创造了条件。

附图说明

图1是血管内超声示意图；

图2是血管壁应变分析系统框图；

图3是本发明方法中主程序流程图；

图4是本发明方法中血管壁标记点示意图；

图中红色十字表示原始点坐标，绿色表示与红色原始点匹配的点坐标；

图5a是本发明方法中血管壁单层划分示意图；

图5b是本发明方法中血管壁三层划分示意图；

图6a是本发明中的一幅彩色血管壁水平方向应变图像；

图6b是本发明中的一幅彩色血管壁垂直方向应变示意图；

图6c是本发明中的一幅彩色血管壁剪切方向应变示意图。

具体实施方式

本发明一种基于有限元形变理论的血管壁弹性分析方法，其基本在于以下几点：

根据血管内超声图像特征的统计分析，标记出血管壁的内外膜点。标记点时需要选择同一心动周期内两幅同一位置不同压力下的图像，分别定义为原始图像和模板图像；

采用模板匹配算法，在模板图像上搜索与原始图像中标记出来的点相匹配的点位置；

提出一种基于有限元的血管壁弹性分析方法；

提出基于线段微元应变的三角形应变计算方法，将血管壁划分为若干三角形微元；

将计算后的结果通过颜色填充标记在原始图像上，显示血管壁弹性图。

本发明的技术方案参见图2、图3所示，基于有限元形变理论的血管壁弹性分析方法，是由血管内超声机采集图像，并将超声信号转换为图像电信号输入到计算机进行处理等操作，该方法特征在于计算机处理主要是通过超声机接口软件对血管内超声图像进行读/写处理，对处理后的图像进行应变计算，并经显示器显示结果，它依次包括下述步骤：

a、计算机从超声机接口读入血管内超声图像信号，并保存在内存中；

b、在保存在内存中的图像中，选择同一心动周期内不同血压下的两幅图像。一幅定义为原始图像，另一幅为模板图像；

c、统计血管内超声图像中血管壁内外膜的特征，在原始图像中分别标记出内外膜点。要求内膜点的个数与外膜点的个数相同，且疏密程度大致相同。

d、在模板图像中，使用模板匹配法搜索与原始图像中标记的点相匹配的点，模板匹配法的公式为：

$C(i,j)=\frac{\underset{m=1}{\overset{M}{\Sigma }}\underset{n=1}{\overset{N}{\Sigma }}{S}^{\mathrm{ij}}(m,n)\times T(m,n)}{\sqrt{\left\{\underset{m=1}{\overset{M}{\Sigma }}\underset{n=1}{\overset{N}{\Sigma }}\left[S^{\mathrm{ij}}{(m,n)}^2\right]\right\}}\times \sqrt{\left\{\underset{m=1}{\overset{M}{\Sigma }}\underset{n=1}{\overset{N}{\Sigma }}\left[T{(m,n)}^2\right]\right\}}}$

式中T(m，n)表示模板中的某个像素值；S(m，n)表示模板覆盖下的那块区域中的某个像素值；C(i，j)表示模板与模板覆盖区域的相似度；

上述步骤a进入血管壁应变计算模块，分为以下几步：

①根据步骤c标记的点，将血管壁标记点中临近的内膜与外膜点划分为一组。按照插值公式将血管壁划分问为若干层，插值公式为：

$\mathrm{arrsum}(x,y)=\mathrm{arrInner}(x,y)+\frac{\mathrm{num}-1}{\mathrm{nLayer}}[\mathrm{arrOuter}(x,y)-\mathrm{arrInner}(x,y)]$

式中nLayer代表血管壁划分的层数；arrsum表示包括内外膜在内的所有点序列；arrInner表示内膜点序列；arrOuter表示外膜点序列；num表示点的序列号；

②根据内外膜点的组合，将血管壁划分为相连的三角形，每个三角形与两个三角形相连；

③根据微元线段的形变公式：

ds^′2-ds²＝2ε_xxdx²+4ε_xydxdy+2ε_yydy²

式中ds、ds′表示微元线段变化前后的平方值。分量ε_xx、ε_yy、ε_xy通称为Green应变张量ε_ij。其中，ε_xx和ε_yy为正应变，分别表示变形前后沿x轴和y轴的微元线段的形态变化，ε_xy为角应变，表示微元线段的方向变化，这三个量代表了血管壁在压力作用下不同方向上产生的应变。设原始三角形ABC，形变后的匹配三角形为A′B′C′，它们对应点的坐标分别为(x_a，y_a)，(x_b，y_b)，(x_c，y_c)和(x_a′，y_a′)，(x_b′，y_b′)，(x_c′，y_c′)。计算出整个三角形的应变值，联立单个三角形三条边的应变公式如下：

$\left(\begin{array}{c}\frac{1}{2}({\overline{A\prime B\prime }}^2-{\overline{\mathrm{AB}}}^2)={ϵ}_{\mathrm{xx}}{(x_b-x_a)}^2+2{\gamma }_{\mathrm{xy}}(x_b-x_a)(y_b-y_a)+{ϵ}_{\mathrm{yy}}{(y_b-y_a)}^2\\ \frac{1}{2}({\overline{B\prime C\prime }}^2-{\overline{\mathrm{BC}}}^2)={ϵ}_{\mathrm{xx}}{(x_c-x_b)}^2+2{\gamma }_{\mathrm{xy}}(x_c-x_b)(y_c-y_b)+{ϵ}_{\mathrm{yy}}{(y_c-y_b)}^2\\ \frac{1}{2}({\overline{A\prime C\prime }}^2-{\overline{\mathrm{AC}}}^2)={ϵ}_{\mathrm{xx}}{(x_c-x_a)}^2+2{\gamma }_{\mathrm{xy}}(x_c-x_a)(y_c-y_a)+{ϵ}_{\mathrm{xy}}{(y_c-y_a)}^2\end{array}\right)$

6)血管壁微元颜色填充。根据计算出每个应变分量的大小，按照颜色由红到黄到绿的顺序填充到每个单元内，称为应变图，如图6a、图6b、图6c所示。其中，考虑到应变图的效果，将颜色填充分为两中：静态填充、动态填充。

①静态填充是在程序中设定256个颜色等级，无应变的三角形单元的应变值为0，颜色为黄色；大于0逐渐变红；小于0逐渐变绿。

②动态填充是找出应变的最值，根据最大与最小值的差值将颜色分为256个等级，从大到小的颜色变化是由红到黄、由黄再到绿，没有固定的0值颜色。

本发明的效果见图5a、图5b、图6a、图6b图6c。可以看出采用基于有限元形变计算出来的血管壁应变图像，可以准确的刻画血管壁病变的特征，为病变分析创造了条件。

最后应说明的是：以上实施例仅用以说明本发明而并非限制本发明所描述的技术方案；因此，尽管本说明书参照上述的各个实施例对本发明已进行了详细的说明，但是，本领域的普通技术人员应当理解，仍然可以对本发明进行修改或等同替换；而一切不脱离发明的精神和范围的技术方案及其改进，其均应涵盖在本发明的权利要求范围当中。