1 引言

放射性疗法是治疗恶性肿瘤最有效的方法之一。据统计,欧洲恶性肿瘤的病人数量将在2025年突破450万^[1],其中超过50%的病人可以通过放疗得到治愈^[2-3]。如何提高放疗的精准度是放疗技术一直以来的研究热点,实现精准放疗需要在确保对肿瘤部位提供足够的照射剂量的同时确保照射束的准确性,既要对目标区域提供足够的辐照剂量又要尽可能地减少放疗照射束对周围正常组织及器官的伤害^[4-5]。

当前国内临床的放疗手段多为在术前根据肿瘤的形态特征确定病人的放疗辐照区域并勾勒标记于病人体表、在手术中直接对指定区域进行照射治疗。但人体是一个非刚性结构,不同的组织器官存在非同步的运动状态,并且人体的心跳和呼吸运动也具有一定的运动复杂度,两者都不遵循严格的周期性和线性。这就会导致人体的组织器官在手术过程中会随着患者的心跳和呼吸产生难以预计的形变和位移^[6-10]。实验发现,在人体呼吸运动过程中,肺部会在上下、左右、前后三个方向分别移动(2.4±1.3) mm、(2.4±1.4) mm 和(3.9±2.6) mm^[11],HOISAK团队利用透视法得到病人在静息态下肺下叶肿瘤的上下运动幅度为14.3mm^[12],这就可能导致治疗过程中的放疗射线会照射到治疗部位以外的正常组织和器官,引发器官组织炎症并诱发其他并发症,不利于病人的治疗和康复。呼吸运动还会导致肿瘤溢出目标范围,使得针对肿瘤区域的辐照剂量不足,影响治疗效果。在这种情况下,控制呼吸运动造成的误差以提高放疗的精准程度就成为了一个亟待解决的问题。

目前解决这个问题的方案主要有两种。一种是在放疗过程中对肿瘤进行实时监测和追踪,并根据肿瘤的实时位置对射线进行调整^[13];该方案准确率高,但技术实现难度大,对硬件设备要求高,难以在短时间内实现该方案在医疗场所的普及。另一种方案是在手术过程中对病人的呼吸运动进行量化,同时获取病人的实时呼吸运动信息,以此来指导手术进行^[14]。针对呼吸运动检测问题,当前国内外在临床上主要的技术方案有屏气或腹部压迫、呼吸门控、四维放射性疗法及实时追踪技术等。其中,通过主动或被动屏气技术控制患者的呼吸或压迫腹部以减少肿瘤及周围器官的运动的方法是直接且可行的,但这种方法需要患者自身肺部功能基本正常且接受了一定时间的呼吸训练,对急症或肺部损伤的病人并不适用,因此该方案的耐受性较差^[15]。呼吸门控技术可以调整放射束的辐照时间与病人呼吸周期的特定时间相同步,以此来减少辐照期间肿瘤在窗口中的移动。但线性加速器的选择性打开延长了加速器的工作周期,使得治疗时长成倍增加^[16]。三维立体放疗技术能够获取人体内部结构及器官的实时位置信息,为肿瘤定位、治疗方案的制定提供技术支持。但该技术的实施需要持续的透视图像采集,过长的辐照时间不但会对患者造成二次伤害,而且容易在纵方向上产生伪影。同时术前对患者的平稳呼吸培训也增加了技术实施难度。近几年也有一些新技术的出现,但这些技术也存在其他的缺点和问题^[17]。综上所述,手术中对病人的呼吸运动进行直观而准确的定量分析,对临床上指导精准放疗具有现实意义,通过对患者呼吸运动的实时表征,可以灵活调整射线和治疗床的位置,确保射线和目标肿瘤一直处于稳定的相对位置上。

当前大多数的实时追踪算法仅在追踪刚体时有出色的表现,但人体表面属于非刚性结构且会随着呼吸运动产生不规则的波动及形变。为解决这个问题,实验中采用一个刚性体表标记物作为追踪目标进而实现对体表的实时追踪。实验结果表明,利用追踪体表标记物可以实现对人体体表呼吸运动的表征,并获取实时呼吸运动曲线以及位移和角度变化等参数。通过对比实验测试结果和实际位移值,验证了实验算法的可靠性。

2 实验

2.1 实验场景搭建

实验1采集了实验对象在5个自然呼吸周期下的图像信息,并通过激光位移传感器测得人体体表的实时位移变化量。将体表的实际位移与算法预测的感兴趣区域(ROI)的左下角点帧间像素差进行对比,验证了实验方案的可行性。实验场景如图1所示,实验过程中令实验对象保持水平平躺状态,将一个木制正方体标志物紧密地连接在实验对象体表肚脐与剑突中间的位置(经验证该部位体表随呼吸运动的波动起伏最明显),将一台采集频率为645Hz的激光位移传感器(Micro-Epsilon optoNCDTBL)安装于支架上,其与实验对象体表的距离为50cm,并调整射线使其照射光斑贴近标记物左下角点附近的体表上。同时架设一台距离实验对象50cm的摄像头,使其正对体表标记物中心点并且使二者保持水平,采集整个呼吸过程。同时记录激光位移传感器反馈的体表实时位移及算法追踪得到的标记物左下方角点的像素值变化。

图 1. 实验1场景示意图

Fig. 1. Schematic diagram of experiment 1

下载图片 查看所有图片

为了准确地控制实验变量,实验2将实验对象替换为仿真体模,在中国科学院近代物理研究所深层肿瘤治疗实验室进行实验。实验场景图如图2所示,包括一张放射治疗床、一台参数可调的仿真人体模型、一个摄像头、一台激光位移传感器,以及一个与实验1完全相同的正方体标记物。图像采集摄像头(罗技科技,瑞士)的动态分辨率为1280pixel×720pixel,最大帧频为30frame·s^-1;正方体标记物为木制的,其尺寸为3cm×3cm×3cm;放射治疗床(兰州科近泰基有限公司)有x,y,z三个方向自由度,可根据临床上的不同需要调整旋转角度;蓝光位移传感器(optoNCDT2300-50BL,米铱,德国)的绝对误差为40μm, 全量程光学误差(FSO)≤±0.08%,分辨率为7.5μm,FSO为0.015%;仿真医学体模内部由PLC(Programmable Logic Controller)构成,可通过调节电机升降参数来模拟人体不同状态的呼吸运动。由可编程逻辑控制器来控制仿真体模内部电机的工作参数。通过调节参数来模拟不同振幅和频率的呼吸模式。摄像头用来采集图像信息,通过图像处理算法实现体表标记物的追踪。仿真模型和摄像头都通过USB 3.0接口连接到计算机上用于输入和输出图像。

图 2. 实验2示意图

Fig. 2. Schematic diagram of experiment 2

下载图片 查看所有图片

2.2 实验过程

为验证最佳追踪效果的拍摄角度,实验2通过调整摄像头的位置,分别从三个角度(即水平侧视、从顶部向下俯视、从底部向上仰视)对体模进行图像采集。其中,水平侧视采用与实验1相同的位置参数:摄像头与标记物保持水平,距离标记物50cm,激光位移传感器距离标记物上表面50cm;从顶部向下俯视采集时,将摄像头置于体模头部上方,并使其与标记物保持水平,且距离标记物50cm,激光位移传感器距离标记物上表面50cm;对于从底部向上采集的仰视视角,将摄像头置于体模脚部下方,使其与标记物保持水平,其距离标记物50cm,激光位移传感器距离标记物表面50cm。

在确定图像采集视角后,考虑到放疗过程中的患者可能会由于各种突发状况出现不同的呼吸表现,设定三组参数来模拟病人最常见的三种呼吸模式:正常呼吸(振幅15mm,频率20次/min),微弱呼吸(振幅10mm,频率10 次/min),急促呼吸(振幅8mm,频率23次/min)。同时将激光位移传感器安装在固定支架上,调整激光射线使光斑对准标记物上表面的中心点位置,图像采集视角都为水平侧视。最后将激光位移传感器反馈的位移量与标记物角点的帧间像素差进行比较,得到实验结论。该算法的应用及验证都是在MATLAB 6.0 (MathWorks Inc.)上实现的。

2.3 图像预处理

为了减少背景中其他运动物体对标记物追踪的干扰,采用帧间处理的方法来提升追踪效果,并将处理后的图像降采样到每帧图像大小为540pixel×960pixel(每像素大小为0.58mm×0.58mm)。对降采样后的图像进行两个阶段的预处理:边缘增强和边缘检测。考虑到传统的均值滤波器有模糊边缘的缺点,为了保留标记物完整的边缘信息,依次选用一个模板为15×15、噪声方差为155的Lee滤波器,一个模板为13×13的Kuwahara滤波器,以及一个方差为1的局部零均值高斯滤波器对图像进行去噪,再使用“5抽头系数”的方法计算图像的梯度幅度和方向,得到边缘增强结果。图像的边缘信息使用Canny边缘检测算子来检测得到。预处理流程图如图3所示。

图 3. 图像预处理流程图

Fig. 3. Flow chart of image pre-processing

下载图片 查看所有图片

2.4 Kanade-Lucas-Tomasi算法

Kanade-Lucas-Tomasi(KLT)算法一种是基于特征点的目标追踪算法,是光流法的一种。它的提出主要基于三个假设条件:恒定的亮度、一致的空间以及微小的运动。以上三个假设是为了确保:1)目标不受亮度影响;2)可以在目标区域找到特征点;3)在同一窗口中所有特征点的位移大小都相等。该算法最先由Lucas和Kanade提出^[18],他们将图像配准问题从传统的滑动窗口搜索方式转变为计算像素点偏移量的过程。后来,Shi和Tomasi^[19]对KLT算法进行了进一步完善,给定图像 Fx和Gx,图像配准问题就转换为寻找一个视差向量 h以使 Fx与Gx的差最小的问题,第一个图像 Gx中的点 x移动到第二个图像 Fx中的点 Ax+d。那么,相邻两帧图像的变化可以表示为

式中: x为位移向量的分量; A是变形矩阵; d是平移向量。对于给定的两幅图 Gx和Fx,就是通过计算 A和d来实现对目标的追踪。但由于图像存在噪声以及仿射运动模型并不完善,因此对 A和d的求解需要对(2)式进行最小化实现:

式中: Σ为最小化数值; W是一个给定的特征窗口; wx是加权函数,在最简单的情况下, wx=1。

KLT算法依靠特征点匹配实现追踪,因此选取的特征点的质量直接影响到最后的追踪结果,而角点由于其丰富的形状信息,尤其适合作为特征点。将物体角点作为特征点的Harris角点检测算法具有良好的鲁棒性和定位能力,此外还能够很好地适应目标的旋转和缩放,因此本研究选取Harris算法作为角点检测算法。

2.5 感兴趣区域追踪

首先利用摄像头将采集到的图像上传到计算机,此时需要手动标记第一帧图像中标记物的4个角点以获取追踪的初始帧图像。ROI区域的宽度 w和高度 h的计算公式为

将4个角点横坐标的最大值 xmax减去最小值 xmin得到变量 x,而将4个角点纵坐标的最大值 ymax减去最小值 ymin得到的差 y即为变量 h的值。另外,此过程中还可以获取ROI左上角点的横坐标 hL及纵坐标 vL,即

式中: hm是ROI 4个角点横坐标的最小值; vm是ROI上方左、右两个角点纵坐标的最小值。根据上述计算得到的坐标点位置,该算法可自动获取ROI,使用基于最小特征值原理的Harris特征点检测得到角点位置,KLT算法通过追踪特征点来实现对目标区域的追踪。随机样本共识(RANSAC)算法用于构造相邻帧图像间的变换仿射模型,估计除第一帧外每一帧的角点位置及ROI区域。同时,算法会自动计算被追踪标记物的位移信息、翻转角度变化等相关参数,追踪流程图如图4所示。手动标定第一帧图像目标标记物的4个角点,将其余帧图像依次经过Harris角点检测、KLT追踪及仿射变换后,得到预测的ROI位置,循环上述步骤实现目标追踪。

图 4. 追踪算法流程图

Fig. 4. Flow chart of tracking algorithm

下载图片 查看所有图片

3 实验结果

3.1 图像预处理结果

图5为实验2中自头部向下拍摄的视角对应的预处理结果,由图可见标记物的边缘信息被明显增强,而背景中其他边缘信息被平滑。该过程有效地将标记物从背景中分割出来,减小了背景中可能存在的其他运动物体对目标追踪造成的影响,有效地提高了算法的追踪精度。其中,Lee滤波器利用图像局部统计特性进行图像斑点滤波,针对斑点噪声有较好的去噪效果。Kuwahara滤波器通过计算图像模板中邻域内的均值和方差,选择图像灰度值较为均匀区域的均值替代模板中心像素灰度值,达到平滑图像信息的目的。

图 5. 预处理结果示意图。(a)原始图像;(b)降采样后的图像;(c)经过Lee滤波器后的图像;(d)经过Kuwahara滤波器后的图像;(e)经过高斯滤波器后的图像;(f)采用Canny算子进行图像边缘检测后的结果

Fig. 5. Results of pre-processing. (a) Original image; (b) down-sampled image; (c) image filtered by Lee filter; (d) image filtered by Kuwahara filter; (e) image denoised by Gaussian filter; (f) edge detection result with Canny operator

下载图片 查看所有图片

3.2 实验1结果

实验采集受试者25s的呼吸过程进行实验。实验1的结果如图6所示。体表实际位移及标记物左上角点的纵坐标数值如表1所示。对两组数值进行相同的归一化处理以进行比较,以实际位移量为例进行如下数据处理。

式中: Nx为归一化后的数值; dx、d1、dmax、dmin分别为位移传感器得到的帧间实际位移差数组的第x个值、第一个值、最大值和最小值。通过上述归一化公式计算得到每组实际位移的对比值。如图6(c)所示,算法预测的ROI与标记物的实际位置基本吻合,为了进一步验证实验方案的追踪精度,将体表的实时位移与ROI左上角点纵坐标帧间像素差的归一化值进行对比,如表1所示。

图 6. 实验1追踪结果。(a)手动标定第一帧图像中标记物的4个角点;(b)算法根据手动标定的角点位置计算出的ROI; (c)KLT算法对第100帧图像的追踪结果

Fig. 6. Tracking results of experiment 1. (a) Mark four corner points of marker in first frame manually; (b) ROI calculated by algorithm according to marked coordinates of corner points; (c) tracking result of 100^th frame with KLT algorithm

下载图片 查看所有图片

实验中激光位移传感器的采集频率为624Hz,远大于算法追踪频率。为了使实验与算法处理图像的帧数相一致,将传感器采集的数据从18027降采样到398。另外,考虑到激光位移传感器与摄像头没有统一的触发开关,传感器的数据采集时长稍长于摄像头,因此仅保留了前370帧的数据用于绘制对比曲线,结果如图7所示。圆圈为终点的曲线表示呼吸运动的实际位移曲线,三角形为终点的曲线表示ROI角点纵坐标的帧间像素差曲线。由图可见,体表的位移变化曲线与标记物角点坐标曲线基本重合。人体体表呼吸运动的频率、幅值及波动状况可以通过算法追踪的角点信息得到表征,这证实了实验方案的可行性。

图 7. 体表实际位移与ROI角点纵坐标帧间像素差的对比曲线

Fig. 7. Curves of actual displacement of body surface and pixel difference between ordinate frames of ROI corner points

下载图片 查看所有图片

3.3 实验2结果

图8所示为三种视角的追踪实验结果。由图可见:在自头部向下拍摄和自脚部向上拍摄时,ROI位置与标记物的实际位置存在明显偏差,而水平侧视下ROI与标记物基本重合。分析其可能的原因是:KLT算法是一种基于特征点检测的追踪算法,实验中可能出现特征点因为被遮挡而丢失的情况,进而导致追踪效果的稳定性变差,甚至可能丢失追踪目标。

图 8. 三种视角下的追踪结果

Fig. 8. Tracking results from three perspectives

下载图片 查看所有图片

实验1已经验证可以通过一个体表标记物来表征体表的呼吸运动,而算法在图像中框选的就是算法预测的标记物的实时位置,以及ROI角点及中心点的实时位置坐标、旋转角度及帧间像素差等信息。为进一步验证算法表征呼吸运动的准确性,通过调整仿真体模的相关参数,模拟人体在三种不同呼吸状态下的运动信息(正常呼吸:振幅15mm,频率20次/min;微弱呼吸:振幅10mm,频率10次/min;急促呼吸:振幅8mm,频率23 次/min)进行追踪实验。

图9为在水平侧视视角下,三种不同呼吸模式[图9(a)]下的体模的运动曲线及追踪实验结果,其中图9(b)为模型PLC界面显示的电机运动轨迹,图9(c)为ROI中心点纵坐标在前1200帧图像的数值变化。将两组曲线进行对比,可以看出基于特征点的追踪曲线与模型运动曲线基本吻合,同时呼吸运动的频率、幅值以及其他关键参数信息都可通过计算获取,这证实了该算法表征人体呼吸运动的可行性。

图 9. 三种呼吸状态下体模运动曲线及算法追踪曲线对比。(a)三种呼吸状态;(b)体模运动曲线;(c)算法追踪曲线

Fig. 9. Comparison of phantom motion curves and algorithm tracking curves under three breathing conditions. (a) Three breathing conditions; (b) phantom motion curves; (c) algorithm tracking curves

下载图片 查看所有图片

为进一步验证追踪的准确性,将算法预测的ROI中心点的帧间坐标像素差与标记物的实际位移进行对比。虽然人体的呼吸运动不遵循严格的周期性,但体表在呼吸运动过程中的上下起伏存在一定的规律,经观察发现,每40帧的图像间隔可以有效表述一个呼吸运动周期。表2是以正常呼吸(振幅15mm,频率20 次/min)为例,将前360帧图像中每40帧作为一组,计算每组数据中ROI中心点的像素差值及标记物的实际位移,并对数值进行归一化处理,其中A代表帧间实际位移,B代表帧间像素差,C代表归一化处理后的帧间实际位移,D代表归一化处理后的帧间像素差。从表中归一化的结果可得:两组数据误差最大不超过0.03。

为更清晰地量化两组数据的误差,选取前800帧图像,将其分为20组,对数值进行归一化处理,依据表2将两组数据绘制成如图10所示的折线图,三角形标记点表示实际位移变化曲线,正方形标记点表示ROI特征点的像素值变化量,由图可见,两条折线的变化趋势基本一致,并且选取的定点特征值误差范围在0.03以内。

图 10. 标记物实际位移与特征点像素差对比曲线

Fig. 10. Actual displacement of marker and pixel difference between feature points

下载图片 查看所有图片

追踪算法常见的一个问题是追踪效果会随实验时长的增加而变差。而追踪效果的持续性和稳定性被认为是追踪任务的关键指标。为探讨本实验算法在完成呼吸追踪任务时的稳定性表现,在正常呼吸模式下重复实验3次,依次选取实验结果的第2帧、第100帧、第150帧、第200帧及第250帧图像,将选取图像的ROI特征点像素差与标记物的实际位移相比较,绘制如图11所示的误差棒。由图11中可以清晰地看到:追踪误差并没有随着追踪时长的增加而增大,未出现随着时间的推移而丢失目标的情况。

图 11. 追踪实验结果误差棒

Fig. 11. Error bar of tracking experimental results

下载图片 查看所有图片

4 结论

针对放疗场景下人体呼吸运动实时监测的问题,提出了一种基于KLT算法的通过追踪体表标记物表征体表呼吸运动的方法。经实验验证,该方法能在实时表征人体呼吸运动的同时,提供呼吸频率、呼吸幅值、呼吸状态变化等关键信息,为放疗手术提供数据参考。实验结果表明,算法追踪得到的呼吸运动曲线与实际呼吸运动曲线基本重合,归一化后实际误差小于0.03。

相比于其他技术难度大、实现复杂的呼吸监测方法,所提出的方法仅需要一个摄像头,通过图像处理算法即可有效采集人体呼吸运动的实时信息,这不仅降低了操作难度,更减少了传统磁共振成像呼吸监测方案对病人的辐照伤害,所提算法能够以更低的成本、更高的追踪精度实现患者的呼吸运动管理。

实验2确定了水平侧视为最佳图像采集视角,在其他两个视角下,ROI与目标的实际位置存在偏差的可能原因是:随着呼吸运动的进行,标记物的左下和右下两个角点被模型的主体部分遮挡,在角点检测中仅剩余两个角点被测得,而ROI的尺寸信息是基于第一帧图像计算得到的。因此在后续的ROI位置预测时,算法会按照已确定的ROI大小,基于检测到的两角点来绘制ROI,这就可能造成较大的追踪误差。虽然Harris角点检测是逐帧进行,当丢失的特征角点再次出现时算法就会自动校正,只要采样频率足够高,实际误差不会对真实场景的检测造成影响。但在某些视角下,丢失的角点可能在长时间内不会再次出现,算法就失去了自我校正的能力。事实上,在设计算法时,可以按照图像情况减少选择的特征点个数并绘制ROI,但这会降低ROI生成和追踪的准确率,本研究中不能以牺牲准确率为代价提升追踪的稳定性。因此,应通过选择最佳的拍摄角度来解决特征角点被遮挡的问题。

此外,在验证实验中,特征点的帧间像素差及标记物的实际位移皆为取自竖直方向上的量,其原因是算法得到的特征点纵坐标帧间平均差值为0.48pixel,而横坐标帧间平均差值为0.02pixel。因此呼吸运动在竖直方向上的运动范围更大,运动趋势也更加明显,进而能够更准确地表征呼吸运动。

实验算法的计算复杂度较高,如何在不降低追踪精度的前提下提高计算速度将是下一步工作的重点。