首页> 中国专利> 基于单演信号特征距离和互相关变换光流算法的电影核磁共振图像序列运动估计方法

基于单演信号特征距离和互相关变换光流算法的电影核磁共振图像序列运动估计方法

摘要

基于单演信号特征距离和互相关变换光流算法的电影核磁共振图像序列运动估计方法,本发明涉及电影核磁共振图像序列运动估计方法。本发明是要解决电影核磁共振成像技术难以找到密集的对应特征点,对于电影核磁共振图像的运动估计比对于加标记的核磁共振图像的运动估计更具有难度的问题,该方法是通过一、将三种特征联合构造三维单演信号特征矩阵;二、利用基于零均值归一化互相关系数的公式计算匹配特征E

著录项

  • 公开/公告号CN105631897A

    专利类型发明专利

  • 公开/公告日2016-06-01

    原文格式PDF

  • 申请/专利权人 哈尔滨工业大学;

    申请/专利号CN201510974766.4

  • 申请日2015-12-22

  • 分类号G06T7/20;

  • 代理机构哈尔滨市松花江专利商标事务所;

  • 代理人杨立超

  • 地址 150001 黑龙江省哈尔滨市南岗区西大直街92号

  • 入库时间 2023-12-18 15:29:29

法律信息

  • 法律状态公告日

    法律状态信息

    法律状态

  • 2018-07-03

    授权

    授权

  • 2016-06-29

    实质审查的生效 IPC(主分类):G06T7/20 申请日:20151222

    实质审查的生效

  • 2016-06-01

    公开

    公开

说明书

技术领域

本发明涉及电影核磁共振图像序列运动估计方法,特别涉及基于单演信号特征距离 和互相关变换光流算法的电影核磁共振图像序列运动估计方法。

背景技术

心功能不全是一种临床常见综合症,病情严重可导致心力衰竭症状的出现,早期评估 预测患者的心肌运动特点对心脏疾病的诊断和治疗有重大意义。心肌运动的运动估计与跟 踪对于临床心功能量化评估、图像引导的手术导航具有重要的意义,对于医学工程领域中 关于非刚体建模和运动仿真也具有重要的指导作用。心脏疾病特别是心血管疾病也会导致 异常的心脏运动。心肌梗死患者的左心室就会出现运动功能不全的现象。如果心脏疾病阻 碍人体局部心肌血液灌注,则会出现运动机能减退、非同步收缩运动、反常运动甚至不能 运动等四种异常的心肌运动形式。心脏是一个始终处于运动状态中的器官,对左心室进行 运动分析目的在于通过提取心脏序列图像中左心室的特征信息,跟踪心肌在一个心动周期 中的运动轨迹。

核磁共振成像(MagneticResonanceImaging,MRI)技术已经成为心脏疾病临床诊断的 重要辅助手段。该技术能够无侵入地检测人体的组织和器宫,而且其成像机理使得该方法 对生物体内如心脏这样的软组织特别有效。在心脏疾病诊断和心功能评估等领域中,基于 核磁共振成像序列的心脏运动和形变估计是一个不可忽视的研究方向。由于电影核磁共振 图像的灰度非常相近,难以找到密集的对应特征点,因而对于电影核磁共振图像的运动估 计比对于加标记的核磁共振图像的运动估计更具有难度,目前国内外研究较少。而针对电 影核磁共振图像的灰度相近,特征稀少的挑战,引入更多的特征用于运动估计就显得很重 要了。

发明内容

本发明的目的是为了解决核磁共振成像技术难以找到密集的对应特征点,对于电影核 磁共振图像的运动估计比对于加标记的核磁共振图像的运动估计更具有难度的问题,而提 出的基于单演信号特征距离和互相关变换光流算法的电影核磁共振图像序列运动估计方 法。

上述的发明目的是通过以下技术方案实现的:

步骤一、利用空间正交滤波器在电影核磁共振图像中提取单演信号的局部相位、局部 方位和局部振幅三种特征,将三种特征联合构造三维单演信号特征矩阵;其中,设置电影 核磁共振图像序列中任意两幅相邻的图像为I1和I2

步骤二、用三维单演信号特征矩阵分别代替彩色RGB图像的3种颜色通道信息作为 光流算法的输入图像;利用基于零均值归一化互相关系数的公式计算光流算法的输入图像 局部特征的匹配特征项Ed(V)作为光流算法的匹配特征项;其中,光流算法由数据的匹配 公式和平滑公式两部分构成;

步骤三、利用双边滤波器根据基于零均值归一化互相关系数的光流算法的平滑公式对 两幅相邻的图像I1和I2进行滤波,限定电影核磁共振图像的点的位移V的最小化匹配误 差即电影核磁共振图像局部特征的平滑特征项Es(V);

步骤四、将步骤二计算的匹配特征项Ed(V)和步骤三计算的电影核磁共振图像局部特 征的平滑特征项Es(V)表示为联合的能量公式:

E(V)=ε·Ed(V)+Es(V)(1)

其中,ε为平衡参数,在光流算法中,有两个变量需要估计,分别是电影核磁共振图 像像素点的水平方向位移u和电影核磁共振图像像素点的垂直方向的位移v;V=(u,v); Ed(V)为电影核磁共振图像的最小化匹配误差即图像局部特征的匹配特征项;Es(V)为电 影核磁共振图像局部特征的平滑特征项。

发明效果

本发明的目的是为了提高电影核磁共振序列图像运动估计的精度。而提出基于单演信 号距离和互相关变换光流法的电影核磁共振图像序列运动估计方法。

在本发明中,我们利用ASSESS工具箱生成模拟心肌运动图像序列和位移场来确定 新型光流算法的精度。通过ASSESS软件生成模拟的电影核磁共振图像序列 D30R20P0F20,D30R20P3F20,D30R20P0F34,D30R20P3F34,这些序列的运动结果是 提前设定好的,通过这些值与估计出的位移值进行比对,就可以算出运动估计算法精度。 D30表示收缩或舒张30%,R20表示旋转角度20度,P0表示健康的序列,P3表示有疾 病的序列,F20表示序列中有20帧,F34表示序列中有34帧。如图2(a)中的电影核磁共 振图像和图2(b)中的加标记的核磁共振图像比起来特征稀少,所以将图2(a)中的图像分解 成单演信号的振幅信息如图2(c),相位信息如图2(d),方位信息如图2(e),然后合成伪彩色图 像如图2(f)。图2(f)比图2(a)具有更多的特征信息。

本发明采用平均角度误差AAE(AverageAngleError)衡量运动估计效果

角度误差常用于直观的给出光流估计结果的好或者差。其计算方法为,对图像每一点 计算其估计的速度与光流真值速度的夹角,夹角越小表示角度结果越准确。平均角度误差 AAE(AverageAngleError),平均角度误差是一种重要的量化评估光流估计结果的方法。设 图像任意点p的角度误差为EAE(p))是两个流向量v0(p)=(u0,v0)和v1(p)=(u1,v1)在p点在 2D空间的角度(u0,v0)和(u1,v1).我们可以得到:

EAE=arccos(v0,v1)

角度误差适合用来评价大速度和小速度的位移。如果位移的结果给定了,那么EAE(p)是 估计值和真值之间的角度误差。AAE是EAE(p)平均值。

并且本发明采用平均终止点误差AEP(averageend-point)衡量运动估计效果 终止点误差(EndPointError)用来衡量两个光流终止点之间的距离,公式是

(u0-u1)2+(u1-v1)2

如果真实结果给定,那么v0(p)=(u0,v0)表示真值,v1(p)=(u1,v1)表示估计值。AEE (AverageEnd-pointError)是终止点误差的平均值。标准差(standarddeviation,STD)用来 评价运动估计结果的稳定性。

本发明主要采用光流算法,所以和目前效果较好的Sun光流算法,还有经典的LK光 流算法进行比较。比较效果如表1所示。

表1平均角度误差和平均终止点误差(均值±标准差)(角度误差:度,终止点误差:像素)

附图说明

图1为具体实施三提出的单演信号分解示意图;其中,为单演相位,θ为单演方位, A为单演振幅,r表示单演信号相位向量,p为p(x),q1为q1(x),q2为q2(x),q为q(x));

图2(a)为具体实施方式一提出的电影核磁共振图像;

图2(b)为具体实施方式一提出的加标记的核磁共振图像;

图2(c)为具体实施方式一提出的局部振幅图像;

图2(d)为具体实施方式一提出的局部相位图像;

图2(e)为具体实施方式一提出的局部方位图像;

图2(f)为具体实施方式一提出的合成伪彩色图像;

图3为具体实施方式三提出的基于单演信号特征的3维矩阵构建示意图。

具体实施方式

具体实施方式一:本实施方式的基于单演信号特征距离和互相关变换光流算法的电影 核磁共振图像序列运动估计方法,具体是按照以下步骤制备的:

步骤一、利用空间正交滤波器(sphericalquadraturefilters,SQFs)在电影核磁共振图像中 提取单演信号的局部相位、局部方位和局部振幅三种特征,将三种特征联合构造三维单演 信号特征矩阵;并且三种特征是互相独立的;其中,设置电影核磁共振图像序列中任意两 幅相邻的图像为I1和I2;电影核磁共振图像具体为医学诊断中最常用、最普通的磁共振检 查方式。

二维解析信号是一维解析信号的一个扩展,原始的二维图像信息通过具有旋转不变性 的广义Hilbert变换,以一种非线性的方式被映射到了虚平面;与一维解析信号类似,空 间正交滤波器涵盖了原始图像在实平面和虚平面上的不同信息,因而可以据此提取出相应 的单演幅度,单演相位和单演方向特征信息;与其他提取图像相位信息的方法如基于 Gabor滤波器的相位计算方法相比,该模型的相位计算不需要对方向进行采样,也不需要 根据计算的方向值来调整Hilbert变换;

步骤二、用三维单演信号特征矩阵分别代替彩色RGB图像(红绿蓝图像)的3种颜 色通道信息作为光流算法的输入图像;利用基于零均值归一化互相关系数的公式计算光流 算法的输入图像局部特征的匹配特征项Ed(V)作为光流算法的匹配特征项;其中,光流算 法已经广泛的应用于医学图像处理中;由数据的匹配公式和平滑公式两部分构成;

步骤三、利用双边滤波器根据基于零均值归一化互相关系数的光流算法的平滑公式对 两幅相邻的图像I1和I2进行滤波,限定电影核磁共振图像的点的位移为V的最小化匹配 误差即电影核磁共振图像局部特征的平滑特征项Es(V);

步骤四、将步骤二计算的匹配特征项Ed(V)和步骤三计算的电影核磁共振图像局部特 征的平滑特征项Es(V)表示为联合的能量公式:

E(V)=ε·Ed(V)+Es(V)(1)

其中,ε为平衡参数,用来调节Ed(V)的占比;在光流算法中,有两个变量需要估计, 分别是电影核磁共振图像像素点的水平方向位移u和电影核磁共振图像像素点的垂直方向 的位移v;V=(u,v);Ed(V)为电影核磁共振图像的最小化匹配误差即图像局部特征的匹配 特征项;Es(V)为电影核磁共振图像局部特征的平滑特征项;

在光流算法中,匹配条件其实变化不大,一般是像素和像素的匹配或者像素梯度的匹 配,为了获得更好的运动估计精度和鲁棒性,更多的匹配条件应该被引入到光流算法;

数字图像相关性(Digitalimagecorrelation,DIC)已经被广泛的应用于图像表面变形的 度量,求两个特定窗口大小的信号之间的平方和之差(sumofsquareddifferences,SSD)的最 小值是一个基本的方法,SSD特征的匹配增加了光流算法的抗噪性;零均值归一化互相 关系数(Zero-meanNormalizedCrossCorrelation,ZNCC)实际上就是给定信号相关变换的 SSD距离;ZNCC特性在纹理稀疏区域表现出了更好的可靠性;采用零均值归一化互相关 系数ZNCC特征作为两个图像区域之间相似性的度量公式。

本实施方式效果:

本实施方式的目的是为了提高电影核磁共振序列图像运动估计的精度。而提出基于单 演信号距离和互相关变换光流法的电影核磁共振图像序列运动估计方法。

在本实施方式中,我们利用ASSESS工具箱生成模拟心肌运动图像序列和位移场来 确定新型光流算法的精度。通过ASSESS软件生成模拟的电影核磁共振图像序列 D30R20P0F20,D30R20P3F20,D30R20P0F34,D30R20P3F34,这些序列的运动结果是 提前设定好的,通过这些值与估计出的位移值进行比对,就可以算出运动估计算法精度。 D30表示收缩或舒张30%,R20表示旋转角度20度,P0表示健康的序列,P3表示有疾 病的序列,F20表示序列中有20帧,F34表示序列中有34帧。如图2(a)中的电影核磁共 振图像和图2(b)中的加标记的核磁共振图像比起来特征稀少,所以将图2(a)中的图像分解 成单演信号的振幅信息如图2(c),相位信息如图2(d),方位信息如图2(e),然后合成伪彩色图 像如图2(f)。图2(f)比图2(a)具有更多的特征信息。

本实施方式采用平均角度误差AAE(AverageAngleError)衡量运动估计效果

角度误差常用于直观的给出光流估计结果的好或者差。其计算方法为,对图像每一点 计算其估计的速度与光流真值速度的夹角,夹角越小表示角度结果越准确。平均角度误差 AAE(AverageAngleError),平均角度误差是一种重要的量化评估光流估计结果的方法。设 图像任意点p的角度误差为EAE(p))是两个流向量v0(p)=(u0,v0)和v1(p)=(u1,v1)在p点在 2D空间的角度(u0,v0)和(u1,v1).我们可以得到:

EAE=arccos(v0,v1)

角度误差适合用来评价大速度和小速度的位移。如果位移的结果给定了,那么EAE(p)是 估计值和真值之间的角度误差。AAE是EAE(p)平均值。

并且本实施方式采用平均终止点误差AEP(averageend-point)衡量运动估计效果 终止点误差(EndPointError)用来衡量两个光流终止点之间的距离,公式是

(u0-u1)2+(u1-v1)2

如果真实结果给定,那么v0(p)=(u0,v0)表示真值,v1(p)=(u1,v1)表示估计值。AEE (AverageEnd-pointError)是终止点误差的平均值。标准差(standarddeviation,STD)用来 评价运动估计结果的稳定性。

本实施方式主要采用光流算法,所以和目前效果较好的Sun光流算法,还有经典的 LK光流算法进行比较。比较效果如表1所示。

表1平均角度误差和平均终止点误差(均值±标准差)(角度误差:度,终止点误差:像素)

具体实施方式二:本实施方式与具体实施方式一不同的是:步骤一中空间正交滤波器 由1个偶数阶滤波器he(x)和2个奇数阶滤波器ho1(x)和ho2(x)组成;其中,x=(x,y)是电 影核磁共振图像中任意一点像素的坐标。其它步骤及参数与具体实施方式一相同。

具体实施方式三:本实施方式与具体实施方式一或二不同的是:步骤一中利用空间正 交滤波器(sphericalquadraturefilters,SQFs)提取单演信号的局部相位、局部方位和局部振幅 三种特征,将三种特征联合构造三维单演信号特征矩阵具体过程为如图3:

(1)、因为高频信息能够反映图像的细节,采用了巴特沃斯高通滤波器;在频率域中 公式如下:

He(ω)=11+(ωc/ω)2n---(2)

其中,ωc是滤波器的截止频率;n为巴特沃斯高通滤波器的阶数;He(·)是经傅里叶 变换后得到的频域下的偶数阶滤波器;

(2)、奇滤波器由偶滤波器计算获得,在频域里计算公式如下:

Ho1(ω)=-x|ω|·He(ω),Ho2(ω)=-y|ω|·He(ω)---(3)

ω=[ωxy]T是正则化的角频率,ωx为图像在x轴方向的正则化的角频率;ωy为图 像在y轴方向的正则化的角频率;j表示虚部;Ho1(ω)为ho1(x)经傅里叶变换后得到的频 域下的奇数阶滤波器;Ho2(ω)为ho2(x)经傅里叶变换后得到的频域下的奇数阶滤波器

(3)、将3个滤波器用于计算图像I的单演相位单演方位θ(x)和单演振幅A(x),计算 公式如下(见图1):

p(x)=(I*he)(x),q1(x)=(I*ho1)(x),q2(x)=(I*ho2)(x),q(x)=[q1(x),q2(x)]TA(x)=q1(x)2+q2(x)2+p(x)2;θ(x)=arctan(q2(x)q1(x))---(4)

其中,p(x)是图像I经过偶数阶滤波器变换的结果,q1(x)为图像I经过奇数阶滤波 器ho1变换的结果,q2(x)为图像I经过奇数阶滤波器ho2变换的结果,p(x),q1(x),q2(x) 构成正交三维向量空间,q(x)是q1(x)和q2(x)构成的向量;A(x)与p(x)的夹角为单演相 位q(x)与q1(x)的夹角为单演方位θ(x);*表示2维卷积;局部相位,局部方位和 局部振幅是互相独立的;

(4)、将局部相位如图2(b)、局部方位如图2(c)和局部振幅如图2(d)正则化为0到255 之间;然后将这3个特征单演相位单演方位θ(x)和单演振幅A(x)利用公式(5)构造成 3维的矩阵作为三维单演信号特征矩阵,计算公式如下:

其中,Im(x)是3种特征合成后的单演信号特征矩阵;其它步骤及参数与具体实施方 式一或二相同。

具体实施方式四:本实施方式与具体实施方式一至三之一不同的是:步骤二中利用基 于零均值归一化互相关系数的公式计算光流算法的输入图像局部特征的匹配特征项Ed(V) 作为光流算法的匹配特征项的具体过程:

数据的匹配公式是衡量两个像素或者区域的相似性;

步骤二一、根据三维单演信号特征矩阵,设R和T是具有相同维数的图像区域;SSD(·) 是一个常用的R和T相似度衡量标准,计算公式如下:

SSD(R,T)=ΣsN(Rs-Ts)2---(6)

其中,N是图像区域R或T所有像素的集合;Rs为图像区域R中的s点;Ts为图像区域T 中的s点;

步骤二二、在三维单演信号特征矩阵中引入ZNCC作为匹配特征项描述:

C(i):=[I(s)-μ(i)σ(i)],sNi---(7)

其中,Ni为图像区域中i点的邻域;C(i)为i像素的描述符;I(s)为图像区域Ni中s点的灰 度;μ(i)是图像区域所有像素的集合的灰度均值;σ(i)为图像区域所有像素的集合的灰度 方差;i∈Ω;其中,Ω为电影核磁共振图像的整个区域;

ZNCC(R,T)为相同维数的图像区域R和T的零均值归一化互相关系数;

步骤二三、根据公式(6)和(7)定义描述符的相似度;由像素i及其周围区域的相关变换求 得:

ZNCC(R,T)=1|N|·<R-μR,T-μT>σR·σT---(8)

|N|表示图像区域的大小,μR表示图像区域R的灰度平均值;μT表示图像区域T的灰 度平均值;σR为表示图像区域R的灰度标准差;σT表示图像区域T的灰度标准差;〈,〉表示 标准点乘;

步骤二四、如果R和T完全相同,那么ZNCC(R,T)度量值为1;最好的匹配是表达式 ZNCC(R,T)趋近于1;则对于图像I1和图像I2,不是之前点对点的匹配灰度,而是根据公 式(6)将I1的描述符C1和I2的描述符C2进行匹配;

步骤二五、针对每个点i定义的位移量dV=V-V0应该满足向量等式:

C2(i+dV(i))=C1(i);(9)

dV(i)表示V(i)的导数;V(i)为图像中第i像素点从图像I1到图像I2之间的位移距离;

步骤二六、基于这个等式C2(i+dV(i))=C1(i),定义匹配特征项Ed(V)为两个描述符距离 之和:

Ed(V)=ΣiΩ1|Ni|||C2(j+dV(i))-C1(i)||2---(10)

设I1,I2是序列中两幅相邻的图像;而且V=(u,v)是I1和I2之间的运动场; Ω={(x,y)|1≤x≤N,1≤y≤M}是标准的二维笛卡尔坐标,(N,M)表示图像的大小(图像的 长为N,宽为M);

每个像素有一个二维的索引i∈Ω和一个邻域命名为Ni;假设Vi是i点的位移向量;

以上的公式在位置i和索引k的匹配能量误差表示为凸函数如下:

Ed(V)=ΣiΩΣk(Ci(i,k)+C(i,k)·(Vi-V0,i))2

其中,Vi为像素点i的位移;k为RGB图像的3种颜色通道k=1,2,3;V0,i像素点i 的初始位置;Ci(i,k)表示在像素i点在k通道的描述符;为C的梯度;

总之,有两步建立两个像素/区域/特征之间的相关性,首先是建立描述符,然后是 求他们之间的距离。其它步骤及参数与具体实施方式一至三之一相同。

具体实施方式五:本实施方式与具体实施方式一至四之一不同的是:步骤三中限定电 影核磁共振图像的点的位移V的最小化匹配误差即电影核磁共振图像局部特征的平滑特 征项Es(V)公式为:

双边滤波对图像进行规范化;双边滤波是一种鲁棒性高的边缘保持滤波器,它是非线 性的并在平滑图像的同时保持边缘信息;它已经被用在很多图像处理和计算机视觉的应用 中;在本方法中,双边滤波被用来减少边缘的过度平滑和图像的噪声;双边滤波结合了 两种滤波器:空间滤波和域滤波用来分别用来衡量中心点和其相邻点的空间距离和域距 离;这两个滤波器通常用高斯分布;随着离中心点距离的增加,点的比重会下降;另外, 它们的单演信号特征差别越大,该点所占的比重就会越低;双边滤波简单的表示了像素点 空间上的差别和图像函数之间的差别;只要它们之间的距离是可以计算的,我们可以定义 多种函数;所以单演信号特征距离也是可以被应用于双边滤波就像应用于彩色图像一样;

光流算法公式的第二部分是正则化的部分,正则化部分是为了反映边缘向平缓区域传 播的过程;这个过程是基于图像的空间一致性的;也就是像素应该和周围的像素有着相同 的运动流;

根据光流的灰度不变假设平滑公式如下:

Es(j)=ΣsNiBFi,s·||Vs-Vi||1=ΣsNi{BFi,s·|us-ui|+BFi,s·|vs-vi|}

BFi,s为针对i像素或s像素的双边滤波;

Vs为图像中像素点s的位移;Vi图像中像素点i的位移;us为图像中像素点s的水平 方向位移;vs图像为I1或I2中像素点s的垂直方向的位移;ui为图像中像素点i的水平位 移;vi图像中像素点i的垂直方向的位移;

这里双边滤波BF度量同属一个目标的两个像素i和s相似度;如果s是像素i的邻 近像素而且它们的单演特征一致那么i和s很可能属于同一个目标,而且BFi,s接近1; 否则,BFi,s接近0,忽略这个像素;如果在理想状态,s和i同属于一个目标,那么它 们的运动也应该是一致的(Vs-Vi=0);

整体的平滑公式具体为:

Es(V)=ΣiΩΣsNiBFi,s·||Vs-Vi||1---(11)

其中,Vi为像素点i的位移,Vs为图像中像素点s的位移;

针对电影核磁共振图像,基于单演信号特征矩阵的双边滤波器定义如下:

BFi,s=e-(Δm(i,s)2σm2+Δd(i,s)2σd2)---(12)

其中,△m(i,s)表示i和s之间的单演信号特征距离,△d(i,s)=||i,s||2是i和s位置之 间的距离;参数σd和σm控制着相似性的测量;对于每个点i,s都表示它邻域内的所有 的点;σd为空间距离的标准差;σm单演信号特征距离的标准差。其它步骤及参数与具体 实施方式一至四之一相同。

去获取专利,查看全文>

相似文献

  • 专利
  • 中文文献
  • 外文文献
获取专利

客服邮箱:kefu@zhangqiaokeyan.com

京公网安备:11010802029741号 ICP备案号:京ICP备15016152号-6 六维联合信息科技 (北京) 有限公司©版权所有
  • 客服微信

  • 服务号