CN103190932B - 一种冠状动脉血管壁应力和应变的估算方法 - Google Patents

Abstract

一种冠状动脉血管壁应力和应变的估算方法，其技术方案是，所述方法首先从IVUS灰阶图像序列中获取血管壁组织随血流施压变化产生的运动场，然后在此基础上求得各帧血管壁横截面图像上各点的应力和应变值。本发明直接从常规IVUS灰阶图像序列中获得血管壁的应力和应变值，算法非常简单，降低了本方法的应用条件。与基于射频信号分析的方法相比，本发明能够计算每帧图像中血管横截面上的各点在任意方向上的应变值（包括两个主应变和一个剪应变），能更全面地反映血管壁和斑块组织的弹性特征，为分析冠状动脉粥样硬化的发病机理、评价粥样硬化斑块的易损性和预测其发展趋势等提供了更为可靠的依据。

Description

一种冠状动脉血管壁应力和应变的估算方法

技术领域

本发明涉及一种根据血管内超声图像序列估算由血流引入的冠状动脉血管壁应力和应变的方法，属于医学成像技术领域。

背景技术

冠状动脉粥样硬化斑块破裂和血栓形成是急性冠状动脉综合征的主要发病机制。斑块破裂的位置与斑块内及其附近区域的应力有关。冠状动脉血管壁和斑块的相关力学性质，可作为早期检测动脉粥样硬化和筛查不稳定斑块的评价指标。

血管内超声（intravascular ultrasound，IVUS）成像是目前临床广泛采用的诊断血管病变，尤其是冠心病以及检测易损斑块的微创介入影像技术。它提供了活体中的血管壁、管腔及斑块的横截面图像，根据动态图像序列还可获得心动周期中冠状动脉血管壁及斑块的变形情况。

目前，利用血管内超声测量动脉血管壁的应力和应变以及管壁弹性模量，多是通过对超声信号的深入分析（包括射频信号、背向散射信号及弹性特征等），了解组织特征。由于在将原始射频信号转化为灰阶图像的过程中损失了部分信息，所以理论上讲此类方法较之基于灰度图像处理的方法具有更高的分辨率和精度。但它计算的是一维径向应变，而实际中管壁和斑块组织的变形是三维的，包括径向、周向和纵向三个分量，所以无法全面显示管壁和斑块组织的弹性信息，这种局限性是该技术难以克服的。此外，只有在组织运动方向与射频信号方向一致时才能获得可靠的应变估计。由于射频信号低信噪比的影响，很难获得对较短的舒张末期区间中较小应变的估计。同时，IVUS基本的信号被转化为标准显示的常规灰度图像，并非所有IVUS成像系统都允许对原始射频信号进行采集，因而限制了此类方法的临床广泛应用。

发明内容

本发明的目的在于针对现有技术之弊端，提供一种冠状动脉血管壁应力和应变的估算方法，为冠心病的计算机辅助诊断、斑块易损性的评估及其发展趋势的预测提供依据。

本发明所述问题是以下述技术方案实现的：

一种冠状动脉血管壁应力和应变的估算方法，所述方法首先从IVUS灰阶图像序列中获取血管壁组织随血流施压变化产生的运动场，然后在此基础上求得各帧血管壁横截面图像上各点的应力和应变值，具体步骤如下：

a、从IVUS灰度图像序列中获取血管壁组织随血流施压变化产生的运动场：

对于IVUS灰度图像序列中一帧图像中的像素(x, y)，以其为中心的足够小的正方形邻域内的像素数为N，则该像素的位移(V _x , V _y )由以下方程组求得：

，

式中，I(x, y, t)是像素(x, y)在时刻t的灰度值；I _x 、I _y 和I _t 分别是I对x、y和t的偏导数：

，

，

，

如此遍历整幅图像，得到各像素的位移；

b、根据各帧图像中血管壁横截面上各点的位移得到其应力和应变：

一帧图像中像素(x, y)处的应变值由下式求得：

，

式中，为平面应变矢量，ε _x 、ε _y 为的两个主应变，γ _xy 为的剪应变；

像素(x, y)处的三个应力分量σ _x 、σ _y 和σ _xy 由下式求得：

，

式中，υ为血管壁的泊松比，E为血管壁的弹性模量。

上述冠状动脉血管壁应力和应变的估算方法，在从IVUS灰度图像序列中获取血管壁组织随血流施压变化产生的运动场之前，应先去除各帧IVUS图像中的背景区域，具体方法如下：

首先，将各帧横向视图转换到极坐标系中，得到极坐标视图；然后按照下式将极坐标视图中的背景灰度设置为255：

，

其中，r和θ分别是像素的极径和极角，I(r, θ)和I'(r, θ)分别是去除背景前后的极坐标视图中像素(r, θ)处的灰度值，ImageHeight为横向视图的高度（单位：像素），之后，再经极坐标逆变换，得到去除背景后的横向视图。

上述冠状动脉血管壁应力和应变的估算方法，计算像素(x, y)的位移(V _x , V _y )时，以像素(x, y)为中心的正方形邻域的大小为3×3（单位：像素），即正方形邻域内的像素数N=9。

本发明中的IVUS灰阶图像序列是采用目前商售的IVUS成像系统（例如Jomed Endosonic超声成像仪、Volcano血管内超声成像仪等）在进行常规心导管检查的过程中获取的。其成像原理是在X射线透视下将0.014英寸导引钢丝插至靶血管远端并固定之，然后沿导引钢丝将高频（20～50MHz）超声探头导管插入靶血管段远端，在马达的驱动下缓慢回撤导管至近端，同时对血管横截面进行360°成像，获得血管壁横截面的一系列超声图像。

本发明直接从常规IVUS灰阶图像序列中获得血管壁的应力和应变值，无需采集IVUS设备的原始射频信号和背向散射信号，无需其它影像资料（如X射线血管造影图像或CT血管造影图像）的辅助，亦无需对血管进行三维重建，因而算法非常简单，降低了本方法的应用条件。与基于射频信号分析的方法相比，本发明能够计算每帧图像中血管横截面上的各点在任意方向上的应变值（包括两个主应变和一个剪应变），能更全面地反映血管壁和斑块组织的弹性特征，为分析冠状动脉粥样硬化的发病机理、评价粥样硬化斑块的易损性和预测其发展趋势等提供了更为可靠的依据。

附图说明

下面结合附图对本发明作进一步说明。

图1是去除IVUS图像的背景区域示意图。其中(a)是原始横向视图；(b)是去除背景前的极坐标视图；(c)是去除背景后的极坐标视图；(d)是去除背景后的横向视图。

图2是连续四帧IVUS图像的运动场估计结果图。为了清晰起见，每隔四个像素显示一个运动向量，并将位移扩大2倍。

图3是一帧含钙化斑块的IVUS图像的应变分布图，采用伪彩编码表示估计出的应变分量值。其中(a)是原始图像；(b)是主应变ε _x 的分布图；(c)是主应变ε _y 的分布图；(d) 是剪应变γ _xy 的分布图。

文中所用符号清单为：r、θ、像素的极径和极角；I(r, θ)、I'(r, θ)、去除背景前、后的极坐标视图中像素(r, θ)处的灰度值；ImageHeight、横向视图的高度（单位：像素）；N、邻域内的像素数；V _x 、V _y 、邻域内的像素在x方向和y方向上的位移；I(x, y, t)、像素(x, y)在时刻t的灰度值；I _x 、I _y 、I _t 、I对x、y和t的偏导数；、平面应变矢量；ε _x 、ε _y 、的两个主应变；γ _xy 、的剪应变；σ _x 、σ _y 、σ _xy 、血管壁受到的三个应力分量；υ、血管壁的泊松比；E、血管壁的弹性模量。

具体实施方式

下面详细说明本发明的处理步骤：

(1)  去除图像背景：

IVUS图像序列的前后帧之间非常相似，背景区域基本相同，如附图1(a)所示。为了减少后续步骤的计算量，需将其去除。首先，将各帧横向视图转换到极坐标系中，得到极坐标视图（附图1(b)）；然后，由于背景区域位于极坐标视图的下方，因此按照下式将极坐标视图中的背景灰度设置为255（附图1(c)）：

                                (1)

其中，r和θ分别是像素的极径和极角；I(r, θ)和I'(r, θ)分别是去除背景前后的极坐标视图中像素(r, θ)处的灰度值；ImageHeight为横向视图的高度（单位：像素）。之后，再经极坐标逆变换，得到去除背景后的横向视图（附图1(d)）。

(2)  估计相邻帧之间血管壁的运动场：

IVUS图像采集速率的常用值为30帧/秒，因此对于一帧图像中的每个像素，假设以其为中心的足够小的正方形邻域内，相邻两帧图像之间的运动近似为线性的，即邻域内全部N个像素的位移是相同的，全部为(V _x , V _y )；将其代入光流约束方程：

                                     (2)

得到由N个方程组成的方程组：

                                      (3)

式中，I(x, y, t)是像素(x, y)在时刻t的灰度值；I _x 、I _y 和I _t 分别是I对x、y和t的偏导数：

     (4)

   (5)

  (6)

采用最小二乘法求解方程组(3)，得到邻域中心像素的位移(V _x , V _y )。如此遍历整幅图像，得到各像素的位移。

正方形邻域的大小是影响运动场估计结果的重要参数，采用较小的邻域可获得较高的估计精度，但运算量较大；反之，采用较大邻域，则所得结果的精度不及较小邻域，但运算量较小。本发明方法为了保证估计精度，采用像素级的最小邻域，即3×3（单位：像素）邻域，则N=9。

(3)  估算血管壁横截面上各点的二维应力和应变值：

相邻帧IVUS图像之间的冠状动脉血管壁位移主要由血管腔的不规则几何形态、管腔内复杂的血液流场以及周期性心脏运动三个因素共同作用产生。但是，在回撤超声导管采集图像序列的过程中，常用的导管回撤速度是0.5mm/s，帧采样率是30帧/秒，确保了高速的图像采集速率（即1mm的长度内采集60帧）。因而，由血管腔的不规则几何形态变化引起的血管壁位移远小于由血液流动及心脏运动引起的位移，可将其忽略不计。基于上述原因，本发明方法不再对步骤(2)中得到的血管壁位移加以区分，直接用于计算由血流引入的血管壁二维应力和应变值。

根据弹性力学知识，平面应变矢量包含两个主应变ε _x 和ε _y 以及一个剪应变γ _xy ，它们与位移之间的关系是：

                                                 (7)

将步骤(2)得到的血管壁横截面上各点的位移，即u=V _x 和v=V _y 代入到式(7)中，即得到其应变值。

弹性模量、应力和应变之间的关系为

              (8)

式中，σ _x 、σ _y 和σ _xy 分别是血管壁受到的三个应力分量，υ为血管壁的泊松比，E为血管壁的弹性模量。对于固定的血管材料，υ和E可视为常量。根据式(8)，即可由平面应变分量(ε _x , ε _y , γ _xy )得到三个应力分量(σ _x , σ _y , σ _xy )。

附图3是一帧含钙化斑块的IVUS图像以及用伪彩编码表示的三个应变分量的分布图，钙化斑块内可见到强回声，且其后有负性声影。由图可知，整幅图像上的应变分布不均匀，且以剪应变即压缩应变为主，靠近血管壁内膜的组织应变大于外膜区域，斑块部位产生的应变大于正常组织，在斑块集中区域出现应力集中现象，这与临床的观察结果是相符的。

Claims (3)

1.一种冠状动脉血管壁应力和应变的估算方法，其特征是，该方法首先从IVUS灰阶图像序列中获取血管壁组织随血流施压变化产生的运动场，然后在此基础上求得各帧血管壁横截面图像上各点的应力和应变值，具体步骤如下：

a、从IVUS灰度图像序列中获取血管壁组织随血流施压变化产生的运动场：

对于IVUS灰度图像序列中一帧图像中的像素(x, y)，以其为中心的足够小的正方形邻域内的像素数为N，则该像素的位移(V _x , V _y )由以下方程组求得：

，

式中，I(x, y, t)是像素(x, y)在时刻t的灰度值；I _x 、I _y 和I _t 分别是I对x、y和t的偏导数：

，

，

，

如此遍历整幅图像，得到各像素的位移；

b、根据各帧图像中血管壁横截面上各点的位移得到其应力和应变：

一帧图像中像素(x, y)处的应变值由下式求得：

，

式中，为平面应变矢量，ε _x 、ε _y 为的两个主应变，γ _xy 为的剪应变；

像素(x, y)处的三个应力分量σ _x 、σ _y 和σ _xy 由下式求得：

，

式中，υ为血管壁的泊松比，E为血管壁的弹性模量。

2.根据权利要求1所述的冠状动脉血管壁应力和应变的估算方法，其特征是，在从IVUS灰度图像序列中获取血管壁组织随血流施压变化产生的运动场之前，应先去除各帧IVUS图像中的背景区域，具体方法如下：

首先，将各帧血管横截面图像转换到极坐标系中，得到极坐标视图；然后按照下式将极坐标视图中的背景灰度设置为255：

，

其中，r和θ分别是像素的极径和极角，I(r, θ)和I'(r, θ)分别是去除背景前后的极坐标视图中像素(r, θ)处的灰度值，ImageHeight为血管横截面图像的高度（单位：像素），之后，再经极坐标逆变换，得到去除背景后的血管横截面图像。

3.根据权利要求1或2所述的冠状动脉血管壁应力和应变的估算方法，其特征是，计算像素(x, y)的位移(V _x , V _y )时，以像素(x, y)为中心的正方形邻域的大小为3×3，单位为像素，即正方形邻域内的像素数N=9。