CN106373089B - B型与弹性双模态超声跟腱图像拼接方法 - Google Patents

Abstract

本发明公开一种B型与弹性双模态超声跟腱图像拼接方法。首先，基于边界矢量场提出一种刚性配准的新相似性测度，通过粒子群优化算法最大化该测度实现两幅B型超声跟腱图像的粗配准。其次，对B型图像进行非刚性配准以实现细配准，运用的技术为B样条自由形态形变模型或指数微分同胚模型。然后，利用与B型图像配准时同样的空间变换关系，对两幅弹性超声跟腱图像进行配准。接着，利用加权平均法分别对两幅B型与两幅弹性图像进行图像融合。最后，将多个跟腱节段的多幅超声图进行拼接，分别生成B型和弹性超声图像的全景图。该方法可准确拼接多幅B型与弹性超声图以提供跟腱的双模态全景图，将有助于跟腱损伤的诊断与治疗评估。

Description

B型与弹性双模态超声跟腱图像拼接方法

所属技术领域

本发明涉及一种图像拼接方法，尤其涉及一种B型与弹性双模态超声跟腱图像拼接方法，属于数字图像处理技术领域。

背景技术

跟腱(Achilles tendon,AT)是人体最长最强的肌腱。跟腱的主要功能是屈小腿和足跖屈，奔跑、跳跃的实现都依赖于跟腱。对跟腱损伤的位置和程度做准确诊断是制定有效的治疗方案及康复评价的重要依据。B型超声具有实时便捷、无电离辐射、成本低等优点，因此在跟腱检查中颇具应用价值。弹性超声(sonoelastography,SE)是利用超声测量组织弹性的新技术，它通过测量组织的应变或剪切波速度来反映组织弹性。与跟腱的B型超声相比，弹性超声提供了跟腱生物力学特性的补充信息。因此，关于跟腱SE成像的研究渐趋热门。

由于跟腱的长度较长，在超声检查时通常将其分为多个连续重叠的节段进行扫描。因此，医生难以快捷、准确的对跟腱疾病的严重性做出全面评价，尤其是在有多个病灶存在的情况下。为了准确判断跟腱疾病病情，医生希望能够在一张图上得到连续的跟腱图，即跟腱的全景图。本发明提出的拼接方法可以提供跟腱B型与弹性超声图像的全景图，完整展现跟腱全程的结构及弹性功能。

发明内容

本实用新型要解决的技术问题是提供一种B型与弹性双模态超声跟腱图像拼接方法。

为解决上述问题，本实用新型采用的技术方案是：

一种B型与弹性双模态超声跟腱图像拼接方法，包括以下步骤：

步骤1：B型超声图像粗配准：基于边界矢量场计算刚性配准的相似性测度，并通过粒子群优化算法最大化所述刚性配准的相似性测度，实现两幅B型超声跟腱图像的粗配准；

步骤2：B型超声图像细配准：运用B样条自由形态形变模型或微分同胚模型，对经过粗配准的两幅B型超声跟腱图像进行非刚性配准，以实现细配准；

步骤3：弹性超声图像配准：利用与B型超声跟腱图像配准相同的空间变换关系，对两幅弹性超声跟腱图像进行配准；

步骤4：两幅B型图融合与两幅弹性图融合：利用加权平均法分别实现两幅B型跟腱图像的融合和两幅弹性跟腱图像的融合；

步骤5：将多幅B型图中相邻的两幅逐一进行如步骤1、2、4的拼接，形成B型超声跟腱图像的全景图；将多幅弹性图中相邻的两幅逐一进行如步骤3、4的拼接，形成弹性超声跟腱图像的全景图。

步骤1中所述的B型超声图像粗配准方法包括以下具体步骤：

步骤1-1：建立特征空间：两幅B型超声跟腱图像分别生成各自的边缘增强图像，作为特征图，特征图的所有像素均作为刚性配准的图像特征点；

步骤1-2：计算相似性测度：所述特征图的归一化互信息作为相似性测度NMI_BVF：

NMI_BVF＝[H(u₁)+H(u₂)]/H(u₁,u₂)+[H(v₁)+H(v₂)]/H(v₁,v₂) (1)

式中，H(u_i)和H(v_i)分别是图像I_i的边界矢量场的水平和垂直分量的香农熵；i＝1,2分别代表参考图像I₁和目标图像I₂；H(u₁,u₂)和H(v₁,v₂)分别是边界矢量场的水平和垂直分量的联合熵；

步骤1-3：搜索策略：设定目标图像I₂平移、旋转(t_x,t_y,θ)量为粒子，t_x,t_y分别为水平平移量和垂直平移量，θ为旋转角度；在每次迭代后更新粒子位置时，对粒子进行两两交叉操作；计算参考图像I₁和变换后的目标图像I₃之间的归一化互信息NMI_BVF(t_x,t_y,θ)；当NMI_BVF(t_x,t_y,θ)达到最大值时，得到最优的变换矢量(t_x,t_y,θ)；

步骤1-4：公共区域定位：定位并截取参考图像I₁和目标图像I₂中相互重叠的部分作为参考图像I₁和目标图像I₂间的公共区域，参考图像I₁中与公共区域对应的部分图像为第一公共区域I₁’，目标图像I₂中与公共区域对应的部分图像为第二公共区域I₂’；

步骤1-5：对第一公共区域和第二公共区域进行刚性配准，获取两者间的刚性变换关系；并使用所述刚性变换对第二公共区域进行变换，得到变换后的第二公共区域I₃’。

步骤1-1中所述的边缘增强图像生成方法，包括以下具体步骤：

步骤1-1-1：对B型超声跟腱图像I(x,y)进行高斯滤波得到边界图；

步骤1-1-2：将边界图转换为二值图；

步骤1-1-3：分别沿着二值图的水平和垂直方向逐行逐列扫描计算水平势能函数和垂直势能函数；由水平势能函数和垂直势能函数的梯度计算BVF场v_BVF(x,y)＝[u(x,y),v(x,y)]，其中水平分量u(x,y)和垂直分量v(x,y)分别用于强调纵向和横向的组织边缘。

所述步骤1-4中以下具体步骤组成：

步骤1-4-1：设定比例参数n与m的初始值；比例参数n与m均为正整数且n<m；

步骤1-4-2：截取参考图像I₁的一侧n/m部分图像作为第一公共区域I₁’，第一公共区域I₁’长度为d₁＝floor(d×n/m)；相应的目标图像I₂的对应侧长度为d₁的部分作为第二公共区域I₂’；

步骤1-4-3：计算n₁＝floor(t_x×m/d)，若n₁不等于0，则n＝n+n₁,并跳转到步骤1-4-2,若n₁等于0则执行步骤1-4-4；横向位移t_x通过粒子群优化算法计算得到；

步骤1-4-4：由比例参数n与m得到第一公共区域I₁’与第二公共区域I₂’。

所述步骤2中设定B样条FFD非刚性配准中第一公共区域I₁’与第二公共区域I₃’之间的非刚性位移场P＝(p,q)的代价函数为：

位移场(p,q)通过均匀分布的控制点插值计算得到。

所述步骤2中，使用对数微分同胚模型实现细配准的方法，由以下具体步骤组成：

步骤2-1：对数微分同胚的算法的目标能量函数为：

式中，u为更新场，s为形变；

步骤2-2：设定初始化位移场v；

步骤2-3：计算形变s＝exp(v)，通过最小化目标能量函数得到对应的u；

步骤2-4：更新位移场v，v←u+v；

步骤2-5：对更新场u做一次高斯卷积：u←K_diff*u；K_diff是高斯卷积核；

步骤2-6：判断是否满足收敛阈值，如果是，输出当前形变s，否则，转向步骤2-3；

步骤2-7：目标图像I₂在公共区域经粗配准与细配准后变换成公共区域目标图像I₄’，在非公共区域经同样的刚性变换及最近邻插值后的非刚性变换，得到非公共区域目标图像I₄”。

步骤4中，第一公共区域I₁’与公共区域目标图像I₄’进行线性加权融合，得到融合图像I_s＝αI₁’+(1-α)I₄’；渐变系数α在第一公共区域最左端和最右端分别取1和0；中间位置则从1线性渐变到0。

采用上述技术方案，本发明产生的有益效果在于：

本发明可准确拼接多幅B型与弹性超声图以提供跟腱的双模态全景图，将有助于跟腱损伤的诊断与治疗评估。

附图说明

图1是本发明超声跟腱图像拼接的流程图；

图2是本发明B型超声跟腱图像的第一公共区域；

图3是本发明B型超声跟腱图像的第二公共区域；

图4是本发明B型超声跟腱图像经粒子群优化的配准结果图；

图5是本发明B型超声跟腱图像经B样条FFD非刚性配准的结果图；

图6是本发明中的图像融合示意图；

图7是本发明B型超声跟腱图像经对数微分同胚非刚性配准的结果图。

具体实施方式

下面结合附图和具体实施例对本发明进行详细地描述，

实施例一：

如图1所示，一种B型与弹性双模态超声跟腱图像拼接方法，包括以下步骤：

步骤1：B型超声图像粗配准：基于边界矢量场计算刚性配准的相似性测度，并通过粒子群优化算法最大化所述刚性配准的相似性测度，实现两幅B型超声跟腱图像的粗配准；

步骤2：B型超声图像细配准：运用B样条自由形态形变模型或微分同胚模型，对经过粗配准的两幅B型超声跟腱图像进行非刚性配准，以实现细配准；

步骤3：弹性超声图像配准：利用与B型超声跟腱图像配准相同的空间变换关系，对两幅弹性超声跟腱图像进行配准；

步骤4：两幅B型图融合与两幅弹性图融合：利用加权平均法分别实现两幅B型跟腱图像的融合和两幅弹性跟腱图像的融合；

步骤5：将多幅B型图中相邻的两幅逐一进行如步骤1、2、4的拼接，形成B型超声跟腱图像的全景图；将多幅弹性图中相邻的两幅逐一进行如步骤3、4的拼接，形成弹性超声跟腱图像的全景图。

步骤1中所述的B型超声图像粗配准方法包括以下具体步骤：

步骤1-1：建立特征空间：两幅B型超声跟腱图像分别生成各自的边缘增强图像，作为特征图，特征图的所有像素均作为刚性配准的图像特征点；

步骤1-2：计算相似性测度：所述特征图的归一化互信息作为相似性测度NMI_BVF：

NMI_BVF＝[H(u₁)+H(u₂)]/H(u₁,u₂)+[H(v₁)+H(v₂)]/H(v₁,v₂) (1)

式中，H(u_i)和H(v_i)分别是图像I_i的边界矢量场的水平和垂直分量的香农熵；i＝1,2分别代表参考图像I₁和目标图像I₂；H(u₁,u₂)和H(v₁,v₂)分别是边界矢量场的水平和垂直分量的联合熵；

步骤1-3：搜索策略：设定目标图像I₂平移、旋转(t_x,t_y,θ)量为粒子，t_x,t_y分别为水平平移量和垂直平移量，θ为旋转角度；在每次迭代后更新粒子位置时，对粒子进行两两交叉操作；计算参考图像I₁和变换后的目标图像I₃之间的归一化互信息NMI_BVF(t_x,t_y,θ)；当NMI_BVF(t_x,t_y,θ)达到最大值时，得到最优的变换矢量(t_x,t_y,θ)；

步骤1-4：公共区域定位：定位并截取参考图像I₁和目标图像I₂中相互重叠的部分作为参考图像I₁和目标图像I₂间的公共区域，参考图像I₁中与公共区域对应的部分图像的为第一公共区域I₁’，目标图像I₂中与公共区域对应的部分图像为第二公共区域I₂’；

步骤1-5：对第一公共区域和第二公共区域进行刚性配准，获取两者间的刚性变换关系；并使用所述刚性变换对第二公共区域变换，得到变换后的第二公共区域I₃’。

步骤1-1中所述的边缘增强图像生成方法，包括以下具体步骤：

步骤1-1-1：对B型超声跟腱图像I(x,y)进行高斯滤波得到边界图；

步骤1-1-2：将边界图转换为二值图；

步骤1-1-3：分别沿着二值图的水平和垂直方向逐行逐列扫描计算水平势能函数和垂直势能函数；由水平势能函数和垂直势能函数的梯度计算BVF场v_BVF(x,y)＝[u(x,y),v(x,y)]，其中水平分量u(x,y)和垂直分量v(x,y)分别用于强调纵向和横向的组织边缘。

所述步骤1-4中以下具体步骤组成：

步骤1-4-1：设定比例参数n与m的初始值；比例参数n与m均为正整数且n<m；

步骤1-4-2：截取参考图像I₁的一侧n/m部分图像作为第一公共区域I₁’，第一公共区域I₁’长度为d₁＝floor(d×n/m)；相应的目标图像I₂的对应侧长度为d₁的部分作为第二公共区域I₂’；

步骤1-4-3：计算n₁＝floor(t_x×m/d)，若n₁不等于0，则n＝n+n₁,并跳转到步骤1-4-2,若n₁等于0则执行步骤1-4-4；横向位移t_x通过粒子群优化算法计算得到；

步骤1-4-4：由比例参数n与m得到第一公共区域I₁’与第二公共区域I₂’。

所述步骤2中设定B样条FFD非刚性配准中第一公共区域I₁’与变换后的第二公共区域I₃’之间的非刚性位移场P＝(p,q)的代价函数为：

位移场(p,q)通过均匀分布的控制点插值计算得到。

步骤4中，第一公共区域I₁’与公共区域目标图像I₄’进行线性加权融合，得到融合图像I_s＝αI₁’+(1-α)I₄’；渐变系数α在第一公共区域最左端和最右端分别取1和0；中间位置则从1线性渐变到0。

拼接顺序可以“从左到右”，也可以全部调换为“从右到左”。

实施例二：

实施例二共有5个步骤，其中步骤1、3、4、5与实施例一的相应步骤相同，步骤2为：对B型图像进行非刚性配准以实现细配准，运用的技术为对数微分同胚模型。

所述步骤2中，使用对数微分同胚模型实现细配准的方法，由以下具体步骤组成：

步骤2-1：对数微分同胚的算法的目标能量函数为：

式中，u为更新场，s为形变；

步骤2-2：设定初始化位移场v；

步骤2-3：计算形变s＝exp(v)，通过最小化目标能量函数得到对应的更新场u；

步骤2-4：更新位移场v，v←u+v；

步骤2-5：对更新场u做一次高斯卷积：u←K_diff*u；K_diff是高斯卷积核；

步骤2-6：判断是否满足收敛阈值，如果是，输出当前形变s，否则，转向步骤2-3；

步骤2-7：目标图像I₂在公共区域经粗配准与细配准后变换成公共区域目标图像I₄’，在非公共区域经同样的刚性变换及最近邻插值后的非刚性变换，得到非公共区域目标图像I₄”。

最后应说明的是：以上各实施例仅用以说明本发明的技术方案，而非对其限制；尽管参照前述各实施例对本发明进行了详细的说明，本领域的普通技术人员应当理解：其对前述各实施例所记载的技术方案进行修改，或者对其中部分或者全部技术特征进行等同替换，并不使相应技术方案的本质脱离本发明各实施例技术方案的范围。

Claims (5)

1.一种B型与弹性双模态超声跟腱图像拼接方法，其特征在于：包括以下步骤：

步骤1：B型超声图像粗配准：基于边界矢量场计算刚性配准的相似性测度，并通过粒子群优化算法最大化所述刚性配准的相似性测度，实现两幅B型超声跟腱图像的粗配准，包括五个子步骤，即：

步骤1-1：建立特征空间：两幅待拼接双模态超声跟腱图像的B型图像I(x,y)分别生成边缘增强图像，其中x，y表示横纵坐标，包括对I(x,y)进行高斯滤波得到边界图，将边界图转换为二值图，分别沿着二值图的水平和垂直方向逐行逐列扫描计算水平势能函数和垂直势能函数，由水平势能函数和垂直势能函数的梯度计算边界矢量场v_BVF(x,y)＝[u(x,y),v(x,y)]，其中水平分量u(x,y)和垂直分量v(x,y)分别用于强调纵向和横向的组织边缘，即为边缘增强图像；将所述边缘增强图像作为特征图，特征图的所有像素均作为刚性配准的图像特征点；

步骤1-2：计算相似性测度：所述特征图的归一化互信息作为相似性测度NMI_BVF：

NMI_BVF＝[H(u₁)+H(u₂)]/H(u₁,u₂)+[H(v₁)+H(v₂)]/H(v₁,v₂) (1)

式中，H(u_i)和H(v_i)分别是图像I_i的边界矢量场的水平和垂直分量的香农熵，i＝1,2分别代表参考图像I₁和目标图像I₂；H(u₁,u₂)和H(v₁,v₂)分别是边界矢量场的水平和垂直分量的联合熵；

步骤1-3：建立搜索策略：设定目标图像I₂平移、旋转(t_x,t_y,θ)为粒子，t_x,t_y分别为水平平移量和垂直平移量，θ为旋转角度；在每次迭代后更新粒子位置时，对粒子进行两两交叉操作；计算参考图像I₁和变换后的目标图像I₃之间的归一化互信息NMI_BVF(t_x,t_y,θ)；当NMI_BVF(t_x,t_y,θ)达到最大值时，得到最优的变换矢量(t_x,t_y,θ)；

步骤1-4：公共区域定位：定位并截取参考图像I₁和目标图像I₂中相互重叠的部分作为参考图像I₁和目标图像I₂间的公共区域，参考图像I₁中与公共区域对应的部分图像为第一公共区域I₁’，目标图像I₂中与公共区域对应的部分图像为第二公共区域I₂’；

步骤1-5：公共区域变换：对第一公共区域和第二公共区域进行刚性配准，获取两者间的刚性变换关系；并使用所述刚性变换对第二公共区域变换，得到变换后的第二公共区域I₃’；

步骤2：B型超声图像细配准：运用B样条自由形态形变模型或微分同胚模型，对经过粗配准的两幅B型超声跟腱图像进行非刚性配准，以实现细配准；目标图像I₂在公共区域经粗配准与细配准后变换成公共区域目标图像I₄’，在非公共区域经同样的刚性变换及最近邻插值后的非刚性变换，得到非公共区域目标图像I₄”；

步骤3：弹性超声图像配准：利用与B型超声跟腱图像配准相同的空间变换关系，对两幅弹性超声跟腱图像进行配准；

步骤4：两幅B型图融合与两幅弹性图融合：利用加权平均法分别实现两幅B型跟腱图像的融合和两幅弹性跟腱图像的融合；

步骤5：将多幅B型图中相邻的两幅逐一进行如步骤1、2、4的拼接，形成B型超声跟腱图像的全景图；将多幅弹性图中相邻的两幅逐一进行如步骤3、4的拼接，形成弹性超声跟腱图像的全景图。

2.根据权利要求1所述的B型与弹性双模态超声跟腱图像拼接方法，其特征在于：

所述步骤1-4中以下具体步骤组成：

步骤1-4-1：设定比例参数n与m的初始值；比例参数n与m均为正整数且n<m；

步骤1-4-2：截取参考图像I₁的一侧n/m部分图像作为预设的第一公共区域I₁’，I₁’长度为d₁＝floor(d×n/m)，其中d为I₁矩阵的列数；I₂的对应侧长度为d₁的部分作为预设的第二公共区域I₂’；

步骤1-4-3：通过粒子群优化算法计算得到横向位移t_x；计算n₁＝floor(t_x×m/d)，若n₁不等于0，则n＝n+n₁,并跳转到步骤1-4-2,若n₁等于0则执行步骤1-4-4；

步骤1-4-4：由比例参数n与m得到第一公共区域I₁’与第二公共区域I₂’。

3.根据权利要求1所述的B型与弹性双模态超声跟腱图像拼接方法，其特征在于：

所述步骤2中设定B样条FFD非刚性配准中第一公共区域I₁’与刚性变换后的第二公共区域I₃’之间的非刚性位移场P＝(p,q)的代价函数为：

位移场(p,q)通过均匀分布的控制点插值计算得到。

4.根据权利要求1所述的B型与弹性双模态超声跟腱图像拼接方法，其特征在于：

所述步骤2中，使用对数微分同胚模型实现细配准的方法，由以下具体步骤组成：

步骤2-1：对数微分同胚的算法的目标能量函数为：

式中，u为更新场，s为形变；

步骤2-2：设定初始化位移场v；

步骤2-3：计算形变s＝exp(v)，通过最小化目标能量函数得到对应的更新场u；

步骤2-4：更新位移场v，v←u+v；

步骤2-5：对更新场u做一次高斯卷积：u←K_diff*u；K_diff是高斯卷积核；

步骤2-6：判断是否满足收敛阈值，如果是，输出当前形变s，否则，转向步骤2-3。

5.根据权利要求1所述的B型与弹性双模态超声跟腱图像拼接方法，其特征在于：

步骤4中，第一公共区域I₁’与公共区域目标图像I₄’进行线性加权融合，得到融合图像I_s＝αI₁’+(1-α)I₄’；渐变系数α在第一公共区域最左端和最右端分别取1和0；中间位置则从1线性渐变到0。