CN1818973A - 基于hl一致性条件的ct投影数据射束硬化效应校正方法 - Google Patents

Abstract

本发明公开了一种基于数据相关性的CT投影数据射束硬化效应的校正方法，该方法根据医用诊断X线物理学成像模型为基础，构造一个通用的CT投影数据的射束硬化效应校正模型；然后根据H－L一致性条件，构造CT投影数据之间的关系矩阵方程；结合上述校正模型和关系矩阵，得到校正模型的参数的求解方程；为了实现更加精确的校正结果，根据初次校正重建结果和CT重投影，估计投影数据中的高密度物质的衰减比例。本发明适合于各类X线CT设备的射束硬化校正功能的实现。和现有的只采用模体校正的方法相比较，能够分别考虑不同密度组织的射束硬化效应的差异，自动适应成像对象的不同，具有更好的校正精度。

Description

基于HL一致性条件的CT投影数据射束硬化效应校正方法

                          技术领域

本发明属于X线CT重建成像技术领域，涉及基于等效非线性模型和H-L一致性条件的X线序列投影图像数据的射束硬化矫正方法。

                          背景技术

CT重建的原理基于RADON变换及其反变换理论[1]，其基础是投影数据来自于线性投影模型，比如，对于二维断层成像来讲，有如下的投影公式：

$p(\mathrm{\θ},l)={\∫}_{L(\mathrm{\θ},l)}\mathrm{\μ}(x,y)\mathrm{dxdy}$

其中，L(θ，l)为射线穿透路径，θ为投影角度，l为投影数据的采样点；

虽然CT技术已经从第一代发展到当今的第五代，扫描轨迹从最早的直线，一直到现在的圆形和螺旋轨迹、乃至任意轨迹，投影方式有扇形和锥形束，其投影和重建理论的基础一直基于线性投影模型。所进行的重建理论研究总是围绕着不同的扫描轨迹和投影方式进行的。

理论上，如果CT的X线源是单能的，那么上述的线性投影模型的要求则是满足的，这是因为物质对特定能量X线光子的衰减成指数衰减形式，取对数后即可满足所谓的线性投影模型，如下式：

$I(\mathrm{\θ},l)=I_0\exp \{-{\∫}_{L(\mathrm{\θ},l)}u(x,y)\mathrm{dxdy}\}$

上式中，I₀为球管发出X线光子强度，u(x，y)为成像对象的质量衰减系数乘以物质密度。

事实上，现代CT的球管发出的X线光谱是具有一定能谱分布的宽谱函数。以头部CT常用的120kVp电压下，球管发出的X线光谱为例，其X线光子能量从30keV到120keV按照一定的规律连续分布，而且其分布函数在58-68keV处还存在4根不连续的特征辐射光谱，如图(1)所示，图中CT球管发出的X线光谱函数图形(上部)，下部图形是经过35mm铝滤过(根据CT设备通常的做法)后的光谱函数(放大了5倍)，所形成的投影成像数据如下：

$I(\mathrm{\θ},l)=I_0{\∫}_0^{V}P(V,E)\exp \{-{\∫}_{L(\mathrm{\θ},l)}u(E,x,y)\mathrm{dxdy}\}\mathrm{dE}$

上式中，P(V，E)为球管电源取VkVp时，发出的X线光子谱乘以探测器的吸收谱的归一化函数，u(E，x，y)仍为成像对象的质量衰减系数乘以物质密度，但是考虑到该参数随X线光子能量E变化的关系。上式取对数后，不再满足上述的线性投影模型。以人体常见组织-骨皮质、肌肉和脂肪为例，不同厚度的组织对120kVp电压下，球管发出的宽谱X线光子的衰减数据取对数后函数关系如图(2)所示(假设X线探测器是Na:CsI晶体，下同)，图(2)中的骨组织、肌肉和脂肪在120kVp电压球管发出的光谱中，衰减量和厚度的关系，其中，粗线代表骨组织，下面两条细虚线代表肌肉和脂肪，拱形的细实线表示骨组织衰减与线性关系的差异(放大了10倍)，另外两条细虚线则分别反映了肌肉和脂肪的情形(放大了10倍)。

线性投影模型要求上述函数是直线形式，但实际上，该函数并不满足其要求。这样直接运用重建理论得到的结果是：对于同样的人体组织，不同空间区域的重建数值(所谓CT值)会出现不同，由于CT设备提供主要的信息反映在CT值上，因此这会给诊断带来至关重要的误差，这就是所谓的“射束硬化效应(Beam Hardening)”所形成的CT值误差被称为“射束硬化伪影”，是CT设备的主要重建伪影之一。

现有CT设备中，实现所谓射束硬化校正(Beam Hardening Correction)的方法来自于工程实践。观察图(2)，在宽谱条件下，虽然投影数据与组织的厚度不是线性关系，但与直线相差不大，因此，可以考虑采用一个多项式来描述该函数，并且根据该多项式系数来实现线性化校正，考虑到校正精度，CT设备在球管的出线处加上一定厚度金属材料(一般等效为35mm厚度的铝)，将光谱中低能光子尽可能地滤除，使穿透人体的X线投影数据的线性度尽可能地好些，再进行多项式校正，因此，得到现有CT设备的射束硬化校正措施1。

校正措施1：将投影数据取对数后，利用一个已知的多项式进行线性化校正，而该多项式的系数是通过事先的模体(水模)实验来获得，由于存在球管老化等问题，需要定期校正。

注意到图(2)所示，不同的人体组织(骨皮组织、软组织)具有不同的校正多项式系数，而在投影数据校正时，并不知道人体中不同组织成份的多少，因此，为了避免人体不同结构带来的校正误差，现有的CT设备采用了第二种校正措施。

校正措施2：在球管到探测器的空间中，除了人体组织所占的空间外，其它空间填满一种类似于人体软组织衰减特征的填充物(树脂)，这样，不同角度投影时，射线穿透的物质基本上都类似人体软组织(血管、血液和肌肉等组织的衰减性质与软组织相当，脂肪组织的衰减特性基本上和软组织相当，只是密度不同)，而高密度骨组织所占的比例很小，因此，可以假设校正方程系数与投影角度和纵向坐标(longitudinal)基本上无关，故可以采用统一的多项式进行射束硬化校正，产生的误差控制在允许的范围内。

事实上，现有医学CT的射束硬化校正方法是上述两种措施的组合。显然，这样的方法存在以下问题。

现有CT设备射束硬化校正方法存在的缺点：

1.每一台CT需要进行定期的水模实验，以获取校正多项式系数，增加使用难度和成本。

2.为了适应人体各部分的差异，需要按照部位分别设计模体，进行实验，来获得不同的校正多项式系数。

3.对于同一人体部位，存在个体差异、以及纵向(longitudinal)差异，因此用统一系数的多项式校正，可能会存在问题，这里最典型的是头部CT中的硬膜外伪影，目前尚无理想解决方案，在大部分的CT设备中，需要靠医生的临床经验来判断区分，而在其它的一些CT设备中，有采用后处理方法校正的方案，但存在重建数据不真实的风险。

4.为了采用统一系数的多项式校正，而且实现一定的校正精度，CT设备不得不采用填充物使不同角度的投影数据均匀，这实际上带来两个缺点：A、增加病人的辐射剂量和球管的热容量；B、均匀的投影数据意味着较小的动态范围，因此，重建的噪声相应增加，降低了重建的灰度精度，为了降低噪声，不得不增加X线剂量，这又增加了病人的辐射伤害风险。

CT设备的核心虽然是重建算法，但为了实现准确的重建效果，众多的校正算法一直是各类CT设备关键核心技术。除了射束硬化校正，还有各种基准校正、探测器通道非线性校正等等。在这些校正中，只有射束硬化校正与扫描的成像对象有关，其它则只与设备相关。与设备相关的校正无论复杂与否，均可以通过一个特定的程序精确地进行。而与成像对象有关的射束硬化校正，如果只采用现有的模体校正方法，那么在理论上不可能做到精确。也就是说，该方法只是一种工程化的解决问题方案，缺乏严格的理论支持。

申请人在对CT设备射束硬化校正的深入研究过程中，注意到这样的事实：在不同角度的投影成像中，成像对象中的每一个坐标点被多个射线穿过，也就是说各个不同角度的投影数据均包含了该点的信息，这意味着不同角度的投影数据是相关的。如果存在一个具有有限参数的校正方程，那么，可能可以通过这种相关性解出校正方程的系数，从而实现校正的目标，而且，这样的校正来自于对象的投影数据本身，是自适应的，理论上会有更好的精度，有可能避免上述缺点。

描述成像对象空间信息与投影数据的相关性可以用H-L一致性条件表述。但该条件描述的是二维平行束投影的情形，不符合X线CT(XCT)的投影情况，申请人构造了二维扇形束的H-L一致性条件公式，并且在此基础上，求解预先设定的非线性方程参数，从而达到了线性校正的目的。

考虑到现有的研究结果表明，XCT的投影数据校正可以由多项式完成，因此，理论上该方法可以用于宽谱XCT的投影数据校正。以上工作无疑给基于投影数据自身相关性进行CT数据线性化校正的研究工作，提供了可靠的理论基础。

                          发明内容

本发明的目的在于，提供一种基于HL一致性条件的CT投影数据射束硬化效应校正方法，该方法适用于各类CT设备重建的数据校正。

实现上述发明目的的技术解决方案是：

一种基于HL一致性条件的CT投影数据射束硬化效应校正方法，其特征在于，该方法根据等效射束硬化失真校正模型和H-L一致性条件，基于不同角度的投影图像数据，通过对等效射束硬化失真校正模型参数的求解来实现对投影数据的校正，具体包括以下步骤：

步骤1：构建投影数据的等效非线性失真校正模型

对于医学CT设备，根据成像物理学模型，将观测到的投影数据分解为高密度组织、低密度组织分别衰减而成的具体模型，构造等效射束硬化失真校正模型：

$g(t,\mathrm{\β})=P({\stackrel{\→}{\mathrm{\α}}}_k,r(t,\mathrm{\β})f(t,\mathrm{\β}))+P({\stackrel{\→}{\mathrm{\η}}}_k,(1-r(t,\mathrm{\β}))f(t,\mathrm{\β}))\·\·\·\left(1\right)$

上式中，f(t，β)为探测器获得的投影数据，β为投影角度，t为投影数据的采样点，g(t，β)为校正后的投影数据，r(t，β)表示每一条投影射线的路径上，高密度骨组织对投影数值的贡献，取值范围是(0，1)，P()为线性校正的多项式方程，

为多项式方程的系数， 为系数的多项式形式如下：

$P({\stackrel{\→}{\mathrm{\η}}}_k,x)={\mathrm{\η}}_{1}x+{\mathrm{\η}}_2{x}^2+{\mathrm{\η}}_{1/2}{x}^{1/2}+{\mathrm{\η}}_3{x}^3+{\mathrm{\η}}_{1/3}{x}^{1/3}+\·\·\·\·\·\·+{\mathrm{\η}}_L{x}^L+{\mathrm{\η}}_{1/L}{x}^{1/L}\·\·\·\left(2\right)$

其中，

为系数的多项式形式与 的相同；

步骤2：构造符合CT投影模式的H-L一致性条件的表达公式

根据已知平行束投影的HL一致性性条件，推广至圆形扫描轨迹、扇形投影方式的表达公式是：

$m_{i,k}=\∫{\∫}_C{x}^i{y}^{k}u(x,y)\mathrm{dxdy},i\≥0,k\≥0\·\·\·\left(3\right)$

$v_d\left(\mathrm{\β}\right)={\∫}_{-\∞}^{+\∞}{g}_f(t,\mathrm{\β}-\mathrm{arctg}\frac{t}{D}){\left(\frac{\mathrm{tD}}{\sqrt{D^2+t^2}}\right)}^d\left(\frac{D^3}{{(D^2+t^2)}^{3/2}}\right)\mathrm{dt}\·\·\·\left(4\right)$

$Q_d\left(\mathrm{\β}\right)=\underset{r=0}{\overset{d}{\mathrm{\Σ}}}\left(\begin{array}{c}d\\ r\end{array}\right)m_{r,d-r}{\cos }^r\left(\mathrm{\β}\right){\sin }^{d-r}\left(\mathrm{\β}\right)\·\·\·\left(5\right)$

        v_d(β)＝Q_d(β)                               (6)

其中，0≤d≤N-1，而N为投影角度的个数；

上列公式中，原始的断层图像为u(x，y)，R为其支撑半径，C＝{(x，y)|x²+y²≤R²}，β为扇型束的投影角度；

步骤3：根据等效非线性失真校正模型和H-L一致性条件，设计模型参数的求解矩阵，通过计算得到模型参数的求解，从而完成对投影数据的校正过程。

本发明的方法适合于各类X线CT设备的射束硬化校正功能的实现。射束硬化现象对于许多细小组织病变的诊断和工业CT等应用场合意义至关重要。由于射束硬化校正与扫描的成像对象有关，本发明的方法与现有各类CT中，只采用模体校正的方法相比较，能够分别考虑不同密度组织的射束硬化效应的差异，自动适应成像对象的不同，具有更好的校正精度。

                          附图说明

图1是CT球管发出的X线光谱函数图形(上部)，下部图形是经过35mm铝滤过(根据CT设备通常的做法)后的光谱函数(放大了5倍)；

图2是骨组织、肌肉和脂肪在120kVp电压球管发出的光谱中，衰减量和厚度的关系，其中，粗线代表骨组织，下面两条细虚线代表肌肉和脂肪，拱形的细实线表示骨组织衰减与线性关系的差异(放大了10倍)，另外两条细虚线则分别反映了肌肉和脂肪的情形(放大了10倍)。

图3是圆形扫描轨迹、扇形投影束的示意图，β为投影角度，S为射线源，g_f(t，β)为投影数据，O为旋转中心，D为S到O的距离；

图4是实验用的FORBILD头部模型，用于射束硬化校正的说明，其中，亮的部分是骨皮组织，灰色部分是软组织，黑色部分是空气。投影参数为：球管电压120kVp，附加滤过为35mm厚度的铝。图中显示利用所提出的基于H-L一致性条件校正方法获得的重建图像，重建方法采用FBP。为了显示校正前后的重建效果，在后续的图中，给出图中亮线所指示位置的剖线数据。

图5(a)是水平剖线位置上的校正前后重建数据的对比。三条曲线分别表示经过HL一致性校正方程校正后重建的结果、未经数据校正的重建结果、校正前后数据的差值，其中校正前后数据的差值放大了5倍。

图5(b)是垂直剖线位置上的校正前后重建数据的对比。三条曲线分别表示经过HL一致性校正方程校正后重建的结果、未经数据校正的重建结果、校正前后数据的差值，其中校正前后数据的差值放大了5倍。

图6是利用所提出的基于分割后的H-L一致性条件校正方法获得的重建图像，重建方法采用FBP。

图7(a)是水平剖线位置上的校正前后重建数据的对比。三条曲线分别表示经过公式(22)校正方程校正后重建的结果、未经数据校正的重建结果、校正前后数据的差值。其中，r(t，β)是通过预校正重建图像分割后计算而获得的，校正前后数据的差值放大了5倍。

图7(b)垂直剖线位置上的校正前后重建数据的对比。三条曲线分别表示经过公式(22)校正方程校正后重建的结果、未经数据校正的重建结果、校正前后数据的差值。其中，r(t，β)是通过预校正重建图像分割后计算而获得的，校正前后数据的差值放大了5倍。

以下结合附图对本发明作进一步的详细说明。

                          具体实施方式

1.构建投影数据的等效非线性失真校正模型

对于医学CT设备，根据成像物理学模型，将观测到的投影数据分解为高密度组织、低密度组织分别衰减而成的具体模型，构造等效射束硬化失真校正模型：

$g(t,\mathrm{\β})=P({\stackrel{\→}{\mathrm{\α}}}_k,r(t,\mathrm{\β})f(t,\mathrm{\β}))+P({\stackrel{\→}{\mathrm{\η}}}_k,(1-r(t,\mathrm{\β}))f(t,\mathrm{\β}))\·\·\·\left(1\right)$

上式中，f(t，β)为探测器获得的投影数据，β为投影角度，t为投影数据的采样点，g(t，β)为校正后的投影数据，r(t，β)表示每一条投影射线的路径上，高密度骨组织对投影数值的贡献，取值范围是(0，1)，P()为线性校正的多项式方程，

为多项式方程的系数。

为系数的多项式形式如下：

$P({\stackrel{\→}{\mathrm{\η}}}_k,x)={\mathrm{\η}}_{1}x+{\mathrm{\η}}_2{x}^2+{\mathrm{\η}}_{1/2}{x}^{1/2}+{\mathrm{\η}}_3{x}^3+{\mathrm{\η}}_{1/3}{x}^{1/3}+\·\·\·\·\·\·+{\mathrm{\η}}_L{x}^L+{\mathrm{\η}}_{1/L}{x}^{1/L}\·\·\·\left(2\right)$

而 为系数的多项式形式与

的类似：

$P({\stackrel{\→}{\mathrm{\α}}}_k,x)={\mathrm{\α}}_{1}x+{\mathrm{\α}}_2{x}^2+{\mathrm{\α}}_{1/2}{x}^{1/2}+{\mathrm{\α}}_3{x}^3+{\mathrm{\α}}_{1/3}{x}^{1/3}+\·\·\·\·\·\·+{\mathrm{\α}}_L{x}^L+{\mathrm{\α}}_{1/L}{x}^{1/L}\·\·\·\left(3\right)$

2.圆形扫描轨迹、扇形投影束下的H-L一致性条件表达式

投影示意图请见图(3)，推导的公式如下：

$m_{i,k}=\∫{\∫}_C{x}^i{y}^{k}u(x,y)\mathrm{dxdy},i\≥0,k\≥0\·\·\·\left(4\right)$

$v_d\left(\mathrm{\β}\right)={\∫}_{-\∞}^{+\∞}{g}_f(t,\mathrm{\β}-\mathrm{arctg}\frac{t}{D}){\left(\frac{\mathrm{tD}}{\sqrt{D^2+t^2}}\right)}^d\left(\frac{D^3}{{(D^2+t^2)}^{3/2}}\right)\mathrm{dt}\·\·\·\left(5\right)$

$Q_d\left(\mathrm{\β}\right)=\underset{r=0}{\overset{d}{\mathrm{\Σ}}}\left(\begin{array}{c}d\\ r\end{array}\right)m_{r,d-r}{\cos }^r\left(\mathrm{\β}\right){\sin }^{d-r}\left(\mathrm{\β}\right)\·\·\·\left(6\right)$

则有推广的H-L一致性条件：v_d(β)＝Q_d(β)                  (7)

其中，0≤d≤N-1，而N为投影角度的个数。

3.构造一种通用的投影数据非线性失真的校正方法

假定成像过程中，X光机绕原点旋转一周，共获得N幅锥束图像，图像序列经过预处理，可以得到扇束投影

其中β取值分别为β₁，β₂，…β_N。假定β下的扇束投影有M(M＞N)个均匀采样数据，

的离散化形式可以用矩阵表示为：

这里 是非线性变换后的结果，记真实的投影结果为：

不考虑噪声的影响，则有：

${\hat{g}}_{i,j}=T(\mathrm{\Θ},g_{i,j})(1\≤i\≤M,1\≤j\≤N)\·\·\·\left(10\right)$

其中Θ为非线性变换T的K维参数向量。

考察H-L一致性条件，有：

$\underset{i=1}{\overset{M}{\mathrm{\Σ}}}{g}_f(t_i,{\mathrm{\β}}_j-\mathrm{arctg}\frac{t_i}{D}){\left(\frac{t_{i}D}{\sqrt{D^2+{t_i}^2}}\right)}^d\left(\frac{D^3}{{(D^2+{t_i}^2)}^{3/2}}\right)$

$=\underset{r=0}{\overset{d}{\mathrm{\Σ}}}\left(\begin{array}{c}d\\ r\end{array}\right)m_{r,d-r}{\cos }^r\left({\mathrm{\β}}_j\right){\sin }^{d-r}\left({\mathrm{\β}}_j\mathrm{\β}\right)(1\≤j\≤N,0\≤d\≤N-1)\·\·\·\left(11\right)$

如果把G，Θ，m_r，d-r(d≤N-1)作为未知数，共有：MN+K+N(N+1)/2个变量，公式(10)与(11)联立则有MN+N²个方程，只要K≤N(N-1)/2，则方程组有解，也就是说可以求得Θ，完成旋转图像序列的非线性校正。

对于二维切片图像，假定从旋转序列图像中获得的扇束投影为 而真实的投影应该是G。由于非线性变换T是单调增函数，故T可逆，有：

$g_{i,j}=T^{-1}\left({\hat{g}}_{i,j}\right)(1\≤i\≤M,1\≤j\≤N)\·\·\·\left(12\right)$

这里T^-1表示T的逆变换，它也是一个单调增的非线性变换。只要我们能够确定T^-1，就可以完成矫正。然而T^-1是一个形式未知的函数，由函数逼近理论，它可以用某个函数空间上的基函数的组合进行逼近，只要阶数足够高，逼近的效果就足够好。

设S_k是 所张成的函数空间，即：

$S_k=\mathrm{span}\{1,x,x^2,\sqrt[2]{x},\·\·\·x^k,\sqrt[k]{x}\}\·\·\·\left(13\right)$

则有：S₂S₃…S_k…S_∞                      (14)

那么对于任何一个非线性函数T^-1，都可以用S_k中的函数来逼近，k越大，逼近效果越好。也就是说T^-1(x)可以写作：

$T^{-1}\left(x\right)=a_0+a_{1}x+a_2{x}^2+b_2\sqrt{x}+\·\·\·a_k{x}^k+b_k\sqrt[k]{x}+O\left(\min (x^k,\sqrt[k]{x})\right)\·\·\·\left(15\right)$

如果能够确定系数，便可以达到矫正的目的。由于T的定义知道，这里的a₀＝0。当2k＜N-1时，令：

$p(j,\mathrm{\γ})=\underset{i=1}{\overset{M}{\mathrm{\Σ}}}[{\hat{g}}_f(t_i,{\mathrm{\β}}_j-\mathrm{arctg}\frac{t_i}{D}){]}^{\mathrm{\γ}}\left(\frac{D^3}{{(D^2+{t_i}^2)}^{3/2}}\right)\·\·\·\left(16\right)$

X＝[a₁ a₂ b₂ … a_k b_k]^T                            (18)

 E＝[1 1 1 … 1 1]^T                                (19)

其中

可以通过插值得到，根据d＝0时的H-L一致性条件可得：

AX＝m_0，0E                                         (20)

公式(20)是一个2k元一次的线性方程组，由最小二乘法，可以得到其最优解：

                 X＝m_0，0(A^TA)^-1A^TE                (21)

当2k＞N-1，可以利用d＝0，1条件下的H-L约束条件，进行类似的推导。实际上，一般k很小(k≤4)，2k≤N-1的条件很容易满足。只要确定了参数，矫正工作便可以通过式(15)完成。必须指出系数m_0，0虽然是不确定的，但它可被灰度映射过程所吸收。

4.构造一种基于H-L一致性条件的射束硬化校正方法

为了获得准确的射束硬化效应校正方法，需要考虑在同一断面投影时，不同投影角度下校正方程系数的不同，以及纵向(longitudinal)不同断面的校正方程的不同。根据成像的物理模型，在医学X线诊断光谱范围内，人体物质对射线的衰减主要取决于康普顿散射衰减和光电效应衰减，其中的自由度只有二维。这意味着，任何物质对X线光子的衰减函数可以用两种不同物质衰减的线性组合来代替。考虑到每一种物质的射束硬化效应的校正可以用一个多项式函数来校正，因此，我们首先构造了一个可以精确校正射束硬化效应的校正方程(该方程针对扇束投影，锥束投影的校正方程也可以根据同样的原理来构造)，如公式(1)。所谓的校正过程就是通过H-L一致性条件来求解上述方程的未知量：

r(t，β)。

在图4所示的模体的基础上，分两种情况进行讨论。

1).r(t，β)＝0时的校正方法

这里对应的是现有CT设备相类似的校正方案，即在同一断层的投影数据中，认为不同角度的投影数据可以用同一系数的多项式校正。图4是利用我们提出的基于H-L一致性条件校正方法获得的图像，采用的重建方法是FBP，并采用R-L型RAMP滤波函数。为了说明校正效果，我们选择显示两条纵、横的剖线下，校正前后的数据对比，其中的两条亮线用于指示显示的剖线位置，结果如图5所示。

图5给出了水平剖线的校正前后重建数据的对比。图中的红色实线表示未经射束硬化校正的重建结果，兰色实线表示经过射束硬化校正的重建结果，兰色虚线表示校正前后数据的差值。从校正前后数据的差值可见对中间区域软组织的校正效果明显高于两边区域的软组织，这一方面反映了校正前重建图像中存在典型的射束硬化失真现象，另一方面表明H-L一致性条件校正方法的确对软组织区域的射束硬化失真起到校正效果。

上述r(t，β)＝0时的校正方案相当于现有CT设备校正方法的改进，其主要优点如下：

A.校正方程系数自适应获得，准确度更高，而且无需事先的模体实验，也就不需要定期的CT参数校正；

B.校正方程系数在纵向坐标(longitudinal)上自适应变化，基于人体对称的特征，这一点可以克服传统方法的不足，将具有较好的应用价值。

2).估计r(t，β)，求

的校正方法

该方法的关键在于需要首先求出r(t，β)，在此基础上，通过我们已经构造好的校正方法直接求得

估计r(t，β)可以这样进行：先假设r(t，β)＝0，进行校正后重建；然后在重建图像中，通过阈值分割的方法将高密度的骨性组织区域同低密度的软组织区域区分开；根据分割的结果和不同组织的物理特征来估计r(t，β)；在此基础上，进行校正和重建。

图6是利用申请人提出的基于图像分割和H-L一致性条件校正方法获得的重建图像，采用的重建方法是FBP，并采用R-L型RAMP滤波函数。为了说明校正效果，我们选择显示两条纵、横的剖线下，校正前后的数据对比，其中的两条亮线用于指示显示的剖线位置，结果如图7所示。

从图7(a)可以看出，在左右骨组织的内侧软组织密度值得到了明显改善，这一点优于图5(a)中的校正效果，这说明该方法对硬膜外伪影这一CT设备难以克服的问题，将会有较好的解决方案途径。

Claims (2)

1.一种基于HL一致性条件的CT投影数据射束硬化效应校正方法，其特征在于，该方法通过构建投影数据的等效非线性失真校正模型，构造符合CT投影模式的H-L一致性条件的表达公式，根据等效射束硬化失真校正模型和H-L一致性条件，基于不同角度的投影图像数据，通过对等效射束硬化失真校正模型参数的求解来实现对投影数据的校正，具体包括以下步骤：

步骤1：构建投影数据的等效非线性失真校正模型

对于医学CT设备，根据成像物理学模型，将观测到的投影数据分解为高密度组织、低密度组织分别衰减而成的具体模型，构造等效射束硬化失真校正模型：

$g(t,\mathrm{\β})=P({\stackrel{\→}{\mathrm{\α}}}_k,r(t,\mathrm{\β})f(t,\mathrm{\β}))+P({\stackrel{\→}{\mathrm{\η}}}_k,(1-r(t,\mathrm{\β}))f(t,\mathrm{\β}))---\left(1\right)$

上式中，f(t，β)为探测器获得的投影数据，β为投影角度，t为投影数据的采样点，g(t，β)为校正后的投影数据，r(t，β)表示每一条投影射线的路径上，高密度骨组织对投影数值的贡献，取值范围是(0，1)，P()为线性校正的多项式方程， 为多项式方程的系数，

为系数的多项式形式如下：

$P({\stackrel{\→}{\mathrm{\η}}}_k,x)={\mathrm{\η}}_{1}x+{\mathrm{\η}}_2{x}^2+{\mathrm{\η}}_{1/2}{x}^{1/2}+{\mathrm{\η}}_3{x}^3+{\mathrm{\η}}_{1/3}{x}^{1/3}+\·\·\·\·\·\·+{\mathrm{\η}}_L{x}^L+{\mathrm{\η}}_{1/L}{x}^{1/L}---\left(2\right)$

其中，

为系数的多项式形式与 的相同；

步骤2：构造符合CT投影模式的H-L一致性条件的表达公式

根据已知平行束投影的HL一致性性条件，推广至圆形扫描轨迹、扇形投影方式的表达公式是：

m_i，k＝∫∫_Cxⁱy^ku(x，y)dxdy i≥0，k≥0      (3)

$v_d\left(\mathrm{\β}\right)={\∫}_{-\∞}^{+\∞}(t,\mathrm{\β}-\mathrm{arctg}\frac{t}{D}){\left(\frac{\mathrm{tD}}{\sqrt{D^2+t^2}}\right)}^d\left(\frac{D^3}{{(D^2+t^2)}^{3/2}}\right)\mathrm{dt}---\left(4\right)$

$Q_d\left(\mathrm{\β}\right)=\underset{r=0}{\overset{d}{\mathrm{\Σ}}}\left(\begin{array}{c}d\\ r\end{array}\right)m_{r,d-r}{\cos }^r\left(\mathrm{\β}\right){\sin }^{d-r}\left(\mathrm{\β}\right)---\left(5\right)$

v_d(β)＝Q_d(β)                              (6)

其中，0≤d≤N-1，而N为投影角度的个数。

上列公式中，原始的断层图像为u(x，y)，R为其支撑半径，C＝{(x，y)|x²+y²≤R²}，β为扇型束的投影角度；

步骤3：根据等效非线性失真校正模型和H-L一致性条件，设计模型参数的求解矩阵，通过计算得到模型参数的求解，从而完成对投影数据的校正过程。

2.如权利要求1所述的方法，其特征在于：求解等效非线性失真校正模型的未知参数

r(t，β)，能够形成以下两种校正方法，具体方法如下：

2.1首先假设r(t，β)＝0，即假设射线穿透路径上绝大大部分物质均与人体的软组织相当，通过下列矩阵和运算方法求解 得到方法I的射束硬化校正结果：

在CT成像过程中，X光机绕原点旋转一周可以得到扇束投影

经过离散化后得到

其中1≤i≤M，1≤j≤N；

$p(j,\mathrm{\γ})=\underset{i=1}{\overset{M}{\mathrm{\Σ}}}{\left[{\hat{g}}_f(t_i,{\mathrm{\β}}_j-\mathrm{arctg}\frac{t_i}{D})\right]}^{\mathrm{\γ}}\left(\frac{D^3}{{(D^2+{t_i}^2)}^{3/2}}\right)---\left(7\right)$

X＝[a₁ a₂ b₂ … a_k b_k]^T            (9)

E＝[1 1 1 … 1 1]^T                 (10)

其中 可以通过插值得到，根据d＝0时的H-L一致性条件有：

AX＝m_0，0E                         (11)

X为未知量，可以根据许多数值求解方法，求得X的最优解；

X即

只要确定了该参数，矫正工作便能够通过式(2)完成，必须指出，系数m_0，0虽然是不确定的，但它可被灰度映射过程所吸收；

2.2在校正方法I得到的结果基础上完成CT图像的重建工作，再运用以下方法得到r(t，β)的计算结果：

假定CT图像射束硬化失***要是由骨性物质引起的，对已经得到的CT图像：

首先，经过简单的阈值图像分割(或者其他经过精心设计的图像分割方法)，将CT图像分割为骨性区域和非骨性区域；

然后，沿着(t，β)方向将非骨性区域的数值求和得Φ(t，β)，该过程也常称为重投影，计算公式为：

$\mathrm{\Φ}(t,\mathrm{\β})=\underset{i=1}{\overset{N_I}{\mathrm{\Σ}}}{f}_I(x_i,y_i){\mathrm{\Ψ}}_{(t,\mathrm{\β})}(x_i,y_i)---\left(12\right)$

其中N_I为整个CT重建图像的象素总数，f_I为校正方法I所得到的CT重建结果，Ψ_(t，β)为空间隶属函数，即

${\mathrm{\Ψ}}_{(t,\mathrm{\β})}(x,y)=\left\{\begin{array}{c}1,(x,y)\∈(M\∩l(t,\mathrm{\β}))\\ 0,(x,y)\∉(M\∩l(t,\mathrm{\β}))\end{array}\right.---\left(13\right)$

其中M表示分割后的软组织区域，l(t，β)表示特定X射线路径；

最后，计算

$r(t,\mathrm{\β})=\{{\hat{g}}_f(t,\mathrm{\β})-\mathrm{\Φ}(t,\mathrm{\β})\}/{\hat{g}}_f(t,\mathrm{\β});---\left(14\right)$

2.3在已知r(t，β)的基础上，通过下列运算方法，求解 和 得到公式(1)描述的方法II校正结果：

首先，根据步骤2.1中所描述的方法，以 为投影求解

以Φ(t，β)为投影求解

然后根据公式(1)，将

的矫正结果和Φ(t，β)的矫正结果相加，得到经过矫正后的投影；

对经过矫正后的投影进行重建得到经过矫正的CT图像；

最后，采用步骤2.2中的方法更新r(t，β)；

迭代完成以上步骤，直到结果收敛到一个确定值。

Images