基于双目视觉的呼吸运动实时跟踪方法研究

朱超华 陈武凡 徐子海 陈超敏 贺志强*（南方医科大学生物医学工程学院，广州 5055）

2（解放军303医院放射治疗中心，南宁 530021）

引言

对肿瘤靶区的精确性定位，是体外分次放疗的关键问题。使用立体定向技术，可以为头颈部的肿瘤靶区制订出高精度的治疗计划［1］。但是，在呼吸等生理因素的影响下，位于胸腹部的肿瘤及其周围组织动态地变化。相关研究表明，胸隔膜在呼吸运动的影响下活动范围可达0.5～2.0cm，并且个体的差异性很大［2-3］。为了减少呼吸运动造成的不利后果，传统解决方案一般采用等中心位移和呼吸门控等手段［4］，这些方法对患者的呼吸和状态要求严格，治疗效率低，正常组织仍然受到照射伤害。采取扩大边界的方式往往会危及到正常器官，因此外扩靶区的治疗手段极不安全。

实时跟踪肿瘤靶区的关键就是计算出肿瘤靶区在空间中随时间变化的具体位置。有研究者采用平板kV级X射线成像系统［5］、红外线跟踪体外标记加正交X线成像（Cyberknife-赛博刀）定位系统［6］和4DCT成像系统［7］等引导放疗，在一定程度上解决了肿瘤靶区精确定位的问题。但是，成人的正常呼吸周期为4～5s，由呼吸运动导致的肿瘤靶区运动时刻都在发生，而一般的X线机的系统成像周期为1min［6］，加上图像配准及电机驱动等延时，因此并不能做到准确实时，并且不停的X线扫描具有一定的副作用。基于4D-CT成像的基本原理，它需要病人严格配合并进行呼吸训练，复杂性较高，由于不是在治疗的同时获得图像，因此也不是真实意义上的实时肿瘤靶区跟踪。研究和探索一种能够实时跟踪靶区、提高治疗精度和效果且给患者伤害最小的方法，仍然是众多科研工作者的共同努力方向。

下面介绍一种基于图像尺度不变特征变换（SIFT）的双目视觉实时跟踪腹部表面运动的测量方法，对动态自适应中的跟踪模型、时间响应和跟踪精度等几个关键技术问题进行了研究和实验。大家都知道，在双目视觉的三维坐标计算过程中，准确匹配出左右视觉图片中的同一物体是计算精度最重要的保证。由于图像的SIFT特征对尺度、旋转、亮度、仿射、噪声等都具有良好的不变性，在混乱和遮挡情况下也可鲁棒地识别目标，所以本系统采用图像的SIFT特征进行配准。Mikolajczyk K比较了各种局部特征检测和描述符算法，指出SIFT特征算法具有最好的匹配效果［8］，但是计算时间长，一般情况下难以满足大尺寸图片配准的实时性要求。因此，本系统采用该算法进行配准时，在配准图像的选取过程中采用动态采集左右视频图像感兴趣区域的方法，减少了提取整张视频图片SIFT特征所需要的大量时间，较好地满足了双目视觉跟踪的实时性要求。

1 方法

双目视觉测量跟踪的基本原理是：通过已经标定好的两个摄像头，对以不同角度同时拍摄到的两幅平面图像进行校准，再利用已知被跟踪目标图像的SIFT特征，对所得到的两幅图像进行配准，计算出被跟踪物体分别在两幅图像质心的坐标。然后，依据双目视觉成像的基本模型，计算出物体在三维空间中的实际位置。通过引入时间变量，可以得到物体在空间的实时运动情况而进行运动跟踪。

1.1 跟踪目标图像SIFT特征的配准

视频跟踪的效果取决于视频图像的可靠性和图像配准的速度。基于特征的匹配稳定可靠，是较适合于视频图像中对特定目标配准的方法。DG Lowe教授总结了现有的、基于不变量技术的特征检测方法，提出了一种基于尺度空间的、对图像缩放和旋转甚至仿射变换保持不变的图像局部特征描述算子—SIF算子，这是一种多尺度技术［9］。SIFT特征匹配算法思想是：首先在尺度空间进行特征检测，确定关键点（key points）的位置，然后使用关键点邻域梯度的主方向作为该点的方向特征，以实现SIFT特征对尺度和方向的无关性。SIFT特征匹配算法分为生成SIFT特征向量和特征向量匹配两个过程。

1.1.1 生成SIFT特征向量

生成SIFT特征向量包括以下4步：检测尺度空间极值点、精确定位极值点、为每个关键点指定方向参数和生成关键点的特征向量。

尺度空间理论是模拟图像数据的多尺度特征，尺度空间极值检测的主要过程是：在尺度空间内利用唯一的线性核——高斯核，建立高斯金字塔［10］。为了有效地在尺度空间检测到稳定的关键点，进一步提出了高斯差分尺度空间（DOG scale-space），利用不同尺度的高斯差分核与图像卷积生成，有

在DOG金字塔内进行极值检测.最后可以初步确定特征点的位置及所在尺度。为了寻找尺度空间的极值点，每一个采样点要和它所有的相邻点比较，看其是否比它的图像域和尺度域的相邻点大或小。

精确定位极值点的作用是消除低对比度极值点和不稳定的边缘响应点，获得图像局部特征点。通过拟和三维二次函数，以精确地确定关键点的位置和尺度（达到亚像素精度）；同时删掉低对比度的关键点和不稳定的边缘响应点（因为DOG算子会产生较强的边缘响应），以增强匹配稳定性、提高抗噪声能力。

对于选定的图像局部特征点，利用特征点邻域像素的梯度方向分布特性，为每个特征点指定其方向参数，使得SIFT算子具有旋转不变性。

下式为（x，y）处梯度的模和方向公式，有式中，L所用的尺度为每个关键点各自所在的尺度。

以特征点为中心取8×8的窗口，在每4×4的小块上计算8个方向的梯度方向累加值，绘制每个梯度方向的直方图，即可形成一个种子点，每个种子点有8个方向的向量信息。Lowe建议对每个特征点使用4×4共16个种子点来描述，这样对于一个特征点就可以产生128个数据，最终形成128维的SIFT特征向量。此时，SIFT特征向量已经去除了尺度变化、旋转等几何变形因素的影响，再继续将特征向量的长度归一化，则可以进一步去除光照变化的影响。

1.1.2 SIFT特征向量的匹配

当两幅图像的SIFT特征向量生成后，进一步采用关键点特征向量的欧式距离作为两幅图像中关键点的相似性判定度量。取图像1中的某个关键点，并找出其与图像2中欧式距离最近的前两个关键点。Lowe采用的方法是使用欧氏距离作为特征点的相似性度量，特征点a、b间的欧氏距离Uab表示为式中，n为特征向量的维数。

为了排除因为图像遮挡和背景混乱而产生的无匹配关系的特征点，Lowe提出通过比较最邻近距离和次邻近距离来消除错配，有

式中，Umin为最邻近距离，U1为次邻近距离，当它们的比值小于距离比例阈值R时判定为正确匹配，否则为错误匹配。降低比例阈值R，会使SIFT匹配点数目减少，但更加稳定。

1.1.3 SIFT算法的优化

虽然SIFT特征匹配算法具有高度的鲁棒性，但是需要大量的计算时间。为了提高特征可靠性和配准速度，在具体情况下，对SIFT算法做以下优化。

1）SIFT算法的计算量主要耗费在建立尺度空间和生成特征向量，图像越大，生成的特征向量越多。为此，在第一次完成标记物图像和左右视频图像的匹配后，在后续的左右视频图像中，以标记物所在的上一帧图像位置为中心，动态选取匹配区域来进行下一步的配准，然后再把所获得的准确位置映射到原图像中，所选动态区域的大小决定了视频匹配的效率。

2）由于所需配准的标记物的图片并不在复杂环境中，因此需在图像配准过程中，进一步减小最邻近和次邻近距离的比值R，使配准更加稳定、运行速度更快。

1.2 双摄像机测量原理

物体P在双摄像头像平面空间中的成像关系如图2所示。

依据上述分析，得到物体P在左像平面上的坐标（xPl，yPl，zPl），可以写为

同理，物体P在右像平面上的坐标（xPr，yPr，zPr）可以写为

图1 平行双目视觉原理Fig.1 Parallel binocular stereo vision schematic diagram

图2 立体坐标Fig.2 Stereo coordinate system

式中，（xl，yl）和（xr，yr）分别为物体在两个像平面上的坐标，可以分别由在两副图像中的像素坐标、通过已知的摄像机内部参数计算得到（内部参数通过双摄像头标定得到）。

依据式（5）和式（6），可以通过获取目标点在两幅图像中的位置，在摄像机坐标系中标定目标点的三维坐标［11］。

1.3 双摄像机立体标定

摄像机标定是通过所摄取的图像，求解摄像机的内外参数。摄像机的外部参数表示摄像机的位置和方位相对于一个实际坐标系的坐标变换，内部参数表示摄像机的光学本质特性。在本系统中，对摄像机的标定方案是先分别对两个摄像机进行单个标定，再进行立体标定。单个摄像机的标定是分别求出它们的内参数矩阵、畸变系数、旋转矩阵R和平移向量T，前两者构成了摄像机的内参数，R和T构成了物体位置和方向的摄像机外参数。在立体标定中，利用旋转矩阵R和平移向量T来联系左右摄像机。对于给定物体坐标系中的任意3D点P，分别代入到左右摄像机的坐标系中，有

相应地，P点在两台摄像机的关联关系为

式中，R和T分别为两个摄像机之间的旋转矩阵和平移向量。

联立等式（7）和式（8），可以得到

即得到了两台摄像机之间的内参数，可以用来帮助计算物体在空间中的位置。

2 实验设计

2.1 双摄像机跟踪系统的设计

由上述分析可知，可以计算静态状况下选定标记物的三维坐标。如果要进行三维视觉跟踪，则需要对时间序列的视频图像进行计算。为了实现上述功能，将整个工程细分成几个模块实现，包括系统初始化模块、视频图像校准模块、标记物图像SIFT特征提取模块、立体匹配模块、三维坐标计算模块和运动跟踪实现模块，具体流程如图3所示。

图3 标记物跟踪处理流程Fig.3 Markers tracking process

1）初始化模块。该模块的主要功能是对摄像机进行立体标定，先分别获取左右摄像机的内外参数和畸变向量，再进行立体标定，计算得到两摄像机之间的内参数（R，T）。其中，单个摄像机的畸变系数输入到视频图像校准模块，对视频图像进行矫正，两摄像机之间的内参数（R，T）则用来计算三维坐标。具体方法是：两台摄像机同时获取标定棋盘板的系列图片（共15副7×8的棋盘板图片），对于每张图片，提取棋盘角点、内参数和畸变系数的计算。

2）视频图像校准模块。用标定的内参数矩阵，对获取的每一帧视频图像进行校准。具体过程是：先将原图像乘以畸变系数，得到非畸变的图像，再对其校正和裁剪成两幅图像相重叠区域的图片。

3）标记物图像SIFT特征提取模块和立体匹配模块。获取标记物图像和左右摄像机图像匹配点。具体方法如下［12］：首先分别建立标记物图像Image_object的SIFT特征点集合Sb，左摄像机图像Image_left的SIFT特征点集合Sl和右摄像机图像Image_right的SIFT特征点集合Sr。用集合Sb对集合Sl进行搜索，提取具有相同特征的元素，建立集合Sbl；同理，建立集合Sbr。

下面计算配准标记物的坐标。设空间匹配点集中元素总数为n，（xi，yi）为点集中的坐标，则该点集的中心坐标为

将中心坐标取整后作为标定点图像坐标，分别记为（xl，yl）和（xr，yr）。

4）三维坐标计算模块。依据本文1.2节所述的原理，结合已经求得的（xl，yl）和（xr，yr），可以计算出物体三维坐标。

5）三维运动计算模块。根据物理学中速度的定义，要实现对标记物的实际速度测量，需要获得两个信息量，一个是标记物的运动时间间隔Δt，另一个是在这段时间内标记物的运动路程S。对于第k帧到第l帧的时间间隔tkl，可以通过计算目标运动期间所经历的视频序列帧数和采集频率得到，即tkl=|l-k|Δt。对于标记物的云动路程Skl，可以用空间欧氏几何知识来计算，有

式中，（xk，yk，zk）为第k帧目标质心的三维坐标。

标记物的平均运动速度为：

标记物在x、y、z轴方向的运动速度，则可以计算为

6）后续帧的配准图像选择。为了减少计算量、提高处理速度，在提取第一帧视频图像的SIFT特征后，其后续的图像的配准采用以下方式：以前一帧图片中计算出的标记物的坐标（xl，yl）为中心，选定一个配准区域图像imge_block进行匹配，依据人体腹部随呼吸运动的状况和所采用的视频采集频率，选定配准区域为

式中，（vx，vy）为上一帧所计算出来的标记物的运动速度，t为视频的采样周期，（x，y）为上一帧中标记物的坐标，w为标记物图像的宽；系数1.5是为了提高匹配的准确性，如果过大会增大所选图像的配准区域，降低效率，过小则有可能丢失标记物图像，导致匹配失败。

如此循环执行，如果标记物丢失，则再到第3步重新执行。

2.2 实验实施方法

实验的硬件由2个平行并排置于治疗床上方的摄像头和1台PC组成，PC的配置为4GB RAM，2.66 GHz*4 CPU，Nvidia Geforce 9800 GT显卡。为了方便实验进行，作者基于C++平台开发了一套实验软件系统，用来完成视频图像采集、摄像机标定、标记物的匹配和三维运动计算等工作，最后将标记物的运动结果动态地显示出来（摄像机的采样频率为15Hz）。依据文中第2.1节所述的工作流程，软件系统先采集标定棋盘图片，完成立体标定。如图4中的窗口所示，选定一个条形码作为标记物，将其固定于一个平躺人的腹部。在图4中，条形码实际位置的边界由系统自动实时界定。

3 结果

根据实验记录，前15帧的数据如表1所示，表中记录了每一帧图像中提取的SIFT特征点的个数和时间，以及标记物（条形码）在各帧图像中的二维坐标（x，y），三维坐标（x，y，z）为实时计算出的标记物在三维空间中的具体坐标。其中，左右图像中坐标的单位为标记物在图像中的像素点位置。

4 讨论

4.1 时间与匹配准确性

从表1中可以发现，在提取第一帧左右图像的SIFT特征向量时，在预设的条件下，左右图像的特征个数分别为590和610，耗费时间分别达到了350.56和430.36ms；如果再加上匹配和坐标计算时间，完成第一帧图像的三维坐标的计算将接近1s，远不能满足实时性。但是，在第二帧以后，由于采用了文中第2.1节6）所述的优化方法，在左右图像中提取的SIFT特征数分别为35和30，耗费时间分别降为13.01和14.12ms，加上匹配和坐标计算时间，每一帧的处理时间可以在50ms内完成。在试验中还发现，当标记物随腹部表面运动时会发生旋转，角度变换和小部分遮挡导致左右图像中的特征向量数目（13）小于标记物图像特征向量数（25），但系统依然能够标记出标记物在图像中的位置，说明采用此方法的准确性能够满足试验的设计要求。

图1 腹部标记物跟踪。（a）左摄像头；（b）右摄像头Fig.1 Abdominal markers tracking figure.（a）the picture from the left camera；（b）the picture from the right camera

表1 实验数据Tab.1 Experimental data

Alper Yilmaz指出，基于图像的目标跟踪存在一些复杂因素，主要在于以下几个方面［13］：由2D图像投射到3D图像过程中信息丢失、图像噪声、目标的复杂运动、目标的部分遮挡、复杂的目标形状、光照变化和实时性的要求。在本实验中，先分别匹配出2D图像中目标的坐标位置，避免了图像直接投射到3D场景中再计算三维坐标而导致的信息丢失。由于采用SIFT算法作为配准的方法，能很好地解决图像的光照和仿射变换等问题。上述的实验结果也表明，系统能够解决目标发生部分遮挡、旋转和角度变换的情况。

4.2 测量精度

将各轴向的坐标变化在时间（帧数）轴上表示出来，可得到如图5所示的运动时间曲线。三维坐标中x表示标记物的左右运动方向，如曲线（a）所示；y表示标记物在人体头脚方向运动的情况，如曲线（b）所示；z表示标记物的运动深度方向，即为标记物随腹部上下运动的情况，如曲线（c）所示。

图5 位置-帧数（时间）跟踪曲线。（a）对应于x轴，人体左右方向；（b）对应于y轴，人体的头脚方向；（c）对应于z轴，腹部的上下运动方向Fig.5 The Position-Time curve of the tracking.（a）xcoordinate axis for the left and right direction of the body；（b）y coordinate axis for the head and foot direction of the body；（c）z coordinate axis for the up and down direction of the abdomen

由图5可见，此呼吸周期约为3.5s，标记物的运动范围分别为前后0.34cm、左右0.25cm、上下1.43cm，与实际监测值的上下运动距离1.60cm的差为0.17cm。由于标记物放置于腹部表面，所以位移变化主要集中在z轴方向，符合人体腹部随呼吸运动变化的规律。

计算中的误差主要来源于以下方面：一是摄像机的内部参数和实际的误差，虽然摄像机系统经过立体标定，但是标定板尺寸和两摄像间的固定距离存在一定的误差；二是摄像机自身的像素限制。

4.3 多点跟踪与呼吸运动模型

前面所讨论的是基于单点跟踪呼吸运动的模型，但在实际的放疗过程中，位于胸腹部的不同位置、具有不同的体积大小、处于不同分期的肿瘤有着不同的运动模式，单独考虑体表标记点的呼吸运动情况，将其等效于肿瘤靶区的呼吸运动情况有一定的误差性。因此，在下一步的工作中，要对单点跟踪模型进行改进，拟采用多点跟踪呼吸运动跟踪法，并结合4DCT图像，为个体构建腹部标记物-肿瘤相关性模型。在实际应用过程中，将建立的肿瘤靶区实时运动轨迹输入到反馈跟踪系统中，对治疗过程中的肿瘤靶区运动进行逆向补偿，以期得到更好的胸腹部肿瘤放射治疗效果。

5 结论

本研究采用提取动态范围图像和标记物图像的SIFT特征进行匹配的方法，实现了双目视觉的三维跟踪。试验结果表明，通过上述方法测得的标记物三维运动轨迹精确度高、实时性好，并且具有良好的稳定性。该方法是对腹部表面单点进行运动跟踪，如果要建立整个腹部表面的运动模型，则还需要对腹部多特征的运动进行跟踪分析。同时，该方法仅仅是对位于腹部表面的标记物跟踪，在实际临床应用中，还需要研究胸腹部肿瘤靶区与标记物之间的对应关系，这一点还需要做很多的工作。

［1］Karger C，J kel O，Debus J，et al.Three-dimensional accuracy and interfractional reproducibility of patient fixation and positioning using a stereotactic head mask system［J］.International Journal of Radiation Oncology Biology Physics，2001，49（5）：1493-1504.

［2］Wu Huanmei，Sharp G，Zhao Qingya，et al.Statistical analysis and correlation discovery of tumor respiratory motion［J］.Physics in Medicine and Biology，2007，52：4761-4774.

［3］张丹丹，张红志，韩伟，等.呼吸运动对肺部肿瘤靶区受照剂量分布影响的研究［J］.中华放射肿瘤学杂志，2009，（003）：191-196.

［4］Minohara S，Kanai T，Endo M，et al.Respiratory gated irradiation system for heavy-ion radiotherapy［J］.International Journal of Radiation Oncology Biology Physics.2000，47（4）：1097-1103.

［5］Berbeco RI，Jiang Steve B，Sharp GC，et al.Integrated radiotherapy imaging system（IRIS）：design considerations of tumourtracking with linacgantry-mounted diagnosticx-ray systems with flat-panel detectors［J］.Phys Med Biol，2004，49（2）：243-255.

［6］Chang S，Main W，Martin D，Gibbs I，et al.An analysis of the accuracy of the CyberKnife：a robotic frameless stereotactic radiosurgical system［J］.Neurosurgery，2003，52（1）：140-147.

［7］Weiss E，Wijesooriya K，Dill S，et al.Tumor and normal tissue motion in the thorax during respiration：Analysis of volumetric and positional variations using 4D CT［J］.International Journal of Radiation Oncology，Biology and Physics，2007，67（1）：296-307.

［8］Mikolajczyk K，Leibe B，Schiele B.Local features for object class recognition［J］.Tenth Ieee International Conference on Computer Vision，Vols 1 and 2，Proceedings，2005：1792-1799.

［9］Lowe DG.Distinctive image features from scale-invariant keypoints［J］.International Journal of Computer Vision，2004，60（2）：91-110.

［10］Babaud J，Witkin AP，Baudin M，et al.Uniqueness of the Gaussian Kernel for Scale-Space Filtering［J］.IEEE Transactions on Pattern Analysis and Machine Intelligence，1986，8（1）：26-33.

［11］吴福朝.计算机视觉中的数学方法［M］.北京：科学出版社，2008：94-98.

［12］孟浩，程康.基于SIFT特征点的双目视觉定位［J］.哈尔滨工程大学学报，2009，（006）：649-652.

［13］Yilmaz A，Javed O，Shah M.Object tracking：A survey［J］.Acm Computing Surveys.2006，38（4）：Article No.13.