CN102422326B - 用于生成断层摄影重建滤波器的系统和方法 - Google Patents

Abstract

一种用于生成用于成像扫描器的重建滤波器的系统包括：模型库(124)，其包括用于生成所述重建滤波器的模型；滤波器准则库(122)，其包括用于生成所述重建滤波器的滤波器准则；以及滤波器生成器(120)，其基于所述滤波器模型和所述滤波器准则生成所述重建滤波器。在一个非限制性例子中，该模型是基于将代价函数最小化的，该代价函数包括诸如图像对比度的预先确定的滤波器准则。在另一个非限制性例子中，所述代价函数包括与空间分辨率、噪声和噪声存在下的信号视觉感知相关的项。

Description

用于生成断层摄影重建滤波器的系统和方法

技术领域

下文总体而言涉及重建滤波器，具体而言可应用于计算机断层摄影（CT）。然而，其还可以用于其他医学成像应用和非医学成像应用。

背景技术

计算机断层摄影（CT）扫描器通常包括安装在旋转机架上的X射线管，其围绕纵轴或Z轴在检查区域周围旋转，并且发射辐射。探测器阵列成一个角度弧，其在X射线管对面与检查区域相对。探测器阵列探测穿过检查区域（以及其中的受检者或对象）的辐射，并且生成表示受检者或对象的投影数据。重建器基于例如滤波反投影重建算法的重建算法对投影数据进行重建，并且生成表示所扫描的受检者或对象的体积图像数据。该体积图像数据可以被处理以生成一个或多个表示所扫描的受检者或对象的图像。

锥形束CT的常规滤波反投影重建算法是基于对数据进行二维（2D）图像傅里叶分解且然后进行傅里叶合成的。通常，具有灰度级的CT图像可以被认为是在受检者或对象的切片内的X射线衰减系数的2D函数。根据傅里叶定理，这一函数可以被分解为在横向平面中传播的波（傅里叶分量）。因此，这一函数可以被写为波的线性组合或叠加。在傅里叶域中通过波的频率k的绝对值并通过由极角θ给出的波的方向来参数表示所述波。在重建期间，基于Radon切片定理使用处理步骤来估计波的振幅和相位。然后对这些波进行求和以形成图像。

在对波进行求和之前，有意地对在重建期间所估计的波振幅进行修改。由重建滤波器控制这一修改。通过示例的方式，由于X射线束的有限宽度，从测量结果估计的高频波的振幅小于它们在X射线衰减校正图内的真实值。通过预先获知振幅的估计值和真实值之间的频率依赖比（frequencydependent ratio）（系统的调制传递函数（MTF）），可以强调高频波的所估计的振幅以补偿这些频率处的MTF下降。然而，降低对所估计的振幅的值的强调能够降低重建的图像中的噪声功率谱（NPS）。这一噪声的主要来源是要到达探测器的光子的泊松概率分布。

对用于特定扫描器的重建滤波器进行修改通常包括对某些频率内容进行强调或降低强调。照惯例，基于用户（例如放射学家、应用专家等）反馈来执行对用于特定应用的特定滤波器的优化。遗憾的是，基于用户反馈根据经验来对重建滤波器进行优化可能使优化过程在相对长的时间段内延长。此外，这样的优化通常限于将滤波器频率依赖性调整到取值范围为若干预先定义的参数的一小组函数。

发明内容

本申请的各方面解决上述问题及其他问题。

根据一方面，一种用于生成用于成像扫描器的重建滤波器的系统包括：模型库，其包括用于生成重建滤波器的模型；滤波器准则库，其包括用于生成重建滤波器的滤波器准则；以及滤波器生成器，其基于滤波器模块和滤波器准则生成重建滤波器。

在另一个实施例中，一种方法包括：将表示为两项之和的代价函数最小化，其中第一项对归因于成像系统的空间分辨率损失进行处罚，并且第二项对由于图像噪声的信号探测损失进行处罚。

在另一个实施例中，一种计算机可读存储介质包括指令，当由计算机执行时，所述指令使计算机执行以下动作：将表示为第一项和第二项之和的代价函数最小化，其中第一项对归因于成像系统的空间分辨率损失进行处罚，并且第二项对由于图像噪声的信号探测损失进行处罚。

在另一个实施例中，一种方法包括：通过将基于预先确定的对比度准则的代价函数最小化来生成滤波器，该代价函数包括与空间分辨率损失相关的项和与噪声存在下的视觉信号感知损失相关的项。

附图说明

本发明可以体现为各种部件和部件的布置，并且可以体现为各种步骤和步骤的布置。附图仅用于举例说明优选实施例的目的，并且不应被解释为限制本发明。

图1举例说明了示例性成像系统；

图2举例说明了示例性方法。

具体实施方式

图1举例说明了一种成像系统或CT扫描器100。成像系统100包括固定机架102和旋转机架104，旋转机架104由固定机架102可旋转地支撑，并且围绕沿着纵轴或Z轴的旋转框架旋转中心在检查区域106周围旋转。例如X射线管的辐射源108由旋转机架104支撑。辐射源108发射辐射。源准直器将发射的辐射进行准直以形成大致锥形、扇形、楔形或其他形状的射束，该射束穿过检查区域106。辐射敏感探测器阵列110成一角度弧，其跨过检查区域106与辐射源108相对，并且探测穿过检查区域106的辐射。探测器阵列110生成表示所探测的辐射的投影数据或信号。例如躺椅的患者支架支撑检查区域106中的对象或受检者。

重建器112基于所选择的重建算法对来自探测器阵列110的信号进行重建，并且生成表示该信号的体积图像数据。一般用途的计算机系统作为操作者控制台114。驻留在控制台114中的软件使操作者能够通过选择扫描协议、接受或改变用于所选择的协议的重建滤波器等来控制系统100的操作。滤波器识别器116在滤波器库118中识别所选择的滤波器，所述滤波器库118可以包括针对不同应用进行各种调整的多个不同的滤波器（例如，尖锐的、平滑的等）。

重建滤波器生成器120生成一个或多个可以由重建器112使用的重建滤波器（例如，尖锐的、平滑的等）。滤波器生成器120可以为系统100的部分，或者远离系统100，例如作为分离的计算机、工作站等的部分。滤波器生成器120可以基于存储在滤波器准则库122中的各种预先确定的准则、存储在滤波器模型库124中的滤波器模型和/或其他信息来生成重建滤波器。通过非限制性示例的方式，在一个例子中，滤波器生成器120基于预先确定的图像对比度准则和代价函数来生成重建滤波器（不具有用户干涉，例如不具有对任何参数的调整）。在其他实施例中，以其他方式生成滤波器。

在傅里叶域中，代价函数可以被表示为两个在k上的积分的和，两个积分依据系统的MTF和NPS来表达。在不对所估计的波振幅进行强调或降低强调的情况下，MTF和NPS取决于扫描器的几何结构——例如探测器像素尺寸和X射线管到探测器的距离，并且可以被测量、模拟或近似。这可以通过针对低于频率v_c的所有频率将重建滤波器值选择为1来实现，所述频率v_c例如被定义为MTF的截止频率。通常，代价函数的第一个积分对应由于空间拖尾效应（smearing）的信号探测能力损失，并且代价函数的第二个积分对应由于图像噪声的探测信号能力损失。两个积分可以包含类似或不同的要素（integrant），并且可以考虑先验信息。

如下文更加详细地描述的那样，在一个非限制性的例子中，通过将代价函数最小化来生成滤波器。在这一示例中，代价函数依据信号功率谱（SPS）来表示，SPS表示在傅里叶域中信号的统计学能量分布。在CT中，可以基于衰减校正图的总体来定义SPS，并且SPS对波频率的依赖性可以被表示为如等式1所示的：

等式1：SPS = SPS(k, r)，

其中k表示频率的绝对值，并且r表示图像像素和旋转框架旋转中心之间的距离。假设对r的依赖性是不重要的，SPS可以被表示为如等式2所示的：

等式2： $\mathrm{SPS}\left(k\right)\∝{\left(\frac{k_0}{k}\right)}^p,p~1.9,$

在所举例说明的实施例中，重建滤波器针对在一（1）到一又十分之九（1.9）的范围内的p进行微小变化。

在代价函数中的分辨率损失可以通过考虑作为针对给定重建滤波器F(k)获得的图像（排除噪声和混叠伪影）和表示辐射衰减系数图的图像之间的差别的误差图像的功率谱而进行处罚。基于上文，代价函数的第一个积分的要素可以被表示为振幅和功率谱之间的二次关系，如等式3所示：

等式3：INT_s(k)∝SPS(k)·(1-MTF(k)·F(k))²。

在这一示例中，代价函数的第二个积分考虑了函数I(k`/k)，其中k`和k分别表示信号波的频率绝对值和噪声波的频率绝对值。函数I(k`/k)描述了频率比k`/k如何对人类观测者在噪声成分存在的情况下对信号进行探测的能力施加影响。

各种定义可以被用于I(k`/k)。下文描述了一个适当的示例。假定：信号图像仅由一个正弦曲线组成，所述正弦曲线例如以频率绝对值k`在该图像内围绕零振荡的平面波。还假定：噪声图像是通过使用以零为中心的随机高斯生成器建立白噪声图像（在临近像素的值之间不具有相关性）而形成的。可以在白噪声图像上执行2D傅里叶分解，并且针对除了频率绝对值落入k-δ和k+δ之间的那些波之外的所有波，振幅可以被设置为零。然后可以在分解的信号上执行傅里叶合成。

注意：由上述步骤获得的噪声图像包含不同方向（θ）的波。换言之，这些波在频域中占据了由半径k-δ和k+δ界定的圆环。上文描述的噪声图像和信号图像的方差可以被表达为如等式4所示的：

等式4： $\mathrm{\σ}=\underset{N\→\∞}{\lim }\frac{\underset{S=1}{\overset{N}{\mathrm{\Σ}}}{(i\left(s\right)-\underset{q=1}{\overset{N}{\mathrm{\Σ}}}i\left(q\right)/N)}^2}{N-1},$

其中，N为图像像素的数量，并且i(s)\i(q)为图像在像素s\q中的值。可以增大信号图像的方差，以使得当噪声图像被加入信号图像时人类观测者将能够探测到信号。或者，可以减小信号图像的方差，以使得在噪声图像被加入信号图像之后人类观测者将不再探测到信号。

因而，假定信号图像方差和噪声图像方差之间的最小比，其中观测者仍然探测到信号。允许观测者对图像进行放大和缩小，这一最小方差比仅取决于在k`和k之间的比。在一个实施例中，I(k`/k)被认为是这一最小比。在一个例子中，I(k`/k)可以根据经验进行估计。之后经验测量点可以被拟合到I对log(k`/k)的高斯依赖性。

下文讨论了如何将信息信号由于位于具有以k为中心的绝对值的圆环形频率区域的噪声波而被退化（degrade）的量包括到代价函数中。这一量与NPS(k)和在k处的重建滤波器（由F(k)表示）的平方的乘积成比例。因此，第二个积分的要素通常遵循等式5：

等式5：INT_噪声(k)∝NPS(k)·F(k)²，

其中NPS(k)表示针对低于v_c的所有频率使用F(k)=1而获得的在k处的噪声功率谱。因此，等式5可以被用于确定在k处的噪声功率密度。

第二个积分的要素还可以与这一功率密度可以使对由SPS(k)分布的信号的探测退化的效果成比例。在一个实施例中，这一效果可以被估计为I(k`/k)在信号频率k`上的积分，这一积分根据在这一频率处的信号功率谱的值进行加权，如等式6所示：

等式6： $e_n\left(k\right)=\underset{0}{\overset{v_c}{\∫}}\mathrm{dk}`k`\·\mathrm{SPS}(k`)\·I(k`,k),$

其中k`为从笛卡尔到极坐标的转换矩阵的雅可比行列式。如上文所述，SPS(k`)表示衰减系数图的总体。在一个例子中，例如可以通过将等式6中的要素乘以MTF(k`)<sup>2</sup>·F(k`)²来将该总体转换为重建的CT图像的总体。

根据等式5和等式6，代价函数的要素包括F(k)²·F(k`)<sup>2</sup>的乘积。然而，在另一个实施例中，作为近似值，将等式6中的F(k`)²替换成1，这将简化滤波器F(k)的求解，并且为代价函数提供全局最小值。等式6中的e_n(k)还可以被表达为依据I(k`/k)的加权平均值的无量纲函数，如等式7所示：

等式7： $e_n^A\left(k\right)=\frac{\underset{0}{\overset{v_c}{\&Integral;}}\mathrm{dk}`k`\&CenterDot;\mathrm{SPS}(k`)\&CenterDot;{\mathrm{MTF}(k`)}^2\&CenterDot;I(k`,k)}{\underset{0}{\overset{v_c}{\&Integral;}}\mathrm{dk}`k`\&CenterDot;\mathrm{SPS}(k`)\&CenterDot;\mathrm{MTF}{(k`)}^2},$

其中，A表示近似值或k`上的平均值。

根据等式1、3、5和7，代价函数可以被表示为如等式8所示的：

等式8：

$E_{\mathrm{\&alpha;}}\left(F\left(k\right)\right)=\underset{0}{\overset{v_c}{\&Integral;}}\mathrm{dkk}\&CenterDot;\mathrm{SPS}\left(k\right)\&CenterDot;{(1-\mathrm{MTF}\left(k\right)\&CenterDot;F\left(k\right))}^2+\mathrm{\&alpha;}\&CenterDot;\underset{0}{\overset{v_c}{\&Integral;}}\mathrm{dkk}\&CenterDot;\mathrm{NPS}\left(k\right)\&CenterDot;e_n^A\left(k\right)\&CenterDot;F{\left(k\right)}^2,$

其中E_α(F(k))表示取决于重建滤波器F(k)的泛函（functional）值或标量值。乘以代价函数的第二项的参数α被用于在由空间分辨率损失引起的信号退化和由噪声引起的信号退化之间进行平衡。在一个实施例中，如接下来描述的，这一参数与重建滤波器F(k)一起被自动计算出。

点(x, y)处的图像噪声方差可以被表达为如等式9所示的：

等式9：

$\mathrm{\&sigma;}(x,y)={\left(\frac{\mathrm{\&pi;}}{M_{\mathrm{proj}}}\right)}^2\&CenterDot;2\&CenterDot;\mathrm{\&tau;}\&CenterDot;\left(\underset{i=1}{\overset{M_{\mathrm{proj}}}{\mathrm{\&Sigma;}}}\frac{1}{\stackrel{\&OverBar;}{N}(x,y,\mathrm{\&theta;}\left(i\right))}\right)\&CenterDot;{\&Integral;}_0^{v_c}\mathrm{dk}{k}^{2}F{\left(k\right)}^2,$

其中M_proj表示反投影到像素上的投影的数量，τ表示读数之间的径向增量，并且表示到达属于以θ角与点(x, y)相交的投影i的读数的探测器的光子的数量的平均值。

使用重建滤波器F(k)获得的图像噪声方差和通过将其所有值都用1来替换获得的图像噪声方差之间的比可以被表示为如等式10所示的：

等式10： $\mathrm{rv}\left(F\left(k\right)\right)\&equiv;\mathrm{\&sigma;}\left(F\left(k\right)\right)/\mathrm{\&sigma;}(F\left(k\right)\&DoubleRightArrow;1)=\frac{\underset{0}{\overset{v_c}{\&Integral;}}\mathrm{dk}{k}^2{F\left(k\right)}^2}{\underset{0}{\overset{v_c}{\&Integral;}}\mathrm{dk}{k}^2}\&CenterDot;$

由于不同的效应，如对投影的扇形到平行的重排、探测器像素之间的串扰和其他，rv对F(k)的依赖性可能不同于在等式10中的近似值。然而，其仍然可以基于所测量或模拟的NPS(k)而被计算为如等式11所示的：

等式11： $\mathrm{rv}\left(F\left(k\right)\right)=\frac{\underset{0}{\overset{v_c}{\&Integral;}}\mathrm{dkkNPS}\left(k\right){F\left(k\right)}^2}{\underset{0}{\overset{v_c}{\&Integral;}}\mathrm{dkkNPS}\left(k\right)}\&CenterDot;$

在一个例子中，对于由rv_g表示的rv(F(k))的给定值，通过以下快速收敛迭代程序可以自动同时计算出参数α和滤波器F(k)。通过示例的方式，α被设置为α₀=1，并且F₀(k)被计算为将代价函数E_α0(F(k))最小化的滤波器。这可以被表达为如等式12所示的：

等式12： $F_0\left(k\right)=\underset{\left\{F\left(k\right)\right\}}{\min }{E}_{{\mathrm{\&alpha;}}_0}\left(F\left(k\right)\right)\&CenterDot;$

如等式13所示来执行n次迭代：

等式13：

$\begin{array}{c}{\mathrm{\&alpha;}}_1=\frac{{\mathrm{\&alpha;}}_0\&CenterDot;\mathrm{rv}\left(F_0\left(k\right)\right)}{{\mathrm{rv}}_g}{F}_1({\mathrm{rv}}_g,k)=\underset{\left\{F\left(k\right)\right\}}{\min }{E}_{{\mathrm{\&alpha;}}_1}\left(F\left(k\right)\right)\\ {\mathrm{\&alpha;}}_2=\frac{{\mathrm{\&alpha;}}_1\&CenterDot;\mathrm{rv}\left(F_1\left(k\right)\right)}{{\mathrm{rv}}_g}{F}_2({\mathrm{rv}}_g,k)=\underset{\left\{F\left(k\right)\right\}}{\min }{E}_{{\mathrm{\&alpha;}}_2}\left(F\left(k\right)\right)...\\ \&CenterDot;\&CenterDot;\&CenterDot;\\ {\mathrm{\&alpha;}}_n=\frac{{\mathrm{\&alpha;}}_{n-1}\&CenterDot;\mathrm{rv}\left(F_{n-1}\left(k\right)\right)}{{\mathrm{rv}}_g}{F}_n({\mathrm{rv}}_g,k)=\underset{\left\{F\left(k\right)\right\}}{\min }{E}_{{\mathrm{\&alpha;}}_n}\left(F\left(k\right)\right),\end{array}$

如等式14所示，具有为F_n(rv_g,k)设置的停止准则。

等式14：

$\mathrm{abs}(\frac{\mathrm{rv}\left(F_n({\mathrm{rv}}_g,k)\right)}{{\mathrm{rv}}_g}-1)\&le;\mathrm{\&epsiv;},$

其中，ε的典型值为大约10^-4，并且重建滤波器可以被表达为如等式15所示的：

等式15：F(rv_g,k)=F_n(rv_g,k)。

现有市面上的CT扫描器包括一组重建滤波器，其为rv(F(k))提供一组不同的值。对于本文的实施例，可以可选地根据这一组不同的值设置rv_g的值。通过例如减小像素尺寸来改变扫描器几何结构可以激活对rv_g的新值（在这一情况下更高）的包括。在任何情况下，甚至当针对新滤波器取的rv_g的值等于旧滤波器的rv_g的值时，由旧滤波器和由通过等式8和等式12-14获得的新的最优滤波器不同地对频率进行强调或降低强调。

基于等式7，将等式8中的代价函数最小化的重建滤波器可以被表示为如等式16所示的：

等式16： $F_i\left(k\right)=\frac{\mathrm{MTF}\left(k\right)}{\mathrm{MTF}{\left(k\right)}^2+{\mathrm{\&alpha;}}_i\&CenterDot;{(k/k_0)}^p\&CenterDot;\mathrm{NPS}\left(k\right)\&CenterDot;e_m^A\left(k\right)}\&CenterDot;$

在NPS(k～0)=0并且MTF(k～0)=1的例子中，等式16在k=0时得出F_i(k)的零阶导数。由于射束的有限宽度，MTF(k)通常在NPS(k)之前下降。基于等式16，在的区域内，F(k)可能快速下降。

图2举例说明了一种用于自动生成用于成像系统的重建滤波器的方法。在步骤202，识别预先确定的滤波器准则。在这一示例中，该准则为图像对比度，包括在噪声存在下人类观测者探测信号的能力，并且生成滤波器以使预先确定的图像对比度准则最优化。在其他实施例中，可以基于其他准则生成滤波器。在步骤204，识别系统的成像特性，例如MTF、NPS、SPS和/或系统的其他特性。如上文所述，可以对MTF和NPS进行测量、模拟或近似。在步骤206，识别依据成像特性表示的代价函数。在步骤208，基于预先确定的滤波器准则、MTF和NPS，通过将代价函数最小化来生成滤波器。如本文描述的，可以在傅里叶域中在波的频率（k）上表示代价函数。这样的代价函数可以考虑由于空间拖尾效应的信号探测能力损失和由于图像噪声的探测信号能力损失。在步骤210，采用该滤波器以对由成像系统生成的投影数据进行重建。在步骤212，可以针对一个或多个滤波器或者一个或多个类似或不同的成像系统重复上述操作。

上述操作可以作为扫描器的控制台应用程序和/或图像处理或规划工作站来实施。作为示例，上述操作可以以计算机可读指令的方式来实施，当被一个或多个计算机处理器（控制台或工作站的处理器）执行时，所述计算机可读指令使处理器执行上述动作。在这样的情况下，在与相关计算机关联或者对于相关计算机可访问的计算机可读存储介质中存储这些指令。

在本文已经参考各种实施例描述了本发明。在阅读了本文的描述的情况下，其他人员可以想到各种修改和替代。意在将本发明解释为包括所有落在随附的权利要求及其等价要件的范围内的这样的修改和替代。

Claims (13)

1.一种用于生成用于成像扫描器的重建滤波器的系统，所述系统包括：

模型库（124），所述模型库（124）包括用于生成所述重建滤波器的模型；

其中，所述模型基于被表示为两项的和的代价函数，其中，第一项对归因于所述成像扫描器的空间分辨率损失进行处罚，并且第二项对由于图像噪声的信号探测损失进行处罚，

滤波器准则库（122），所述滤波器准则库（122）包括用于生成所述重建滤波器的滤波器准则；以及

滤波器生成器（120），所述滤波器生成器（120）基于用于生成所述重建滤波器的所述模型和所述滤波器准则来生成所述重建滤波器。

2.根据权利要求1所述的系统，其中，所述代价函数对所述滤波器准则的退化进行建模。

3.根据权利要求1所述的系统，其中，所述滤波器生成器（120）通过将所述代价函数最小化来生成所述滤波器。

4.根据权利要求1所述的系统，其中，所述滤波器准则包括使图像对比度最优化。

5.根据权利要求1到4中任意一项所述的系统，其中，所述第一和第二项取决于所述扫描器的调制传递函数、所述扫描器的噪声功率谱以及所扫描的受检者或对象的统计学平均信号功率谱。

6.根据权利要求1到4中任意一项所述的系统，其中，所述第二项基于预先确定的信噪比。

7.根据权利要求6所述的系统，其中，所述预先确定的信噪比对应于最小信噪方差比的近似值，在该最小信噪方差比，人类观测者能够从噪声中辨别出已知的信号。

8.根据权利要求7所述的系统，其中，所述信噪方差比为在特定频率的单个正弦信号的方差与包含噪声功率谱的各向同性噪声图像的方差之间的比，该噪声功率谱仅在傅里叶域中的有限宽度的圆环上有限。

9.根据权利要求7所述的系统，其中，所述最小信噪方差比是根据经验进行测量，并且使用高斯拟合进行平滑。

10.根据权利要求1到4中任意一项所述的系统，其中，所述滤波器生成器（120）通过对所述第一和第二项进行平衡的平衡参数将所述代价函数最小化。

11.根据权利要求10所述的系统，其中，所述平衡参数通过至少一个预先定义的参数来确定。

12.根据权利要求11所述的系统，其中，所述预先定义的参数为对应于最优滤波器的图像噪声方差与对应于预先确定的参考滤波器的图像噪声方差之间的比。

13.一种用于生成成像扫描器的重建滤波器的方法，其包括：通过将基于预先确定的对比度准则的代价函数最小化来生成重建滤波器，其中，所述代价函数包括与空间分辨率、噪声和噪声存在下的信号视觉感知相关的项。