基于四阶全变分流的CT图像归一化的金属伪影校正法

摘要

基于四阶全变分流的CT图像归一化的金属伪影校正的方法，其采用基于凸性分裂方法进行求解，以实现金属伪影的校正。具体实现步骤包括：采集数据、重建图像、分割金属区域、计算先验图像、获取金属区域和先验图像的重投影、归一化投影数据、四阶全变分方程校正、反归一化校正后的数据、再次重建图像、对图像金属区域恢复。其具有能够有效的去除金属伪影，并很好的保留金属及其周围的结构信息，最大程度抑制了二次伪影的出现等优点。

权利要求

1.基于四阶全变分流的CT图像归一化的金属伪影校正法，其特征在于，包括以下步骤：(1)原始扫描数据采集

利用微计算机断层成像系统采集360°的生物体投影数据，得到扇束投影数据g(m,n),m＝1,…M,n＝1,…N，其中M是探测器总个数，N是投影角度总数；

(2)确定金属区域

利用步骤(1)中得到的扇束投影数据g(m,n)，对投影数据用FBP重建，得到含有金属伪影的CT图像fΩ(x,y)，并选择阈值hm提取金属区域：式中 表示金属区域的像素值；

(3)计算先验图像

对步骤(2)中重建的图像进行高斯平滑滤波，然后通过阈值分割出空气区域、软组织区域和骨骼区域，将空气区域的CT值设定为-1000HU，软组织区域的CT值设定为0HU，骨骼区域的CT值保持不变，金属区域的CT值用其周围的值来填充；

(4)对金属区域和先验图像进行重投影

利用Joseph投影法分别对金属区域和先验图像进行前向投影，得到每个投影角度下的重投影即正弦图，确定金属投影数据的范围gmetal和先验图像的前向投影gprior；

gmetal(n)＝{[s(i,n),e(i,n)]},i＝1,…l(n)其中l(n)是金属区域在角度n下总的投影点数，s(i,n),e(i,n)是金属投影数据在每个角度下的起始和结束位置，先验图像的前向投影gprior，用来对原始扫描数据归一化；

(5)原始扫描数据的归一化

原始投影数据g除以先验图像的前向投影，得到归一化后的数据gnorm：其中teps为正数以避免除数为零；

(6)对金属区域投影数据进行校正

在步骤(5)中得到了归一化的投影数据gnorm，根据步骤(4)中确定的金属投影数据的范围gmetal，采用四阶TV方程对归一化数据gnorm在金属投影数据的范围gmetal进行数据校正，得到校正后的投影数据gTV；

(7)进行反归一化，其包括：

gdenorm＝gTV*(gprior+teps)其中teps取值和步骤(5)中保持一致；

(8)重建图像并恢复图像金属信息步骤，其包括：对非金属部分的投影数据，用原投影数据g对gdenorm进行替代：并使用FBP重建，得到图像 将原图像fΩ(x,y)中的金属区域取代 金属区：

图像f′(x,y)即为校正后的图像。

2.根据权利要求1所述的基于四阶全变分流的CT图像归一化的金属伪影校正法，其特征在于：步骤(6)中对金属投影区域校正时利用了步骤(3)中所获得的原重建图像计算得到的先验信息。

3.如权利要求1所述的基于四阶全变分流的CT图像归一化的金属伪影校正法，其特征在于：所述四阶TV方程校正具体包括：且

f为需要校正的图像，即为gnorm，需要校正的区域为gmetal，u为迭代过程中校正图像f的变化结果，δ为平滑参数0＜δ＜＜/1。

4.如权利要求3所述的基于四阶全变分流的CT图像归一化的金属伪影校正法，其特征在于：运用凸面分裂的思想计算公式(1)，其具体包括：当λ＝0时，公式(1)为梯度流用H-1范数表示其能量：2

公式(1)中的后项也是梯度流，用L范数表示其能量：然后对E1进行分裂：E1＝E11-E12同理，E2也可以分裂为：E2＝E21-E22根据以上讨论的分裂思想可得到以下时步方法：式中C1＞1/δ，C2＞λ0，Uk表示经过上式第k次迭代后的结果，即图像gnorm在时间kΔt后的结果，即可得到校正后的投影数据gTV。

说明书

基于四阶全变分流的CT图像归一化的金属伪影校正法

技术领域

[0001] 本发明属于医学影像处理领域，涉及CT图像处理领域，特别是一种金属伪影校正方法。

背景技术

[0002] 目前，CT的应用越来越广泛，随之而来的是一系列的新问题，重建图像中的金属伪影就是其中之一。金属伪影在图像中表现为明暗相间的条状伪影，这些伪影降低图像质量，严重影响图像的后续处理和分析。金属伪影的形成主要是由于射束硬化。CT实际产生的X射线覆盖了一个很宽的频谱，而大多数材料的线性衰减系数又是和能量相关的，当具有多能谱的X射线射向一软组织时，能量较低的光子比能量较高的光子被吸收得更多，导致高能光子的比率变大，X射线变“硬”。射线通过金属物质这种高吸收系数的物质后，能量大大降低，这时射束硬化效果会越来越明显，最后的透射强度位于数据电子设备动态范围的底部，导致投影数出现跳跃变化，经过滤波处理后，这种弱连续性进一步扩大，就在图像中形成明暗相间的“条状”伪影。

[0003] 在此之前，人们也提出了很多方法来校正金属伪影。这些方法基本上可以分为以下几类：正弦图修补法、迭代法、统计学方法和滤波法。最为常用的是正弦图修补法，首先确定金属投影区域，将金属影响的数值视为丢失数据，然后对它们进行修补，最后重建图像。这类方法较为简单，便于实际操作，修补的方法不同，也会带来二次的伪影。迭代法能有效去除金属伪影和抑制噪声，而且能很好的呈现金属物体的结构，但其运算量非常大，速度很慢，难以实用化。统计学方法和滤波法对噪声不太敏感，可很好呈现金属物体的结构，但操作复杂，在现实医用CT中有很大局限性。

发明内容

[0004] 本发明的目的是针对当前主流算法存在的不足，提出一种新的金属伪影校正方法，可以保留金属结构的边缘信息，使正弦图修补过程中原始数据和修补数据之间足够平滑，最大程度抑制重建后金属结构周围的二次伪影。

[0005] 其包括，基于四阶全变分流的CT图像归一化的金属伪影校正法，基于四阶全变分流的CT图像归一化的金属伪影校正法，包括以下步骤：

[0006] （1）原始扫描数据采集

[0007] 利用微计算机断层成像系统采集360°的生物体投影数据，得到正弦图g(m,n),m=1,…M,n=1,…N，其中M是探测器总个数，N是投影角度总数；

[0008] （2）确定金属区域

[0009] 利用步骤（1）中得到的扇束投影数据g(m,n)，对投影数据用FBP重建，得到含有金属伪影的CT图像fΩ(x,y)，并选择阈值hm提取金属区域：

[0010]

[0011] 式中 (x,y)表示金属区域的像素值；

[0012] （3）计算先验图像

[0013] 对步骤（2）中重建的图像进行高斯平滑滤波，然后通过阈值分割出空气区域、软组织区域和骨骼区域，将空气区域的CT值设定为-1000HU，软组织区域的CT值设定为0HU，骨骼区域的CT值保持不变，金属区域的CT值用其周围的值来填充；

[0014] （4）对金属区域和先验图像进行重投影

[0015] 利用Joseph投影法分别对金属区域和先验图像进行前向投影，得到每个投影角度下的重投影即正弦图，确定金属投影数据的范围gmetal和先验图像的前向投影gprior；

[0016] gmetal(n)={[s(i,n),e(i,n)]},i=1,…l(n)

[0017] 其中l(n)是金属区域在角度n下总的投影点数，s(i,n),e(i,n)是金属投影数据在每个角度下的起始和结束位置，先验图像的前向投影gprior(m,n),m=1,…M,n=1,…N，用来对原始扫描数据归一化；

[0018] （5）原始扫描数据的归一化

[0019] 原始扫描数据g除以先验图像的前向投影gprior，得到归一化后的数据gnorm：

[0020]

[0021] 其中teps为正数以避免除数为零；

[0022] （6）对金属区域投影数据进行校正

[0023] 在步骤（5）中得到了归一化的投影数据gnorm，根据步骤（4）中确定的金属重投影数据区域gmetal，采用四阶TV方程对归一化数据gnorm在金属投影区域gmetal进行数据校正，得到校正后的投影数据gTV；

[0024] (7)进行反归一化，其包括：

[0025] gdenorm=gTV*(gprior+teps)

[0026] 其中teps取值和步骤（5）中保持一致;

[0027] (8)重建图像并恢复图像金属信息步骤，其包括：

[0028] 对非金属部分的投影数据，用原投影数据g对校正后的投影数据gdenorm进行替代：

[0029]

[0030] 并使用FBP重建，得到图像 将原图像fΩ(x,y)中的金属区域取代金属区：

[0031]

[0032] 图像f′(x,y)即为校正后的图像。

[0033] 在上述技术方案耳朵基础上，步骤（6）中对金属投影区域校正时利用了步骤（3）中所获得的原重建图像计算得到的先验信息。

[0034] 在上述技术方案耳朵基础上，所述四阶TV方程校正具体包括：

[0035]

[0036] 且

[0037]

[0038] f为需要校正的图像，即为gnorm，需要校正的区域为gmetal，u为迭代过程中校正图像f的变化结果，δ为平滑参数0<δ＜＜1。

[0039] 在上述技术方案耳朵基础上，运用凸面分裂的思想计算公式(1)，其具体包括：

[0040] 当λ=0时，公式(1)为梯度流用H-1范数表示其能量：

[0041]2

[0042] 公式(1)中的后项也是梯度流，用L范数表示其能量：

[0043] 然后对E1进行分裂：E1=E11-E12

[0044]

[0045]

[0046] 同理，E2也可以分裂为：E2=E21-E22

[0047]

[0048]

[0049] 根据以上讨论的分裂思想可得到以下时步方法：

[0050]

[0051] 式中C1>1/δ，C2>λ0，Uk表示经过上式第k次迭代后的结果，即图像gnorm在时间k△t后的结果，即可得到校正后的投影数据gTV。

[0052] 相对于现有技术，本发明采用归一化对扫描的原始投影数据预处理，可有效保留金属结构的边缘信息，尤其对骨骼中含有金属的情况，本发明处理效果更优于其他方法。同时，本发明采用四阶TV方程对金属投影区域数据进行校正，可有效的去除金属伪影，并很好的保留金属及周围的结构信息，可有效的抑制二次伪影。

附图说明

[0053] 图1为本发明的流程图；

[0054] 图2为需要扫描重建的理性模型原图；

[0055] 图3为未经金属伪影校正的仿真重建结果；

[0056] 图4为本发明实施案例的仿真重建结果。

具体实施方式

[0057] 下面结合附图具体实施例来详细描述本发明的技术方案。

[0058] 参照图1介绍一下本发明的主要步骤：

[0059] 步骤1，原始扫描数据采集

[0060] 通过软件控制旋转台旋转小角度，一般为不大于1°，微计算机断层成像系统采集360°的生物体投影数据。利用预处理软件对计算机断层成像投影数据进行补偿坏点坏线、亮场暗场校正、几何校正预处理，即得到正弦图g(m,n),m=1,…M,n=1,…N，M是探测器总个数，N是投影角度总数。

[0061] 步骤2，确定金属区域

[0062] 利用步骤1中得到的扇束投影数据g(m,n)，对投影数据直接用FBP重建，得到含有金属伪影的CT图像fΩ(x,y)，选择适当的阈值hm提取金属区域：

[0063]

[0064] 式中 表示金属区域的像素值。

[0065] 步骤3，计算先验图像

[0066] 先验图像是一个不包含伪影而又特别接近重建图像的模型图像。为了获取先验图像，要在步骤2中重建的图像g(m,n)上划分出空气区域、软组织区域和骨骼区域。先对重建图像g(m,n)进行高斯平滑滤波，然后通过适当的阈值分割出空气区域、软组织区域和骨骼区域。将空气区域的CT值设定为-1000HU，软组织区域的CT值设定为0HU，骨骼区域的CT值保持不变，而金属区域的CT值用其周围的值来填充，因为之后正弦图修补利用的是金属周边的数据，且金属区域的数据是用不到的，同时这样也不会影响之后的归一化和正弦图的修补。

[0067] 步骤4，对金属区域和先验图像进行重投影

[0068] 对通过阈值分割出来的金属区域进行重投影，确定金属投影数据的区域，这样才能对该范围内的投影数据进校正。利用Joseph投影法或距离驱动投影法分别对金属区域和先验图像进行前向投影，得到每个投影角度下的重投影即正弦图，确定金属投影数据的范围gmetal和先验图像的前向投影gprior。

[0069] gmetal(n)={[s(i,n),e(i,n)]},i=1,…l(n)

[0070] l(n)是金属区域在角度n下总的投影点数，s(i,n),e(i,n)是金属投影数据在每个角度下的起始和结束位置。先验图像的前向投影gprior(m,n),m=1,…M,n=1,…N，用来对原始扫描数据归一化。

[0071] 步骤5，原始扫描数据的归一化

[0072] 为了使得原扫描数据更为平坦，需要对其归一化。原始扫描数据g除以先验图像的前向投影gprior，其中除法是针对每个像素而言的，即可得到归一化后的数据gnorm：

[0073]

[0074] 这里选用一个较小的正数teps来避免除数为零。

[0075] 步骤6，对金属区域投影数据进行校正

[0076] 步骤5中得到了归一化的投影数据gnorm，根据步骤4中确定的金属重投影数据区域gmetal，对归一化数据gnorm在金属投影区域gmetal进行数据校正。本发明中运用的是四阶TV方程校正：

[0077]

[0078] 且

[0079]

[0080] f为需要校正的图像，本发明中即为gnorm，需要校正的区域为gmetal，u为迭代过程中校正图像f的变化结果，δ为平滑参数0<δ＜＜1。运用凸面分裂的思想可快速计算公式-1(2)，当λ=0时，公式(1)为梯度流用H 范数表示其能量：

[0081]

[0082] 公式(2)中的后项也是梯度流，用L2范数表示其能量：

[0083]

[0084] 然后对E1进行分裂：E1=E11-E12

[0085]

[0086]

[0087] 同理，E2也可以分裂为：E2=E21-E22

[0088]

[0089]

[0090] 根据分裂思想可得到以下时步方法：

[0091]

[0092] 式中C1>1/δ，C2>λ0，Uk表示经过上式第k次迭代后的结果，即图像gnorm在时间k△t后的结果。可得到校正后的投影数据gTV。

[0093] 步骤7，对校正后的投影数据进行反归一化

[0094] 步骤6中的校正过程是在归一化的数据上进行的，为了得到实际的投影数据，需要对校正后的投影数据gTV进行反归一化，得到gdenorm：

[0095] gdenorm=gTV*(gprior+teps)

[0096] 这里的数teps取值要和步骤5中保持一致。

[0097] 步骤8，重建图像并恢复图像金属信息

[0098] 由步骤7得到的投影数据gdenorm虽然对金属区域的投影数据进行了校正，同时也对非金属区域的数据进行了一定程度的平滑。对于非金属部分的投影数据，要用原投影数据g对校正后的投影数据gdenorm进行替代：

[0099]

[0100] 然后直接使用FBP重建，得到图像 并恢复图像金属部分信息。在图像中，金属伪影得到基本消除，但是却丢失了金属的信息，所以需要对图像金属信息进行恢复。只需要将原图像fΩ(x,y)中的金属区域取代 金属区域即可：

[0101]

[0102] 图像f′(x,y)即为校正后的图像。

[0103] 下面结合附图2、附图3和附图4对本发明的重建结果做进一步的描述。

[0104] 本案例中，采用扇束解析法模拟投影数据，投影角度为0.72°，共500个投影度，探测器数目为730个，重建图像大小为512*512。

[0105] 附图2需要扫描重建的理性模型原图。其中，白色部分为金属物体。

[0106] 附图3为未经金属伪影校正的仿真重建结果。其中，金属伪影明显，图像质量非常差，而且金属结构不清晰。旋转角度为360°，采集的投影数据的个数是500个。

[0107] 图4为本发明实施案例的仿真重建结果。其中，金属伪影消失，图像质量提高，而且金属结构比较清晰。旋转角度为360°，采集的投影数据的个数是500个。

[0108] 将本发明的重建效果附图4和未经金属伪影校正的重建结果附图3、需要扫描重建的理性模型附图2进行比较，可以看出重建结果中的金属伪影得到了校正，可以清晰的得到了实验模型的重建结果，并且金属及其周围的结构变得清晰，说明本发明克服了现有技术校正后丢失金属周围结构信息的缺点，有效的完成计算机断层成像的金属伪影校正。