文丽丽 罗洪艳* 张绍祥 郑小林 吴 毅 李 敏

1(生物流变科学与技术教育部重点实验室(重庆大学)，重庆大学生物工程学院，重庆 400044)

2(第三军医大学基础医学部解剖学教研室，重庆市计算医学研究所，重庆 400038)

基于运动跟踪获取人体呼吸曲线的初步研究

文丽丽1罗洪艳1*张绍祥2郑小林1吴 毅2李 敏1

1(生物流变科学与技术教育部重点实验室(重庆大学)，重庆大学生物工程学院，重庆 400044)

2(第三军医大学基础医学部解剖学教研室，重庆市计算医学研究所，重庆 400038)

4D-CT的实现需要在采集CT图像的同时，同步获取患者的呼吸曲线，再据此对CT图像进行排序和三维重建。本研究提出一种基于运动跟踪获取人体呼吸曲线的方法。即在人体腹部贴上标记物，用两台摄像机在CT扫描过程中成正交位同步跟踪呼吸运动，将得到的视频提取为单帧图像，利用区域生长算法提取标记物的质心，并跟踪质心的运动，描绘坐标的任意维度与时间的关系曲线，即得到的呼吸曲线。通过此方法获得了正常人在深呼吸时的呼吸曲线，并在该曲线上确定了屏气模式下7个时相点的位置，腹部表面在z轴方向上的运动幅度最大，约为3.5 cm。实验结果表明该方法切实可行，并且相比于其他的呼吸信号提取方法，具有简单易行的优点。

呼吸曲线;运动跟踪;4D-CT;呼吸运动;放疗

引言

基于CT、MRI的图像引导放射治疗和调强放射治疗基本解决了静止靶区的精确放疗问题［1］，但对于运动的靶区，却达不到很好的解决效果［2］。这就需要将时间纳入三维适形放疗，产生了4D-CT的概念［3］。将4D-CT应用于放疗的一种思路是采集呼吸功能正常的健康人的4D-CT数据，据此建立相应的有限个时相点的胸腔模型，再结合呼吸运动规律，采用空间配准或插值等方法获得一个呼吸周期中任意时相点的模型，即构建出胸腔动态仿真模型作为标准模型，然后运用图像融合等方法，将临床获取的病人的CT图像中反映出的个体特征融入到该标准模型中，从而获得个性化的胸腔动态仿真模型，以模拟、预测呼吸运动对患者胸腔内放疗靶区的影响。因此，整个工作包括4D-CT数据获取、标准胸腔动态仿真模型构建、个性化胸腔动态仿真模型构建、临床应用探讨等4个方面。其中，4D-CT数据获取是前提和基础，也是一个难点问题，其精度决定了动态模型的质量。它的基本思路［4］是同步采集受试者在自由呼吸或深呼吸时的胸腔CT图像和呼吸曲线，再将两者关联，确定所获取的CT图像序列所对应的呼吸时相。

整个4D-CT数据的获取包含两个关键环节，一是CT图像的获取，目前普遍使用的是电影扫描和螺旋扫描两种扫描方式。电影扫描［3］分床位扫描，每个床位至少扫描一个呼吸周期，所获取的CT图像在同一床位内的时相关系是确定的，但在相邻床位间的时相关系存在不确定性，需将扫描图像重新排序后拼接［5］。螺旋扫描是对整个目标区域连续扫描多个呼吸周期，所获取的序列CT图像有确定的时间关系，无需配准［6］;二是呼吸曲线的获取。现有方法大多是通过测量记录呼吸过程中产生的一些物理量的变化提取呼吸曲线，比如肺量测定法、压力差测定法等［7－8］。还可以通过测量体表随呼吸起伏的高度差［9］来反映呼吸运动，较为成熟的是Varian公司已经商业化的实时位置跟踪系统(RPM)。这些方法存在各自的利弊，改进的空间很大。

本研究针对现有呼吸信号获取方法存在的问题与不足，提出一种基于运动跟踪获取人体呼吸曲线，从呼吸周期中抽取7个具有代表性的时相点，采用螺旋扫描方式获取各时相点的胸腔CT数据以获取4D-CT数据集的方法，旨在提高4D-CT模型精确度的同时，简化其数据获取，以便于医院临床实践，推动4D-CT相关研究在我国的开展。

1 材料与方法

采用双摄像机跟踪记录受试者腹部标记物随呼吸的运动，并以螺旋扫描方式获取呼吸周期中具有代表性的7个时相点(包括吸气早期、吸气中期、吸气末、呼气早期、呼气中期、呼气晚期、呼气末)的CT图像，根据获得的呼吸曲线确定各时相点在呼吸曲线上的具体位置，从而实现了CT图像序列与呼吸曲线的关联，为后续构建标准胸腔动态仿真模型打下基础。整个方法的流程如图1所示，算法用Matlab编写实现。

图1 基于运动跟踪获取人体呼吸曲线方法流程图Fig.1 Flow chart of acquiring the respiration curve based on motion tracking

1.1 实验准备

实验场景如图2所示。受试者为20岁呼吸功能正常的健康男性，取仰卧位，双手上举抱头，在其体侧和脚部位置分别放置1台摄像机。在受试者体表的肋骨尖端、乳头纵线、前正中线贴有浸染显影剂的棉签(能够在CT图像上清晰显现，用于后续动态仿真模型的构建)，同时在其胸腹部随呼吸运动起伏明显处贴有5个边长为2 cm的正方形黑色标记物，后期将通过跟踪该标记物的运动来获得呼吸曲线。此外，分别在 CT床位正前方、侧边、床面贴校准参照图若干。

图2 实验场景图Fig.2 Experimental scene picture

1.2 获取视频数据和CT图像

采用GE Lightspeed VCT 64层螺旋CT获得CT数据，参数设置为管电压120 kV，管电流50 mAs，扫描层厚为2.5 mm，扫描范围包括整个肺组织及主要肋骨。

考虑到呼吸运动是半自主的运动，存在着重复性差［10］的特点，而对受试者的呼吸进行言语指导能够有效地提高呼吸的规律性、可重复性和稳定性［11］。因此在本实验中，先对受试者进行了呼吸训练，即让受试者进行连续深呼吸，并根据其呼吸情况制定一个呼吸节奏，通过言语指导让受试者按照该节奏进行多次呼吸，同时调整两个摄像机的位置以获取最佳拍摄效果;待受试者深呼吸平稳后，摄像机成正交位拍摄约2 min以提取深呼吸曲线。随后，训练受试者在深呼吸过程中达到指定的7个呼吸时相时分别屏气。待受试者已熟悉整个过程后，正式采集屏气7个时相的视频信号和CT数据。根据受试者的手势确定各时相的采集起始点，触发CT扫描，并同步视频拍摄该过程。

1.3 提取呼吸曲线及确定时相点位置

在标记物平面建立世界坐标系，x轴与CT床的宽度方向一致，y轴为标记物的长度方向，即与受试者由脚到头的方向一致，z轴则垂直于标记物平面。位于受试者脚部位置的摄像机拍摄的视频能够反映标记物在x、z轴的运动，而位于其体侧的摄像机拍的视频信号反映标记物在y、z轴的运动。根据该坐标系，对采集获得的视频信号进行处理以提取出呼吸曲线，整个过程如图3所示，包括序列单帧图像截取、标记物定位及质心提取、标记物运动跟踪、绘制呼吸曲线等4个步骤。

1.3.1 序列单帧图像截取

将用摄像机拍摄得到的深呼吸模式下和屏气模式7个时相点的视频信号导入视频播放软件(如KMPLAYER)，每隔0.2 s捕获一帧图像，将视频信号截取为序列单帧图像。

1.3.2 标记物定位及质心提取

这一过程包括确定标记物的位置、区域生长和质心提取等3个步骤。首先在单帧图像中确定标记物的位置，截取包含标记物的一个矩形子图，见图4(a)。然后对该子图采用区域生长算法去除背景，并转换为二值图像，从中提取出5个标记物的轮廓，见图4(b)。区域生长算法的基本思路是将具有相似性质的像素集合起来构成区域［12］，具体处理方法是首先选取一组种子点，并预先定义一个生长准则，然后将类似于种子点的邻域像素附加到每个种子点，最后形成区域。种子点选择灰度值为30的像素点，生长准则定义为邻域像素与种子点的灰度差小于40。该阈值通过分析矩形子图的灰度直方图确定。然后再根据如下公式确定标记物的质心坐标(xc，yc)［13］:

式中，m00、m10与m01为图像f(i，j)的零阶矩和1阶矩，可根据图像矩的定义(见式(2))计算。

式中，p、q为矩的阶数，f(i，j)为M ×N的图像。

质心提取结果如图4(c)所示，将跟踪标记物的运动转化为跟踪标记物质心的运动，以提高跟踪的准确性。

图4 标记物定位及质心提取结果图。(a)原图(黑色方框为标记物);(b)去背景后的标记物轮廓图;(c)标记物的质心图(黑色圆点为质心位置)Fig.4 Results of locating markers and extracting centroids.(a)Originalimage (black rectangle denotes markers);(b)Contour image of the markers after removing background;(c)Image of the markers’centroids(marked with black dots)

1.3.3 标记物的运动跟踪

分别选取5个标记物在呼气末时的质心位置为参考点，得到标记物在图像中相对于参考点的运动变化轨迹。为了进一步获得标记物随呼吸变化的实际运动轨迹，还需要根据图像在x轴方向上的单位距离与实际空间距离间的比例尺，在y轴方向上的拉伸比例，对图像进行标定，得到标记物相对于参考点的运动位移。图5是从放置在脚部位置的摄像机拍摄的视频中截取的单帧图像。以最靠近肚脐处的标记物为例，首先在床的边沿找到与该标记物的质心Q处在同一垂直平面上的E、F点，得到这两点在图像中的距离像素，已知CT床的实际宽度EF=58 cm，可计算出该方向的比例尺rx=而在CT机洞口处放置的校准参照图中A点与B点、C点与D点间的实际距离为AB=10.9 cm、CD=9.45 cm，它们在视频图像中的距离分别为像素像素，由此可计算出y轴方向的拉伸比例标定后得到该标记物此时沿x轴和 z轴的呼吸运动幅度 (xQ，zQ)。再跟踪从体侧位拍摄的视频图像，就可得到该标记物在y轴的运动幅度yQ。按上述方法分别对深呼吸模式时两台摄像机拍摄的视频信号中截取的序列单帧图像中的5个标记物进行处理即可获得其实际的运动轨迹。

图5 放置在脚部处的摄像机拍摄的视频中截取的单帧图像Fig.5 A frame of the video sequence acquired from the camera placed near the feet

1.3.4 绘制呼吸曲线及确定时相点的位置

分别以质心坐标的 x、y、z坐标为纵轴，图像采集的时刻为横轴，描绘出多个呼吸周期的连续呼吸曲线，再进行叠加平均，获得一个周期的平均呼吸曲线，并在 Matlab中对该曲线进行拟合，以得到反映呼吸曲线变化规律的数学表达式。对与CT图像同步获取的屏气7个时相点的视频信号进行上述处理，分别得到标记点在7个时相点相对于同一参考点的位移信息，据此在平均呼吸曲线上确定这7个时相点的位置。

2 结果

利用上述方法跟踪标记物的运动，发现标记物在x、y轴方向上运动的幅度远小于在z轴的运动幅度，例如图5中紧靠肚脐的标记物在x轴方向相对于原点的最大运动位移为0.266 cm，而在 z轴方向上的最大运动位移为3.9078 cm，且贴在腹部的5个标记物中，紧靠肚脐的标记物运动幅度最大，因此本研究重点关注该标记物在z轴方向的运动。

从受试者深呼吸模式下的视频数据中获取了12个周期的连续呼吸信号，一个周期的平均呼吸曲线及相应的拟合曲线分别如图6和图7所示。图中虚线为实测的呼吸曲线，实线为拟合曲线，其表达式为

式中各参数的取值如表1所示，对连续呼吸曲线及平均呼吸曲线的拟合度分别达到0.973 0和0.999 9。而在平均呼吸曲线上确定出的屏气7个时相点的位置如图7所示。

表1 连续呼吸曲线和平均呼吸曲线拟合公式中各参数的取值Tab.1 Parameter values in the curve－ fitting formula of the continuous and average respiratory curve

3 讨论

图6 z轴方向的连续呼吸曲线Fig.6 The continuous respiratory curve in z axis

图7 z轴方向的平均呼吸曲线及屏气7个时相点的位置Fig.7 The average respiratory curve in z axis and the locations ofseven discrete breathing phases

从实验结果可以看出，采用所提出的方法提取的呼吸曲线类似于正弦曲线，符合Lujan等提出的呼吸运动符合高阶余弦函数的理论模型［14］。该曲线反映了受试者的呼吸运动规律和呼吸的不可重现性。一个深呼吸模式下的呼吸周期约6 s左右，腹部标记物在垂直于腹部平面的方向上的运动幅度在3.5 cm左右，与文献报道［15］相符，说明了该方法的有效性。因呼吸引起的胸腹部器官的运动是一种三维运动，本实验中采用双摄像头获取视频数据，可同时得到标记物在3个方向上的位置变化信息，从而对胸腔的运动情况进行更为全面的分析。分析结果表明在z轴方向上的运动幅度最大，因此在以后获取呼吸曲线的实验中，可以采用单摄像头，仅对z轴方向上的运动情况进行分析获取呼吸曲线。同时通过跟踪标记物的运动发现呼吸运动，引起器官的运动主要是在垂直于腹部平面的方向上，而在其它方向的运动幅度很小。屏气7个时相点在平均呼吸曲线上均找到了对应的位置，据此得到基于CT图像序列构建的模型分别所处的时间点。但因呼吸的幅值和周期在任意两周期间存在差异，使得吸气最大时的幅值与平均呼吸曲线的最大幅值存在0.04 cm的偏差，且在训练受试者到达指定的呼吸时相屏气的过程中，呼吸的深度由受试者自主确定，以致这7个时相点在呼吸曲线上的分布不够均匀，因此在保证呼吸的平稳性及控制时相点的精确性方面还有待改进。

现有的获取呼吸曲线的方法，会给病人带来诸多不适。如以呼吸量的周期性变化来表征呼吸运动变化规律的肺量测定法［7］，是用数字肺活量计测量进出入肺的气流，在一定程度上反映了呼吸运动的规律。但重建的4D-CT质量与肺活量计的测量精度相关并存在基线漂移。通过测量呼吸的压力差来反映呼吸运动规律的压力测定法［8］，是将一个圆柱形弹性气囊用皮带捆扎在人体腹部，吸气时，腹部扩张，气囊内的气压增加，呼气时与上述过程相反。该系统可视为一个无漂移的系统。但是它会受到压力传感器的摆放位置的影响。高度差测定法通过测量胸腹部的高度差变化来表征呼吸运动，较为成熟的RPM系统是在人体腹部放置一个塑料块，体外通过红外跟踪此塑料块上两个红外反射标记物的运动，国外大多有关4D-CT的研究及呼吸门控均基于此系统，较为方便，但该系统需要在CT机上安装经过标定的专用接口及配套软硬件才能使用，且价格较昂贵。而且因整体封装，其CT图像的获取和门控治疗的过程都完全是自动的，不利于改进和控制。同时塑料块的压迫可能导致受试者或患者呼吸的变化，尤其对于腹部肿瘤、腹部的障碍物还会影响CT成像和治疗［16］。另据报道该系统在一些情况下识别吸气末时相的准确性较低［17］。鉴于上述原因，呼吸门控在临床上应用较少，因此国内绝大部分的CT机未配备该系统。与这些提取呼吸信号的方法相比，本研究的方法具有简单易行、经济适用、便于移植，针对性强等优点，且不需要对受试者施加会影响呼吸运动的检测物，克服了RPM系统在这方面的不足。此外，标记物紧贴于腹部表面，能很好的跟随腹部器官的运动，反映呼吸运动器官的运动规律。

根据CT图像采集所涉及的呼吸相位的不同，有前瞻性与回顾性4D-CT之分［17］。在现有的采用螺旋扫描方式进行的回顾性4D-CT的研究中，通常是根据呼吸曲线划分时相，对CT图像进行分组排序，因此呼吸信号在一个周期内划分的时相数越多，4D-CT体现器官随呼吸运动的精确度越高，但每个时相所占有的 CT图像越少，模型的空间分辨率降低。因此大多存在着提高4D-CT精确度和增加三维模型空间分辨率难以平衡的问题［8］。本研究结合回顾性重建和前瞻性重建的优点，选取7个具有代表性的时相点，让受试者在指定的时相点处于屏气状态下采集 CT图像，与自由呼吸状态下采集的CT图像相比，因呼吸引起的运动伪影和失真会大大减少，且每个时相点均完整扫描整个待重建区域，有利于提高所构建4D-CT的空间分辨率。但该方法对一些患者(如呼吸功能障碍患者)可能不适用，且受试者在屏气状态和自由呼吸状态下的胸腔活动也可能存在差异，因此尚需这方面的对比分析研究。

4 结论

采用摄像机跟踪记录人体腹部标记物在呼吸过程中位置变化、再结合图像处理的方法，成功获取了健康人的呼吸曲线，并确定了在屏气模式下采集的CT图像序列分别对应的时相点位置，为用于构建胸腔动态仿真模型的4D-CT数据集获取、呼吸信号提取及实现其与CT图像序列的关联，提供了一种经济、方便、实用的方法。在后续研究中将考虑采用一定的视听装置辅助测试，以获得更为平稳的呼吸曲线，实现CT数据多个采集时相点在呼吸曲线上的均匀分布，并对屏气模式与自由呼吸模式下获取的4D-CT数据进行对比分析研究。

［1］ 张书旭，周凌宏，陈光杰，等.4D-CT重建及其在肺癌放疗中的应用研究进展［J］.中国辐射卫生，2008，17(3):375－377.

［2］ 张书旭，陈光杰，周凌宏，等.运动靶区CT三维重建对其体积与形状的影响［J］.华南国防医学杂志，2007，21(5):14－17.

［3］ Lu Wei，Parikh PJ，Hubenschmidt JP，et al.A comparison between amplitude sorting and phase-angle sorting using external respiratory measurement for 4D CT ［J］.Medical Physics，2006，33(8):2964－2974.

［4］ Wink NM，Panknin C，Solberg TD，et al.Phase versus amplitude sorting of 4D-CT data［J］.Journal of Applied Clinical Medical Physics，2006，7(1):77－85.

［5］ 张书旭，周凌宏，徐海荣，等.基于相邻图像最相似原理的4D-CT图像重建研究［J］.中国医学物理学杂志，2009，26(5):1409－1414.

［6］ 杨克柽.基于单螺旋CT的四维成像［J］.实用放射学杂志，2008，24(9):1269－1272.

［7］ Low DA，Nystrom M，Kalinin E，et al.A method for the reconstruction of four－dimensional synchronized CT scans acquired during free breathing［J］.Medical Physics，2003，30(6):1254－1263.

［8］ Vedam SS，Keall PJ，Kini VR， et al. Acquiring a fourdimensional computed tomography dataset using an external respiratory signal［J］.Physical in Medicine and Biology，2003，48(1):45－62.

［9］ Keall PJ，Starkschall G，Shukla H，et al.Acquiring 4D thoracic CT scans using a multislice helical method［J］.Physical in Medicine and Biology，2004，49(10):2053－2067.

［10］ Keall PJ，Kini VR，Vedam SS，et al.Potential radiotherapy improvements with respiratory gating［J］.Australas Phys Eng Sci Med，2002，25(1):1－6.

［11］ Mageras GS，Yorke E， Rosenzweig K， et al.Fluoroscopic evaluation of diaphragmatic motion reduction with a respiratory gated radiotherapy system ［J］. Applied ClinicalMedical Physics，2001，2(4):191－200.

［12］ 罗洪艳，李敏，谭立文，等.一种数字人脑部切片图像分割新方法［J］.中国医学影像技术，2009，25(8):1488－1491.

［13］ 王冰，职秦川，张仲选，等.灰度图像质心快速算法［J］.计算机辅助设计与图形学学报，2004，16(10):1360－1365.

［14］ Lujan AE， Larsen EW， Balter JM， et al. A method for incorporating organ due to breathing into 3D dose calculations［J］.Medical Physics，1999，26(5)，715 －720.

［15］ Shi Chengyu，Papanikolaou N.Tracking versus gating in the treatment of moving targets［J］.European Oncological Disease，2007，1(1):83－87.

［16］ 宁徐彷，姚毅，全红，等.双摄像机获取人体呼吸信号的探讨［J］.武汉大学学报，2007，28(3):319－ 321.

［17］ Pan T.Comparison of helical and cine acquisitions for 4D-CT imaging with multislice CT ［J］.Medical Physics，2005，32(2):627－634.

A Preliminary Study of Acquiring the Respiration Curve Based on Motion Tracking

WEN Li-Li1LUO Hong-Yan1*ZHANG Shao-Xiang2ZHENG Xiao-Lin1WU Yi2LI Min1
1(Key Laboratory of Biorheological Science and Technology(Chongqing University)，Ministry of Education，Bioengineering College，Chongqing University，Chongqing 400044，China)
2(Department of Anatomy，College of Medicine，The Third Military Medical University，Lab of Compute Medical Science of Chongqing，Chongqing 400038，China)

The realization of 4D-CT typically requires the synchronous recording of the respiratory curve of the subject while CT scanning，and then sort and three-dimensional reconstruction of CT images accordingly.This paper presents a method based on motion tracking for human respiratory curve.Initially，a few markers were attached to the abdomen of the subject.Two cameras were placed orthogonally and used to monitor the respiratory motion synchronously in the CT scanning process.Then the single frames were extracted from the videos.The centroid of each marker in these images was determined using region growing algorithm and its movement was tracked.Finally，the arbitrary coordinate dimension of the centroid vs time curve，i.e.the respiration curve，was depicted.Through this method，the respiration curve of a normal person in deep breathing was acquired，and the locations of seven discrete breathing phases in breath－holding mode were determined.The motion of abdominal surface in the z axis direction has the maximum amplitude of 3.5 cm.The experimental results indicated that this method was feasible and easy to realize，compared with some other respiratory signal extraction methods.

respiration curve;motion tracking;4D-CT;respiratory movement;radiotherapy

R318

A

0258-8021(2012)04-0512-06

10.3969/j.issn.0258-8021.2012.04.00

2012-01-04，录用日期:2012-05-29

国家自然科学基金(60871099);中央高校基本科研业务费(CDJXS10231122)

*通信作者。E-mail:cqu_lhy@163.com