CN101625757A - 一种抑制冠状动脉内超声图像序列中刚性运动伪影的方法 - Google Patents

Abstract

一种抑制冠状动脉内超声图像序列中刚性运动伪影的方法，用于提高纵向视图的视觉效果、同时保证图像数据集合的完整性。其技术方案是：它是在建立血管壁刚性运动模型的基础上，首先采用三维自动分割技术从ICUS图像序列的各帧中分割出管腔轮廓，并计算相邻帧之间管腔轮廓的刚性运动参数(包括平移和旋转)，然后从中分离出刚性运动分量，最后通过对原始ICUS图像中的各像素点进行反向平移和旋转，使其刚性运动伪影得到补偿。同现有方法相比较，本发明不仅应用成本低、简单易行，而且可保证图像数据集合的完整性，可为冠状动脉血管病变的诊断和治疗提供更为丰富的信息。

Description

一种抑制冠状动脉内超声图像序列中刚性运动伪影的方法

技术领域

本发明涉及一种能够有效抑制冠状动脉内超声图像序列中由心脏运动和搏动的血流所造成的刚性运动伪影，同时保证图像数据集合完整性的方法，属医学成像技术领域。

背景技术

在进行冠状动脉内超声(Intracoronary Ultrasound，ICUS)成像时，由心脏运动和搏动的血流所造成的刚性运动伪影是影响ICUS图像纵向视图的视觉效果、准确定量分析和三维重建的主要因素。

目前，抑制ICUS成像运动伪影的方法主要有两种：一种方法是采用专用的ECG门控图像采集装置，只在每个心动周期中的相同相位处(一般是R波)采集ICUS图像，达到抑制运动伪影的目的。这种方法的缺点在于需要专门的ECG门控图像采集装置，并且与连续回撤超声导管采集图像相比，这种方式需要更长的采集时间(每个心动周期只采集一帧图像)，从而延长测试时间。另一种方法不采用ECG门控设备，而是连续回撤超声导管采集ICUS图像序列，同时记录ECG信号。待检查结束后，由医生根据同步记录的ECG信号从图像序列中选取在相同心脏相位处采集的图像。这种方式虽然不会延长测试的时间，但是其结果的客观性和可重复性差，受医生的专业知识和主观因素影响较大，而且受病人本身和不同病人之间的差异影响，很难选取出心动周期中的最佳采样时刻。同时，这两种门控方法都是每个心动周期只采集或选择一帧图像，因而可能会丢失很多对病变诊断和治疗有价值的信息。

发明内容

本发明的目的是克服现有技术的不足、提供一种既不需要专门的ECG门控图像采集装置，其结果也不受医生主观因素的影响，并且能保证图像数据集合完整性的抑制冠状动脉内超声图像序列中刚性运动伪影的方法。

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

一种抑制冠状动脉内超声图像序列中刚性运动伪影的方法，它是在建立血管壁刚性运动模型的基础上，首先采用三维自动分割技术从ICUS图像序列的各帧中分割出管腔轮廓，并计算相邻帧之间管腔轮廓的刚性运动参数(包括平移和旋转)，然后从中分离出刚性运动分量，最后通过对原始ICUS图像中的各像素点进行反向平移和旋转，使其刚性运动伪影得到补偿，具体步骤如下：

a、建立ICUS图像序列中血管壁的刚性运动模型：

ICUS图像序列中，相邻帧之间管腔轮廓的旋转和平移由刚性运动分量和几何分量构成：

$\left\{\begin{array}{c}\mathrm{\Δx}={\mathrm{\Δx}}_d+{\mathrm{\Δx}}_g\\ \mathrm{\Δy}={\mathrm{\Δy}}_d+{\mathrm{\Δy}}_g\\ \mathrm{\Δ\α}=\mathrm{\Δ}{\mathrm{\α}}_d+\mathrm{\Δ}{\mathrm{\α}}_g\end{array}\right.$

式中(Δx，Δy)是相邻帧之间管腔边界曲线重心之间的位移，Δα是相邻帧之间管腔边界曲线之间的旋转角，脚标d表示刚性运动分量，g表示几何分量；

b、采用三维自动分割技术从各帧ICUS图像中分割出管腔轮廓：

首先根据ICUS图像中管腔形状和灰度特征的先验知识，从ICUS纵向视图中分割出管腔的左右边界，并将其映射到各帧横截面ICUS图像中，得到每一帧中的初始管腔轮廓，然后采用快速推进法对初始管腔轮廓进行演化变形，得到最终的管腔轮廓；

c、计算ICUS图像序列中相邻帧之间管腔轮廓的刚性运动参数，包括位移和旋转角：

对于ICUS图像序列中的第k-1帧和第k帧图像(k＝2，3，...，M，M为ICUS图像序列的总帧数)，由分割出的管腔边界曲线γ_k-1和γ_k的重心C_k-1＝(xc_k-1，yc_k-1)和C_k＝(xc_k，yc_k)计算出其位移量(Δx_k，Δy_k)和旋转角Δα_k：

$\left\{\begin{array}{c}{\mathrm{\Δx}}_k=x{c}_k-x{c}_{k-1}\\ {\mathrm{\Δy}}_k=y{c}_k-y{c}_{k-1}\end{array}\right.$

$\mathrm{\Δ}{\mathrm{\α}}_k=\mathrm{arctg}\left(\frac{y{c}_k}{x{c}_k}\right)-\mathrm{arctg}\left(\frac{y{c}_{k-1}}{x{c}_{k-1}}\right)$

$C_{k-1}=(x{c}_{k-1},y{c}_{k-1})=(\frac{1}{N_{k-1}}\underset{i=1}{\overset{N_{k-1}}{\mathrm{\Σ}}}{x}_{k-1,i},\frac{1}{N_{k-1}}\underset{i=1}{\overset{N_{k-1}}{\mathrm{\Σ}}}{y}_{k-1,i})$

$C_k=({\mathrm{xc}}_k,y{c}_k)=(\frac{1}{N_k}\underset{i=1}{\overset{N_k}{\mathrm{\Σ}}}{x}_{k,i},\frac{1}{N_k}\underset{i=1}{\overset{N_k}{\mathrm{\Σ}}}{y}_{k,i})$

其中Δx_k、Δy_k和Δα_k均包含刚性运动分量(Δx_k，d、Δy_k，d和Δα_k，d)和几何分量(Δx_k，g、Δy_k，g和Δα_k，g)，(x_k-1，i，y_k-1，i)、γ_k-1上第i个点在以导管中心为坐标原点的坐标系中的坐标；(x_k，i，y_k，i)、γ_k上第i个点在以导管中心为坐标原点的坐标系中的坐标；N_k-1、γ_k-1上点的总数；N_k、γ_k上点的总数。

d、分离ICUS图像序列中相邻帧之间管腔轮廓刚性运动参数的刚性运动分量和几何分量：

首先对计算出的Δx_k、Δy_k和Δα_k(k＝2，3，...，M)分别进行傅立叶变换，得到其幅度谱ΔX_k、ΔY_k和ΔA_k，然后采用高通滤波器分别对ΔX_k、ΔY_k和ΔA_k进行滤波，得到刚性运动分量的幅度谱ΔX_k，d、ΔY_k，d和ΔA_k，d，对其进行逆傅立叶变换即可得到Δx_k，d、Δy_k，d和Δα_k，d的估计值；

e、补偿原始ICUS图像序列中的刚性运动伪影：

对于ICUS图像序列中的第k帧图像I_k(x，y)(k＝2，3，...，M)，根据估计出的刚性运动分量：位移(Δx_k，d，Δy_k，d)和旋转角Δα_k，d，将I_k(x，y)的所有像素点进行反向平移和旋转，得到消除刚性运动伪影后的图像I_k(x′，y′)，坐标变换关系如下：

$\left[\begin{array}{c}{x}^{\′}\\ y^{\′}\end{array}\right]=\left[\begin{array}{cc}\cos (-\underset{i=2}{\overset{k}{\mathrm{\Σ}}}\mathrm{\Δ}{\mathrm{\α}}_{i,d})& \sin (-\underset{i=2}{\overset{k}{\mathrm{\Σ}}}\mathrm{\Δ}{\mathrm{\α}}_{i,d})\\ -\sin (-\underset{i=2}{\overset{k}{\mathrm{\Σ}}}\mathrm{\Δ}{\mathrm{\α}}_{i,d})& \cos (-\underset{i=2}{\overset{k}{\mathrm{\Σ}}}{\mathrm{\Δ\α}}_{i,d})\end{array}\right]\left[\begin{array}{c}x-\underset{i=2}{\overset{k}{\mathrm{\Σ}}}{\mathrm{\Δx}}_{i,d}\\ y-\underset{i=2}{\overset{k}{\mathrm{\Σ}}}{\mathrm{\Δy}}_{i,d}\end{array}\right],$

上述抑制冠状动脉内超声图像序列中刚性运动伪影的方法，所述对ΔX_k、ΔY_k和ΔA_k进行滤波的高通滤波器的通带截止频率设定为病人的心率值(次/秒)。

上述抑制冠状动脉内超声图像序列中刚性运动伪影的方法，所述病人的心率值按照如下步骤从原始ICUS图像序列中估计：

①对ICUS图像序列进行逐帧比较，按下式计算ICUS图像序列中的第i帧图像I_i和第j帧图像I_j之间灰度特征的差异值d_i，j：

d_i，j＝1-C(i，j)，i＝1，2，...，M-1，j＝i+1，i+2，...，M

$C(i,j)=\left\{\begin{array}{ccc}\mathrm{NCC}(I_i,I_j),& \mathrm{if}& \mathrm{NCC}(I_i,I_j)\&GreaterEqual;0\\ 0,& \mathrm{if}& \mathrm{NCC}(I_i,I_j)<0\end{array}\right.$

$\mathrm{NCC}(I_i,I_j)=\frac{\underset{x=0}{\overset{N_1-1}{\mathrm{\&Sigma;}}}\underset{y=0}{\overset{N_2-1}{\mathrm{\&Sigma;}}}[I_i(x,y)-{\mathrm{\&mu;}}_i][I_j(x,y)-{\mathrm{\&mu;}}_j]}{\sqrt{\left[\underset{x=0}{\overset{N_1-1}{\mathrm{\&Sigma;}}}\underset{y=0}{\overset{N_2-1}{\mathrm{\&Sigma;}}}{[I_i(x,y)-{\mathrm{\&mu;}}_i]}^2\right]\left[\underset{x=0}{\overset{N_1-1}{\mathrm{\&Sigma;}}}\underset{y=0}{\overset{N_2-1}{\mathrm{\&Sigma;}}}{[I_j(x,y)-{\mathrm{\&mu;}}_j]}^2\right]}}$

其中，N₁、N₂分别为ICUS图像像素的总行数和总列数，μ_i为第i帧图像I_i的平均灰度值，μ_j为第j帧图像I_j的平均灰度值；M是ICUS图像序列的总帧数；NCC是两帧图像之间灰度值的归一化互相关。

②计算ICUS图像序列中间隔为i-1帧的两帧图像I_j和I_j+(i-1)差异度的平均值c(i)：

$c\left(i\right)=\frac{1}{M-(i-1)}\underset{j=1}{\overset{M-(i-1)}{\mathrm{\&Sigma;}}}{d}_{j+(i-1),j}$

当i＝1，2，...，M时，则得到一条曲线c(i)～i，并且c(1)＝0。

③计算曲线c(i)～i的第一个谷点(极小值)所对应的i值，记为R，则病人心率的近似值为：F/R次/秒(Hz)

其中，F为ICUS图像的采集速率，单位是帧/秒。

本发明既不需要专门的心电门控图像采集装置，也不需要同步记录心电信号，而是运用图像分析技术，通过对连续回撤超声导管采集到的、覆盖多个心动周期的ICUS图像序列数据进行处理，实现对原始ICUS图像序列中刚性运动伪影的补偿。与现有的心电门控方法相比较，本发明具有应用成本低和简单易行的优点。同时，本发明维持了连续回撤导管采集的ICUS图像数据集合的完整性，从而为冠状动脉血管病变的诊断和治疗提供更为丰富的信息。

附图说明

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

图1是ICUS图像纵向视图中的运动伪影；

图2是ICUS图像横向视图中的运动伪影；

图3是ICUS图像刚性运动参数示意图；

图4是ICUS图像序列三维自动分割方法流程图；

图5是一个ICUS图像序列中的管腔边界重心在x方向上的位移(Δx)变化曲线。

图中各符号为：γ₁、γ₂、时刻t₁和t₂的管腔边界曲线；C₁、C₂、γ₁和γ₂的重心；Δx、Δy、C₁和C₂之间分别在x和y方向上的位移；Δα、γ₁和γ₂之间的旋转角；(x₀，y₀)、导管中心的坐标。

文中所用符号：γ₁(x₁，y₁)、时刻t₁的管腔边界曲线；γ₂(x₂，y₂)、时刻t₂的管腔边界曲线；C₁、γ₁的重心；(xc₁，yc₁)、C₁的坐标；C₂、γ₂的重心；(xc₂，yc₂)、C₂的坐标；(Δx，Δy)、C₁和C₂之间的位移；Δx_d、Δx的刚性运动分量；Δx_g、Δx的几何分量；Δy_d、Δy的刚性运动分量；Δy_g、Δy的几何分量；Δα、γ₁和γ₂之间的旋转角；Δα_d、Δα的刚性运动分量；Δα_g、Δα的几何分量；M、ICUS图像序列的总帧数；γ_k、从第k帧图像中分割出的管腔边界曲线；γ_k+1、从第k+1帧图像中分割出的管腔边界曲线；C_k、γ_k的重心；C_k+1、γ_k+1的重心；(Δx_k，Δy_k)、γ_k和γ_k+1之间管腔边界重心之间的位移；Δα_k、γ_k和γ_k+1之间的旋转角；ΔX_k、Δx_k的幅度谱；ΔY_k、Δy_k的幅度谱；ΔA_k、Δα_k的幅度谱；ΔX_k，d、Δx_k的刚性运动分量的幅度谱；ΔY_k，d、Δy_k的刚性运动分量的幅度谱；ΔA_k，d、Δα_k的刚性运动分量的幅度谱；d_i，j、第i帧图像I_i和第j帧图像I_j之间灰度特征的差异值；N₁、ICUS图像像素的总行数；N₂、ICUS图像像素的总列数；μ_i、第i帧图像I_i的平均灰度值，μ_j、第j帧图像I_j的平均灰度值；NCC、两帧图像之间灰度值的归一化互相关；c(i)、第j帧图像I_j和第j+(i-1)帧图像I_j+(i-1)之间差异度的平均值；F、ICUS图像的采集速率(单位：帧/秒)；R、曲线c(i)～i的第一个谷点所对应的i值；(x_k-1，i，y_k-1，i)、γ_k-1上第i个点在以导管中心为坐标原点的坐标系中的坐标；(x_k，i，y_k，i)、γ_k上第i个点在以导管中心为坐标原点的坐标系中的坐标；N_k-1、γ_k-1上点的总数；N_k、γ_k上点的总数；I_k(x，y)、ICUS图像序列中的第k帧图像；I_k(x′，y′)、消除刚性运动伪影后的第k帧ICUS图像。

具体实施方式

(1)依照已有技术建立ICUS图像序列中血管壁的刚性运动模型：

(1.1)ICUS图像序列中运动伪影的表现形式：

ICUS图像序列中的运动伪影表现为：①序列图像中血管横截面的平移和旋转，即帧间的错位，为刚性运动；②管腔的径向变形，即沿血管半径方向的膨胀/收缩，为非刚性运动。

本发明仅对刚性运动进行定量估计和补偿。由于单纯的变形不会改变弹性体的重心和方向，即径向变形不会影响刚性运动参数的估计结果，因此本发明建立的运动模型中不包含径向变形。

(1.2)产生刚性运动的原因：

由于存在运动伪影，使得管腔长轴方向的视图中呈现出锯齿状的管壁边缘(参看图1)，在横向视图中表现为沿血管切线方向的灰度偏移(参看图2)。在连续回撤超声导管采集到的ICUS图像序列中，相邻帧之间管腔横截面空间方向的改变和其重心的位移主要由运动和几何两方面的因素造成：

①运动：即外部因素，指由周期性心脏运动和搏动的血流造成的导管相对于管腔的运动和血管本身的运动所引起；

②几何：即内部因素，指由血管本身不规则的几何形状所引起。

本发明用从各帧ICUS图像中分割出的管腔边界曲线重心之间的位移(Δx，Δy)和边界曲线之间的旋转角Δα来分别表示相邻帧之间管腔横截面的相对平移和旋转，它们分别由两部分组成：

$\left\{\begin{array}{c}\mathrm{\&Delta;x}={\mathrm{\&Delta;x}}_d+{\mathrm{\&Delta;x}}_g\\ \mathrm{\&Delta;y}={\mathrm{\&Delta;y}}_d+{\mathrm{\&Delta;y}}_g\\ \mathrm{\&Delta;\&alpha;}=\mathrm{\&Delta;}{\mathrm{\&alpha;}}_d+\mathrm{\&Delta;}{\mathrm{\&alpha;}}_g\end{array}\right.---\left(1\right)$

式中脚标d表示刚性运动分量，g表示几何分量。本发明通过将二者分开，仅补偿其中的刚性运动分量，完成对刚性运动伪影的抑制。

(1.3)建立血管壁的刚性运动模型：

参看图3，设时刻t₁的ICUS图像中管腔边界曲线为γ₁(x₁，y₁)，其重心为C₁＝(xc₁，yc₁)；时刻t₂的ICUS图像中管腔边界曲线为γ₂(x₂，y₂)，其重心为C₂＝(xc₂，yc₂)。C₁和C₂之间的位移为(Δx，Δy)：

(xc₂，yc₂)＝(xc₁+Δx，yc₁+Δy)         (2)

其中Δx和Δy如式(1)所示分别包含刚性运动和几何两个分量。设γ₁和γ₂之间的旋转角为Δα，旋转中心为导管中心(x₀，y₀)：

$\mathrm{\&Delta;}\mathrm{\&alpha;}=\mathrm{arctg}\left(\frac{y{c}_2}{x{c}_2}\right)-\mathrm{arctg}\left(\frac{y{c}_1}{x{c}_1}\right)---\left(3\right)$

那么γ₁(x₁，y₁)与γ₂(x₂，y₂)之间的几何变换关系为：

$\left[\begin{array}{c}{x}_2\\ y_2\end{array}\right]=\left[\begin{array}{cc}\cos \mathrm{\&Delta;\&alpha;}& \sin \mathrm{\&Delta;\&alpha;}\\ -\sin \mathrm{\&Delta;\&alpha;}& \cos \mathrm{\&Delta;\&alpha;}\end{array}\right]\left[\begin{array}{c}{x}_1+\mathrm{\&Delta;x}\\ y_1+\mathrm{\&Delta;y}\end{array}\right]---\left(4\right)$

对于ICUS图像序列中的各相邻帧，分别计算出位移(Δx，Δy)和旋转角Δα，并将其刚性运动分量(Δx_d、Δy_d和Δα_d)和几何分量(Δx_g、Δy_g和Δα_g)分离开，即可完成对刚性运动的定量估计。将原始图像按照估计出的旋转角Δα_d和位移(Δx_d，Δy_d)分别进行反向旋转和平移即可实现对刚性运动伪影的补偿。

(2)提取各帧ICUS图像中的管腔轮廓：

参看图4，本发明采用三维自动分割方法，对ICUS图像序列进行分割，提取出各帧图像中的管腔轮廓，包括如下步骤：

步骤1获取ICUS纵向视图：按照采集顺序以空间体素的形式保存图像，获得ICUS图像序列的纵向视图。

步骤2利用管腔形状和灰度特征的先验知识(在纵向视图中，内膜位于无回声的内腔两侧)，自动提取出ICUS纵向视图中管腔的左右边界。

步骤3获取横截面ICUS图像中的初始管腔轮廓：将纵向视图中的边界线映射到各帧横截面ICUS图像上，得到每一帧中的初始管腔轮廓。

步骤4提取横截面ICUS图像中的管腔轮廓：利用快速推进法(M.R.Cardinal，J.Meunier，G.Soulez，R.L.Maurice，E.Therasse，G.Cloutier.Intravascular ultrasound imagesegmentation：A three-dimensional fast-marching method based on gray leveldistributions.IEEE Transactions on Medical Imaging.2006，vol.25，no.5，pp.590-601)对初始管腔轮廓进行演化变形，得到最终的管腔轮廓。

该方法可实现对ICUS图像序列的自动分割，无须操作者的手动参与，避免操作者技术和主观因素的影响。同时，可实现对ICUS图像序列的并行处理，大大提高计算效率，缩短处理时间。

从第k帧(k＝1，2，...，M，M为ICUS图像序列的总帧数)ICUS图像中分割出管腔边界曲线γ_k之后，计算γ_k的几何中心C_k作为其重心的近似：

$C_k=(x{c}_k,y{c}_k)=(\frac{1}{N_k}\underset{i=1}{\overset{N_k}{\mathrm{\&Sigma;}}}{x}_{k,i},\frac{1}{N_k}\underset{i=1}{\overset{N_k}{\mathrm{\&Sigma;}}}{y}_{k,i})---\left(5\right)$

其中(x_k，i，y_k，i)(i＝1，2，...，N_k)是γ_k上第i个点在以导管中心为坐标原点的坐标系中的坐标，N_k是γ_k上点的总数。

(3)计算ICUS图像序列中相邻帧之间管腔轮廓的刚性运动参数：

对于ICUS图像序列中的第k-1帧和第k帧图像(k＝2，3，...，M)，由分割出的管腔边界曲线γ_k-1和γ_k的重心C_k-1＝(xc_k-1，yc_k-1)和C_k＝(xc_k，yc_k)计算出其位移量和旋转角：

$\left\{\begin{array}{c}{\mathrm{\&Delta;x}}_k=x{c}_k-x{c}_{k-1}\\ {\mathrm{\&Delta;y}}_k=y{c}_k-y{c}_{k-1}\end{array}\right.---\left(6\right)$

$\mathrm{\&Delta;}{\mathrm{\&alpha;}}_k=\mathrm{arctg}\left(\frac{y{c}_k}{x{c}_k}\right)-\mathrm{arctg}\left(\frac{y{c}_{k-1}}{x{c}_{k-1}}\right)---\left(7\right)$

如式(1)所示，Δx_k、Δy_k和Δα_k均包含刚性运动(Δx_k，d、Δy_k，d和Δα_k，d)和几何(Δx_k，g、Δy_k，g和Δα_k，g)两个分量，需将二者分离开。

(4)分离ICUS图像序列中相邻帧之间管腔轮廓刚性运动参数的刚性运动分量和几何分量：

本发明采用谱分析的方法完成刚性运动分量和几何分量的分离，即对计算出的Δx_k、Δy_k和Δα_k(k＝2，3，...，M)分别进行傅立叶变换，得到其幅度谱ΔX_k、ΔY_k和ΔA_k。由于由血管本身几何形状的变化所引起的相邻帧之间管腔横截面空间方向和重心位置的变化速度，应远小于由周期性心脏运动所致的变化速度，因此在ΔX_k、ΔY_k和ΔA_k中，高频分量对应于刚性运动分量Δx_k，d、Δy_k，d和Δα_k，d，而低频分量对应于几何分量Δx_k，g、Δy_k，g和Δα_k，g，如附图5所示。

本发明采用高通滤波器分别对ΔX_k、ΔY_k和ΔA_k进行滤波，则滤波器的输出即是刚性运动分量的幅度谱ΔX_k，d、ΔY_k，d和ΔA_k，d，对其进行逆傅立叶变换即可得到Δx_k，d、Δy_k，d和Δα_k，d的估计值。

由于刚性运动分量Δx_k，d、Δy_k，d和Δα_k，d主要由心脏的周期性运动引起，因此本发明将高通滤波器的通带截止频率设定为病人的心率值(单位：次/秒，即Hz)。

本发明按照以下步骤从原始ICUS图像序列中估计病人的心率值(单位：次/秒，即Hz)：

步骤1对ICUS图像序列进行逐帧比较，计算各帧之间灰度特征的差异值d_i，j：

d_i，j＝1-C(i，j)，i＝1，2，...，M-1，j＝i+1，i+2，...，M    (8)

$C(i,j)=\left\{\begin{array}{ccc}\mathrm{NCC}(I_i,I_j),& \mathrm{if}& \mathrm{NCC}(I_i,I_j)\&GreaterEqual;0\\ 0,& \mathrm{if}& \mathrm{NCC}(I_i,I_j)<0\end{array}\right.---\left(9\right)$

$\mathrm{NCC}(I_i,I_j)=\frac{\underset{x=0}{\overset{N_1-1}{\mathrm{\&Sigma;}}}\underset{y=0}{\overset{N_2-1}{\mathrm{\&Sigma;}}}[I_i(x,y)-{\mathrm{\&mu;}}_i][I_j(x,y)-{\mathrm{\&mu;}}_j]}{\sqrt{\left[\underset{x=0}{\overset{N_1-1}{\mathrm{\&Sigma;}}}\underset{y=0}{\overset{N_2-1}{\mathrm{\&Sigma;}}}{[I_i(x,y)-{\mathrm{\&mu;}}_i]}^2\right]\left[\underset{x=0}{\overset{N_1-1}{\mathrm{\&Sigma;}}}\underset{y=0}{\overset{N_2-1}{\mathrm{\&Sigma;}}}{[I_j(x,y)-{\mathrm{\&mu;}}_j]}^2\right]}}---\left(10\right)$

其中I_i和I_j分别是ICUS图像序列中的第i帧和第j帧图像，其大小均为N₁×N₂像素，其平均灰度值分别为μ_i和μ_j；M是ICUS图像序列的总帧数；NCC是两帧图像之间灰度值的归一化互相关。

步骤2计算ICUS图像序列中间隔为i-1帧的两帧图像I_j和I_j+(l-1)差异度的平均值c(i)：

$c\left(i\right)=\frac{1}{M-(i-1)}\underset{j=1}{\overset{M-(i-1)}{\mathrm{\&Sigma;}}}{d}_{j+(i-1),j}---\left(11\right)$

当i＝1，2，...，M时，则得到一条曲线c(i)～i，并且c(1)＝0。

步骤3计算曲线c(i)～i的第一个谷点(极小值)所对应的i值，作为以帧为单位的心动周期长度的估计值，记为R。根据已知的ICUS图像采集速率F(单位：帧/秒)(典型值为F＝30帧/秒)，可得到病人心率的近似值：F/R次/秒(Hz)。

(5)补偿原始ICUS图像序列中的刚性运动伪影：

估计出ICUS图像序列中的第k帧图像I_k(x，y)(k＝2，3，...，M)的刚性运动分量：位移(Δx_k，d，Δy_k，d)和旋转角Δα_k，d之后，即可对其各像素点的刚性运动进行补偿。方法是将I_k(x，y)中、基于以导管中心为坐标原点的坐标系的各像素点的坐标先反向平移

再反向旋转即得到消除刚性运动伪影后的图像I_k(x′，y′)：

$\left[\begin{array}{c}{x}^{\&prime;}\\ y^{\&prime;}\end{array}\right]=\left[\begin{array}{cc}\cos (-\underset{i=2}{\overset{k}{\mathrm{\&Sigma;}}}\mathrm{\&Delta;}{\mathrm{\&alpha;}}_{i,d})& \sin (-\underset{i=2}{\overset{k}{\mathrm{\&Sigma;}}}\mathrm{\&Delta;}{\mathrm{\&alpha;}}_{i,d})\\ -\sin (-\underset{i=2}{\overset{k}{\mathrm{\&Sigma;}}}\mathrm{\&Delta;}{\mathrm{\&alpha;}}_{i,d})& \cos (-\underset{i=2}{\overset{k}{\mathrm{\&Sigma;}}}{\mathrm{\&Delta;\&alpha;}}_{i,d})\end{array}\right]\left[\begin{array}{c}x-\underset{i=2}{\overset{k}{\mathrm{\&Sigma;}}}{\mathrm{\&Delta;x}}_{i,d}\\ y-\underset{i=2}{\overset{k}{\mathrm{\&Sigma;}}}{\mathrm{\&Delta;y}}_{i,d}\end{array}\right]---\left(12\right).$

Claims (3)

1、一种抑制冠状动脉内超声图像序列中刚性运动伪影的方法，其特征是，它是在建立血管壁刚性运动模型的基础上，首先采用三维自动分割技术从ICUS图像序列的各帧中分割出管腔轮廓，并计算相邻帧之间管腔轮廓的刚性运动参数，所述刚性运动包括平移和旋转，然后从中分离出刚性运动分量，最后通过对原始ICUS图像中的各像素点进行反向平移和旋转，使其刚性运动伪影得到补偿，具体步骤如下：

a、建立ICUS图像序列中血管壁的刚性运动模型：

ICUS图像序列中，相邻帧之间管腔轮廓的旋转和平移由刚性运动分量和几何分量构成：

$\left\{\begin{array}{c}\mathrm{\&Delta;x}={\mathrm{\&Delta;x}}_d+{\mathrm{\&Delta;x}}_g\\ \mathrm{\&Delta;y}={\mathrm{\&Delta;y}}_d+{\mathrm{\&Delta;y}}_g\\ \mathrm{\&Delta;\&alpha;}={\mathrm{\&Delta;\&alpha;}}_d+{\mathrm{\&Delta;\&alpha;}}_g\end{array}\right.$

式中(Δx，Δy)是相邻帧之间管腔边界重心之间的位移，Δα是相邻帧之间管腔边界曲线之间的旋转角，脚标d表示刚性运动分量，g表示几何分量；

b、采用三维自动分割技术从各帧ICUS图像中分割出管腔轮廓：

首先根据ICUS图像中管腔形状和灰度特征的先验知识，从ICUS纵向视图中分割出管腔的左右边界，并将其映射到各帧横截面ICUS图像中，得到每一帧中的初始管腔轮廓，然后采用快速推进法对初始管腔轮廓进行演化变形，得到最终的管腔轮廓；

c、计算ICUS图像序列中相邻帧之间管腔轮廓的刚性运动参数，包括位移和旋转角：

对于ICUS图像序列中的第k-1帧和第k帧图像(k＝2，3，...，M，M为ICUS图像序列的总帧数)，由分割出的管腔边界曲线γ_k-1和γ_k的重心C_k-1＝(xc_k-1，yc_k-1)和C_k＝(xc_k，yc_k)计算出其位移量(Δx_k，Δy_k)和旋转角Δα_k：

$\left\{\begin{array}{c}{\mathrm{\&Delta;x}}_k={\mathrm{xc}}_k-{\mathrm{xc}}_{k-1}\\ {\mathrm{\&Delta;y}}_k={\mathrm{yc}}_k-{\mathrm{yc}}_{k-1}\end{array}\right.$

${\mathrm{\&Delta;\&alpha;}}_k=\mathrm{arctg}\left(\frac{{\mathrm{yc}}_k}{{\mathrm{xc}}_k}\right)-\mathrm{arctg}\left(\frac{{\mathrm{yc}}_{k-1}}{{\mathrm{xc}}_{k-1}}\right)$

$C_{k-1}=({\mathrm{xc}}_{k-1},{\mathrm{yc}}_{k-1})=(\frac{1}{N_{k-1}}\underset{i=1}{\overset{N_{k-1}}{\mathrm{\&Sigma;}}}{x}_{k-1,i},\frac{1}{N_{k-1}}\underset{i=1}{\overset{N_{k-1}}{\mathrm{\&Sigma;}}}{y}_{k-1,i})$

$C_k=({\mathrm{xc}}_k,{\mathrm{yc}}_k)=(\frac{1}{N_k}\underset{i=1}{\overset{N_k}{\mathrm{\&Sigma;}}}{x}_{k,i},\frac{1}{N_k}\underset{i=1}{\overset{N_k}{\mathrm{\&Sigma;}}}{y}_{k,i})$

其中Δx_k、Δy_k和Δα_k均包含刚性运动分量(Δx_k，d、Δy_k，d和Δα_k，d)和几何分量(Δx_k，g、Δy_k，g和Δα_k，g)，(x_k-1，i，y_k-1，i)、γ_k-1上第i个点在以导管中心为坐标原点的坐标系中的坐标；(x_k，i，y_k，i)、γ_k上第i个点在以导管中心为坐标原点的坐标系中的坐标；N_k-1、γ_k-1上点的总数；N_k、γ_k上点的总数。

d、分离ICUS图像序列中相邻帧之间管腔轮廓刚性运动参数的刚性运动分量和几何分量：

首先对计算出的Δx_k、Δy_k和Δα_k(k＝2，3，...，M)分别进行傅立叶变换，得到其幅度谱ΔX_k、ΔY_k和ΔA_k，然后采用高通滤波器分别对ΔX_k、ΔY_k和ΔA_k进行滤波，得到刚性运动分量的幅度谱ΔX_k，d、ΔY_k，d和ΔA_k，d，对其进行逆傅立叶变换即可得到Δx_k，d、Δy_k，d和Δα_k，d的估计值；

e、补偿原始ICUS图像序列中的刚性运动伪影：

对ICUS图像序列中的第k帧图像I_k(x，y)(k＝2，3，...，M)，根据估计出的刚性运动分量：位移(Δx_k，d，Δy_k，d)和旋转角Δα_k，d，将I_k(x，y)的所有像素点进行反向平移和旋转，得到消除刚性运动伪影后的图像I_k(x′，y′)，坐标变换关系如下：

$\begin{array}{}\end{array}\left[\begin{array}{c}{x}^{\&prime;}\\ y^{\&prime;}\end{array}\right]=\left[\begin{array}{cc}\cos (-\underset{i=2}{\overset{k}{\mathrm{\&Sigma;}}}{\mathrm{\&Delta;\&alpha;}}_{i,d})& \sin (-\underset{i=2}{\overset{k}{\mathrm{\&Sigma;}}}{\mathrm{\&Delta;\&alpha;}}_{i,d})\\ -\sin (-\underset{i=2}{\overset{k}{\mathrm{\&Sigma;}}}{\mathrm{\&Delta;\&alpha;}}_{i,d})& \cos (-\underset{i=2}{\overset{k}{\mathrm{\&Sigma;}}}{\mathrm{\&Delta;\&alpha;}}_{i,d})\end{array}\right]\left[\begin{array}{c}x-\underset{i=2}{\overset{k}{\mathrm{\&Sigma;}}}{\mathrm{\&Delta;x}}_{i,d}\\ y-\underset{i=2}{\overset{k}{\mathrm{\&Sigma;}}}{\mathrm{\&Delta;y}}_{i,d}\end{array}\right].$

2、根据权利要求1所述抑制冠状动脉内超声图像序列中刚性运动伪影的方法，其特征是，所述对ΔX_k、ΔY_k和ΔA_k进行滤波的高通滤波器的通带截止频率设定为病人的心率值次/秒。

3、根据权利要求2所述抑制冠状动脉内超声图像序列中刚性运动伪影的方法，其特征是，所述病人的心率值按照如下步骤从原始ICUS图像序列中估计：

①对ICUS图像序列进行逐帧比较，按下式计算ICUS图像序列中的第i帧图像I_i和第j帧图像I_j之间灰度特征的差异值d_i，j：

d_i，j＝1-C(i，j)，i＝1，2，...，M-1，j＝i+1，i+2，...，M

$C(i,j)=\left\{\begin{array}{ccc}\mathrm{NCC}(I_i,I_j),& \mathrm{if}& \mathrm{NCC}(I_i,I_j)\&GreaterEqual;0\\ 0,& \mathrm{if}& \mathrm{NCC}(I_i,I_j)<0\end{array}\right.$

$\mathrm{NCC}(I_i,I_j)=\frac{\underset{x=0}{\overset{N_1-1}{\mathrm{\&Sigma;}}}\underset{y=0}{\overset{N_2-1}{\mathrm{\&Sigma;}}}[I_i(x,y)-{\mathrm{\&mu;}}_i][I_j(x,y)-{\mathrm{\&mu;}}_j]}{\sqrt{\left[\underset{x=0}{\overset{N_1-1}{\mathrm{\&Sigma;}}}\underset{y=0}{\overset{N_2-1}{\mathrm{\&Sigma;}}}{[I_i(x,y)-{\mathrm{\&mu;}}_i]}^2\right]\left[\underset{x=0}{\overset{N_1-1}{\mathrm{\&Sigma;}}}\underset{y=0}{\overset{N_2-1}{\mathrm{\&Sigma;}}}{[I_j(x,y)-{\mathrm{\&mu;}}_j]}^2\right]}}$

其中，N₁、N₂分别为ICUS图像像素的总行数和总列数，μ_i为第i帧图像I_i的平均灰度值，μ_j为第j帧图像I_j的平均灰度值；M是ICUS图像序列的总帧数；NCC是两帧图像之间灰度值的归一化互相关。

②计算ICUS图像序列中间隔为i-1帧的两帧图像I_j和I_j+(i-1)差异度的平均值c(i)：

$c\left(i\right)=\frac{1}{M-(i-1)}\underset{j=1}{\overset{M-(i-1)}{\mathrm{\&Sigma;}}}{d}_{j+(i-1),j}$

当i＝1，2，...，M时，则得到一条曲线c(i)～i，并且c(1)＝0。

③计算曲线c(i)～i的第一个谷点所对应的i值，记为R，则病人心率的近似值为：F/R次/秒，

其中，F为ICUS图像的采集速率，单位是帧/秒。

Images