基于Python程序的量子探测效率分析处理系统的研制

徐 桓 刘 敏 赵庆军 刘 洋 李晓亮

量子探测效率(detection quantum efficiency，DQE)是描述数字X射线影像探测器成像性能的一项主要参数，可测量探测器对入射到探测器表面的X射线的吸收能力[1]。DQE值高则表明X射线的转化效率越高，系统能够在低入射剂量条件下获得更高的图像质量，其涉及到数字X射线影像探测器的灵敏度、线性、噪声剂量、调制传递函数(modulation transfer function，MTF)等多项参数，能够实现对影像探测器的全面性能评价[2]。为了规范DQE的检验，国家出台了医药行业标准《医用电气设备 数字X射线成像装置特性 第1部分：量子探测效率的测定》(YY/T0590.1)和《医用电气设备 数字X射线成像装置特性第1部分：量子探测效率的测定 乳腺X射线摄影用探测器》(YY/T0590.2)。基于DQE检测的研究很多，但该项目作为日常检测项目开展仍存在一定难度，主要是由于DQE检测涉及调制传递函数、噪声功率谱等一系列的数学运算，要求测试人员掌握一定的数字图像处理的专业知识，并编制相应的图像处理程序方能完成，检测难度较大[3-6]。美国物理师协会(American Association of Physicists in Medicine，AAPM)第162号报告为了规范化DQE的测试，发布了一套软件系统，但该系统只支持Linux与Mac系统，同时，无图形用户界面，操作并不十分方便[7]。为此，本研究针对上述问题设计开发出一套DQE分析处理系统，采用图形化用户界面，全面规范DQE测量过程，使DQE的测试成为一种可操作性强的标准化测试方法，能有效解决从事DQE检测的工程技术人员开展DQE检测难度大的问题。

1 DQE分析处理系统设计

DQE分析处理系统(DQE analyzer)的目的是解决DQE检测工作实际开展难度大的问题。系统设计涵盖了DQE测试所需的关键流程，分别为：①获取转换函数，实现数字X射线图像的线性化处理；②MTF的分析、噪声功率谱(noise-power spectrum，NPS)分析、DQE的分析计算等功能，为了方便系统的使用，设计实现了包括图像的感兴趣区域(region of interest，ROI)选取、图像翻转等基本的图像处理工具，为用户进行DQE分析提供支持。系统设计了可视化的图形用户界面，提供灵活的参数输入和设置，全面实现了规范化的DQE分析过程，具有自动化程度高、用户操作简洁高效、结果呈现图表化等特点，能够极大提高DQE测试的准确性和效率。软件主界面和工作流程见图1。

图1 DQE分析处理系统主界面和工作流程图

2 DQE分析处理系统开发工具

本研究利用Python程序语言开发实现。Python是一个高层次的结合了解释性、编译性、互动性和面向对象的脚本语言，并且可跨平台、能够兼容UNIX，Windows和Macintosh等系统[8]。Python程序的最大优势是其丰富的扩展库，其扩展库包括Pydicom、Numpy、Matplotlib和Pyqt。Pydicom库是用于处理医学数字成像及通信(digital imaging and communication of medicine，DICOM)格式的医学图像库，利用Pydicom实现了DICOM格式的图像数据的读取[9]；NumPy系统是Python库的一种开源数值计算扩展库，用于存储和处理大型矩阵，本研究主要利用Numpy库实现图像数据的矩阵运算；Matplotlib是一个强大的Python数据可视化库，能够实现散点图、折线图、直方图、柱状图、箱线图等多种样式图标的绘制，利用Matplotlib库实现系统中各种分析计算结果的图表化显示；Pyqt5是QT开发框架的Python语言实现，利用Pyqt5开发了DQE分析处理系统的图形用户界面(graphical user interface，GUI)，为用户使用提供了便利[10]。

3 DQE分析处理系统实现

DQE分析处理系统严格按照国家医药行业标准《医用电气设备数字X射线成像装置特性 第1部分：量子探测效率的测定》YY/T0590.1-2018(IEC 62220-1-1:2015)规定的要求实现。相比之前的YY/T0590.1-2005，2018版的标准做了一些修订，包括标准中关于射线质量(radiation quality)的模型有了更新，对应的输入噪声功率谱的值也发生了变化；对于MTF的测量，边缘器件需要连续旋转4个90°来完成测量，同时增加了对45°对角线的MTF测量要求等。

3.1 数据线性化

进行DQE分析之前，必须确保图像数据与曝光剂量的线性关系。数据线性化是通过获取转换函数实现，按照国家医药行业标准YY/T0590.1的要求设置测试装置，并用剂量计测量不同曝光水平下的空气比释动能。空气比释动能的测量范围应包括从0到基准空气比释动能的4倍。不同的曝光次数取决于转换函数的形式，如果只是检查转换函数的线性，只需在期望的测定范围内均匀的划分5次曝光点即可；如果是测定完整的转换函数，空气比释动能的最大增量按对数计算(以10为底)≤0.1。

将不同曝光水平下得到的图像保存，推荐按照曝光顺序对图像进行重命名。通过DQE分析处理系统打开需要进行线性化分析的图像数据，并录入不同图像对应的空气比释动能，系统会自动选取图像中心区域的100×100像素区域进行分析，并拟合出转换函数，并得到拟合结果的相关系数R。系统按照线性拟合、对数线性拟合和多项式拟合3种不同的方式进行，如果线性拟合结果R2(R为相关系数)＞99%，则接受该拟合结果；如果拟合结果为非线性，则自动通过拟合得到的转换函数对图像进行线性化处理，然后才能将线性化的图像用于MTF和NPS分析。

3.2 MTF分析

对于数字化X射线设备，MTF可以采用狭缝或者边缘器件来实现。由于狭缝测量对于摆位精度要求过高，国家医药行业标准YY/T0590.1中规定MTF测试用边缘器件来实现。基于边缘器件进行MTF测试的研究有很多，具体的实现有一定区别[11-13]。本系统的实现严格按照国家医药行业标准YY/T0590.1的要求进行，实现MTF的自动分析。①测量时的试验器件按照YY/T0590.1的要求实施，得到一张基准水平的黑白边缘图像，边缘应与探测器的像素矩阵的行或列的排列有一个1.5°～3°的夹角；②确定MTF需要全长度的边缘扩展函数(edge spread function，ESF)，需选取跨越边缘选取100 mm×50 mm的ROI进行分析，在该ROI中，提取跨越边缘的N行(或列)像素数据，用于产生过采样的ESF，得到过采样ESF的方式有很多种，本系统严格按照国家医药行业标准YY/T0590.1实施；③对ROI中的数据重复上述运算，得到多组过采样ESF并进行平均，对平均处理后的过采样ESF进行差分运算，得到线扩展函数；④对线扩展函数进行傅里叶变换，即可得到MTF。

对于MTF分析需要获知一些关键参数，包括像素尺寸、ROI尺寸、边缘与图像矩阵的行或列的夹角角度α、过采样ESF的N值、对有限差分算法导致的频率平滑效应进行校正的窗函数等。像素尺寸和ROI尺寸用于确定进行MTF计算的数据范围，也是确定空间频率范围的关键参数。像素尺寸直接由DICOM图像的文件信息中提取；ROI尺寸则在设置中输入。本系统提供了ROI选择工具，利用该工具在图像上选取跨越边缘的ROI，系统会自动在其中定位边缘中心点，并根据中心点提取由输入的ROI尺寸确定的区域，一般为100 mm×50 mm；边缘夹角α的准确性对于确定过采样ESF至关重要，由于过采样的N值计算是由该夹角α确定，即N=1/tanα，一般取最接近的整数。为了得到准确的N值，本系统对黑白边缘图像进行滤波、边缘提取后，通过霍夫变换计算得到边缘的直线方程，进而确定夹角α。为了增加系统的灵活性，提供手动输入N值的方式，某些情况下通过适当调整N值，可以得到更为理想的MTF结果。

3.3 NPS分析

功率谱是功率谱密度的简称，定义为单位频带内的信号功率，表示信号功率在频域的分布状态。NPS即为图像中的噪声信号的功率谱，由图像噪声协方差的傅里叶变换后的模计算得到，其计算为公式1：式中ΔxΔy为水平和垂直方向的像素间距的乘积；M为ROI的数量；I(xi,yi)为—线性化数据；S(xi,yi)为可选的拟合多项式。

NPS由将若干幅均匀图像分割成的多个大小相同的ROI计算得到。为了精确计算功率谱，需要至少4百万个独立的像素数据[14]。而国家医药行业标准YY0590.1中提到，为了保证测量精度，应至少计算960个256×256的ROI，即包含1 600万个独立像素，因此需要从多幅图像中得到这些数据。

本系统实现了计算NPS的全部数学计算过程，只需要设置数个基本参数，并在一幅图像的均匀区域选取好一个矩形的ROI区域。点击计算功率谱按钮即可得到结果，其计算过程如下：①选择若干幅均匀图像，并对图像数据进行线性化处理；②利用一个二维二阶多项式去除图像中的趋势影像，即去掉图像非均匀性的影响。系统自动对均匀区域的图像数据进行分析，包括自动定位256×256的ROI，并按照标准要求的方法在水平和垂直方向移位。获取全部ROI中的数据后，进行二维傅里叶变换，得到二维的NPS，在二维NPS的图像上，在X轴和Y轴两个中心轴两侧各选取7行数据进行平均，计算出X轴和Y轴方向的一维NPS。计算一维NPS时中心轴该步骤涉及的关键参数包括空气比释动能、计算所用的图像数、像素尺寸以及拟合用的多项式阶数等。空气比释动能需要根据实际的辐照数据输入，用于确定输入的NPS；图像数、像素尺寸则根据实际选取的图像计算得到。

3.4 DQE分析

DQE可以认为是探测器利用输入的量子数产生图像信噪比的相对效率，是对数字X射线成像系统的性能整体评价，定义计算为公式2：

式中MTF(u,v)为数字X射线成像装置的预采样调制传递函数；Win(u,v)为在探测器表面辐射野的噪声功率谱；Wout(u,v)为数字X射线成像装置输出的噪声功率谱。

通过之前的步骤，计算DQE的参数大部分已分析完成，而Win(u,v)为输入的噪声功率谱，数值根据测量所选择的射线辐射质量查表得到。系统按照公式(2)进行计算分析，并按照标准中规定的频率间隔进行插值和平均处理后，即可得到X轴和Y轴两个方向的DQE数值。同时，本系统还引入了有效的量子探测效率(effective DQE，eDQE)的分析计算，用于评估临床相关的系统效率。

3.5 非标准射线质Win(u,v)计算

DQE的测量通常是在规定的辐射质量下进行，如常见的RQA3～RQA9的标准辐射条件。国家医药行业标准YY0590.1中也规定了上述标准辐射条件下的空气比释动能信噪比的平方(SNR2)值，该值是计算Win(u,v)所必须。对于非标准辐射条件，由于光谱特性未知，无法获取SNR2值，从而无法得到Win(u,v)，则无法进行DQE分析计算。为了弥补这一问题，本系统集成了自定义的射束条件仿真程序，能够实现对常见的X射线常规摄影和乳腺摄影的射束进行仿真，计算得到光谱曲线和SNR2值，用于计算Win(u,v)。自定义射线质Win(u,v)计算见图2。

图2 自定义射线质计算

4 不同剂量X射线数字图像DQE分析

4.1 DQE分析方法

对3种不同剂量水平的X射线数字图像进行DQE分析，并对结果进行比较。具体分析方法如下。

(1)对X射线数字图像进行数据线性化分析。根据6个不同的剂量水平的图像与实际测试得到的空气比释动能进行了线性化拟合，其数据见表1。

表1 数据线性化的实测值

(2)对图像分别沿X轴和Y轴两个方向进行MTF分析。对3个剂量水平(分别为3.25 μGy、10.8 μGy、34.2 μGy)的图像进行了NPS分析。测量MTF和NPS的参数设置分别见表2和表3。

(3)根据上述步骤得到的MTF和NPS计算得到DQE。

表2 调制传递函数分析的参数设置

表3 噪声功率谱分析的参数设置

4.2 DQE分析结果

经拟合后证实图像数据相对于输入的空气比释动能为对数关系，相应的转换函数及拟合的曲线见图3。

图3 拟合的转换函数

对X轴和Y轴两个方向进行MTF分析。由于进行了过采样处理，MTF的结果会根据过采样的N值扩展至N倍的图像截止频率，实际的结果只截取了低于截止频率的数据。MTF结果见图4。

图4 调制传递函数分析结果

噪声功率谱在不同的剂量水平下的结果有明显差异，对3个剂量水平的图像进行了分析，每个剂量水平选取10幅图像，按照256×256的ROI进行划分，其结果见图5。最终计算得到的DQE结果见图6。

4.3 结果验证

图5 噪声功率谱结果

图6 量子探测效率测量结果

MTF和NPS是计算DQE的两个关键步骤，其准确与否直接影响DQE的结果准确性。为了保证DQE分析结果的准确性，分别对MTF、NPS和DQE的分析结果进行了验证。首先分别建立了MTF和NPS的理论模型，并与实际的计算结果进行比较。为了验证MTF，根据设置边缘的角度生成一张理想的无噪声的黑白边缘，并将该图像输入本系统，计算得到MTF。同时，由于理想边缘函数ESF得到的线扩展函数(line spread function，LSF)应为矩形函数，进行快速傅里叶变换(fast Fourier transform，FFT)后应得到sinc函数:

理论计算的结果与实际计算结果进行了对比，结果具有很好的一致性，见图7。

图7 sinc函数与理想边缘的结果对比

对于NPS的验证，根据利用程序自动生成了8幅均值为10 000的均匀图像，图像大小与NPS计算时选取的大的ROI尺寸相同。然后在图像上添加了方差为100的高斯白噪声。高斯白噪声的功率谱在不同频率上是相同的，该幅值计算为公式3：

式中μ为白噪声图像的平均信号值；σ2为白噪声图像的方差；du和dv分别为水平和垂直方向上的频率增量；N为白噪声图像的全部像素数量。

实际计算得到的NNPS幅值为3.99×10-8。然后用本系统的计算过程对8幅高斯白噪声图像进行了分析，与理论计算值进行了比较，显示其具有很好的一致性，见图8。

图8 噪声功率谱的验证

最后对DQE分析结果进行了验证，对与已发表文献[9]中的图像数据进行了分析，并与文献中的结果进行对比，二者在整个频率范围上均具备很好的一致性，截止频率均为3 mm-1。对比结果见图9。

图9 量子探测效率的结果验证

5 结论

本研究利用Python程序语言开发实现了DQE分析处理系统，该软件系统能够实现DQE分析处理的全流程，为数字X射线成像装置的质量控制提供了有效的工具。系统的设计实现完全依据国家医药行业标准YY/T0590.1实施，能够保证DQE测试流程的规范性和测试结果的准确性。为了进一步证实测试结果的准确性，本研究对软件系统的关键计算步骤进行了验证，通过与理论模型的对比，证实了理论结果与模拟测试结果具有很好的一致性。相比已有的程序，本软件系统的优势在于具有良好的用户交互界面，除必须的输入参数外，每个单独的分析计算基本都通过点击按钮即可完成，且系统具有良好的灵活性和扩展性。对于标准中未明确的步骤，本系统既提供了默认的参考值，让不熟悉计算原理的用户能够得到可接受的分析结果，又可根据实际测试实验开展情况，提供该参数的录入，保证得到最优的结果。

本研究开发实现的DQE分析处理系统为边缘角度的计算提供了不同的方法，既可以通过自动分析计算得到边缘角度，用户也可根据实际测试的情况对角度值进行调整，确保能够得到最优的结果。系统还为非标准的射线辐射条件提供了仿真计算，为实现非标准辐射条件下的DQE计算提供了可能。软件系统能够提高DQE测试的规范性和测试效率，为从事该项工作的同行提供参考。