CN103310432B - 基于四阶全变分流的ct图像归一化的金属伪影校正法 - Google Patents

Abstract

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

Description

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

技术领域

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

背景技术

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

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

发明内容

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

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

（1）原始扫描数据采集

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

（2）确定金属区域

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

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

（3）计算先验图像

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

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

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

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

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

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

原始扫描数据g除以先验图像的前向投影g_prior，得到归一化后的数据g_norm：

$g_{\mathrm{norm}}=\frac{g}{g_{\mathrm{prior}}+t_{\mathrm{eps}}}$

其中t_eps为正数以避免除数为零；

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

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

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

g_denorm=g_TV*(g_prior+t_eps)

其中t_eps取值和步骤（5）中保持一致;

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

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

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

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

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

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

$u_t=-\mathrm{\Δ}\▿\·\left(\frac{\▿u}{{\sqrt{|\▿u|}}^2+{\mathrm{\δ}}^2}\right)+\mathrm{\λ}(f-u)---\left(1\right)$

且

$\mathrm{\λ}(x,y)=\left\{\begin{array}{cc}0& \mathrm{\λ}(x,y)\∈g_{\mathrm{metal}}\\ {\mathrm{\λ}}_0& \mathrm{\λ}(x,y)\∈g_{\mathrm{norm}}\backslash g_{\mathrm{metal}}\end{array}\right.$

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

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

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

$E_1={\&Integral;}_{g_{\mathrm{prior}}}|\&dtri;u|\mathrm{dx}$

公式(1)中的后项也是梯度流，用L²范数表示其能量： $E_2=\mathrm{\&lambda;}{\&Integral;}_{g_{\mathrm{prior}}}{(f-u)}^2\mathrm{dx}$

然后对E₁进行分裂：E₁=E₁₁-E₁₂

$E_{11}={\&Integral;}_{g_{\mathrm{prior}}}\frac{C_1}{2}{|\&dtri;u|}^2\mathrm{dx}$

$E_{12}={\&Integral;}_{p_{\mathrm{prior}}}-|\&dtri;u|+\frac{C_1}{2}{|\&dtri;u|}^2\mathrm{dx}$

同理，E₂也可以分裂为：E₂=E₂₁-E₂₂

$E_{12}={\&Integral;}_{g_{\mathrm{prior}}}\frac{C_2}{2}{|\&dtri;u|}^2\mathrm{dx}$

$E_{22}=\frac{1}{2}{\&Integral;}_{g_{\mathrm{prior}}}-\mathrm{\&lambda;}{(f-u)}^2+C_2{|\&dtri;u|}^2\mathrm{dx}$

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

$\frac{U_{k+1}-U_k}{\mathrm{\&Delta;t}}+C_1\mathrm{\&Delta;\&Delta;}{U}_{k+1}+C_2{U}_{k+1}=C_1\mathrm{\&Delta;\&Delta;}{U}_k-\mathrm{\&Delta;}(\&dtri;\&CenterDot;\left(\frac{{\&dtri;U}_k}{\left|{\&dtri;U}_k\right|}\right))+C_2{U}_k+\mathrm{\&lambda;}(f-U_k)$

式中C₁>1/δ，C₂>λ₀，U_k表示经过上式第k次迭代后的结果，即图像g_norm在时间k△t后的结果，即可得到校正后的投影数据g_TV。

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

附图说明

图1为本发明的流程图；

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

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

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

具体实施方式

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

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

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

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

步骤2，确定金属区域

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

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

步骤3，计算先验图像

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

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

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

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

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

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

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

$g_{\mathrm{norm}}=\frac{g}{g_{\mathrm{prior}}+t_{\mathrm{eps}}}$

这里选用一个较小的正数t_eps来避免除数为零。

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

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

$u_t=-\mathrm{\&Delta;}\&dtri;\&CenterDot;\left(\frac{\&dtri;u}{{\sqrt{|\&dtri;u|}}^2+{\mathrm{\&delta;}}^2}\right)+\mathrm{\&lambda;}(f-u)---\left(2\right)$

且

$\mathrm{\&lambda;}(x,y)=\left\{\begin{array}{cc}0& \mathrm{\&lambda;}(x,y)\&Element;g_{\mathrm{metal}}\\ {\mathrm{\&lambda;}}_0& \mathrm{\&lambda;}(x,y)\&Element;g_{\mathrm{norm}}\backslash g_{\mathrm{metal}}\end{array}\right.$

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

$E_1={\&Integral;}_{g_{\mathrm{prior}}}|\&dtri;u|\mathrm{dx}$

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

$E_2=\mathrm{\&lambda;}{\&Integral;}_{g_{\mathrm{prior}}}{(f-u)}^2\mathrm{dx}$

然后对E₁进行分裂：E₁=E₁₁-E₁₂

$E_{11}={\&Integral;}_{g_{\mathrm{prior}}}\frac{C_1}{2}{|\&dtri;u|}^2\mathrm{dx}$

$E_{12}={\&Integral;}_{p_{\mathrm{prior}}}-|\&dtri;u|+\frac{C_1}{2}{|\&dtri;u|}^2\mathrm{dx}$

同理，E₂也可以分裂为：E₂=E₂₁-E₂₂

$E_{21}={\&Integral;}_{g_{\mathrm{prior}}}\frac{C_2}{2}{|\&dtri;u|}^2\mathrm{dx}$

$E_{22}=\frac{1}{2}{\&Integral;}_{g_{\mathrm{prior}}}-\mathrm{\&lambda;}{(f-u)}^2+C_2{|\&dtri;u|}^2\mathrm{dx}$

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

$\frac{U_{k+1}-U_k}{\mathrm{\&Delta;t}}+C_1\mathrm{\&Delta;\&Delta;}{U}_{k+1}+C_2{U}_{k+1}=C_1\mathrm{\&Delta;\&Delta;}{U}_k-\mathrm{\&Delta;}(\&dtri;\&CenterDot;\left(\frac{{\&dtri;U}_k}{\left|{\&dtri;U}_k\right|}\right))+C_2{U}_k+\mathrm{\&lambda;}(f-U_k)$

式中C₁>1/δ，C₂>λ₀，U_k表示经过上式第k次迭代后的结果，即图像g_norm在时间k△t后的结果。可得到校正后的投影数据g_TV。

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

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

g_denorm=g_TV*(g_prior+t_eps)

这里的数t_eps取值要和步骤5中保持一致。

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

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

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

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

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

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

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

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

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

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

Claims (4)

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)，并选择阈值h_m提取金属区域：

式中表示金属区域的像素值；

(3)计算先验图像

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

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

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

g_metal(n)＝{[s(i,n),e(i,n)]},i＝1,…l(n)

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

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

原始投影数据g除以先验图像的前向投影，得到归一化后的数据g_norm：

$g_{norm}=\frac{g}{g_{prior}+t_{eps}}$

其中t_eps为正数以避免除数为零；

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

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

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

g_denorm＝g_TV*(g_prior+t_eps)

其中t_eps取值和步骤(5)中保持一致；

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

对非金属部分的投影数据，用原投影数据g对g_denorm进行替代：

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

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

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

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

$u_t=-\mathrm{\&Delta;}\&dtri;\&CenterDot;\left(\frac{\&dtri;u}{{\sqrt{|\&dtri;u|}}^2+{\mathrm{\&delta;}}^2}\right)+\mathrm{\&lambda;}(f-u)---\left(1\right)$

且

$\mathrm{\&lambda;}(x,y)=\left\{\begin{array}{cc}0& \mathrm{\&lambda;}(x,y)\&Element;g_{metal}\\ {\mathrm{\&lambda;}}_0& \mathrm{\&lambda;}(x,y)\&Element;g_{norm}\backslash g_{metal}\end{array}\right.$

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

4.如权利要求3所述的基于四阶全变分流的CT图像归一化的金属伪影校正法，其特征在于：运用凸面分裂的思想计算公式(1)，其具体包括：

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

$E_1={\&Integral;}_{g_{prior}}|\&dtri;u|dx$

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

$E_2=\mathrm{\&lambda;}{\&Integral;}_{g_{prior}}{(f-u)}^{2}dx$

然后对E₁进行分裂：E₁＝E₁₁-E₁₂

$E_{11}={\&Integral;}_{g_{prior}}\frac{C_1}{2}|\&dtri;u{|}^{2}dx$

$E_{12}={\&Integral;}_{g_{prior}}-|\&dtri;u|+\frac{C_1}{2}|\&dtri;u{|}^{2}dx$

同理，E₂也可以分裂为：E₂＝E₂₁-E₂₂

$E_{21}={\&Integral;}_{g_{prior}}\frac{C_2}{2}|\&dtri;u{|}^{2}dx$

$E_{22}=\frac{1}{2}{\&Integral;}_{g_{prior}}-\mathrm{\&lambda;}{(f-u)}^2+C_2|\&dtri;u{|}^{2}dx$

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

$\frac{U_{k+1}-U_k}{\mathrm{\&Delta;}t}+C_1{\mathrm{\&Delta;\&Delta;U}}_{k+1}+C_2{U}_{k+1}=C_1{\mathrm{\&Delta;\&Delta;U}}_k-\mathrm{\&Delta;}(\&dtri;\&CenterDot;(\frac{\&dtri;U_k}{|\&dtri;U_k|}\left)\right)+C_2{U}_k+\mathrm{\&lambda;}(f-U_k)$

式中C₁＞1/δ，C₂＞λ₀，U_k表示经过上式第k次迭代后的结果，即图像g_norm在时间kΔt后的结果，即可得到校正后的投影数据g_TV。