心脏超声序列图像配准与跟踪技术研究

兰其斌

(福建师范大学协和学院 信息技术系，福建 福州 350007)

0 引言

心脏的运动是其内部各结构功能性运动(如舒缩变形等)和呼吸、血液涌动等运动引起的非功能性运动(如平移旋转等)的复合运动。国际上的超声心动图检测如M型、组织多普勒及全方向M型心动图都存在这样一个问题：由于检测波束或方向线束是空间固定坐标，所以取得的也只能是复合运动信息，而不是反映该结构机能的功能性运动信息，因此不能准确反映心脏结构某个部位的运动和变形，研究人员已经意识到这个问题的存在，但目前，还没有相关方面的研究报道。本文针对这个问题采用弹性配准技术对心脏超声序列图像区域进行跟踪，从而使得动态信息的检测更为准确和可靠。

心脏是一种弹性变形体，因此文中将基于薄板样条的弹性配准方法[1-2]运用于心脏超声序列图像短轴室壁的轮廓跟踪过程，对其图像进行弹性变换，并使用弹性近似的方法来解决标记点定位误差导致的弹性变换模拟问题。超声序列图像的跟踪步骤为：第一步，对第一帧图像进行轮廓提取，确定轮廓上16点的位置，通过3次B样条曲线拟合得到轮廓线(初始轮廓)；第二步，以16点中的6点作为标记点，在下一帧图像中对应坐标位置周围一定范围内查找匹配点作为标记点的对应点，并运用这6组标记对应点对通过平滑薄板样条弹性配准，得到该帧图像中轮廓线上16点的坐标，并将它们作为下一帧图像的初始轮廓；依据第二步算法进行逐帧计算，直至整个序列图像的轮廓跟踪结束。

1 轮廓跟踪算法

1.1 第一帧图像轮廓提取

1.1.1 预处理

1)运用模糊数学理论对图像进行模糊边缘增强 为有效地提取出超声心动图像心脏短轴室壁的轮廓，我们运用模糊数学理论对图像进行模糊边缘增强[3]，结果如图1(b)，可以看到超声心动图像的心脏短轴内、外膜边缘信息突出，而其他部分信息，尤其是噪声，得到了较好的控制，这为心脏短轴图像边缘的提取奠定了基础。

2)边缘锐化并二值化处理 运用微分算子加权运算[4]计算各像素点的值，即将像素某邻域灰度的平均值的一定比例数，减去这个邻域4个方向上相同大小邻域像素的灰度平均值，即：

(1)

就可以得到该像素增强后的灰度值，在此过程中，平均化相当于突出低频部分，减去平均值相当于减去图像中的低频部分，进而尖锐化边界轮廓。实际上需要的是单像素的边缘轮廓，因此在上述增强后的图像的基础上，必须再进行二值化处理。通过实验选取二值化阈值为128，即任一像素其灰度值若大于128，则令其等于255(白点)；若其灰度值小于128，则令其等于0(黑点)，结果如图1(c)。

图1 图像预处理

1.1.2 轮廓提取

预处理后，即可进行初始轮廓提取，步骤如下：

2)人工干预。由于短轴室壁上乳头肌等的存在，使得所得到的轮廓点可能会受到干扰，而初始轮廓提取的准确度直接影响超声序列图像的跟踪效果，因此在此步骤中，医生可以依据经验进行人工干预，手动修正轮廓点的位置，从而最后得到室壁轮廓上的16个边缘点的位置。

3)曲线拟合。由于提取到的边缘点是不连续的，为了得到连续和光滑的边缘曲线，我们采用3次B样条[6]对前面所得到的边缘点进行拟合，即可得到光滑的边缘曲线。

1.2 标记特征点选取

通过前面的一系列步骤，已经获得了第一帧图像的初始轮廓及16个边缘点，接着提取出下一帧图像，对这两幅图像进行弹性配准。由于是针对心脏短轴室的轮廓进行跟踪的，因此配准得到的是相对应的轮廓线。配准的前提是必须先得到两幅图像中的标记特征点对，因此我们就以第一幅图像(第一帧)上已获取的16个边缘点作为标记特征点，接着从第二幅图像(第二帧)中查找对应点，以第一幅图像中标记点的坐标位置为中心，选取大小为24×24的图像子块作为模板，随后在第二幅图像上以图像子块的中心(即第一幅图像中标记点的坐标位置)为搜索中心，以36×36的窗口为搜索范围，以互相关的相似性测度[7]为匹配测度，在第二幅图像上搜索匹配程度最高的子块，最后以互相关的相似性测度最大的图像子块的中心点作为标记点的对应点，如图2所示。

图2 模板及其搜索图

一个正常人的心动周期约为0.75s～1s，即对于PAL制的超声序列图像来说大约20帧构成一个心动周期。正常人的室壁运动从收缩至最小到舒张到最大移动距离一般约为8 mm，一般幅度从最小到最大大约要8帧左右，即前后帧室壁运动的最大距离大约是1 mm。在屏幕上相邻两个象素点之间的间距约为0.3 mm ～ 0.4 mm，那么两帧图像之间室壁运动大概为3个象素点，因此上述第二幅图像选取36×36的窗口为搜索范围，保证了标记点的对应点落在所取区域内。

1.3 平滑薄板样条弹性配准

通过上述方法得到待配准的两幅图像的对应标记点对，运用匹配运算得到边缘点，再运用3次B样条曲线拟合即可得到轮廓线，但由于得到的边缘点存在误差，因此边缘线就与实际的边缘线存在一定差距。综合考虑轮廓的整体变形，将整个轮廓进行6等分，划分为6个区域，并在这6个区域内查找互相关的相似性测度最大的标记点对，这样就可以得到6组标记点对。鉴于无论是哪一种标记特征点对的选取方法，误差都是不可避免的，因此采用一种平滑的处理方法，对常用的弹性配准的插值方法——薄板样条法进行平滑处理，使得在配准时对标记特征点对的精确性和数量的要求都明显降低，这样通过平滑薄板样条弹性插值得到其余的边缘点，最后通过3次B样条曲线拟合得到轮廓线。依据上述步骤进行逐帧的配准，直到序列结束，从而实现轮廓跟踪，结果如图3所示。

图3 心脏超声序列图像短轴室壁轮廓跟踪

2 结语

如图3所示，对6帧图像进行跟踪，得到的室壁内膜轮廓边缘线与人眼的视觉效果大致吻合，表明文中所使用的算法取得了较好的效果，能够反映出心脏的运动情况。如果以轮廓中心为非功能性运动点，依据轮廓跟踪得到的运动信息，由各帧图像上室壁内膜轮廓上各点的坐标即可得到轮廓中心坐标，进而得到轮廓中心的运动信息，即依据非功能性的运动信息。这样我们可以将轮廓运动中的由于呼吸、血液涌动等运动引起的部分非功能性运动(如平移旋转)去除，得到医生所需要的室壁内膜轮廓上各点的功能性运动(如舒缩变形运动)，从而实现从心脏某结构部位的复合运动中分离出其功能性运动信息。

本文所论述的方法为心脏某结构部位的功能性与非功能性运动的分离提供了基础的研究资料，为心脏某结构部位动态信息的准确检测提供一定的参考。

[参考文献]

[1] Rohr K，Stiehl H S，Sprengel R，et al. Point-based elastic registration of medical image data using approximating thin-plate splines[D]. Germany: SpringerVerlag，1996，1131：297-306.

[2] 罗述谦,阎华.基于薄板样条的MRI图像与脑图谱的配准方法[J].中国生物医学工程学报，2004，23(6): 479-485.

[3] 郭桂蓉.模糊模式识别[M].北京:国防科技出版社,1992：25-36.

[4] 宋光德,陈驷辉,赵伯雷,等. B超序列截面图像轮廓的提取[J].天津大学学报，2000，33(3): 382-383.

[5] 王同，余建国，王威琪.提取心脏截面图像上内腔的边缘曲线[J].复旦学报，2002，41(2)：174-176.

[6] 符祥，郭宝龙. 三次均匀B样条插值曲线和曲面的矩阵形式[J]．光电子·激光，2007，18(8)：992-995.

[7] 杨枝灵，王开.Visual C++数字图像获取、处理及实践应用[M].北京：人民邮电出版社，2003：89-165.