CN104083177A - 数字x射线成像系统的分辨率性能评估方法 - Google Patents

Abstract

本发明提供一种数字X射线成像系统的分辨率性能评估方法，包括：①连续采集多幅刀口图像；②对所采集的多幅刀口图像叠加平均，而后对平均后的刀口图像进行刀口边界检测，并且利用直线拟合获得刀口边界直线；③获得刀口倾斜角度，构建出过采样ESF曲线；④采用单调曲线拟合方法对过采样曲线消噪；⑤选择适当的平移距离，构建真实ESF曲线与其平移距离结果的差值曲线，对其傅里叶变换得到sinc函数；⑥对实际过采样ESF曲线进行适当的平移，而后与原来的过采样曲线相减，得到对称的过采样曲线；⑦利用频域下恶化模型公式求解系统调制传递函数。本发明可以获得较传统刀口测量方法更为准确的MTF曲线。

Description

数字X射线成像系统的分辨率性能评估方法

所属技术领域

本发明属于生物医学工程及计算机领域，涉及一种数字X射线成像系统的分辨率性能评估方法。 

背景技术

调制传递函数(Modulation Transfer Function,MTF)是是调制度的传递函数，是对线性影像系统空间频率传输特性的定量描述，是影像评价方法中的重要进展。此前普遍使用定性描述指标，例如影像密度、对比度、清晰度、分辨率及失真度等来评价成像系统的影像质量，但其结果受个人主观因素影响大。近年来随数字化X射线成像技术迅速发展，MTF作为客观指标已成为放射成像工作者和研究者所关注的重要影像评价手段。同时，MTF也是获得成像系统探测量子效率的必须参数。 

要定量地评价数字X射线成像系统的固有成像质量，只需计算不受个人主观因素影响的系统固有预采样MTF。国际电工委员会推荐采用刀口法作为MTF的测量方法，由于刀口法测量仪器自加工相对容易，在科研实验和常规检测中使用较为广泛。采用刀口法可以获得金属模块刀口的边缘扩散函数(Edge Spread Function,ESF)，它反映了经成像系统后边缘弥散程度，ESF的导数为线扩展函数(Line Spread Function,LSF)，再经过傅立叶变换并进行归一化后即可获得系统MTF。 

为了获得高分辨率的MTF曲线，目前一般采用改进的MTF刀口测量方法，即将刀口方向与图像采样方向成一定角度，以此获得过采样的边缘响应函数ESF，而后通过傅里叶变换得到高分辨MTF曲线。由于LSF曲线的获取是通过对ESF求导实现的，ESF曲线中的噪声被求导过程放大。因此，通过对LSF进行傅里叶变换来获取MTF的方法受到LSF信噪比降低的影响，阻碍了MTF特性的准确获取。 

在常规的系统性能评估报告中，需要给出系统分辨率和噪声特性的综合评估，即采用量子检出效率(detective quantum efficiency,DQE)作为评估标准。DQE是MTF的平方和噪声功率谱倒数的综合作用结果，因此，MTF的准确测量对获取准确的DQE曲线至关重要。而在实际系统评估时，不仅需要给出DQE幅度谱曲线，同时需要给出低频或零频率位置(如0.5mm¹)的DQE值。由于零频率的DQE值会受到系统低频伪差的影响，一般采用低频特性的DQE曲线来拟合得到零频率DQE值。无论是给出低频还是零频率位置的DQE值，低频段MTF的准确测量对该DQE值得获取至关重要。而传统方法通过求导得到的LSF曲线获取MTF，导致LSF信噪比低于原始的ESF曲线，进而阻碍了低频段MTF值无法准确获取。 

因此，探索一种有效的MTF测量方法，避免传统方法中由于求导过程导致的噪声放大，准确实现数字X射线成像系统的分辨率性能评估，已经成为成像系统性能评估的研究重点，这对于成像系统的发展应用具有十分重要的意义。 

发明内容

本发明的主旨是提出一种数字X射线成像系统的分辨率性能评估新方法，以此解决数字放射成像系统性能评估中的关键问题之一：实现数字X射线成像系统的调制传递函数低频特性的准确测量，为进一步全面评估放射成像系统分辨率性能提供有力条件。本发明针对传统刀口法中对过采样ESF曲线求导过程导致的噪声放大引起低频MTF测量精度下降问题，提出一种过采样ESF曲线的解卷积方法，避开获取LSF曲线的求导过程，从而可以获得较传统刀口测量方法更为准确的MTF曲线。该发明为正确获取系统调制传递函数特性，有效评估放射成像系统性能，深入开展放射影像学临床实践和研究提供有力支持。本发明的技术方案如下： 

一种数字X射线成像系统的分辨率性能评估新方法，包括下列步骤： 

①设置数字放射成像系统的曝光参数，放置刀口测试仪器，使其紧贴数字X射线成像系统的探测器，并使刀口与探测器采样方向有一倾斜角度，连续采集多幅刀口图像； 

②对所采集的多幅刀口图像叠加平均，而后对平均后的刀口图像进行刀口边界检测，并且利用直线拟合获得刀口边界直线； 

③对刀口边界直线图进行Hough变换，获得刀口倾斜角度α，而后依据公式N＝round(N_ave)＝round(1/tanα)得到该倾角对应的插值数N，从而构建出过采样ESF曲线g(x)，round符号代表取整函数； 

④计算过采样ESF曲线g(x)的信噪比，当零频率位置对称ESF曲线中噪声与信号频率的幅度谱比值大于0.5％时，采用单调曲线拟合方法对过采样曲线g(x)进行消噪； 

⑤通过对数字放射成像系统的模拟实验得到真实ESF曲线f(x)，根据探测器像素大小对应的分辨率，选择适当的平移距离T，构建真实ESF曲线f(x)与其平移T距离结果f(x-T)的差值曲线，即方波 对其傅里叶变换得到sinc函数T的取值需要保证sinc函数的过零点与所需MTF曲线的取值频率不重合； 

⑥对实际过采样ESF曲线g(x)进行适当的平移T距离得到g(x-T)，而后与原来的过采样曲线g(x)相减，得到对称的过采样曲线对和分别作傅里叶变换得到和

⑦利用频域下恶化模型公式求解系统调制传递函数MTF′(u)，而后 对MTF′(u)采用零频率位置MTF值进行归一化，即MTF(u)＝MTF′(u)/MTF′(u)|_u＝0，最终得到归一化调制传递函数MTF(u)。 

在进行数字X射线成像系统的调制传递函数测量时，通过本发明的基于系统恶化模型的过采样ESF曲线解卷积方法，避免了传统刀口法中对过采样ESF曲线求导过程导致的噪声放大引起低频MTF测量精度下降问题，从而可以获得较传统刀口测量方法更为准确的MTF曲线。该方法的应用，将为准确测量系统调制传递函数，有效评估放射成像系统性能，深入开展放射影像学临床实践和研究提供有力支持。 

附图说明

图1. 8条插值重构过采样ESF曲线。 

图2. 过采样ESF曲线(模拟刀口器件覆盖区域外像素最大值13000，刀口器件覆盖区域内像素最大值300)。 

图3. 图2过采样ESF曲线对应的归一化LSF曲线。 

图4. 对称过采样ESF曲线。 

图5. LSF曲线中信号和噪声的傅里叶幅度谱对比示意图。 

图6. 对称ESF曲线中信号和噪声的傅里叶幅度谱对比示意图。 

图7. 4种方法得到的MTF曲线与真实MTF曲线对比(刀口物体覆盖区域内外的像素最大值分别为300和13000)。 

图8. 4种方法得到的MTF曲线与真实MTF曲线误差对比(刀口物体覆盖区域内外的像素最大值分别为300和13000)。 

图9. 4种方法得到的MTF曲线与真实MTF曲线对比(刀口物体覆盖区域内外的像素最大值分别为150和6500)。 

图10. 4种方法得到的MTF曲线与真实MTF曲线误差对比(刀口物体覆盖区域内外的像素最大值分别为150和6500)。 

图11. 4种方法得到的实际探测器系统的MTF曲线对比。 

图12. 4种方法得到的实际探测器系统的MTF方差曲线对比。 

具体实施方式

下面结合附图和实施例对本发明进行说明。 

针对传统刀口MTF测量方法中对过采样ESF曲线求导过程导致的噪声放大引起低频MTF测量精度下降问题，本发明提出一种数字X射线成像系统的分辨率性能评估新方法，采用过采样ESF曲线的解卷积方法， 避开获取LSF曲线的求导过程，从而可以获得较传统刀口测量方法更为准确的MTF曲线。最后得到的技术方案如下：①设置数字放射成像系统的曝光参数，放置刀口测试仪器紧贴探测器，并使刀口与探测器采样方向有一倾斜角度，连续采集15幅图像；②对15幅刀口图像叠加平均，降低系统噪声，而后对平均后的刀口图像进行刀口边界检测，并且利用直线拟合获得刀口边界直线；③对刀口边界直线图进行Hough变换，获得刀口倾斜角度α，而后依据公式N＝round(N_ave)＝round(1/tanα)得到该倾角对应的插值数N，从而构建出过采样ESF曲线g(x)，round符号代表取整函数。④计算过采样ESF曲线g(x)的信噪比，当信噪比较低时，采用单调非减拟合方法对过采样曲线g(x)进行消噪。⑤根据探测器像素大小对应的分辨率，选择适当的平移距离T，构建理想ESF曲线f(x)与其平移T距离结果f(x-T)的差值曲线，即方波 对其傅里叶变换得到sinc函数T的取值需要保证sinc函数的过零点与所需MTF曲线的取值频率不重合。⑥对实际过采样ESF曲线g(x)进行适当的平移T距离得到g(x-T)，而后与原来的过采样曲线g(x)相减，得到对称的过采样曲线对进行傅里叶变换得到⑦利用恶化模型公式求解系统调制传递函数MTF′(f)，而后对MTF′(f)采用零频率位置MTF值进行归一化，即MTF(f)＝MTF′(f)/MTF′(f)|_f＝0，最终得到归一化调制传递函数MTF(f)。 

1ESF过采样曲线获取及处理 

利用过采样技术重构ESF，需将刀口与探测器矩阵成一定倾角摆放，通常将仪器边缘和数字化探测器阵列方向之间倾角α保持为1.5°～7.1°(对应插值数为8～38条)，α越小(插值数越高)所获MTF曲线的分辨率就越高，但相应的MTF曲线中的噪声也越大。 

采用过采样刀口测量方法获得的刀口示意图如图1所示。图中灰色部分均表示刀口器件覆盖区域，粗黑线表示经同轴准直和理想切割后的刀口边缘。图1中的刀口倾斜角度为7.1°，针对该倾角，要想得到精确的边缘响应函数ESF必须对刀口附近图像进行8条插值。插值数与倾斜角度的关系由下面的公式确定 

N_ave＝1/tanα  (1) 

由于实际刀口倾斜角度对应的插值数一般不是整数，因此之前的测量方法中一般取最接近插值数N_ave的整数N 

N＝round(N_ave)＝round(1/tanα)  (2) 

上式(2)中的round符号代表取整函数。 

采用8条插值来重构过采样ESF曲线可以参照图1中数字顺序，对应插值顺序从1到8，然后从9到16，依此类推完成整条ESF曲线的重构，本发明标记为f(x)。 

一般来说，调制传递函数低频域主要决定图像的对比度，高频域决定图像细节的重现能力及边缘处清晰度。在系统性能评估示，为了给出准确的低频或零频率DQE值，MTF低频段的准确测量尤为重要。 

2基于系统恶化模型的MTF测量方法 

从数字X射线成像系统获取图像的过程是一个包括模糊效应和加性噪声的系统恶化模型。模糊效应是由于探测器有限的分辨率引起的，在频域可以用调制传递函数(MTF)来表征。刀口物体x射线成像过程满足系统恶化模型，而对于从刀口图像中抽取的线扩散函数(ESF)，同样满足该恶化模型。数字X射线成像系统可以看作线性移不变系统。 

对于恶化的倾斜刀口图像，从刀口图像抽取得到的ESF曲线g(x)可以用以下模型表示： 

g(x)＝f(x)⊙h(x)+n(x)  (3) 

这里⊙表示卷积符号，x是过采样空间位置，f(x)是理想的没有模糊效应的ESF曲线，h(x)表征探测器线扩散函数(LSF)，即模糊效应的空间域表达。n(x)代表系统的加性噪声。 

在频域，上述恶化模型，即频率u下的实际过采样ESF曲线G(u)，理想ESF曲线F(u)，以及系统的加性噪声N(u)可以表示为： 

G(u)＝F(u)·MTF(u)+N(u)  (4) 

基于上述的恶化模型，在倾斜刀口图像具备高信噪比的前提下，可以通过频域解卷积方法来获取系统的MTF： 

MTF(u)＝{G(u)-N(u)}/F(u)  (5) 

但是，严格来说，理想ESFf(x)的傅里叶变换F(u)是不存在的，因为f(x)是一个理想的阶跃函数，无法满足傅里叶变换的绝对可积条件。因此，为了满足频域恶化模型的要求，需要对上述建立的恶化模型进行改进。 

基于线性移不变系统理论，移不变表示如果我们给系统施加一个输入，无论该输入是在当前位置，还是在T毫米以后的位置，系统输出是固定不变的，只是输出也平移了对应的T毫米。也就是说，如果系统输出g(x)是系统输入f(x)的响应，那么对于输入f(x-T)所对应的系统输出可以表示为： 

g(x-T)＝f(x-T)⊙h(x)+n₂(x)  (6) 

这里采用一个新的符号n₂(x)来表示噪声，因为考虑到噪声的随机过程特性。 

进一步的，系统的线性表示系统的输入和输出之间满足一个线性映射。即线性加权的输入必然对应线性加权的输出。如果我们对系统输入[f(x)-f(x-T)]，那么对于该线性加权输入的系统输出为： 

[g(x)-g(x-T)]＝[f(x)-f(x-T)]⊙h(x)+[n(x)-n₂(x)]  (7) 

这里 $\hat{g}\left(x\right)=g\left(x\right)-g(x-T),\hat{f}\left(x\right)=f\left(x\right)-f(x-T),\hat{n}\left(x\right)=n\left(x\right)-n_2\left(x\right).$

这时，线性加权的输入是一个方波函数，它能满足傅里叶变换的绝对可积条件。因此，系统MTF可以用下面的公式获取： 

${\mathrm{MTF}}^{\′}\left(u\right)=\{\hat{G}\left(u\right)-\hat{N}\left(u\right)\}/\hat{F}\left(u\right)---\left(9\right)$

这里和分别是和的傅里叶变换。对MTF′(f)采用零频率位置MTF值进行归一化，即MTF(f)＝MTF′(f)/MTF′(f)/_f＝0，最终得到归一化调制传递函数MTF(f)。 

但是，通过上述模型求取MTF的问题仍然没有被完全解决，这是由的特性决定的。对方波函数的傅里叶变换是sinc函数，该函数是周期性的，具有均匀分布的过零点。因此，在sinc函数过零点频率位置的MTF值无法直接通过除以分母的sinc函数来获得(分母不能为零)。在实际中，对于数字X射线成像系统的性能评估，研究人员通常采用离散傅里叶变换来获得在采样频率范围内频域均匀采样的MTF值。因为sinc函数的过零点分布取决于方波函数的宽度，通过采用适当的输入位置平移T，可以使过零点分布偏移所需要的MTF值的频率位置，以此解决上述问题。 

针对所建立的恶化模型公式(9)，分子项表示消噪过程。这一过程需要特别注意，必须权衡消噪后信号的信噪比与消噪引起的信号失真，不恰当的消噪会导致MTF测量误差超出其测量准确性要求的范围。在实际情况下，当过采样ESF曲线具有充分高的信噪比，那么可以忽略模型中的噪声项。同时为了作一比较，本发明也采用了单调非减曲线拟合方法对过采样ESF曲线进行消噪，而后将消噪后的ESF曲线带入系统恶化模型，求解系统MTF。 

本发明将对原始获得的过采样ESF曲线以及的单调非减曲线拟合后的ESF曲线通过公式(9)的系统恶化模型来获取系统MTF，同时将传统方法采用的对ESF求导并傅里叶变换获取MTF的结果给出，作为比较，以此证明本发明的优势。 

本发明分别采用模拟刀口图像和实际刀口图像来验证本发明在MTF测量准确性及不确定性方面的性 能。模拟刀口图像采用目前本领域普遍认可的指数模型刀口图像，即Lorentzian型函数： 

$\mathrm{MTF}\left(u\right)=\frac{r^2}{r^2+{(2\mathrm{\π}\·u)}^2}---\left(10\right)$

这里的r值选取探测器像素的倒数，u为空间频率。具体可以参考Esamei等人发表的MTF测量相关论文。本发明模拟图像采用2度倾角，像素大小为0.1毫米。为了模拟不同X射线曝光剂量对应的信噪比，模拟图像中刀口物体覆盖区域外的像素最大值分别设置为13000和6500，与之对应的，刀口物体覆盖区域内的像素最大值分别设置为300和150。模拟刀口图像中加入泊松噪声，图像像素点的方差等于该点位置的像素值的平方根。 

本发明采用的实际刀口图像采集系统为美国BIOPTICS公司生产的Pixarray100小动物数字放射成像系统。按照国际电气技术委员会(IEC)制定的标准，经微纳设备加工后得到刀口法测试仪器，测试仪器由黄铜板刀口与有机玻璃支架组装而成，刀口测试材料安置固定在支架上。实验中需将精细加工后刀口放置于视场中心，并使刀口与探测器采样方向有一倾斜角度，通过Hough变换即可获得刀口的倾斜角度。设置数字放射成像系统的曝光参数，将刀口测试仪器紧贴探测器，连续采集15幅图像，用于分析本发明方法的不确定性。 

模拟结果如下： 

首先模拟X射线曝光剂量对应刀口物体覆盖区域外像素最大值13000和刀口物体覆盖区域内的像素最大值300的情况，这时候刀口覆盖区域外的信噪比为41.1dB，而刀口覆盖区域内的信噪比为24.5dB。从刀口图像获得的过采样ESF曲线如图1所示，对其求导并归一化后得到对应的LSF如图2所示。 

选取一个适当的平移距离，对过采样ESF曲线平移，而后将与原始的过采样ESF曲线相减，可以得到如图4所示的对称过采样ESF曲线。根据系统恶化模型公式(9)，该对称过采样ESF曲线为对应宽度的方波输入的系统响应。 

为了观察求导过程对信号信噪比的影响。图5和图6分别为LSF和对称ESF曲线中信号和噪声的傅里叶幅度谱的对比示意图。在零频率位置，LSF曲线中噪声与信号频率的幅度谱比值为6.10％，而相应的在对称ESF曲线中噪声与信号频率的幅度谱比值为0.12％。由于低频段对称ESF曲线的信噪比足够高，那么在求解MTF时忽略系统恶化模型中的噪声项是合理的。本发明同时给出了采用单调非减曲线拟合方法对ESF曲线消噪，而后再利用系统恶化模型求解MTF。 

图7给出了采用4种不同的方法获取MTF的结果对比，4种方法包括：本发明提出的基于系统恶化模型的未消噪ESF曲线解卷积方法，本发明提出的基于系统恶化模型的单调非减拟合ESF曲线解卷积方法，对未消噪ESF曲线求导后归一化傅里叶变换的传统方法，对单调非减拟合后ESF曲线求导后归一化傅里叶变换的传统方法。图7同时给出了模拟刀口图像的真实MTF。4种方法获取的MTF与真实MTF的差异如图8所示。在低频 范围0-1mm¹范围，基于系统恶化模型的未消噪ESF曲线解卷积方法获取MTF的平均误差为0.03％，优于基于系统恶化模型的单调非减拟合ESF曲线解卷积方法(平均误差为0.06％)，也优于对未消噪ESF曲线求导后归一化傅里叶变换的传统方法(平均误差为0.19％)以及对单调非减拟合后ESF曲线求导后归一化傅里叶变换的传统方法(平均误差为0.09％)。如果考虑采样频率内的全频范围，本发明提出的基于系统恶化模型的单调非减拟合ESF曲线解卷积方法获取的MTF平均误差最小。 

为了分析本发明在不同x射线曝光剂量的适应性，我们又模拟了X射线曝光剂量对应刀口物体覆盖区域外像素最大值6500和刀口物体覆盖区域内的像素最大值150的情况，这时候刀口覆盖区域外的信噪比为38.1dB，而刀口覆盖区域内的信噪比为21.8dB。图9给出了采用4种不同的方法获取MTF的结果对比，同时给出了模拟刀口图像的真实MTF。4种方法获取的MTF与真实MTF的差异如图10所示。在低频范围0-1mm¹范围，基于系统恶化模型的未消噪ESF曲线解卷积方法获取MTF的平均误差为0.07％，优于基于系统恶化模型的单调非减拟合ESF曲线解卷积方法(平均误差为0.09％)，也优于对未消噪ESF曲线求导后归一化傅里叶变换的传统方法(平均误差为0.59％)以及对单调非减拟合后ESF曲线求导后归一化傅里叶变换的传统方法(平均误差为0.15％)。如果考虑采样频率内的全频范围，本发明提出的基于系统恶化模型的单调非减拟合ESF曲线解卷积方法获取的MTF仍然具有最小的MTF测量平均误差。 

通过已知MTF的模拟图像对本发明方法对低频MTF测量准确性进行了评估，为了进一步评估本发明方法在低频MTF测量不确定性方面的性能，我们采用实际刀口图进行MTF测量对比。由于实际MTF日常评估时，非常需要从单幅刀口图获取数字x射线影像系统的MTF，因此良好的MTF测量的不确定性十分重要。国际电工委员会为MTF测量算法的不确定性制定了标准。本发明采用15幅在同一曝光条件下获得的刀口图像来进行MTF测量算法不确定性的对比。 

图11给出了采用4种方法从15幅刀口图像测量得到的pixarray-100小动物发射影像系统的平均归一化MTF。4种方法得到的归一化MTF的相对方差曲线如图11所示。在低频范围0-2mm¹范围，基于系统恶化模型的未消噪ESF曲线解卷积方法获取MTF的平均相对方差为1.31％，优于基于系统恶化模型的单调非减拟合ESF曲线解卷积方法(平均相对方差为1.86％)，也优于对未消噪ESF曲线求导后归一化傅里叶变换的传统方法(平均相对方差为2.94％)以及对单调非减拟合后ESF曲线求导后归一化傅里叶变换的传统方法(平均相对方差为2.46％)。如果考虑采样频率内的全频范围，本发明提出的基于系统恶化模型的单调非减拟合ESF曲线解卷积方法获取的MTF具有最小的MTF平均相对方差。 

最终结果表明，在进行数字X射线成像系统的调制传递函数测量时，通过本发明的基于系统恶化模型的过采样ESF曲线解卷积方法，避免了传统刀口法中对过采样ESF曲线求导过程导致的噪声放大引起低频MTF测量精度下降问题，从而可以获得较传统刀口测量方法更为准确的MTF曲线。该方法的应用，将为准确测量系统调制传递函数，有效评估放射成像系统性能，深入开展放射影像学临床实践和研究提供有力支持。 

Claims (1)

1.一种数字X射线成像系统的分辨率性能评估新方法，包括下列步骤：

①设置数字放射成像系统的曝光参数，放置刀口测试仪器，使其紧贴数字X射线成像系统的探测器，并使刀口与探测器采样方向有一倾斜角度，连续采集多幅刀口图像。

②对所采集的多幅刀口图像叠加平均，而后对平均后的刀口图像进行刀口边界检测，并且利用直线拟合获得刀口边界直线；

③对刀口边界直线图进行Hough变换，获得刀口倾斜角度α，而后依据公式N＝round(N_ave)＝round(1/tanα)得到该倾角对应的插值数N，从而构建出过采样ESF曲线g(x)，round符号代表取整函数；

④计算过采样ESF曲线g(x)的信噪比，当零频率位置对称ESF曲线中噪声与信号频率的幅度谱比值大于0.5％时，采用单调曲线拟合方法对过采样曲线g(x)进行消噪；

⑤通过对数字放射成像系统的模拟实验得到真实ESF曲线f(x)，根据探测器像素大小对应的分辨率，选择适当的平移距离T，构建真实ESF曲线f(x)与其平移T距离结果f(x-T)的差值曲线，即方波对其傅里叶变换得到sinc函数T的取值需要保证sinc函数的过零点与所需MTF曲线的取值频率不重合；

⑥对实际过采样ESF曲线g(x)进行适当的平移T距离得到g(x-T)，而后与原来的过采样曲线g(x)相减，得到对称的过采样曲线对和分别作傅里叶变换得到和

⑦利用频域下恶化模型公式求解系统调制传递函数MTF′(u)，而后对MTF′(u)采用零频率位置MTF值进行归一化，即MTF(u)＝MTF′(u)/MTF′(u)|_u＝0，最终得到归一化调制传递函数MTF(u)。