CN106991694A - 基于显著区域面积匹配的心脏ct与超声图像配准方法 - Google Patents

Abstract

本发明公开了一种基于显著区域面积匹配的心脏CT与超声图像配准方法，主要解决现有技术配准精度低、速度慢的问题。其实现步骤为：1)将预处理后的CT和超声图像，分别作为参考图像和浮动图像，分别建立这两幅图像的灰度特征金字塔模型和邻域均差特征金字塔模型；2)基于这两种金字塔模型分别生成参考图像和浮动图像的显著图，并对其进行二值化；3)对于二值化结果，根据区域面积特征提取其感兴趣区域，并对感兴趣区域进行基于区域质心距离的粗配准和基于ICP算法的精配准，得到配准参数；5)根据配准参数，对浮动图像进行旋转和平移变换，得到配准结果。本发明配准精度高，速度快，可用于对术前CT图像和术中超声图像的实时配准。

Description

基于显著区域面积匹配的心脏CT与超声图像配准方法

技术领域

本发明属于医学图像处理技术领域，特别是一种图像的配准方法，可用于心脏手术中对术前CT图像和术中超声图像的实时配准。

背景技术

随着现代社会生活节奏的不断加快、生活压力的不断增大，心血管疾病的发病率逐年攀升。据世界卫生组织的调查数据显示，每年心血管疾病的死亡人数占全球死亡总数的30％。预计到2030年，死于心血管疾病的人数将增加至2330万，心血管疾病将继续成为威胁人类健康的头号死因。因此对心血管疾病的诊断研究受到广泛关注，更加准确、高效的医疗诊断成为一种迫切的需要。

目前，对于心脏的诊疗和评估高度依赖于CT、超声和核磁共振等成像技术。其中，心脏超声成像由于成本相对较低、能够实时成像和操作方便等原因，被广泛应用于微创手术中。它可以帮助医生在手术过程中定位病理区域，并为手术工具提供引导。而超声图像固有的斑点噪声、伪影、分辨率差、信噪比较低等局限，使得它难以清晰的显示器官、血管和肿瘤的解剖结构和病理信息。CT图像虽然成像精度和信噪比高，但由于成本高、灵活性差、不能实时成像等缺点，无法反映手术中组织的移位和变形。

将术前CT图像与术中超声图像两种模态图像信息相结合，利用术前图像的信息来引导术中图像，提高手术的精度和准确率，是当前研究的热点。其中术前CT图像和术中超声图像配准技术是图像引导手术中的关键技术。

超声图像与CT图像的配准方法，目前主要包括基于灰度的方法和基于几何或解剖特征的方法。基于灰度的配准方法直接利用整幅图像的灰度信息，建立两幅图像之间的相似性度量，然后使用搜索方法寻找使相似性度量值最大或最小的变化模型的参数值。然而，由于超声图像质量较差、与CT图像灰度差异大，单一的基于灰度的配准精度并不理想且耗时长，不能满足图像引导手术对图像配准精度和速度的要求。基于特征的配准方法将对整个图像的分析转换为对图像特征的分析，所说大大减少了图像处理过程中的计算量，具有较好的实时性，但配准精度的高低主要取决于几何或解剖特征的选择和提取，特征提取的准确程度将对整个配准过程产生很大的影响，由于常用的特征提取主要依赖于人工选取，速度慢且会引入人为偏差，使提取的不同模态下图像特征差异较大。因此迫切需要准确可靠、全自动的特征提取方法以确保图像配准的精度和速度，满足手术要求。

发明内容

本发明的目的在于针对目前已有配准技术的不足，提出一种基于显著区域面积匹配的心脏CT图像与超声图像配准方法，以实现CT图像对超声图像的精确引导，提高手术的精度和成功率。

为实现上述目的，本发明包括如下步骤：

(1)输入待配准心脏CT图像及食道超声心动图TEE，采用双线性插值方法调整CT图像，使其与TEE图像每个像素对应的物理分辨率大小一致；

(2)对CT图像和TEE图像分别进行非局部均值滤波处理，并将滤波后的CT图像作为参考图像r，将滤波后的TEE图像作为浮动图像f；

(3)分别建立参考图像r的灰度特征金字塔模型A_I、邻域均差特征金字塔模型A_M和浮动图像f的灰度特征金字塔模型B_I、邻域均差特征金字塔模型B_M；

(4)根据参考图像r的灰度特征金字塔模型A_I和邻域均差特征金字塔模型A_M生成参考图像的显著图S_r，根据浮动图像f的灰度特征金字塔模型B_I和邻域均差特征金字塔模型B_M生成浮动图像的显著图S_f；

(5)设置灰度阈值T₁＝85，利用灰度阈值T₁分别对参考图像的显著图S_r和浮动图像的显著图S_f进行二值化分割，得到参考图像显著图的分割结果b_r和浮动图像显著图的分割结果b_f；

(6)分别统计参考图像显著图的分割结果b_r和浮动图像显著图的分割结果b_f中每个非连通区域的面积，根据非连通区域的面积特征提取参考图像的感兴趣区域轮廓点集p_r和浮动图像的感兴趣区域轮廓点集p_f；

(7)分别提取参考图像r中感兴趣区域的质心(x_r,y_r)和浮动图像f中感兴趣区域的质心(x_f,y_f)，计算两区域质心间的距离(Δx,Δy)，并将浮动图像f和其感兴趣区域轮廓点集p_f分别横向平移Δx个像素，纵向平移Δy个像素，得到粗配准结果图f'和粗配准结果图f'的感兴趣区域轮廓点集p'_f；

(8)以参考图像r的感兴趣区域轮廓点集p_r和粗配准结果图f'的感兴趣区域轮廓点集p'_f作为输入，使用迭代最近点ICP算法进行精配准，得到两点集之间的旋转矩阵TR和平移矩阵TT；

(9)将粗配准结果图f'按旋转矩阵TR和平移矩阵TT依次进行旋转和平移操作，得到最终配准结果q，输出配准图像q。

本发明与现有技术相比具有以下优点：

1、本发明通过分别建立参考图像和浮动图像的灰度特征金字塔模型和邻域均差特征金字塔模型以提取其显著图，并根据区域面积特征精确地提取了参考图像和浮动图像的感兴趣区域，为后续参考图像和浮动图像的高精度配准奠定了基础。

2、本发明对参考图像和浮动图像进行了基于参考图像和浮动图像感兴趣区域质心距离的粗配准和基于迭代最近点ICP算法的精配准，具有更高的准确性。

3、本发明在整个配准过程中实现了全自动化，无需人机交互，避免了人为因素的干扰，配准效率高，满足图像引导手术对配准速度和配准精度的要求。

附图说明

图1是本发明的实现流程图；

图2是从某医院获取的心脏CT图像和心脏食道超声TEE图像；

图3是用本发明提取CT图像和TEE图像的显著图；

图4是用本发明对CT图像和TEE图像的显著图进行二值化后的结果图；

图5是用本发明对CT图像和TEE图像感兴趣区域的分割结果图；

图6是以CT图为参考图像，TEE图为浮动图像，采用本发明进行配准的结果图；

图7是本发明所使用的同一个人的心脏CT图像和超声TEE图像；

图8是用本发明对图7中CT图像和超声TEE图像进行配准的结果图。

具体实施方法

以下结合附图对本发明的实施例和效果作详细描述。

参照图1，本发明的实现步骤如下：

步骤1.输入心脏CT图像和食道超声TEE图像，调整图像像素大小。

心脏CT图像和食道超声图像是从某医院获取的，如图2所示，其中图2(a)是心脏CT图像，图2(b)是食道超声TEE图像，心脏CT图像的大小是512×512，像素分辨率是0.24×0.24mm，食道超声TEE图像的大小是600×800，像素分辨率是0.17×0.17mm；

采用双线性插值方法调整CT图像，使其与TEE图像像素分辨率大小一致。

步骤2.对待配准图像进行预处理。

将心脏CT图像进行非局部均值滤波，使图像更平滑、清晰，并用滤波后的图像作为参考图像r；

将食道超声TEE图像进行非局部均值滤波，以去除其固有的斑点、伪影这些噪声，使图像更平滑、清晰，并用该滤波后的图像作为浮动图像f。

步骤3.分别建立参考图像r和浮动图像f的灰度特征金字塔模型和邻域均差特征金字塔模型。

(3a)提取参考图像r的灰度特征图I_r和邻域均差特征图M_r：

M_r(i,j)＝|r(i,j)-F(i,j)|

式中I_r(i,j)是参考图像r的灰度特征图I_r在坐标(i,j)处的值，r(i,j)是参考图像r在坐标(i,j)处的灰度值，max()为提取矩阵最大灰度值，min()为提取矩阵最小灰度值，M_r是邻域均差特征图，M_r(i,j)是M_r在坐标(i,j)处的值，F(i,j)是参考图像r在坐标(i,j)处的3×3邻域均值；

(3b)以参考图像r的灰度特征图I_r为第一层，在其基础上进行下采样，得到大小为第一层图像1/2的第二层灰度特征图，再在第二层灰度特征图上进行相同下采样，得到大小为第二层灰度特征图1/2的第三层灰度特征图，并将第一层灰度特征图放置在最下层，第二层灰度特征图置于第一层之上，再将第三层灰度特征图置于第二层之上，构成参考图像r的灰度特征金字塔模型A_I；

(3c)以参考图像r的邻域均差特征图M_r为第一层，在其基础上进行下采样，得到大小为第一层图像1/2的邻域均差特征图，再在第二层邻域均差特征图上进行相同下采样，得到大小为第二层邻域均差特征图1/2的第三层邻域均差特征图，并将第一层邻域均差特征图放置在最下层，第二层邻域均差特征图置于第一层之上，再将第三层邻域均差特征图置于第二层之上，构成参考图像r的邻域均差特征金字塔模型A_M；

(3d)按(3a)相同操作提取浮动图像f的灰度特征图I_f和邻域均差特征图M_f；

(3e)以浮动图像f的灰度特征图I_f为第一层，在其基础上进行下采样，得到大小为第一层图像1/2的第二层灰度特征图，再在第二层灰度特征图上进行相同下采样，得到大小为第二层灰度特征图1/2的第三层灰度特征图，并将第一层灰度特征图放置在最下层，第二层灰度特征图置于第一层之上，再将第三层灰度特征图置于第二层之上，构成浮动图像f的灰度特征金字塔模型B_I；

(3f)以浮动图像f的邻域均差特征图M_f为第一层，在其基础上进行下采样，得到大小为第一层图像1/2的第二层邻域均差特征图，再在第二层邻域均差特征图上进行相同下采样，得到大小为第二层邻域均差特征图1/2的第三层邻域均差特征图，并将第一层邻域均差特征图放置在最下层，第二层邻域均差特征图置于第一层之上，再将第三层邻域均差特征图置于第二层之上，构成浮动图像f的邻域均差特征金字塔模型B_M。

步骤4.生成参考图像和浮动图形的显著图。

目前生成显著图的方法比较成熟，常用的有：基于空间频域分析的FT算法、基于图论的GBVS算法和基于分块的C-A算法。本发明采用但不限于下述方法生成参考图像的显著图S_r和浮动图像的显著图S_f，其步骤如下：

(4a)将参考图像灰度特征金字塔模型A_I中的第一层图像分别与δ₁＝32和δ₂＝64的两个高斯核进行二维卷积，用该层图像减去卷积后图像，再对两个差分结果图求和并归一化得到第一层图像显著图；

(4b)将参考图像灰度特征金字塔模型A_I中的第二、第三层图像分别与δ₁＝32和δ₂＝64的两个高斯核进行二维卷积，用该层图像减去卷积后图像，再对两个差分结果图求和并归一化得到第二、第三层图像显著图；

(4c)对参考图像灰度特征金字塔模型A_I中三层显著图求平均，得到参考图像r的灰度特征显著图S_rI；

(4d)对参考图像r邻域均差特征金字塔模型A_M采用(4a)-(4c)相同操作，得到参考图像r的邻域均差特征显著图S_rM；

(4e)对参考图像r的灰度特征显著图S_rI和邻域均差特征显著图S_rM求和，得到参考图像r的最终显著图S_r；

(4f)对浮动图像f的灰度特征金字塔模型B_I和邻域均差特征金字塔模型B_M采用(4a)-(4e)相同操作，得到浮动图像f的最终显著图S_f，如图3所示，其中图3(a)是提取CT图像的显著图，图3(b)是提取TEE图像的显著图。

步骤5.对参考图像的显著图S_r和浮动图像的显著图S_f进行二值化。

设置灰度阈值T₁＝85，利用灰度阈值T₁分别对参考图像的显著图S_r和浮动图像的显著图S_f进行二值化分割：若图像像素点灰度值大于灰度阈值T₁，则将该像素点灰度值赋值为1，否则赋值为0，得到参考图像显著图的分割结果b_r和浮动图像显著图的分割结果b_f，如图4所示，其中图4(a)是对CT图像的二值化结果，图4(b)是对TEE图像的二值化结果。

步骤6.提取参考图像r和浮动图像f感兴趣区域轮廓点集p_r和p_f。

(6a)统计参考图像最终显著图的分割结果b_r中各个非连通区域的面积，根据面积大小对各连通区域按降序排序；

(6b)按(6a)中排序顺序，从第二个非连通区域开始计算各区域与其前一区域面积的比值，并将该比值与设置的面积比阈值T₂＝0.25进行比较：若该比值大于T₂，则标记该区域属于感兴趣区域，并继续判断下一区域是否属于感兴趣区域，否则，提取所有在迭代过程中标记为感兴趣区域的区域轮廓点，作为参考图像r的感兴趣区域的轮廓点集p_r；

(6c)统计浮动图像最终显著图的分割结果b_f中各非连通区域的面积，根据面积大小对各连通区域按降序排序；

(6d)按(6c)中排序顺序，从第二个非连通区域开始计算各区域与其前一区域面积的比值，并将该比值与设置的面积比阈值T₂＝0.25进行比较：若该比值大于T₂，则标记该区域属于感兴趣区域，并继续判断下一区域是否属于感兴趣区域，否则，提取所有在迭代过程中标记为感兴趣区域的区域轮廓点，作为浮动图像f的感兴趣区域的轮廓点集p_f，如图5所示，其中图5(a)表示提取的参考图像感兴趣区域，图5(b)表示提取的浮动图像感兴趣区域。

步骤7.根据参考图像感兴趣区域和浮动图像感兴趣区域质心间的距离对参考图像和浮动图像进行粗配准。

分别提取参考图像r中感兴趣区域的质心(x_r,y_r)和浮动图像f中感兴趣区域的质心(x_f,y_f)，计算两区域质心间的距离(Δx,Δy)，并将浮动图像f和其感兴趣区域轮廓点集f分别横向平移Δx个像素，纵向平移Δy个像素，得到粗配准结果图f'和粗配准结果图的感兴趣区域轮廓点集p'_f。

步骤8.采用迭代最近点ICP算法进行精配准，得到配准参数。

以参考图像r的感兴趣区域轮廓点集p_r和粗配准结果图f'的感兴趣区域轮廓点集p'_f作为输入，使用迭代最近点ICP算法进行精配准，得到两点集之间的旋转矩阵TR和平移矩阵TT。

步骤9.根据配准参数得到最终配准结果。

将粗配准结果图f'按旋转矩阵TR和平移矩阵TT依次进行旋转和平移操作，得到最终配准结果q，完成配准工作，最终配准结果图如图6所示。

本发明的效果可通过以下仿真实验进一步说明：

实验内容：

为了说明本发明的精确性和高效性，实验所使用的输入图像为来自同一个人的心脏CT图像和超声TEE图像，如图7所示，其中图7(a)是CT图像，图7(b)是TEE图像。

以7(a)所示的CT图像为参考图像，以图7(b)所示的TEE图像为浮动图像，采用本发明对其进行配准，配准结果如图8所示，从图8可以看出，本发明对CT图像和超声图像配准的精确度高。

对图像配准前后目标配准误差和配准时间进行统计，结果如表一。

表一评价指标

由表一可以看出，采用本发明进行配准的目标配准误差与配准前相比有大幅提高，而且运行时间短，满足图像引导手术对配准精度和配准时间的要求。

Claims (5)

1.基于显著区域面积匹配的心脏CT与超声图像配准方法，包括如下步骤：

(1)输入待配准心脏CT图像及食道超声心动图TEE，采用双线性插值方法调整CT图像，使其与TEE图像像素大小一致；

(2)对CT图像和TEE图像分别进行非局部均值滤波处理，并将滤波后的CT图像作为参考图像r，将滤波后的TEE图像作为浮动图像f；

(3)分别建立参考图像r的灰度特征金字塔模型A_I、邻域均差特征金字塔模型A_M和和浮动图像f的灰度特征金字塔模型B_I、邻域均差特征金字塔模型B_M；

(4)根据参考图像r的灰度特征金字塔模型A_I和邻域均差特征金字塔模型A_M生成参考图像的显著图S_r，根据浮动图像f的灰度特征金字塔模型B_I和邻域均差特征金字塔模型B_M生成浮动图像的显著图S_f；

(5)设置灰度阈值T₁＝85，利用灰度阈值T₁分别对参考图像的显著图S_r和浮动图像的显著图S_f进行二值化分割，得到参考图像显著图的分割结果b_r和浮动图像显著图的分割结果b_f；

(6)分别统计参考图像显著图的分割结果b_r和浮动图像显著图的分割结果b_f中各非连通区域的面积，根据非连通区域的面积特征提取参考图像的感兴趣区域轮廓点集p_r和浮动图像的感兴趣区域轮廓点集p_f；

(7)分别提取参考图像r中感兴趣区域的质心(x_r,y_r)和浮动图像f中感兴趣区域的质心(x_f,y_f)，计算两区域质心间的距离(Δx,Δy)，并将浮动图像f和其感兴趣区域轮廓点集p_f分别横向平移Δx个像素，纵向平移Δy个像素，得到粗配准结果图f'和粗配准结果图的感兴趣区域轮廓点集p'_f；

(8)以参考图像r的感兴趣区域轮廓点集p_r和粗配准结果图f'的感兴趣区域轮廓点集p'_f作为输入，使用迭代最近点ICP算法进行精配准，得到两点集之间的旋转矩阵TR和平移矩阵TT；

(9)将粗配准结果图f'按旋转矩阵TR和平移矩阵TT依次进行旋转和平移操作，得到最终配准结果q，完成配准工作。

2.根据权利要求1所述的方法，其中步骤(3)中分别建立参考图像r和浮动图像f的灰度特征金字塔模型和邻域均差特征金字塔模型，按如下步骤进行：

(3a)提取参考图像r的灰度特征图I_r和邻域均差特征图M_r；

(3b)以灰度特征图I_r为第一层，在其基础上进行下采样，得到大小为第一层图像1/2的第二层灰度特征图，再在第二层灰度特征图上进行相同下采样，得到大小为第二层灰度特征图1/2的第三层灰度特征图，以此三层灰度特征图建立参考图像r的灰度特征金字塔模型A_I；

(3c)按照与(3b)相同的操作建立参考图像r的邻域均差特征金字塔模型A_M；

(3d)提取浮动图像f的灰度特征图I_f和邻域均差特征图M_f；

(3e)按照与(3b)相同的操作建立浮动图像f的灰度特征金字塔模型B_I和邻域均差特征金字塔模型B_M。

3.根据权利要求2所述的方法，其中步骤(3a)中提取参考图像r的灰度特征图I_r和邻域均差特征图M_r，通过如下公式进行：

$I_r(i,j)=\frac{r(i,j)-min\left(I\right(i,j\left)\right)}{max\left(I\right(i,j\left)\right)-min\left(I\right(i,j\left)\right)}$

M_r(i,j)＝|r(i,j)-F(i,j)|

式中I_r(i,j)是参考图像r的灰度特征图I_r在坐标(i,j)处的值，r(i,j)是参考图像r在坐标(i,j)处的灰度值，max()为提取矩阵最大灰度值，min()为提取矩阵最小灰度值，M_r是邻域均差特征图，M_r(i,j)是M_r在坐标(i,j)处的值，F(i,j)是参考图像r在坐标(i,j)处的3×3邻域均值。

4.根据权利要求1所述方法，其中步骤(4)中生成参考图像的显著图S_r和浮动图像的显著图S_f，按如下步骤进行：

(4a)将参考图像灰度特征金字塔模型A_I中的第一层图像分别与δ₁＝32和δ₂＝64的两个高斯核进行二维卷积，用该层图像减去卷积后图像，再对两个差分结果图求和并归一化得到第一层图像显著图；

(4b)采用(4a)相同的操作得到参考图像灰度特征金字塔模型A_I中第二、第三层图像显著图；

(4c)对参考图像灰度特征金字塔模型A_I中三层显著图求平均，得到参考图像r的灰度特征显著图S_rI；

(4d)对参考图像r邻域均差特征金字塔模型A_M采用(4a)-(4c)相同操作，得到参考图像r的邻域均差特征显著图S_rM；

(4e)对参考图像r的灰度特征显著图S_rI和邻域均差特征显著图S_rM求和，得到参考图像r的最终显著图S_r；

(4f)对浮动图像f的灰度特征金字塔模型B_I和邻域均差特征金字塔模型B_M采用(4a)-(4e)相同操作，得到浮动图像f的最终显著图S_f。

5.根据权利要求1所述方法，其中步骤(6)中提取参考图像r和浮动图像f感兴趣区域轮廓点集p_r和p_f，按如下步骤进行：

(6a)统计参考图像最终显著图的分割结果b_r中各非连通区域的面积，根据面积大小对各连通区域按降序排序；

(6b)按(6a)中排序顺序，从第二个非连通区域开始计算各区域与其前一区域面积的比值，并将该比值与设定的面积比阈值T₂＝0.25进行比较，若该比值大于T₂，则标记该区域为感兴趣区域，并继续判断下一区域是否为感兴趣区域，否则，提取所有在迭代过程中标记为感兴趣区域的区域轮廓点，作为参考图像r的感兴趣区域的轮廓点集p_r；

(6c)按(6a)-(6b)的相同操作，提取浮动图像f的感兴趣区域的轮廓点集p_f。