CN103528840A - 基于x射线成像系统探测器特性的调制传递函数测量方法 - Google Patents

Abstract

本发明属于生物医学工程及计算机领域，涉及一种基于X射线成像系统探测器特性的调制传递函数测量方法，包括：放置刀口测试仪器紧贴探测器，并使刀口与探测器采样方向有一倾斜角度，连续采集多幅图像；对多幅刀口图像叠加平均，降低系统噪声，而后对平均后的刀口图像进行刀口边界检测，获得刀口边界直线；对刀口边界直线图进行Hough变换，获得刀口倾斜角度α，构建过采样ESF曲线；采用C样条曲线作为回归样条曲线f对过采样ESF曲线进行分段拟合；进行差分运算得到线扩散函数；对线扩散函数进行傅里叶变换再取模，得到调制传递函数；归一化处理。本发明解决了刀口法获取过采样ESF曲线时的拟合误差问题，可以获得较传统刀口测量方法更为准确的MTF曲线。

Description

基于X射线成像系统探测器特性的调制传递函数测量方法

所属技术领域

本发明属于生物医学工程及计算机领域，涉及一种数字X射线成像系统的调制传递函数精确测量新方法。

背景技术

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

要定量地评价数字X射线成像系统的固有成像质量，只需计算不受个人主观因素影响的系统固有预采样MTF。实际中系统MTF常由以下三种扩散函数来计算：点扩散函数、线扩散函数和边缘响应函数，它们分别描述经成像系统后点、线和边缘弥散程度，能够间接地反映系统成像能力。

目前测量MTF的方法中，线对卡方法可以获得高精度的MTF值，但是线对卡方法只能提供有限几个整数空间频率位置的MTF值，无法实现系统传递函数的全面评估。为了获取成像系统空间频率范围下的MTF变化曲线，一般常用的测量方法包括狭缝法和刀口法，狭缝法(Slit Camera)和刀口法(Edge)已被国际放射学会公认为是获得MTF的较好方法，日本将狭缝法定义成测量MTF的标准方法，刀口法也已被国际电气技术委员会(IEC)指定为测量系统MTF的标准方法。将通过狭缝法和刀口法所获得的同一系统下不同MTF曲线进行比较时，可知前者在高频域有较高信噪比，而后者在低频域有较高信噪比。通过狭缝法获得的系统MTF精确，且操作简便、方法成熟。但因其加工难度高（狭缝宽度小于等于10微米，误差在1微米内），导致在实际应用中此法较难被推广。由于刀口法测量仪器自加工相对容易，在科研实验和常规检测中使用较为广泛。

采用刀口法可以获得金属模块刀口的边缘扩散函数(Edge Spread Function，ESF)，它反映了经成像系统后边缘弥散程度，ESF的导数为线扩展函数(Line Spread Function，LSF)，再经过傅立叶变换即可获得系统MTF，目前大部分关于MTF的研究都是基于刀口法技术展开的。

为了获得高分辨率的MTF曲线，目前一般采用改进的MTF刀口测量方法，即将刀口方向与图像采样方向成一定角度，以此获得过采样的边缘响应函数ESF，而后通过傅里叶变换得到高分辨MTF曲线。由于LSF曲线的获取是通过对ESF求导实现的，ESF曲线中的噪声被求导过程放大，阻碍了LSF曲线的准确提取。为了克服这个困难，目前已有学者期望通过对ESF曲线高斯平滑、多阶分段曲线拟合、单调性曲线拟合等方法来提高LSF曲线提取的准确性。但到目前为止，由于上述方法没有充分利用探测器本身性质的先验知识，在MTF评估精度和鲁棒性方面仍存在较大的缺欠。

因此，探索一种有效的MTF刀口测量方法，充分利用探测器本身的先验知识，提高刀口法MTF测量精度，准确实现数字X射线成像系统的调制传递函数特性评估，已经成为成像系统性能评估的研究重点，这对于成像系统的发展应用具有十分重要的意义。

发明内容

本发明的主旨是克服现有技术的上述不足，提供一种数字X射线成像系统的调制传递函数测量新方法，以此解决数字放射成像系统性能评估中的基本问题：实现数字X射线成像系统的调制传递函数的准确测量，为进一步全面评估放射成像系统性能提供有力条件。本发明的技术方案如下：

一种基于X射线成像系统探测器特性的调制传递函数测量方法，包括下列步骤：

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

②对多幅刀口图像叠加平均，降低系统噪声，而后对平均后的刀口图像进行刀口边界检测，并且利用直线拟合获得刀口边界直线；

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

④采用C样条曲线作为回归样条曲线f对ESF曲线进行分段最小二乘拟合，同时满足约束条件，即拟合时保证ESF曲线左半部分为单调非减凹函数，右半部分为单调非减凸函数；

⑤对经过最小二乘拟合的ESF曲线进行差分运算得到线扩散函数LSF(x)；

⑥对线扩散函数LSF(x)进行傅里叶变换再取模，得到调制传递函数MTF′(f)；

⑦对MTF′(f)采用零频率位置MTF值进行归一化，得到归一化的调制传递函数MTF(f)。

作为优选实施方式，步骤④中，采用3阶C样条曲线作为的回归样条曲线f，对ESF曲线进行分段拟合，所满足的约束条件为：f′(x)≥0， $\left\{\begin{array}{cc}{f}^{\″}\left(x\right)\≥0,& x\∈[x_{\min },x_{\mathrm{mid}}]\\ f^{\″}\left(x\right)\≥0,& x\∈[x_{\mathrm{mid}},x_{\max }]\end{array}\right.,$ 这里f′(x)和f′′(x)分别为回归样条曲线f的一阶导数和二阶导数，x_min，x_mid和x_max分别是待拟合ESF曲线的左端点，中心边界和右端点，λ为平滑因子，λ≥0，用于控制曲线的平滑度与误差之间的平衡。

在进行数字X射线成像系统的调制传递函数测量时，通过本发明的基于X射线成像系统探测器特性的调制传递函数测量新方法，解决了刀口法获取过采样ESF曲线时的拟合误差问题，可以获得较传统刀口测量方法更为准确的MTF曲线。该方法的应用，将为准确测量系统调制传递函数，有效评估放射成像系统性能提供有力支持。

附图说明

图1.8条插值重构过采样ESF曲线。

图2.刀口法测试仪器的X光成像结果平均图。

图3.刀口法测试仪器的刀口边界直线拟合结果。

图4.本发明方法与目前3种方法得到的LSF曲线对比。

图5.本发明方法与目前3种方法得到的15条LSF曲线的方差对比。

图6.本发明方法与目前3种方法得到的MTF曲线对比。

图7本发明方法与目前3种方法得到的MTF方差曲线对比。

图8.4种方法得到的LSF曲线与真实LSF曲线对比。

图9.4种方法得到的MTF曲线与真实MTF曲线对比。

图10.4种方法得到的MTF误差曲线。

具体实施方式

针对传统刀口法获取过采样ESF曲线时的拟合误差问题，本发明提出基于X射线成像系统探测器特性的调制传递函数测量新方法，通过融合探测器特性先验知识的C样条曲线拟合方法，以此构建准确的过采样ESF曲线，最终可以获得较传统刀口测量方法更为准确的MTF曲线。下面对本发明进行详细说明。

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

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

对于刀口倾斜角度为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曲线的重构，即得到ESF(x)。

想从边缘响应函数ESF(x)运算得到线扩散函数LSF(x)，对ESF(x)进行差分运算或使用卷积滤波器[-0.5，0，0.5]均可，即

LSF(x)＝dESF(x)/dx      (3)

LSF(x)函数经过傅里叶变换后再取模，就能得到其调制传递函数

MTF′(f)＝|FT(LSF(x))|     (4)

式中FT为傅里叶变换，通常在描述系统调制传递函数时使用归一化处理的结果，其公式为

$\mathrm{MTF}\left(f\right)=\frac{{\mathrm{MTF}}^{\′}\left(f\right)}{{\mathrm{MTF}}^{\′}\left(f\right){|}_{f=0}}---\left(5\right)$

一般来说，调制传递函数低频域主要决定图像的对比度，高频域决定图像细节的重现能力及边缘处清晰度。从MTF曲线的包络面积可判断系统成像质量优劣，通常希望成像系统MTF曲线包括的面积越大越好。

由于MTF测量精度受到过采样ESF曲线中的噪声的制约，需要对刀口法获得的过采样ESF曲线进行处理，目前已有学者提出了几种处理方法。Samei提出了一种高斯加权平滑的方法降低ESF曲线的噪声，之后Samei和Saunders又提出了四阶移动曲线拟合方法提高了MTF高频部分的精度。Maidment and Albert提出了利用ESF曲线自身的单调性来实现曲线拟合和降噪的方法，取得了一定的效果。但到目前为止，由于上述方法没有充分利用探测器本身性质的先验知识，在MTF评估精度和鲁棒性方面还有较大的缺欠。本发明为了提高MTF测量精度，提出了基于X射线成像系统探测器特性的调制传递函数测量新方法，通过3阶C样条分段拟合，融合探测器特性的先验知识，得到准确的ESF拟合结果，而后计算获得MTF曲线。

2基于探测器先验知识融合的ESF曲线拟合

为了提取探测器特性的先验知识，可以从MTF的两种测量方法考虑，即刀口法和狭缝法。通过对刀口设备成像，从刀口设备的黄铜板覆盖区域到空气区域成像区，X射线透射强度逐渐增强，因此，最终得到的刀口图像的横截面曲线（ESF曲线）是一个单调非减曲线。再考察狭缝法，通过对狭缝设备成像，可以得到狭缝图像的横截面为线扩散曲线(LSF)，由于狭缝设备从左到右依次为黄铜板覆盖区域，空气狭缝，黄铜板覆盖区域，因此，X射线透射强度在空气狭缝中央最强，而后向两边逐渐降低，这就形成了LSF曲线的特性：LSF曲线左半部分是单调非减曲线，而其右半部分是单调非增曲线。由于LSF曲线是ESF曲线的一阶导数结果，因此，可以从LSF曲线推导出相应的ESF曲线性质，即以刀口边界为中心，ESF曲线左半部分是凸曲线，而其右半部分是凹曲线。由此得到，与探测器特性直接相关的ESF曲线是由左半部分单调非减的凸曲线和右半部分单调非减的凹曲线构成。

为了将上述探测器先验知识融合到MTF测量方法中，我们采用了约束样条回归方法来实现这一目标。约束样条回归是一种非参数回归方法，它提供了回归函数的有效估计而不需要给出假设参数形式。该方法不仅具备充分的适应性，并且允许局部效应变化，提供了一种融合形状约束的有效框架。

首先给出过采样ESF曲线模型，曲线上在xi位置的点yi可以用如下表达式表达：

y_i＝f(x_i)+ε_i     (5)

曲线位置坐标满足x₁＜x₂＜…x_n，ε_i为测量ESF曲线与真实ESF曲线之间的误差。对于包含较多噪声的ESF曲线，同时考虑形状约束及曲线平滑可以得到更理想的ESF曲线拟合结果，即在回归样条曲线f下最小化下列表达式

$\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{[y_i-f\left(x_i\right)]}^2+\mathrm{\λ}{\∫}_a^b{\left[f^{\left(k\right)}\left(x\right)\right]}^2\mathrm{dx}---\left(6\right)$

同时满足约束条件：f′(x)≥0， $\left\{\begin{array}{cc}{f}^{\″}\left(x\right)\≥0,& x\∈[x_{\min },x_{\mathrm{mid}}]\\ f^{\″}\left(x\right)\≥0,& x\∈[x_{\mathrm{mid}},x_{\max }]\end{array}\right.,$ 这里f′(x)和f′′(x)分别为回归样条曲线f的一阶导数和二阶导数，x_min，x_mid和x_max分别是待拟合ESF曲线的左端点，中心边界和右端点。

上式(6)中的λ为平滑因子，λ≥0，用于控制曲线的平滑度与误差之间的平衡。本发明采用交叉验证方法来获得最优平滑因子。参数k表示样条函数的阶数。考虑到三阶样条具有很好的适应性，并且能够拟合曲线中的快速上升段，因此我们采用3阶样条来实现公式(6)的最小化。

回归样条曲线可以采用几类不同形式的基函数来构建。B样条是一种常用的样条曲线基函数，具备构建的直接性和递归定义。其他基函数包括截断功率函数、M样条、C样条等。本发明为了将ESF曲线的单调性和凹凸性融入MTF测量方法中，我们选择采用C样条作为基函数来构建回归样条曲线。C样条是M样条曲线的二次积分，首先给出M样条曲线表达式。

对于k阶的回归样条曲线，给定待拟合曲线的节点位置为x₁＝t₁＝…t_k＜…＜t_l+k＝…＝t_l+2k＝x_n，在给定曲线段t_i≤x≤t_i+k内，一阶M样条曲线的表达式如下：

从一阶M样条可以获得k阶M样条曲线，表达式如下

通过对M样条曲线函数积分，可以得到k阶I样条曲线

i＝1，…，l+k，x∈[x₁，x_n]     (9)

对k阶I样条曲线积分，得到k阶C样条曲线

i＝1，…，l+k，x∈[x₁，x_n]     (10)

本发明取3阶C样条曲线作为公式(6)中的回归样条曲线f，对ESF曲线进行分段拟合，从而获得高分辨率的精确ESF曲线，由此降低系统MTF评估的误差。

上述过程可以总结为:采用C样条曲线作为回归样条曲线f对过采样ESF曲线进行分段最小二乘拟合，同时满足约束条件，即拟合时保证ESF曲线左半部分为单调非减凹函数，右半部分为单调非减凸函数。

3基于X射线成像系统探测器特性的调制传递函数测量新方法流程

本发明属于生物医学工程及计算机领域，涉及一种数字X射线成像系统的调制传递函数精确测量新方法，其技术流程如下：

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

②对15幅刀口图像叠加平均，降低系统噪声，而后对平均后的刀口图像进行刀口边界检测，并且利用直线拟合获得刀口边界直线；

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

④采用C样条曲线作为回归样条曲线f对过采样ESF曲线进行分段最小二乘拟合，同时满足约束条件，即拟合时保证ESF曲线左半部分为单调非减凹函数，右半部分为单调非减凸函数。

⑤根据公式LSF(x)＝dERF(x)/dx，对ESF曲线进行差分运算得到线扩散函数LSF(x)；

⑥对线扩散函数LSF(x)进行傅里叶变换再取模MTF′(f)＝|FT(LSF(x))|，得到其调制传递函数MTF′(f)，而后对MTF′(f)采用零频率位置MTF值进行归一化，即MTF(f)＝MTF′(f)/MTF′(f)|_f＝0，最终得到归一化调制传递函数MTF(f)。

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

设置数字放射成像系统的曝光参数，将刀口测试仪器紧贴探测器，连续采集15幅图像，对这些图像进行平均，可以获得图2所示的刀口图像。检测刀口边界，并且利用直线拟合获得刀口边界直线，如图3所示，而后通过Hough变换即可获得刀口的倾斜角度。

根据刀口倾斜角度，我们可以获得传统刀口测量方法的整数插值数，再根据5.1节所述方法获得过采样ESF原始曲线，而后采样本发明方法对原始ESF曲线进行处理，为了说明本发明方法的有益结果说明，我们给出了本发明方法与目前采用的三种方法的结果对比。

图4给出了采用高斯加权平滑、4阶移动曲线拟合、单调非减拟合以及本发明的基于X射线成像系统探测器特性的C样条曲线拟合方法得到的15条LSF曲线的平均曲线结果。

为了对LSF结果准确评估，我们给出了4种方法的LSF方差曲线，结果如图5所示，可以看出，本发明方法获得的LSF曲线具有最小的方差，即再4种方法中具有最佳的LSF提取性能。

图6给出了4种方法对应的平均归一化MTF结果。在截止频率以内，4种方法的相对误差小于0.5%，在采样频率以内，4种方法的相对误差小于0.9%。这也说明了目前采样的各种方法在MTF测量中是有效的，差别仅存在与测量精度的高低。

为了对4种方法获得的MTF曲线进行评估对比，图7给出了各方法得到的MTF曲线相对方差曲线。在采样频率0-10mm-1范围内，高斯加权平滑方法的平均相对方差为2.01%，4阶移动曲线拟合方法的平均相对方差为1.65%，单调非减拟合的平均相对方差为4.29%，而本发明方法的平均相对方差仅为0.83%。由此可见本发明方法可以降低MTF测量结果的不确定性（降低MTF曲线方差）。

为了进一步评估本发明方法在MTF测量的准确性能，我们对已知MTF的模拟刀口图像进行了测试，模拟刀口图像采用目前本领域普遍认可的指数模型刀口图像，即Lorentzian型函数：

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

这里的r值选取探测器像素的倒数，u为空间频率。具体可以参考Esamei等人发表的MTF测量相关论文。

图8给出了通过4种方法获得的LSF曲线与真实LSF曲线的对比图。图9给出相应的4种方法获得的MTF曲线与真实MTF曲线的对比图，在零频率到截止频率范围内，采用高斯加权方法获得的MTF测量结果与真实MTF间平均误差为0.65%±0.40%，采用4阶移动平滑曲线拟合，获取MTF的平均误差为0.19%±0.14%，而采用单调拟合方法，MTF误差为0.16%±0.10%。采用本发明的基于探测器特性先验知识融合的C样条曲线拟合方法，MTF误差为0.12%±0.10%。具体的误差曲线如图10所示，可以看出，本发明方法在MTF测量精度上较目前采用的其他3种方法更具有优势。

最终结果表明，在进行数字X射线成像系统的调制传递函数测量时，通过本发明的基于X射线成像系统探测器特性的调制传递函数测量新方法，充分利用了探测器ESF曲线的单调性和凹凸性，可以获得较传统刀口测量方法更为准确的MTF曲线。该方法的应用，将为准确测量系统调制传递函数，有效评估放射成像系统性能提供有力支持。

Claims (2)

1.一种基于X射线成像系统探测器特性的调制传递函数测量方法，包括下列步骤：

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

②对多幅刀口图像叠加平均，降低系统噪声，而后对平均后的刀口图像进行刀口边界检测，并且利用直线拟合获得刀口边界直线；

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

④采用C样条曲线作为回归样条曲线f对ESF曲线进行分段最小二乘拟合，同时满足约束条件，即拟合时保证ESF曲线左半部分为单调非减凹函数，右半部分为单调非减凸函数；

⑤对经过最小二乘拟合的ESF曲线进行差分运算得到线扩散函数LSF(x)；

⑥对线扩散函数LSF(x)进行傅里叶变换再取模，得到调制传递函数MTF′(f)；

⑦对MTF′(f)采用零频率位置MTF值进行归一化，得到归一化的调制传递函数MTF(f)。

2.根据权利要求1所述的测量方法，其特征在于，步骤④中，采用3阶C样条曲线作为的回归样条曲线f，对ESF曲线进行分段拟合，所满足的约束条件为：f′(x)≥0， $\left\{\begin{array}{cc}{f}^{\″}\left(x\right)\≥0,& x\∈[x_{\min },x_{\mathrm{mid}}]\\ f^{\″}\left(x\right)\≥0,& x\∈[x_{\mathrm{mid}},x_{\max }]\end{array}\right.,$ 这里f′(x)和f′′(x)分别为回归样条曲线f的一阶导数和二阶导数，x_min，x_mid和x_max分别是待拟合ESF曲线的左端点，中心边界和右端点，

λ为平滑因子，λ≥0，用于控制曲线的平滑度与误差之间的平衡。

Images