基于DICOM的多平台放射治疗数据可视化与对比分析

于大宇，唐丽玉，李小波，徐本华，张建平

1.福州大学空间数据挖掘与信息共享教育部重点实验室，福建福州350116；2.福州大学地理空间信息技术国家地方工程研究中心，福建福州350116；3.福建医科大学附属协和医院放射治疗科，福建福州350001；4.福建医科大学医学技术与工程学院，福建福州350004；5.福建医科大学临床医学部，福建福州350108

前言

恶性肿瘤已成为威胁生命健康的主要疾病之一。放射治疗、手术治疗和化学治疗作为肿瘤治疗的三大手段，可有效降低肿瘤的威胁。近年来，放疗技术已发展到可根据肿瘤的三维（Three-Dimensional,3D）形状调节辐射强度，使高剂量辐射在分布上适形肿瘤靶区的3D 适形调强放射治疗技术（Three-Dimensional Conformal Radiotherapy,3D-CRT），从而减少靶区周围正常组织的照射剂量，提供更好的治疗效果，3D-CRT现已广泛应用于临床肿瘤治疗［1］。调强放射治疗（Intensity-Modulated Radiotherapy，IMRT）技术在计划制定阶段对目标函数实施逆向优化，可根据不同靶区处方剂量要求得到更加适形的非均匀的剂量分布，相比传统3D-CRT技术更为复杂，对质量保证与质量控制有较高的要求［2］。

IMRT计划的客观目标是根据目标函数寻求最优解，使等剂量曲线达到理想的分布，优化过程依赖于商用3D-CRT治疗计划系统（Treatment Planning System,TPS）中的剂量计算引擎［3］。然而剂量计算引擎在精度上存在局限性，并受时间限制，计算过程中易产生不确定性，可能影响整体计算精度。Jeraj等［4］和李小波等［5］研究展示了IMRT优化过程中使用笔形剂量计算算法存在的系统性误差和收敛误差。Prabhakar等［6］使用60例不同癌肿区域的病人数据研究Eclipse、Pinnacle等不同TPS对肿瘤组织体积的计算差异，结果表明所有TPS对同一肿瘤组织体积的计算结果都存在不同程度的差异。IMRT合作工作组讨论了质量保证的过程和剂量验证的方法［2］。Ayyalasomayajula等［7］将Eclipse TPS计算的3D-CRT 剂量分布与研发的ROPS独立计算的剂量分布作对比，以评估Eclipse 3D剂量分布。Nelson等［8］和Park等［9］分别开发了一个独立的剂量评估系统，用于二次验证及评估TPS计算的复杂的3D适形剂量分布。对IMRT这种高度适形放射治疗技术进行质量控制是十分重要的。

目前，国外的研究组织和医疗机构对放射治疗有相对广泛深入的研究，国内的研究主要集中在TPS的研发，以及放射治疗剂量计算、优化算法等［10-12］，对放射治疗数据的可视化及质量控制的研究相对较少。

DICOM-RT 是DICOM 标准的一个扩展模块，是NEMA协会成立的放射治疗标准工作组专门制定的关于放射治疗的数字影像和传输标准［13］。DICOM-RT已是放射治疗领域应用最广泛的放射治疗标准，主流商用TPS均支持以DICOM-RT格式导出数据。DICOMRT具有良好的可扩展性但部分商用TPS具有封装性，不同商用TPS的部分数据无法互相传输与使用。

本研究为满足研究团队研发的远程放射治疗质量控制数据的集成管理与应用平台的需要，以质量保证/质量控制为重点，并基于DIOCM协议研发放射治疗数据可视化与分析系统（RTViewer）。系统可对TPS 生成的放射治疗计划文件进行独立的剂量体积直方图（DVH）计算，提供不同TPS的DVH、点剂量及等剂量曲线分布图对比，用于在服务器上了解和分析不同放射治疗单位上传的放射治疗数据。

1 DICOM-RT解析及可视化

1.1 DICOM-RT对象及数据结构

DICOM-RT 的数据结构和编码方式与DICOM标准一致。如图1所示，DICOM 文件由文件头和数据集两大部分组成。文件头由导言和前缀组成，主要用来标识此文件是否为DICOM文件；数据集由一个个数据元素组成，数据元素是DICOM的基本组成单元，每个数据元素对应现实世界的一个对象实例。数据元素包含标签、值表示、数据长度和值域这4部分。标签是数据元素的唯一标识，值域是数据元素的值，值的字节长度由值长度定义，值表示描述值域的数据类型。DICOM 文件采用二进制方式编码，符合计算机逻辑运算，易于转换，可靠性高。

图1 DICOM数据结构Fig.1 DICOM data structure

DICOM-RT模块增加了许多与放射治疗相关的数据元素，主要分为5个对象：RT-Image、RT-Structure Set、RT-Plan、RT-Dose、RT-Record（图2）。DICOM标准对每个RT对象都有详尽的定义，不同RT对象含有特定数据元素，负责不同的存储任务。RT-Image存储数字重建透视图、射野验证图像等放射治疗图像信息；医生在TPS中勾画的肿瘤靶区和危及器官等保存在RT-Structure Set对象中；RT-Dose主要存储TPS计算的3D剂量分布和DVH信息；RT-Plan包括体位、处方、射野等照射治疗计划信息；RT-Record为放射治疗记录对象。

图2 DICOM-RT对象关系图Fig.2 DICOM-RT object relation diagram

DICOM-RT的对象之间相互参照、相互关联，对象之间的逻辑关系较强。对DICOM-RT的解析涉及到实体关系模型的分析、元信息解读和数据元素提取等,解析过程较为复杂。根据实际需要，本研究选取在剂量分布评估上重要且必须的RT-Structure Set和RT-Dose对象，论述其解析方法。其中，对DICOM-RT的解析使用DCMTK工具库实现，解析和分析的结果使用Visualize Tookit（VTK）可视化。

1.2 RT-Structure Set对象解析及可视化

RT-Structure Set 文件包含以下数据模块：Patient、Study、Series、Equipment、Structure Set。其中，前4 个为DICOM 文件中的常见数据模块，Structure Set 为RT-Structure Set 对象的特有模块。Structure Set 由Structure Set、ROI Contour、RT ROI Observation 和SOP Conmen 这4 个子数据集组成，每个子数据集嵌套一个或多个数据集。Structure Set中定义了放射治疗中医生勾画的重要组织区域的三维体积轮廓，如身体外轮廓、危及器官、肿瘤靶区（包括肿瘤靶区、临床靶区和计划靶区）及其它感兴趣区（ROI）［14-16］。ROI提取过程见图3。

将ROI 保存到内存后，利用Contour 的Referenced SOP Instance UID 与每一帧CT 影像的SOP Instance UID 做匹配分析，若相等，则认为该Contour 位于该CT 影像上。将Contour 和CT 影像的坐标转换为屏幕坐标系，然后将两个图层（CT影像为第一图层,Contour 为第二图层）进行叠加显示，最终显示结果见图4。

ROI由医生在TPS中根据经验绘制而成，一份完整肿瘤靶区勾画的制作需2～4 h，是放射治疗计划中最为珍贵的数据。将肿瘤靶区勾画从TPS中提取出，导入到Matlab、R 语言等独立的分析环境中，可用于检验自动勾画的准确度、训练神经网络及进行剂量计算等［14］。因此，该系统提供RT-Structure Set输出接口，可以利用ASCⅡ编码方式输出肿瘤靶区结构。

图3 ROI提取Fig.3 Extraction of region of interest

图4 ROI显示结果Fig.4 Representation of region of interest

1.3 RT-Dose对象解析及可视化

TPS 计算的3D 剂量分布数据通常以3D 格网形式表达，在RT-Dose对象中通常用二维剂量平面形式表达3D剂量格网，即用多帧二维剂量图像表示3D剂量分布［12］。

从图5可以看出，完整的3D 剂量格网沿Z 轴按层切分，每一层被等分为体积相同的体素。通常用二维剂量图像表示一层的体素，剂量图像的Pixel Spacing（像素间隔）代表单个体素的长和宽，Slice Thickness（层厚）表示体素的高。体素是3D 剂量分布中最基础的元素，其内部的剂量分布是匀质的，剂量值可通过像素值转换而得，计算公式为：

其中，PV 为图像的像素值，一般为int 类型；DGS 为DoseGridScaling数据元素（Tag：3004,000E）的值域，用以将像素值转为指定单位的剂量值。剂量单位可通过Dose Units数据元素获得。剂量分布的坐标通过Image Position、Pixel Spacing和Grid Frame Offset Vector数据元素隐式地表示，Image Position可认为是3D剂量格网的X、Y、Z 坐标原点的真实世界坐标值。Grid Frame Offset Vector数据元素为存有每一帧剂量图像Z值的一维数组。基于以上数据元素，经计算后，可获得3D剂量格网内任意位置的坐标与剂量值。

除点剂量的显示外，剂量分布常以剂量分布图的形式直观地显示。要获得某一区域的剂量分布图，需先建立一幅对应空间范围的空像素值图像，依次对图像的每个像素填充像素值，像素值通过对应坐标的剂量格网的像素值获取，具体方法见上述分析。填充所有像素后，即可得到灰度剂量分布图。调用VTK库的vtkLookupTable类，创建颜色映射表，可将灰度剂量分布图映射为彩色图像，凸显剂量的差异［1］，结果如图6所示。

图6 剂量分布伪彩显示Fig.6 Dose distribution in pseudo color

DVH 模块是RT-Dose 对象中的一个可选模块，部分TPS导出的RT-Dose内部包含DVH模块。DVH包括有差分型DVH和累积型DVH这两种形式，需先根据DVH Type数据元素判断其类型，然后再根据其类型，使用对应的计算算法将DVH Data数据元素中的数据以DVH 的形式显示。DVH 模块与Structure Set类似，其部分数据元素采用嵌套结构，解析时由外向内的顺序逐层解析。解析完成后，依据Referenced ROI Number匹配到从RT-Structure Set中解析的相对应ROI，获取参照ROI 的名字、颜色等信息。图7为RT-Dose对象中的Eclipse TPS计算的DVH。

图7 Eclipse计算的DVHFig.7 Dose-volume histogram(DVH)calculated with Eclipse

1.4 RTViewer系统的实现

在Windows 10 64-bit操作系统上，采用C++开发语言，集成开发环境为Visual Studio 2015，利用DCMTK 开源库和VTK 可视化开源库，开发形成放射治疗数据可视化与分析系统RTViewer。系统可对TPS 生成的放射治疗计划文件进行独立的DVH 计算，并提供不同TPS的DVH、点剂量及等剂量曲线分布图，用于在服务器上了解和分析不同放射治疗单位上传的放射治疗数据。

2 放射治疗数据分析

2.1 DVH计算算法

DVH表达了3D剂量分布的统计学信息，为放射治疗计划评估及剂量验证提供重要依据。一些TPS（如XIO）并没有在RT-Dose 对象中提供DVH 模块，同时，Ebert等［17］研究也发现不同TPS计算DVH所使用的算法各不一致，计算的DVH 也均有差异，因此，TPS计算的DVH无法用于对比不同TPS剂量分布的差异。目前，已有部分学者针对如何使用剂量分布和Structure Set 计算DVH 进行研究［18-19］。本研究提出一种新的方法，以ROI 的Contour 作为掩膜截取剂量图像，获取ROI内所包含体积对应的剂量，用于计算DVH及其它统计学参数。图8为一个ROI的DVH计算流程图，要获得所有ROI的DVH，遍历所有ROI执行算法即可。

由于CT 影像与剂量图像的帧厚不一致（如Eclipse的CT帧厚为5.0 mm，剂量图像帧厚为2.5 mm）。如果直接将CT的帧厚作为掩膜的厚度，会导致计算结果不精确，故在此算法中，掩膜的厚度x取CT帧厚和剂量图像帧厚的最大公约数。数组A 含有n 个数值，每个数值代表体积为Dpw×Dph×x 的组织的剂量值，其中，Dpw 和Dph 分别为剂量图像的像素宽和像素高。因此，该ROI包含的所有组织体积为：

ROI 在单帧CT 影像上的Contour 数据主要有3种类型（图9）。环状Contour 的ROI 区域为两侧Contour之间部分，然而在计算过程中，会基于内侧和外侧Contour生成两份掩膜，导致提取的剂量比实际情况多，需对其进行特殊处理。故在制作掩膜之前，需先判断该Contour 是否为环状类型，若为环状类型，需用内侧Contour 计算的数组对外侧Contour 的数组取差集。

2.2 数据分析

图8 DVH计算算法流程图Fig.8 Flow chart of DVH calculation algorithm

为满足放射治疗数据控制的需要，实现对不同TPS数据的相互验证和对比分析，该系统对剂量分布进行独立的DVH计算，以提供不同TPS数据的DVH对比分析以及剂量分布的统计信息对比、点剂量对比、剂量分布图对比和等剂量线对比等分析功能。

当前判断放射治疗计划是否达标主要凭借DVH分析。因此，该系统除可显示和对比不同TPS 的DVH 外，还提供DVH 统计指标：最小剂量、最大剂量、平均剂量、V95和均匀性指数。V95是达到95%处方剂量的体积百分比。HI用于分析和量化靶区体积内剂量分布的均匀性，其定义为：

其中，D2%和D98%分别为受照剂量最高的2%和98%，例如，D98%代表至少98%的目标体积达到这个剂量，因此D2%和D98%可认为是最小剂量和最大剂量；DP为处方剂量。文献［20-21］对HI有相同的定义，在理想情况下HI趋向等于0。

图9 Contour的3种类型Fig.9 Three types of contour

对比不同TPS的DVH是评估剂量分布的一个较好方式［8,22-23］。为更好地分析剂量分布，RTViewer 可相互验证和对比分析不同TPS的放射治疗数据，提供基于同一套CT 的不同TPS 的DVH 对比分析。不同TPS 基于相同的序列CT 影像计算生成剂量分布后，导入到RTViewer 中，经系统分别计算DVH 后，对比不同DVH 的差异。利用RTViewer 对比不同TPS 的DVH（图10），对比的数据为XIO 和Eclipse 系统对同一份数据分别使用Superposition算法和AAA算法计算而得的3D-CRT剂量分布。表1为该系统根据不同TPS 中的数据计算的DVH 指标值，数据来源于Eclipse 系统和XIO 系统基于同一份CT、靶区勾画、处方剂量生成的3D-CRT 计划，分别采用AAA 算法和Superposition算法，光子线处方剂量为5 040 cGy。

图10 不同TPS的DVHs对比Fig.10 Comparison of DVH of different treatment planning system(TPS)

DVH分析的缺点在于未考虑剂量分布的空间位置信息，若仅凭借DVH指标判断剂量计划是否达标，会存在安全隐患。因此，为考虑到剂量分布的空间位置信息以全面地分析剂量分布，RTViewer 提供了不同TPS剂量分布的点剂量对比、剂量分布图对比和等剂量线对比等功能。表2给出RTViewer系统计算的点剂量值，可对比任意空间位置的不同TPS计算的剂量值。图11 为等剂量线对比图，可对比不同TPS的等剂量线的差异。红色和粉红色分别为Eclipse和XIO 系统的50 Gy 等剂量线，绿色和浅绿色分别Eclipse和XIO系统的40 Gy等剂量线，蓝色和浅蓝色分别Eclipse和XIO系统的20 Gy等剂量线。图12为不同TPS的剂量分布对比图。

2.3 结果验证

从XIO、Eclipse、Monaco、Pinnacle 及Raystation系统各随机选取10例临床病例数据导入到RTViewer中，将RTViewer 对临床病例解析、显示、分析的结果与原TPS 对比。经过验证，RTViewer 可兼容XIO、Eclipse、Monaco、Pinnacle 及Raystation 输 出 的DICOM 格式数据，对序列CT 的解析及可视化正确，对靶区结构勾画、剂量分布信息的解析及分析与原系统一致。对于DVH 计算，各个TPS 的计算算法不一致，对同一份数据计算的DVH也均有略微差别，根据本研究提出的DVH计算算法计算而得的DVH，与原TPS的差异小于1%。

表1 Eclipse和XIO系统DVH指标对比Tab.1 Comparison of DVH indexes of Eclipse system and XIO system

表2 不同TPS的点剂量对比Tab.2 Comparison of point doses of different TPS

3 结论

图11 不同TPS的等剂量曲线对比Fig.11 Comparison of isodose lines of different TPS

本研究研发的RTViewer能够支持放射治疗数据快速、简便、精准地提取，加载，显示和分析，实现可自由调用主流TPS 的数据并进行相互验证和对比分析。本研究也提供了一种DVH 计算算法，并可根据需要导出RT文件中的靶区勾画、剂量分布等重要数据到Matlab 等环境。商用TPS 的DICOM 格式放射治疗数据，经RTViewer解析和可视化后，系统可对剂量分布进行独立的DVH 计算，并提供不同TPS 数据的DVH对比分析、剂量分布的统计信息对比、点剂量对比、剂量分布图对比、等剂量线对比等功能。经临床病例验证，系统兼容XIO、Eclipse 和Monaco，Pinnacle 和Raystation 等TPS 导出的放射治疗数据，对数据的解析及可视化结果与原系统一致。系统具有友好的用户界面，较强的交互性及操作性，可结合研究团队研发的远程放射治疗质量控制数据的集成管理与应用平台，用于肿瘤放射治疗远程质量控制。

图12 不同TPS的剂量分布图对比Fig.12 Comparison of dose distribution of different TPS