一种使用能量函数方法的锥束CT中杯状伪影的校正方法

摘要

本发明公开了一种使用通过最优化能量函数对锥束CT中杯状伪影进行校正的方法，该方法针对重建后的切片图像，直接面向用户，不对原有锥束CT设备进行任何改动，就可以完成校正工作，能够高效地进行锥束CT的杯状伪影校正的同时还能够提高重建图像的同种物质的CT值均匀性，从而有助于重建图像中，完善的体可视化和基于阈值的可视化技术的发展。该方法应用于锥束CT重建图像(即图像域)校正技术领域。

说明书

技术领域

本发明涉及医学图像处理技术领域，尤其涉及锥束CT图像杯状伪影校正、灰度不均匀校正医学图像处理技术领域。

背景技术

锥束CT作为近年来发展的医疗及工业检测仪器，常用于图像引导治疗、上腹部检查、口腔检查、工业检测等方面。基于平板探测器的锥束CT(CBCT)与传统的二维CT相比，具有突出的优点，主要表现在锥束CT一次圆周扫描周期内，可以得到完成数百甚至上千个断层图像的投影，具有更高的扫描速度和辐射利用率，并有效的减少X射线管的负载输出，降低扫描成本。影响锥束CT重建图像质量的因素有很多，如X线散射、噪声、几何误差、能谱、探测单元响应不一致等。但由于锥束平板CT使用大范围的X射线平板探测器，这使得成像质量与传统CT相比较更易受到X射线散射及射束硬化的影响。因散射和射束硬化而形成的伪影(主要包括杯状伪影和条纹状伪影)严重影响对重建图像的分析与判断。在医疗级别的锥束CT重建图像中，杯状伪影占有很大比重，这些伪影对于基于阈值的可视化显示方面和基于阈值的锥束CT图像分割方面影响非常严重。且杯状伪影的校正可以为其他伪影校正提供反馈参考，为其他无先验的锥束CT散射及射束硬化校正提供验证信息。因此本发明关于锥束CT中的杯状伪影的校正非常有意义。

为了减少杯状伪影(即：CT值不均匀性伪影)的影响，目前，现有技术或文献资料的研究主要集中在对投影图像上的散射校正。早期的伪影校正主要体现在基于硬件的校正，如X射线滤线器、准直器或金属栅格、空气隙方法、扫描狭缝技术和铅条或铅板技术等。最近几年伪影校正研究主要体现在基于蒙特卡洛方法、散射分析估计方法和基于部分散射射线测量的散射校正方法。蒙特卡洛仿真在CBCT散射校正中是很有效的方法，但是计算量巨大。近年来，一些改进的蒙特卡洛仿真算法也被提出来，如使用GPU加速技术,基于模型的体恢复方法等。这些基于蒙特卡洛仿真的思想，都试图在模拟精度和计算代价上建立一个较好的平衡点。但始终局限于计算代价太高，而不易于使用。

基于部分射线遮挡的锥束CT伪影校正的方法，早期的有BSA散射校正板方法，是通过测量在射线遮挡阵列下方的散射量，来插值出到达探测器上的整体散射分布。然后再进行不含射线遮挡阵列的正常扫描，从扫描的投影图像中减去散射分布图像，就得到了校正后的图像。这种方法要进行两次扫描，增加了扫描时间的同时也增加了X射线照射剂量。后来出现了可移动遮挡快的方法，能解决两次扫描问题。

也有对在锥束CT重建后图像进行伪影校正的研究。该研究主要借助CT图像中的解剖结构。这类方法完全依赖变形配准精度，且需要CT图像数据。

现有技术的缺点主要包括：

(1)目前现有技术主要是针对投影图像校正，没有直接针对重建后切片图像的校正，如专利CN104408753A。

(2)现有技术大多数集中在因散射而造成的伪影校正的方法中，并且大多需要添加硬件设备，如：专利200710019084和201310039298，这两个专利都需要在昂贵的锥束CT设备上添加硬件设备，增加了操作的复杂性和对设备造成潜在的安全风险。特别是专利200710019084需要两次扫描被测物体，这样无疑增加的被测物的辐射量。

综上所述，在现有技术或文献资料的方法中，蒙特卡洛模拟方法非常耗费时间，初级射线调制方法中校正结果受限于调制板自身的结构，基于部分散射射线测量方法，需要增加照射剂量，现有方法对散射分布的估计准确度不高。而本发明能够很好地解决上面的问题。

发明内容

本发明目的在于解决了上述现有技术的不足，提出了一种使用通过最优化能量函数对锥束CT中杯状伪影进行校正的方法。本方法针对重建后的切片图像，直接面向用户，不对原有锥束CT设备进行任何改动，就可以完成校正工作，能够高效地进行锥束CT的杯状伪影校正的同时还能够提高重建图像的同种物质的CT值均匀性，从而有助于重建图像中，完善的体可视化和基于阈值的可视化技术的发展。该方法应用于锥束CT重建图像(即图像域)校正技术领域。

本发明解决其技术问题所采取的技术方案是：本发明提供了一种使用锥束CT中重建图像进行杯状伪影校正的最优化能量函数方法，该方法具有很强的鲁棒性，不需要重复扫描被测物体，不增加锥束CT系统的复杂度。

方法流程：

步骤1：获取重建后的切片图像。

步骤2：杯状伪影表示并构建能量函数

$>F(f_s,f_p)=F(u,c,w)=\underset{\Omega }{\int }|f\left(x\right)-\underset{i=1}{\overset{N}{\Sigma }}{c}_i{u}_i\left(x\right)-w^{T}G\left(x\right){|}^{2}dx,$>

其中G(x)＝(g₁(x),…g_M(x))^T为平滑基函数，c_i为常数，满足当x∈Ω_i时，u_i(x)＝1；当时，u_i(x)＝0，w＝(w₁,…w_M)

步骤3：固定的c和u，通过解方程来获得F(u,c,w)的最小值。得到$>\hat{w}=A^{-1}v.$>

步骤4：固定并使用更新的w和u，以为变量的F(u,c,w)最小值解为：

$>{\hat{c}}_i=\frac{\underset{\Omega }{\int }\left(f\right(x)-f_s(x\left)\right)u_i\left(x\right)dx}{\underset{\Omega }{\int }{u}_i^2\left(x\right)dx},i=1,...,N$>

步骤5：固定并使用更新的w和c，以u＝(u₁,…u_N)^T为变量的F(u,c,w)最小值解时，满足如下条件：

$>{\hat{u}}_i=\left(\begin{array}{c}1,i=i_{\min }\left(x\right)\\ 0,i\ne i_{\min }\left(x\right)\end{array}\right)$>

其中，i_min(x)＝argmin{f(x)-c_i-w^TG(x)}。

步骤6：若w稳定或迭代次数超过10次,则执行步骤7，否则回到步骤3。

步骤7：校正后图像为$>f_p=f-{\hat{f}}_s=f-w^{T}G\left(x\right).$>

进一步的，本发明直接面向锥束CT切片数据，不对原有锥束CT原有设备进行任何改动，不需要用户和被测目标的先验信息。

进一步的，本发明证明了杯状伪影是从重建图像分解出来，即：

$>\begin{array}{c}f=\frac{1}{4{\pi }^2}{\int }_0^{2\pi }\frac{{d_{so}}^2}{{\left(d_{so}+r·s\right)}^2}{\int }_{-\infty }^{\infty }\frac{d_{so}}{{\left({d_{so}}^2+t^2+z^2\right)}^{1/2}}·\\ \left(P_{3D}\left(t,z\left(r\right),\phi \right)+S_{3D}\left(t,z\left(r\right),\phi \right)+n\right)·{\int }_{-\infty }^{\infty }\left|\omega \right|e^{j2\pi \omega \left(t\left(r\right)-t\right)}d\omega dtd\phi \\ =f_p+f_s+f_n\end{array},$>

其中f_p表示真实无伪影的切片图像,f_s表示散射和射束硬化造成伪影的切片图像，f_n表示噪声n形成的切片图像。

有益效果：

1、本发明是直接针对重建后的切片图像的杯状伪影校正，该方法计算量相对较小，能够高效的进行锥束CT切片图像杯状伪影校正的同时，提高了同种物质重建图像的CT值均匀性。

2、本发明直接面向CT切片需求用户，不对原有锥束CT原有设备进行任何改动，不需要用户和被测目标的先验信息，很好地完成校正工作。

3、本发明能够增加图像对比度，使校正后图像能更准确的表现被照物体原本信息。

4、本发明很好地改善并且实现了基于阈值的CT图像可视化、分割及病灶检测。

附图说明

图1为本发明的方法流程图。

图2为本发明由重建图像校正前(左)后(右)的人头骨样本轴向视图。

图3为本发明测量的人头骨模体的水平剖面值：纵轴是图像剖面。

图4为本发明人头骨模体选择的感兴趣区域图像示意图。

图5(a)、图5(b)为本发明散射校正前后CatPhan500中CTP486重建图像两个不同切片的样本轴向视图。

图6为本发明中图5所测量模体的水平剖面图。

图7为小鼠骨的选择计算区域图像示意图。

图8(a)、图8(b)分别为散射校正前后小鼠骨重建图像的样本轴向视图。

具体实施方式

下面结合说明书附图对本发明创造作进一步的详细说明。

如图1所示，本发明提供了一种使用能量函数方法的锥束CT中杯状伪影的校正方法，该方法包括如下步骤：

步骤1：根据FDK算法获取重建后的切片图像，且能从理论证明杯状伪影可以从重建后图像分解出来。

步骤2：杯状伪影表示并构建能量函数

$>F(f_s,f_p)=F(u,c,w)=\underset{\Omega }{\int }|f\left(x\right)-\underset{i=1}{\overset{N}{\Sigma }}{c}_i{u}_i\left(x\right)-w^{T}G\left(x\right){|}^{2}dx,$>

其中G(x)＝(g₁(x),…g_M(x))^T为平滑基函数，c_i为常数，满足当x∈Ω_i时，u_i(x)＝1；当时，u_i(x)＝0，w＝(w₁,…w_M)

步骤3：固定的c和u，通过解方程来获得F(u,c,w)的最小值。得到$>\hat{w}=A^{-1}v.$>

步骤4：固定并使用更新的w和u，以为变量的F(u,c,w)最小值解为：

$>{\hat{c}}_i=\frac{\underset{\Omega }{\int }\left(f\right(x)-f_s(x\left)\right)u_i\left(x\right)dx}{\underset{\Omega }{\int }{u}_i^2\left(x\right)dx},i=1,...,N$>

步骤5：固定并使用更新的w和c，以u＝(u₁,…u_N)^T为变量的F(u,c,w)最小值解时，满足如下条件：

$>{\hat{u}}_i=\left(\begin{array}{c}1,i=i_{\min }\left(x\right)\\ 0,i\ne i_{\min }\left(x\right)\end{array}\right)$>

其中，i_min(x)＝argmin{f(x)-c_i-w^TG(x)}。

步骤6：若w稳定或迭代次数超过10次,则执行步骤7，否则回到步骤3。

步骤7：校正后图像为$>f_p=f-{\hat{f}}_s=f-w^{T}G\left(x\right).$>

本发明杯状伪影可以从重建图像分解出来详细过程包括：

锥束CT中重建图像以FDK算法为基础，重建图像集f可以写为如下形式：

$>\begin{array}{c}f=\frac{1}{4{\pi }^2}{\int }_0^{2\pi }\frac{{d_{so}}^2}{{\left(d_{so}+r·s\right)}^2}{\int }_{-\infty }^{\infty }\frac{d_{so}}{{\left({d_{so}}^2+t^2+z^2\right)}^{1/2}}\\ ·I_{3D}\left(t,z\left(r\right),\phi \right)·{\int }_{-\infty }^{\infty }\left|\omega \right|e^{j2\pi \omega \left(t\left(r\right)-t\right)}d\omega dtd\phi \end{array},$>公式1

其中d_so表示射线源到旋转轴的距离，I_3D(t,z(r),φ)表示投影图像的序列。投影图像I_3D可分解为如下形式：

I_3D＝P_3D+S_3D+n，公式2

其中P_3D表示理论真实投影图像，S_3D是由散射和射束硬化造成的伪影成分，n是均值为零的加性噪声。公式1表示为如下形式：

$>\begin{array}{c}f=\frac{1}{4{\pi }^2}{\int }_0^{2\pi }\frac{{d_{so}}^2}{{\left(d_{so}+r·s\right)}^2}{\int }_{-\infty }^{\infty }\frac{d_{so}}{{\left({d_{so}}^2+t^2+z^2\right)}^{1/2}}·\\ \left(P_{3D}\left(t,z\left(r\right),\phi \right)+S_{3D}\left(t,z\left(r\right),\phi \right)+n\right)·{\int }_{-\infty }^{\infty }\left|\omega \right|e^{j2\pi \omega \left(t\left(r\right)-t\right)}d\omega dtd\phi \\ =f_p+f_s+f_n\end{array},$>公式3

其中f_p表示真实无伪影的切片图像,f_s表示散射和射束硬化造成伪影的切片图像，f_n表示噪声n形成的切片图像。

由公式3可知，重建图像可以表示为三个独立成分相加。其中f_s是由S_3D得到，而S_3D是一个光滑的低频投影信号，所有f_s主要区域(即表现为大面积的同种物质区域，杯状伪影是其主要成分)也应该是一个光滑的低频切片图像。为了有效运用f_s和f_p的性质，本发明将f_s表示为一个已知的平滑基函数集g₁,…g_M的线性组合，这适应杯状伪影的平稳变化性质。通过寻找线性组合中的最佳系数w＝(w₁,…w_M)来估计杯状伪影(本发明列举的实验数据中取M＝10)。本发明把f_s(x)表示为f_s(x)＝w^TG(x)的向量形式，其中G(x)＝(g₁(x),…g_M(x))^T。

假设在图像域Ω_i中存在N种组织，则真实切片图像f_p(x)实际上在第i个组织中是关于x的一个常数c_i。每个Ω_i都可以用它的从属函数来u_i表示。在理想条件下，u_i是一个二进制从属函数，满足当x∈Ω_i时，u_i(x)＝1；当时，u_i(x)＝0。当从属函数u_i和常数c_i已知时，f_p可被表示为如下形式：

$>f_p=\underset{i=1}{\overset{N}{\Sigma }}{c}_i{u}_i\left(x\right),$>公式4

能量泛函表达式包括：

在本模型中，本发明考虑如何在重建图像f中，获得使下述能量函数F最小化时的f_s和f_p：

$>F(f_s,f_p)=\underset{\Omega }{\int }|f\left(x\right)-f_p\left(x\right)-f_s\left(x\right){|}^{2}dx,$>公式5

显然如果变量f_s和f_p没有任何约束条件，F的最小化是一个病态问题。实际上，当f_s和f_p为满足条件f_p＝f-f_s的任意值时，能量F(f_s,f_p)都可取得最小值。使用真实图像和杯状伪影图像的表达式，能量F(f_s,f_p)可以被表示为如下形式：

$>F(f_s,f_p)=F(u,c,w)=\underset{\Omega }{\int }|f\left(x\right)-\underset{i=1}{\overset{N}{\Sigma }}{c}_i{u}_i\left(x\right)-w^{T}G\left(x\right){|}^{2}dx,$>公式6

最优化求解包括：

能量F(u,c,w)的所有变量u＝(u₁,…u_N)、c＝(c₁,…c_N)和w都是凸的，这个性质确保了F(u,c,w)对于任何变量都有唯一的最优解。通过交替计算不同变量下F(u,c,w)最小值解，即可达到求解目的。

首先，固定的c和u，本发明可以通过解方程来获得F(u,c,w)的最小值。过程如下：

$>\frac{\partial F}{\partial w}=-2\underset{\Omega }{\int }G\left(x\right)|f\left(x\right)-\underset{i=1}{\overset{N}{\Sigma }}\left(c_i{u}_i\right(x\left)\right)|dx+2w\underset{\Omega }{\int }G\left(x\right)G{\left(x\right)}^{T}dx=0,$>公式7

上述等式可以改写为如下形式：

Aw＝v，公式8

其中$>A=\underset{\Omega }{\int }G\left(x\right)G{\left(x\right)}^{T}dx,v=\underset{\Omega }{\int }G\left(x\right)|f\left(x\right)-\underset{i=1}{\overset{N}{\Sigma }}\left(c_i{u}_i\right(x\left)\right)|dx.$>

容易证明，矩阵A是非奇异矩阵。因此，向量可以被表示为：

$>\hat{w}=A^{-1}v,$>公式9

更新的f_s，可由下式计算获得：

$>{\hat{f}}_s={\hat{w}}^{T}G\left(x\right),$>公式10

固定并使用更新的w和u，以为变量的F(u,c,w)最小值解为：

$>{\hat{c}}_i=\frac{\underset{\Omega }{\int }\left(f\right(x)-f_s(x\left)\right)u_i\left(x\right)dx}{\underset{\Omega }{\int }{u}_i^2\left(x\right)dx},i=1,...,N$>公式11

固定并使用更新的w和c，以u＝(u₁,…u_N)^T为变量的F(u,c,w)最小值解时，满足如下条件：

$>{\hat{u}}_i=\left(\begin{array}{c}1,i=i_{min}\left(x\right)\\ 0,i\ne i_{min}\left(x\right)\end{array}\right)$>公式12

其中，i_min(x)＝argmin{f(x)-c_i-w^TG(x)}。

进行不断的迭代运算更新u、c和w，求出w的稳定解。大量实验表明10次迭代后w的解基本稳定。最后所得到的校正后图像为$>f_p=f-{\hat{f}}_s=f-w^{T}G\left(x\right)$>

实验和结果

定量分析指标定义包括：

定义杯状伪影τ_cup＝100(u_M,edge-u_M,center)/u_M,edge，其中u_M,center和u_M,edge是模体中心和边缘的CT值(HU)。.

均方根对比度表示为$>RMSC=\sqrt{\frac{1}{MN}\underset{i=0}{\overset{N-1}{\Sigma }}\underset{j=0}{\overset{M-1}{\Sigma }}{(I_{ij}-\bar{I})}^2},$>其中I_ij是二维图像中(i,j)位置的像素值。是图像中所有像素平均值。

人头骨的锥束CT切片杯状伪影校正实验包括：

实验图像大小为211×211。如图2和图3所示，校正后图像的CT值均匀性明显提高。分析结果如表1所示。本方法中一个切片的校正过程耗时1.89秒(CPU:i5-2450,RAM:6GB,GPU:NVIDAGeForce610M)。

图3展示了水平剖面校正前后的图像。从中可以看出校正后观察图像的杯状伪影大幅度减少。

CTP486模体锥束CT切片杯状伪影校正实验包括：

使用CTP486模体进行实验。校正图像如图5所示。CTP486模块由CT数为含2％(0-20H)水的均匀材料浇铸而成。图像大小为229×229。一个切片的校正过程耗时1.93秒(CPU:i5-2450,RAM:6GB,GPU:NVIDAGeForce610M)。

水平剖面校正前后的图像如图6所示。从图6中可以看出校正前同种物质切片图CT值不均匀，即(表现为杯状伪影)，而使用本发明的方法后杯状伪影消除到人眼不易觉察的程度。

表1：头骨模体的量化分析。杯状伪影校正前的重建图像(RI_BC)，杯状伪影校正后的重建图像(RI_AC)

小鼠骨锥束CT切片杯状伪影校正实验，具体包括：

小鼠骨散射数据由HiscanM1000(Micro-CT)获得。重建图像的获取包括80kVp,200uA,30ms的360次投影。图8(a)和图8(b)展示了散射校正前后的小鼠骨重建图像的样本轴向视图。图像大小为339×339。一个切片的校正过程耗时2.7秒(CPU:i5-2450,RAM:6GB,GPU:NVIDAGeForce610M)。如图7所示，本方法将杯状伪影由23.8％降低到9.8％，分析结果如表2所示。

表2：小鼠骨的量化分析。杯状伪影校正前的重建图像(RI_BC)，杯状伪影校正后的重建图像(RI_AC)。