CN105160294A - 实时mce序列图像心肌组织感兴趣区域的自动跟踪方法 - Google Patents

Abstract

实时MCE序列图像心肌组织感兴趣区域的自动跟踪方法，属医学超声序列图像分析处理技术领域，其特征在于将基于分级策略的医学图像弹性配准方法与卡尔曼滤波技术相结合，以实现实时MCE序列图像心肌ROI的自动跟踪。本发明方法的步骤为：选取心肌ROI；初始化；采用基于分级策略的医学图像弹性配准方法对序列图像进行配准；利用卡尔曼滤波技术对配准结果进行修正；在序列图像上标注心肌ROI。本发明能够可靠地实现实时MCE序列图像心肌ROI的自动跟踪，且跟踪方法有着较高的精度和准确性，可以有效地避免跟踪误差的帧间积累。本方法也适用于其他医学超声序列图像ROI的自动跟踪。

Description

实时MCE序列图像心肌组织感兴趣区域的自动跟踪方法

技术领域

本发明涉及一种医学超声序列图像感兴趣区域(ROI)的自动跟踪方法，尤其涉及一种实时心肌声学造影(实时MCE)序列图像心肌组织感兴趣区域(ROI)的自动跟踪方法，属医学超声序列图像分析处理技术领域。

背景技术

实时心肌声学造影(实时MCE)作为一种能够无创和定量评估心肌微循环的新技术，正在从基础研究走向临床实践。但在影像医学临床实践中，诸如超声探头位置或方位的细微变化、心脏跳动、呼吸、患者体动、脏器蠕动及其弹性形变等因素是不可能完全避免的，从而导致实时MCE序列图像的帧间失准，使得所选择的心肌ROI在实时MCE序列图像不同帧上的位置发生不同程度的漂移，严重影响了实时MCE定量分析技术对冠心病早期诊断的准确性和客观性。

据查，对于MCE序列图像的帧间失准问题，美国Viginia大学Jayaweera教授基于MCE序列图像像素值的相似性测度提出了一种配准算法(ARJayaweera,etal.,Methodforthequantificationofmyocardialperfusionduringmyocardialcontrasttwodimensionalechocardiography,JAmSocEchocardiogr,1994)，但该算法只适用于对舒张末期的左心室短轴图像进行配准。牛津大学Noble博士基于图像信号间的互相关性提出了一种配准算法(JANoble,etal.Automatednon-rigidalignmentofclinicalmyocardialcontrastechocardiographyimagesequences:comparisonwithmanualalignment.UltrasoundMedBiol,2002)，该算法也只是适用于MCE间歇触发显像方式，且算法的有效性更多地依赖于后续的人工调整。巴西Campinas大学Lope博士综合基于分割与基于像素属性的图像配准技术提出了一种半自动的MCE序列图像的配准算法(MLLopes,etal,Improvingmyocardialcontrastechocardiographyimagesregistrationbyincreasingimagessimilarity.WCU,2003)，该算法的配准精度和有效性取决于分割阶段人工选择的模板，无法实现配准过程的自动化。孙丰荣、王晓婧、贾晓波等人及其研究团队在以往进行了相关的工作(王晓婧，散斑跟踪技术在超声心动图图像分析处理中的应用研究，山东大学硕士学位论文，2007；贾晓波，超声心动图分析处理若干问题的应用研究，山东大学硕士学位论文，2009；孙丰荣等，一种心肌声学造影图像定量分析的方法，ZL200910020745.3,2011；李新彩等，基于散斑跟踪技术的实时心肌声学造影图像定量分析，航天医学与医学工程，2010)，其间分别采用了基于块匹配和光流场的散斑跟踪方法去解决实时MCE定量分析中心肌ROI的自动跟踪问题，但实时MCE序列图像并不能很好地满足块匹配方法的两个基本假设(同一块内的像素具有相同的运动矢量、块内的像素只是平移运动)和光流场方法的三个基本假设(亮度不变、梯度不变和平滑性)，致使方法的跟踪精度和准确性以及鲁棒性不够理想，跟踪误差的帧间积累大，且基于光流场的散斑跟踪方法计算复杂度较高。其后孙丰荣、王丽、高新建等人分别对医学图像弹性配准技术(王丽，医学图像配准技术应用研究，山东大学硕士学位论文，2012)和卡尔曼滤波技术(高新建，实时心肌声学造影定量分析若干图像处理问题的应用研究，山东大学硕士学位论文，2013)在所述心肌ROI自动跟踪问题上的应用可行性进行了初步探索，但两种技术各自的应用性能(尤指两种技术各自的精确性、鲁棒性和实用性)都尚不能满足临床应用需求。

发明内容

针对背景技术所述以往的方法需要人工调整、分析过程缺乏自动化、方法的精度和准确性差、跟踪误差的帧间积累大、鲁棒性和临床实用性低等弊端，本发明综合基于分级策略的医学图像弹性配准方法与卡尔曼滤波技术，提出了一种实时MCE序列图像心肌组织感兴趣区域的自动跟踪方法，以实现实时MCE序列图像心肌ROI的自动跟踪。

本发明的技术方案如下：

一种实时MCE序列图像心肌组织感兴趣区域的自动跟踪方法，由配置一台具有实时MCE定量分析功能的计算机工作站来实现,设实时MCE序列图像由N帧图像构成，分别记作I_k(k＝1,2,…,N)，其中I_k为该序列图像的第k帧图像，初始帧为I₁，该方法步骤如下：

1)由工作站操作医师在初始帧I₁上手动选取心肌ROI

在实时MCE序列图像的初始帧I₁上，由工作站操作医师手动选取椭圆形的心肌组织感兴趣区域，该椭圆形区域的长轴两个端点的位置坐标分别记作(x_a1,y_a1)₁和(x_a2,y_a2)₁，短轴两个端点的位置坐标分别记作(x_b1,y_b1)₁和(x_b2,y_b2)₁；将上述位置坐标组成为向量(X,Y)₁，并称该向量为心肌组织感兴趣区域位置向量，记作R₁，该向量是一个4×2的矩阵，其结构可表示为：

$R_1={(X,Y)}_1=\left(\begin{array}{c}{(x_{a1},y_{a1})}_1\\ {(x_{a2},y_{a2})}_1\\ {(x_{b1},y_{b1})}_1\\ {(x_{b2},y_{b2})}_1\end{array}\right);$

2)初始化

设定基于分级策略的图像配准的分块阈值δ；初始化卡尔曼滤波器的时变的虚拟系统噪声ξ_k和时变的虚拟观测噪声η_k为相互独立的白噪声；令k＝1；

3)采用基于分级策略的医学图像弹性配准方法对I_k与I_k+1进行配准

以I_k为浮动图像、I_k+1为参考图像，采用基于分级策略的医学图像弹性配准方法，对I_k与I_k+1进行配准，获得非线性形变场，进而得到心肌组织感兴趣区域位置向量的配准向量；这里，心肌组织感兴趣区域位置向量(X,Y)_k记作R_k，其配准向量记作且心肌组织感兴趣区域位置向量R_k的结构可表示为：

$R_k={(X,Y)}_k=\left(\begin{array}{c}{(x_{a1},y_{a1})}_k\\ {(x_{a2},y_{a2})}_k\\ {(x_{b1},y_{b1})}_k\\ {(x_{b2},y_{b2})}_k\end{array}\right)$

所述采用基于分级策略的医学图像弹性配准方法对I_k与I_k+1进行配准的具体步骤如下：

(1)对I_k和I_k+1进行全局的刚体配准，得到两幅图像的粗略配准结果；

(2)将参考图像I_k+1划分为多个不重叠的子块，分块后的子图像分别与完整的浮动图像I_k进行局部的刚体配准，配准完成后对参考图像的子图像再分块，如此进行逐步细化的局部刚体配准，当参考图像的子图像的大小小于设定的阈值δ时，停止分块，并返回所述参考图像逐级分块配准过程中每个层级的每个子图像所对应的空间变换的平移和旋转参量；

(3)对返回的一系列空间变换参数进行插值运算，得到平滑的非线性形变场；

(4)将步骤(3)中得到的非线性形变场作用于心肌组织感兴趣区域位置向量R_k得其配准向量其结构可表示为：

${\tilde{R}}_{k+1}={(X,Y)}_{k+1}^R=\left(\begin{array}{c}{(x_{a1},y_{a1})}_{k+1}^R\\ {(x_{a2},y_{a2})}_{k+1}^R\\ {(x_{b1},y_{b1})}_{k+1}^R\\ {(x_{b2},y_{b2})}_{k+1}^R\end{array}\right)$

其中，分别为四个位置坐标(x_a1,y_a1)_k,(x_a2,y_a2)_k,(x_b1,y_b1)_k,(x_b2,y_b2)_k经形变场作用后的配准位置；

4)利用卡尔曼滤波技术对配准结果进行修正

若k＝1，令转到步骤5)；否则，依据如下所示的卡尔曼滤波器模型对配准向量进行修正得

$\left\{\begin{array}{c}{\hat{R}}_{k+1}={\hat{R}}_k+{\mathrm{\ξ}}_{k+1}\\ {\tilde{R}}_k={\hat{R}}_k+{\mathrm{\η}}_{k+1}\end{array}\right.,(k=2,3,4,\mathrm{...})$

其中，卡尔曼滤波器的最优估计为的修正值；为滤波器的观测值；ξ_k和η_k分别为时变的虚拟系统噪声和时变的虚拟观测噪声，这里，向量的结构可表示为：

${\hat{R}}_{k+1}={(\hat{X},\hat{Y})}_{k+1}=\left(\begin{array}{c}{({\hat{x}}_{a1},{\hat{y}}_{a1})}_{k+1}\\ {({\hat{x}}_{a2},{\hat{y}}_{a2})}_{k+1}\\ {({\hat{x}}_{b1},{\hat{y}}_{b1})}_{k+1}\\ {({\hat{x}}_{b2},{\hat{y}}_{b2})}_{k+1}\end{array}\right);$

其中，为配准向量经上述卡尔曼滤波技术修正后的心肌组织感兴趣区域位置向量，分别为的四个位置坐标 经上述卡尔曼滤波技术修正后的位置坐标；

5)若k＜N，令k＝k+1且转到步骤3)，否则转入步骤6)；

6)在实时MCE序列图像上标注心肌组织感兴趣区域

分别以为节点进行抛物拟合，其中n＝1,2,…,N，得拟合曲线分别以节点进行抛物拟合，得拟合曲线由和包围形成的封闭区域即为在I_n帧上得到的心肌组织感兴趣区域，记作ROI_n；由计算机在实时MCE序列图像上自动标注ROI_n，n＝1,2,…,N。

上述的ROI是英文RegionOfInterest的缩写，意指图像处理中的术语“感兴趣区域”。

本发明的有益效果是：可靠地实现了实时MCE序列图像心肌ROI的自动跟踪；且所述方法有着较高的跟踪精度和准确性，可以有效地避免跟踪误差的帧间积累。同时，本发明方法也适用于其他医学超声序列图像ROI的自动跟踪。

附图说明

图1是本发明方法的流程框图，S1～S6为其各个步骤。

具体实施方式

下面结合附图和实施例对本发明作进一步说明，但不限于此。

实施例：

本发明实施例如图1所示，一种实时MCE序列图像心肌组织感兴趣区域的自动跟踪方法，由配置一台具有实时MCE定量分析功能的计算机工作站来实现,设实时MCE序列图像由N帧图像构成，分别记作I_k(k＝1,2,…,N)，其中I_k为该序列图像的第k帧图像，初始帧为I₁，该方法步骤如下：

S1)由工作站操作医师在初始帧I₁上手动选取心肌ROI

在实时MCE序列图像的初始帧I₁上，由工作站操作医师手动选取椭圆形的心肌组织感兴趣区域，该椭圆形区域的长轴两个端点的位置坐标分别记作(x_a1,y_a1)₁和(x_a2,y_a2)₁，短轴两个端点的位置坐标分别记作(x_b1,y_b1)₁和(x_b2,y_b2)₁；将上述位置坐标组成为向量(X,Y)₁，并称该向量为心肌组织感兴趣区域位置向量，记作R₁，该向量是一个4×2的矩阵，其结构可表示为：

$R_1={(X,Y)}_1=\left(\begin{array}{c}{(x_{a1},y_{a1})}_1\\ {(x_{a2},y_{a2})}_1\\ {(x_{b1},y_{b1})}_1\\ {(x_{b2},y_{b2})}_1\end{array}\right);$

S2)初始化

设定基于分级策略的图像配准的分块阈值δ；初始化卡尔曼滤波器的时变的虚拟系统噪声ξ_k和时变的虚拟观测噪声η_k为相互独立的白噪声；令k＝1；

S3)采用基于分级策略的医学图像弹性配准方法对I_k与I_k+1进行配准

以I_k为浮动图像、I_k+1为参考图像，采用基于分级策略的医学图像弹性配准方法，对I_k与I_k+1进行配准，获得非线性形变场，进而得到心肌组织感兴趣区域位置向量的配准向量；这里，心肌组织感兴趣区域位置向量(X,Y)_k记作R_k，其配准向量记作且心肌组织感兴趣区域位置向量R_k的结构可表示为：

$R_k={(X,Y)}_k=\left(\begin{array}{c}{(x_{a1},y_{a1})}_k\\ {(x_{a2},y_{a2})}_k\\ {(x_{b1},y_{b1})}_k\\ {(x_{b2},y_{b2})}_k\end{array}\right)$

所述采用基于分级策略的医学图像弹性配准方法对I_k与I_k+1进行配准的具体步骤如下：

(1)对I_k和I_k+1进行全局的刚体配准，得到两幅图像的粗略配准结果；

(2)将参考图像I_k+1划分为多个不重叠的子块，分块后的子图像分别与完整的浮动图像I_k进行局部的刚体配准，配准完成后对参考图像的子图像再分块，如此进行逐步细化的局部刚体配准，当参考图像的子图像的大小小于设定的阈值δ时，停止分块，并返回所述参考图像逐级分块配准过程中每个层级的每个子图像所对应的空间变换的平移和旋转参量；

(3)对返回的一系列空间变换参数进行插值运算，得到平滑的非线性形变场；

(4)将步骤(3)中得到的非线性形变场作用于心肌组织感兴趣区域位置向量R_k得其配准向量其结构可表示为：

${\tilde{R}}_{k+1}={(X,Y)}_{k+1}^R=\left(\begin{array}{c}{(x_{a1},y_{a1})}_{k+1}^R\\ {(x_{a2},y_{a2})}_{k+1}^R\\ {(x_{b1},y_{b1})}_{k+1}^R\\ {(x_{b2},y_{b2})}_{k+1}^R\end{array}\right)$

其中，分别为四个位置坐标(x_a1,y_a1)_k,(x_a2,y_a2)_k,(x_b1,y_b1)_k,(x_b2,y_b2)_k经形变场作用后的配准位置；

S4)利用卡尔曼滤波技术对配准结果进行修正

若k＝1，令转到步骤S5)；否则，依据如下所示的卡尔曼滤波器模型对配准向量进行修正得

$\left\{\begin{array}{c}{\hat{R}}_{k+1}={\hat{R}}_k+{\mathrm{\ξ}}_{k+1}\\ {\tilde{R}}_k={\hat{R}}_k+{\mathrm{\η}}_{k+1}\end{array}\right.,(k=2,3,4,\mathrm{...})$

其中，卡尔曼滤波器的最优估计为的修正值；为滤波器的观测值；ξ_k和η_k分别为时变的虚拟系统噪声和时变的虚拟观测噪声，这里，向量的结构可表示为：

${\hat{R}}_{k+1}={(\hat{X},\hat{Y})}_{k+1}=\left(\begin{array}{c}{({\hat{x}}_{a1},{\hat{y}}_{a1})}_{k+1}\\ {({\hat{x}}_{a2},{\hat{y}}_{a2})}_{k+1}\\ {({\hat{x}}_{b1},{\hat{y}}_{b1})}_{k+1}\\ {({\hat{x}}_{b2},{\hat{y}}_{b2})}_{k+1}\end{array}\right);$

其中，为配准向量经上述卡尔曼滤波技术修正后的心肌组织感兴趣区域位置向量，分别为的四个位置坐标 经上述卡尔曼滤波技术修正后的位置坐标；

S5)若k＜N，令k＝k+1且转到步骤S3)，否则转入步骤S6)；

S6)在实时MCE序列图像上标注心肌组织感兴趣区域

分别以为节点进行抛物拟合，其中n＝1,2,…,N，得拟合曲线分别以节点进行抛物拟合，得拟合曲线由和包围形成的封闭区域即为在I_n帧上得到的心肌组织感兴趣区域，记作ROI_n；由计算机在实时MCE序列图像上自动标注ROI_n，n＝1,2,…,N。

Claims (1)

1.一种实时MCE序列图像心肌组织感兴趣区域的自动跟踪方法，由配置一台具有实时MCE定量分析功能的计算机工作站来实现,设实时MCE序列图像由N帧图像构成，分别记作I_k(k＝1,2,…,N)，其中I_k为该序列图像的第k帧图像，初始帧为I₁，该方法步骤如下：

1)由工作站操作医师在初始帧I₁上手动选取心肌ROI

在实时MCE序列图像的初始帧I₁上，由工作站操作医师手动选取椭圆形的心肌组织感兴趣区域，该椭圆形区域的长轴两个端点的位置坐标分别记作(x_a1,y_a1)₁和(x_a2,y_a2)₁，短轴两个端点的位置坐标分别记作(x_b1,y_b1)₁和(x_b2,y_b2)₁；将上述位置坐标组成为向量(X,Y)₁，并称该向量为心肌组织感兴趣区域位置向量，记作R₁，该向量是一个4×2的矩阵，其结构可表示为：

$R_1={(X,Y)}_1=\left(\begin{array}{c}{(x_{a1},y_{a1})}_1\\ {(x_{a2},y_{a2})}_1\\ {(x_{b1},y_{b1})}_1\\ {(x_{b2},y_{b2})}_1\end{array}\right);$

2)初始化

设定基于分级策略的图像配准的分块阈值δ；初始化卡尔曼滤波器的时变的虚拟系统噪声ξ_k和时变的虚拟观测噪声η_k为相互独立的白噪声；令k＝1；

3)采用基于分级策略的医学图像弹性配准方法对I_k与I_k+1进行配准

以I_k为浮动图像、I_k+1为参考图像，采用基于分级策略的医学图像弹性配准方法，对I_k与I_k+1进行配准，获得非线性形变场，进而得到心肌组织感兴趣区域位置向量的配准向量；这里，心肌组织感兴趣区域位置向量(X,Y)_k记作R_k，其配准向量记作且心肌组织感兴趣区域位置向量R_k的结构可表示为：

$R_k={(X,Y)}_k=\left(\begin{array}{c}{(x_{a1},y_{a1})}_k\\ {(x_{a2},y_{a2})}_k\\ {(x_{b1},y_{b1})}_k\\ {(x_{b2},y_{b2})}_k\end{array}\right)$

所述采用基于分级策略的医学图像弹性配准方法对I_k与I_k+1进行配准的具体步骤如下：

(1)对I_k和I_k+1进行全局的刚体配准，得到两幅图像的粗略配准结果；

(2)将参考图像I_k+1划分为多个不重叠的子块，分块后的子图像分别与完整的浮动图像I_k进行局部的刚体配准，配准完成后对参考图像的子图像再分块，如此进行逐步细化的局部刚体配准，当参考图像的子图像的大小小于设定的阈值δ时，停止分块，并返回所述参考图像逐级分块配准过程中每个层级的每个子图像所对应的空间变换的平移和旋转参量；

(3)对返回的一系列空间变换参数进行插值运算，得到平滑的非线性形变场；

(4)将步骤(3)中得到的非线性形变场作用于心肌组织感兴趣区域位置向量R_k得其配准向量其结构可表示为：

${\tilde{R}}_{k+1}={(X,Y)}_{k+1}^R=\left(\begin{array}{c}{(x_{a1},y_{a1})}_{k+1}^R\\ {(x_{a2},y_{a2})}_{k+1}^R\\ {(x_{b1},y_{b1})}_{k+1}^R\\ {(x_{b2},y_{b2})}_{k+1}^R\end{array}\right)$

其中，分别为四个位置坐标(x_a1,y_a1)_k,(x_a2,y_a2)_k,(x_b1,y_b1)_k,(x_b2,y_b2)_k经形变场作用后的配准位置；

4)利用卡尔曼滤波技术对配准结果进行修正

若k＝1，令转到步骤5)；否则，依据如下所示的卡尔曼滤波器模型对配准向量进行修正得

$\left\{\begin{array}{c}{\hat{R}}_{k+1}={\hat{R}}_k+{\mathrm{\ξ}}_{k+1}\\ {\tilde{R}}_k={\hat{R}}_k+{\mathrm{\η}}_{k+1}\end{array}\right.,(k=2,3,4,...)$

其中，卡尔曼滤波器的最优估计为的修正值；为滤波器的观测值；ξ_k和η_k分别为时变的虚拟系统噪声和时变的虚拟观测噪声，这里，向量的结构可表示为：

${\hat{R}}_{k+1}={(\hat{X},\hat{Y})}_{k+1}=\left(\begin{array}{c}{({\hat{x}}_{a1},{\hat{y}}_{a1})}_{k+1}\\ {({\hat{x}}_{a2},{\hat{y}}_{a2})}_{k+1}\\ {({\hat{x}}_{b1},{\hat{y}}_{b1})}_{k+1}\\ {({\hat{x}}_{b2},{\hat{y}}_{b2})}_{k+1}\end{array}\right);$

其中，为配准向量经上述卡尔曼滤波技术修正后的心肌组织感兴趣区域位置向量，分别为的四个位置坐标 经上述卡尔曼滤波技术修正后的位置坐标；

5)若k＜N，令k＝k+1且转到步骤3)，否则转入步骤6)；

6)在实时MCE序列图像上标注心肌组织感兴趣区域

分别以为节点进行抛物拟合，其中n＝1,2,…,N，得拟合曲线(n＝1,2,…,N)；分别以节点进行抛物拟合，得拟合曲线(n＝1,2,…,N)；由和包围形成的封闭区域即为在I_n帧上得到的心肌组织感兴趣区域，记作ROI_n；由计算机在实时MCE序列图像上自动标注ROI_n，n＝1,2,…,N。