CN102809494A - 数字x射线成像系统调制传递函数的刀口法测量方法 - Google Patents

Abstract

本发明属于生物医学工程及计算机领域，涉及一种数字X射线成像系统调制传递函数的刀口法测量方法，包括：采集多幅刀口图像；进行刀口边界检测，获得刀口边界直线；对刀口边界直线图进行Hough变换，获得刀口倾斜角度，计算刀口边界附近各点与刀口直线的距离，按照距离从小到大排序，构成非均匀采样的ERF曲线序列；采取基于重心拉格朗日插值方法，获得均匀采样插值位置的ERF新序列；得到线扩散函数；得到归一化调制传递函数。本发明构建准确的过采样ERF曲线，最终可以获得较传统刀口测量方法更为准确的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 Response Function,ERF)，它反映了经成像系统后边缘弥散程度，ERF的导数为线扩展函数(Line Spread Function,LSF)，再经过傅立叶变换即可获得系统MTF，目前大部分关于MTF的研究都是基于刀口法技术展开的。

为了获得高分辨率的MTF曲线，目前一般采用改进的MTF刀口测量方法，即将刀口方向与图像采样方向成一定角度，以此获得过采样的边缘响应函数ERF，而后通过傅里叶变换得到高分辨MTF曲线。过采样过程可以通过插值法实现，插值点数由刀口倾斜角度决定，但是由于实际像素点数只能是整数值，而通过刀口倾斜角度计算得到的插值数一般都是非整数。如果希望摆放合适的刀口倾斜角度，以此获取整数插值数，必须要采用极其严苛的精确定位工作，操作极其困难。

因此，针对刀口倾斜角度造成的插值数误差问题，探索一种有效的MTF测量方法，提高刀口法MTF测量精度，准确实现数字X射线成像系统的调制传递函数特性评估，已经成为成像系统性能评估的研究重点，这对于成像系统的发展应用具有十分重要的意义。

发明内容

针对传统刀口法重构过采样ERF曲线时非均匀采样的插值位置误差的问题，本发明提出一种数字X射线成像系统的调制传递函数精确测量新方法。本发明基于重心拉格朗日插值的刀口法MTF测量方法，通过重心拉格朗日插值获取均匀采样位置的真实像素值，以此构建准确的过采样ERF曲线，最终可以获得较传统刀口测量方法更为准确的MTF曲线。

本发明的技术方案如下：

一种数字X射线成像系统调制传递函数的刀口法测量方法，包括下列步骤：

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

2)对刀口图像叠加平均，对平均后的刀口图像进行刀口边界检测，并且利用直线拟合获得刀口边界直线；

3)设定横坐标x以j表示，代表列数，纵坐标y以i表示，代表行数，像素间隔以p表示，刀口倾斜角度与探测器采样方向的夹角为α，对刀口边界直线图进行Hough变换，获得刀口倾斜角度，依据公式Δd_i，j=jp cosα-ip sinα计算刀口边界附近各点与刀口直线的距离，按照距离从小到大排序，构成非均匀采样的ERF曲线序列；

4)采取基于重心拉格朗日插值方法，获得均匀采样插值位置的ERF新序列；

5)对ERF新序列进行差分运算得到线扩散函数LSF(x)；

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

其中，第4）步中，可依据公式 $L\left(x\right)={\mathrm{\Σ}}_{j=0}^k\frac{{\mathrm{\ω}}_j}{x-x_j}{y}_j/{\mathrm{\Σ}}_{j=0}^k\frac{{\mathrm{\ω}}_j}{x-x_j},$ 计算一系列均匀采样位置x的像素值L(x),获得均匀采样插值位置的ERF新序列，ω_j是重心权重， ${\mathrm{\ω}}_j=\frac{1}{{\mathrm{\Π}}_{i=0,i\≠j}^k{x}_j-x_i}.$

本发明针对传统刀口法重构过采样ERF曲线时非均匀采样的插值位置误差，采用基于重心拉格朗日插值，获得均匀采样位置的真实像素值，从而可以获得较传统刀口测量方法更为准确的MTF曲线。该发明为正确获取系统调制传递函数特性，有效评估放射成像系统性能，深入开展放射影像学临床实践和研究提供有力支持。

附图说明

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

图2.传统刀口测量方法采用的误差插值位置与实际插值位置示意图。

图3.刀口图像中刀口附近像素与刀口边界相对位置示意图。

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

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

图6.传统刀口法与本发明方法获得的ERF曲线对比。

图7.传统刀口法与本发明方法获得的LSF曲线对比。

图8.传统刀口法、本发明方法与线对卡方法获得的MTF结果对比。

具体实施方式

下面结合附图和实施例从几个方面对本发明进行详细说明。

1 MTF过采样刀口测量方法

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

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

N_ave=1/tanα                                            (1)

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

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

LSF(x)=dERF(x)/dx                                       (2)

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

MTF＇(f)=|FT(LSF(x))|                                   (3)

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

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

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

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

N=round(N_ave)=round(1/tanα)                            (5)

由于N_ave是非整数，导致目前传统刀口测量方法采用的插值位置总是会偏离实际的真实插值位置，我们以图2进行说明。图中p为每两个像素的间隔，Δx为当前行像素与下一行对应像素距离刀口距离的差值，由下式决定

Δx=p tanα                                           (6)

传统刀口测量方法采用的插值位置是一些近似插值位置，如图2所示，这些插值点位置是均匀采样的，即传统刀口测量法中认为刀口图像中的相邻像素点到刀口距离的差值相等，则这些插值点序列构成均匀采样。但实际上，由于刀口倾斜角度对应的插值数N_ave为非整数，刀口图像中的像素点实际对应图2中非均匀的实际插值位置（实际插值位置）。因此，为了精确获取MTF曲线，应该根据实际插值位置来构建ERF曲线。

为了构建出均匀采样的插值位置，首先需要计算各像素点距离刀口边界的距离，刀口附近像素与刀口相对位置如图3所示，设定横坐标x以j表示，代表列数，纵坐标y以i表示，代表行数，像素间隔以p表示，纵轴与刀口边界夹角为α，以传统刀口法采用插值的中间行与刀口边界交点处的像素作为坐标原点，这样我们可以计算处于第i行，第j列像素距离刀口边界的距离Δd_i，j

Δd_i，j＝jp cosα-ip sinα                            (7)

根据Δd_i，j的符号及大小对各像素点排序，构成非均匀采样的ERF序列。传统刀口法中，直接以此ERF序列近似构建均匀采样的ERF曲线，从而带来一定的误差。本发明为了避免以上误差，采取基于重心拉格朗日插值方法，获得均匀采样插值位置的ERF新序列，而后计算获得MTF曲线。

2 基于重心拉格朗日插值方法的ERF曲线精确构建

为了精确构建ERF曲线，我们必须根据图2所示的非均匀采样的实际插值位置的像素值来构建出均匀采样位置的像素值，而不是像传统刀口测量方法中，直接以非均匀采样的实际插值位置像素值替代与之有一定距离的均匀采样插值位置的像素值。为了实现该目的，我们采用了重心拉格朗日插值方法。

拉格朗日插值法是以法国十八世纪数学家约瑟夫·路易斯·拉格朗日命名的一种多项式插值方法。许多实际问题中都用函数来表示某种内在联系或规律，而不少函数都只能通过实验和观测来了解。如对实践中的某个物理量进行观测，在若干个不同的地方得到相应的观测值，拉格朗日插值法可以找到一个多项式，其恰好在各个观测的点取到观测到的值。这样的多项式称为拉格朗日（插值）多项式。数学上来说，拉格朗日插值法可以给出一个恰好穿过二维平面上若干个已知点的多项式函数。

对于给定的若n+1个点(x₀,y₀)，(x₁,y₁)，…,(x_n,y_n)，对应于它们的次数不超过n的拉格朗日多项式L只有一个。如果计入次数更高的多项式，则有无穷多个，因为所有与L相差λ(x-x₀)(x-x₁)…(x-x_n)的多项式都满足条件。

对某个多项式函数，已知有给定的k+1个取值点：(x₀,y₀)，(x₁,y₁)，…,(x_k,y_k)

其中x_j对应着自变量的位置，而y_j对应着函数在这个位置的取值。

假设任意两个不同的x_j都互不相同，那么应用拉格朗日插值公式所得到的拉格朗日插值多项式为：

$L\left(x\right):=\underset{j=0}{\overset{k}{\mathrm{\Σ}}}{y}_j{l}_j\left(x\right)---\left(8\right)$

其中每个l_j(x)为拉格朗日基本多项式（或称插值基函数），其表达式为：

$l_j\left(x\right):=\underset{i=0,i\≠j}{\overset{k}{\mathrm{\Π}}}\frac{x-x_i}{x_j-x_i}=\frac{x-x_0}{x_j-x_0}\·\·\·\frac{x-x_{j-1}}{x_j-x_{j-1}}\·\·\·\frac{x-x_k}{x_j-x_k}---\left(9\right)$

拉格朗日基本多项式l_j(x)的特点是在x_j上取值为1，在其它的点x_i,i≠j上取值为0。

拉格朗日插值法的公式结构整齐紧凑，在理论分析中十分方便，然而在计算中，当插值点增加或减少一个时，所对应的基本多项式就需要全部重新计算，于是整个公式都会变化，非常繁琐。这时可以用重心拉格朗日插值法来代替。

重心拉格朗日插值法是拉格朗日插值法的一种改进。在拉格朗日插值法中，运用多项式

l(x)=(x-x₀)(x-x₁)…(x-x_k)                               (10)

可以将拉格朗日基本多项式重新写为：

$l_j\left(x\right):=\frac{l\left(x\right)}{x-x_j}\frac{1}{{\mathrm{\Π}}_{i=0,i\≠j}^k{x}_j-x_i}---\left(11\right)$

定义重心权

${\mathrm{\ω}}_j=\frac{1}{{\mathrm{\Π}}_{i=0,i\≠j}^k{x}_j-x_i}---\left(12\right)$

上面的表达式可以简化为：

$l_j\left(x\right):=l\left(x\right)\frac{{\mathrm{\ω}}_j}{x-x_j}---\left(13\right)$

于是拉格朗日插值多项式变为：

$L\left(x\right)=l\left(x\right)\underset{j=0}{\overset{k}{\mathrm{\Σ}}}\frac{{\mathrm{\ω}}_j}{x-x_j}{y}_j---\left(14\right)$

即所谓的重心拉格朗日插值公式（第一型）或改进拉格朗日插值公式。它的优点是当插值点的个数增加一个时，将每个ω_j都除以x_j-x_k+1，就可以得到新的重心权ω_k+1，计算复杂度为比重新计算每个基本多项式所需要的复杂度降了一个量级。

将以上的拉格朗日插值多项式用来对函数g(x)≡1插值，可以得到，对任意x：

$g\left(x\right)=l\left(x\right)\underset{j=0}{\overset{k}{\mathrm{\Σ}}}\frac{{\mathrm{\ω}}_j}{x-x_j}---\left(15\right)$

因为g(x)≡1是一个多项式。因此，将L(x)除以g(x)后可得到：

$L\left(x\right)=\frac{{\mathrm{\Σ}}_{j=0}^k\frac{{\mathrm{\ω}}_j}{x-x_j}{y}_j}{{\mathrm{\Σ}}_{j=0}^k\frac{{\mathrm{\ω}}_j}{x-x_j}}---\left(16\right)$

有了公式16，我们就可以计算任意位置x_i对应的像素值L(x_i)。这个公式被称为重心拉格朗日插值公式（第二型）或真正的重心拉格朗日插值公式。它继承了（14）式容易计算的特点，并且在代入x值计算L(x)的时候不必计算多项式l(x)。它的另一个优点是，可以很好地模拟给定的函数，使得插值点个数趋于无穷时，最大偏差趋于零。同时，重心拉格朗日插值可以达到极佳的数值稳定性。

因此，采用基于重心拉格朗日插值的方法，我们可以通过非均匀采样插值位置的像素值来构建出均匀采样位置的像素值，从而获得高分辨率的精确ERF曲线，由此降低系统MTF评估的误差。

3 本发明的数字X射线成像系统调制传递函数的刀口法测量新方法技术流程

本发明技术流程如下：

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

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

③对刀口边界直线图进行Hough变换，获得刀口倾斜角度，依据公式Δd_i，j=jp cosα-ip sinα计算刀口边界附近各点与刀口直线的距离，按照距离从小到大排序，构成非均匀采样的ERF曲线序列，α是刀口倾斜角度与探测器采样方向的夹角；

④采取基于重心拉格朗日插值方法，依据公式

计算一系列均匀采样位置x的像素值L(x),获得均匀采样插值位置的ERF新序列，ω_j是重心权重， ${\mathrm{\ω}}_j=\frac{1}{{\mathrm{\Π}}_{i=0,i\≠j}^k{x}_j-x_i};$

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

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

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

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

根据刀口倾斜角度，我们可以获得传统刀口测量方法的整数插值数，但这个插值数与真实插值数（非整数）存在误差。因此，我们计算每一个插值位置距离刀口边界的距离，按照距离从小到大排序，构成非均匀采样的ERF曲线序列(x₀,y₀)，(x₁,y₁)，…,(x_k,y_k)。而后通过基于重心拉格朗日插值方法构建均匀采样插值位置的像素值，而后重构出过采样ERF曲线。图6所示为传统刀口法与本发明方法获得的ERF曲线的对比结果。对ERF曲线进行求导，可以获得各自对应的LSF曲线，图7所示为传统刀口法与本发明方法获得的LSF曲线的对比结果，可以看出，传统刀口法获得的LSF曲线光滑性较差，因此很多研究者建议对曲线进行曲线拟合，但是曲线的选择对系统评估误差影响很大，具体采取什么曲线目前还存在很大争议。而本发明方法获得的LSF曲线已经比较光滑，避免了采取曲线拟合引入误差。

对上述LSF曲线进行傅里叶变换并取模，即可获得系统MTF曲线。为了评估本发明方法获得的MTF曲线的精度，我们采用线对卡方法获得的MTF值作为判别标准，我们采用线对卡方法获取整数位置的MTF值（仅在整数频率位置5，6，7，8，9，10有值)。图8给出了传统刀口法以及本发明方法获得的MTF曲线，以及通过线对卡方法计算获得的MTF值。从图中可以看出，采用传统刀口法测量获得的MTF曲线与线对卡方法获得的MTF准确值有比较明显的差距，而本发明方法获得的MTF曲线更接近于线对卡方法获得的MTF准确值。

最终结果表明，在进行数字X射线成像系统的调制传递函数测量时，通过本发明的基于重心拉格朗日插值的MTF刀口测量方法，解决了重构过采样ERF曲线时非均匀采样的插值位置误差，可以获得较传统刀口测量方法更为准确的MTF曲线。该方法的应用，将为准确测量系统调制传递函数，有效评估放射成像系统性能，深入开展放射影像学临床实践和研究提供有力支持。

Claims (2)

1.一种数字X射线成像系统调制传递函数的刀口法测量方法，包括下列步骤：

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

2)对刀口图像叠加平均，对平均后的刀口图像进行刀口边界检测，并且利用直线拟合获得刀口边界直线；

3)设定横坐标x以j表示，代表列数，纵坐标y以i表示，代表行数，像素间隔以p表示，刀口倾斜角度与探测器采样方向的夹角为α，对刀口边界直线图进行Hough变换，获得刀口倾斜角度，依据公式Δd_i,j=jp cosα-ip sinα计算刀口边界附近各点与刀口直线的距离，按照距离从小到大排序，构成非均匀采样的ERF曲线序列；

4)采取基于重心拉格朗日插值方法，获得均匀采样插值位置的ERF新序列；

5)对ERF新序列进行差分运算得到线扩散函数LSF(x)；

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

2.根据权利要求1所述的数字X射线成像系统调制传递函数的刀口法测量方法，第4）步中，依据公式

计算一系列均匀采样位置x的像素值L(x),获得均匀采样插值位置的ERF新序列，ω_j是重心权重，

Images