基于全广义变分正则化的低剂量x射线ct图像重建方法
【专利摘要】一种基于全广义变分正则化的低剂量X射线CT图像重建方法，包括如下步骤：(1)获取CT设备的系统参数和低剂量扫描协议下对数变换后的投影数据；(2)对投影数据进行Anscombe变换，将服从Poission分布的投影数据转化为近似服从方差为1的Gaussian分布数据u；(3)建立基于全广义变分最小化的理想数据恢复模型，使用Chambolle?Pock算法求解得到恢复后的投影数据f；(4)对步骤(3)得到的恢复后的投影数据f进行Anscombe逆变换，再通过滤波反投影算法得到CT重建图像。该方法可以在投影数据不满足分段常数假设的前提下去除图像的噪声和条形伪影，同时较好地保持图像的空间分辨率。
【专利说明】
基于全广义变分正则化的低剂量X射线CT图像重建方法
技术领域
[0001] 本发明涉及医学影像成像技术领域，特别是涉及一种基于全广义变分正则化的低 剂量X射线CT图像重建方法。
【背景技术】
[0002] X射线计算机断层成像(Computed Tomography，CT)因其在时间、空间分辨率上的 卓越表现，已经广泛应用于临床诊断和治疗。CT图像的质量与X射线辐射剂量密切相关，剂 量越高图像质量越好，然而，过量的X射线照射又会诱发恶性肿瘤、白血病以及其他遗传性 疾病。
[0003] 为了降低X射线的辐射剂量，许多优化的低剂量CT扫描方案以及抑制噪声和伪影 的重建算法相继被提出。目前，直接降低管电流、管电压或者减少扫描角向采样数目是实现 低剂量CT扫描最简单且最有效的方法。对于这一问题有两种传统的解决方法:一是直接对 低剂量CT图像进行滤波，以减少图像的噪声和伪影，属于图像后处理技术;二是根据投影数 据满足的统计学规律，完成基于统计的CT图像迭代重建。
[0004] 第一类方法为后处理技术，直接减少图像的噪声和伪影。图像后处理技术因其简 单且易于操作。由于低剂量CT图像中噪声和伪影分布的复杂性，使得高精度的滤波器设计 困难极大。
[0005] 第二种方法是通过CT系统建模，构建图像重建模型，通过优化目标函数实现图像 重建。相对于经典的FBP算法，迭代重建算法通过系统建模(系统光学模型和系统统计模型） 对CT成像几何、X射线的能谱特性、射束硬化效应、散射和噪声特性进行准确描述，而且易于 加入先验信息约束，因此特别适合低剂量CT图像优质重建。统计迭代重建在抑制图像噪声 和伪影以及提高空间分辨率等方面都有上佳表现，还可以结合其他硬件相关的低剂量成像 技术，进一步降低辐射剂量，提高重建图像的质量。然而，由于统计迭代重建需要反复进行 投影与反投影运算，且CT图像数据量庞大，导致CT图像重建速度特别慢，难以满足临床中实 时交互的需求。
[0006] 不同于上述两种方法的另一种策略是，基于低剂量CT投影数据恢复的图像重建方 法对CT成像系统进行精确建模(包括X线源、探测器、电子噪声、成像物体等），可以实现优质 的低剂量CT图像重建。其目标函数是根据测量数据的噪声特性来构建的，通常包含两项，即 数据保真项和正则化项，前者用于描述投影数据的统计特性，后者用于修正解的特性。然 而，许多临床CT图像并不是完全满足分段常量的假设，在剂量特别低时或者投影角度特别 少时，TV重建的图像会产生阶梯效应和块状伪影。因此TV正则化在基于投影数据恢复的低 剂量CT重建中会失去应有的效力。
[0007] 因此，针对现有技术不足，有必要提供一种基于全广义变分正则化的低剂量X射线 CT图像重建方法，其可以在投影数据不满足分段常数假设的前提下去除图像的噪声和条形 伪影，同时较好地保持图像的空间分辨率。
【发明内容】
[0008] 本发明的目的在于避免现有技术的不足之处而提供一种基于全广义变分正则化 的低剂量X射线CT图像重建方法，该基于全广义变分正则化的低剂量X射线CT图像重建方 法，可以在投影数据不满足分段常数假设的前提下去除图像的噪声和条形伪影，同时较好 地保持图像的空间分辨率。
[0009] 本发明的上述目的通过如下技术手段实现。
[0010] 提供一种基于全广义变分正则化的低剂量X射线CT图像重建方法，包括如下步骤：
[0011] (1)获取CT设备的系统参数和低剂量扫描协议下对数变换后的投影数据q;
[0012] (2)对步骤（1)获取的投影数据q进行Anscombe变换，将服从复合Poission分布的 投影数据q转化为近似服从方差为1的Gaussian分布数据u;
[0013] (3)对步骤(2)获取的数据u建立基于全广义变分最小化的理想数据恢复模型，使 用Chambolle-Pock算法求解得到恢复后的投影数据f;
[0014] (4)对步骤(3)得到的恢复后的投影数据f进行Anscombe逆变换，再通过滤波反投 影算法得到CT重建图像。
[0015]上述步骤(1)中获取的CT设备的系统参数包括X射线入射光子强度1〇。
[0016] 上述步骤(2)中对步骤(1)获取的投影数据q进行Anscombe变换，计算式如下：
[0017]
[0018] U= (U1，U2，…，
UN)1表示投影数据经过Anscombe变化后近似服从方差为1的 Gaussian分布数据，其中，T表示转置运算，ui、U2、"_、un是u的元素，N为元素个数;令f = (fi, 5，~，办)7为待估计投影数据的理想值44、~、&是抒勺分量。
[0019] 上述步骤(3)中对步骤(2)获取的数据u建立全广义变分最小化模型，具体通过如 下方法进行：
[0020] (3-1)投影数据恢复的一般模型通过式(2)表示：
[0021]
[0022]其中，第一项为数据保真项，第二项为正则化项，λ>〇是正则化参数；
[0023] (3-2)通过全广义变分正则化代替式(2)中的正则化项R(f)，得到基于全广义变分 最小化的理想数据恢复模型为：
[0024]
[0025]
[0026]
[0027]
[0041] 其中，τ是步长，τ>〇。
[0042] 本发明的基于全广义变分正则化的低剂量X射线CT图像重建方法，包括如下步骤： (1)获取CT设备的系统参数和低剂量扫描协议下对数变换后的投影数据q; (2)对步骤(1)获 取的投影数据q进行Anscombe变换，将服从复合Poission分布的投影数据q转化为近似服从 方差为1的Gauss ian分布数据u;(3)对步骤(2)获取的数据u建立基于全广义变分最小化的 理想数据恢复模型，使用Chambolle-Pock算法求解得到恢复后的投影数据f;(4)对步骤(3) 得到的恢复后的投影数据f进行Anscombe逆变换，再通过滤波反投影算法得到CT重建图像。 该方法可以在投影数据不满足分段常数假设的前提下去除图像的噪声和条形伪影，同时较 好地保持图像的空间分辨率。
【附图说明】
[0043] 利用附图对本发明作进一步的说明，但附图中的内容不构成对本发明的任何限 制。
[0044] 图1是本发明一种基于全广义变分正则化的低剂量X射线CT图像重建方法的流程 图。
[0045] 图2是本发明实施例2的数值Clock体模及采用不同方法重建的结果。其中，（a)是 真实的Clock体模图像；（b)表示采用斜波滤波反投影方法(FBP算法，Ramp窗）重建的图像； (c)表示采用汉宁滤波反投影方法(FBP算法，Hanning窗）重建的图像；（d)表示采用非单调 最小化方法(NTVM方法)重建的图像;（e)表示采用本发明的方法重建的图像。
[0046] 图3是对应图2数值Clock体模重建结果的局部放大图；每行分别对应一个区域，每 行的第一列分别是真实的Clock体模图像，第二列是采用斜波滤波反投影方法(FBP算法， Ramp窗）重建的图像；第三列是采用汉宁滤波反投影方法(FBP算法，Hanning窗）重建的图 像;第四列是采用非单调最小化方法(NTVM方法)重建的图像;第五列是采用本发明的方法 重建的图像。
[0047] 图4是对应图2中不同算法重建图像的线剖面图。
[0048]图5是本发明实施例2的数值Shepp-Logan体模及采用不同方法重建的结果。其中， (a)是真实的Shepp-Logan体模图像；（b)表示采用斜波滤波反投影方法(FBP算法，Ramp窗） 重建的图像；（c)表示采用汉宁滤波反投影方法(FBP算法，Hanning窗)重建的图像；（d)表示 采用非单调最小化方法(NTVM方法)重建的图像;（e)表示采用本发明的方法重建的图像。 [00 49]图6是对应图5数值Shepp-Logan体模重建结果的局部放大图；每行分别对应一个 区域，每行的第一列分别是真实的Shepp-Logan体模图像，第二列是采用斜波滤波反投影方 法(FBP算法，Ramp窗）重建的图像;第三列是采用汉宁滤波反投影方法(FBP算法，Hann i n g 窗)重建的图像;第四列是采用非单调最小化方法(NTVM方法)重建的图像;第五列是采用本 发明的方法重建的图像。
[0050]图7是对应图5中不同算法重建图像的线剖面图。
【具体实施方式】
[0051 ]结合以下实施例对本发明作进一步描述。
[0052] 实施例1。
[0053]提供一种基于全广义变分正则化的低剂量X射线CT图像重建方法，通过以下步骤 进行的。
[0054] (1)获取CT设备的系统参数和低剂量扫描协议下对数变换后的投影数据q。获取的 CT设备的系统参数包括X射线入射光子强度Io等。
[0055] (2)对步骤（1)获取的投影数据q进行Anscombe变换，将服从复合Poission分布的 投影数据q转化为近似服从方差为1的Gaussian分布数据u。
[0056] 对步骤(1)获取的投影数据q进行Anscombe变换，计算式如下：
[0071] 其中，F = Rnxn，W = R2nxn以及微分算子ε，▽由限差分算子近似得到，Rnxn表示N X N维 的实空间。 L 0082 J 其中，τ是步长，τ>〇。
[0083] (4)对步骤(3)得到的恢复后的投影数据f进行Anscombe逆变换，再通过滤波反投 影算法得到CT重建图像。
[0084] 本发明的一种基于全广义变分正则化的低剂量X射线CT图像重建方法，构造全广 义变分正则化项，可以在投影数据不满足分段常数假设的前提下去除噪声。最后通过 Anscombe逆变换和经典的滤波反投影算法对恢复后的投影数据进行解析重建。数值体模实 验结果表明，本发明可以有效地抑制低剂量CT图像中的噪声和条形伪影，同时可以很好地 保持图像的结构信息和空间分辨率。
[0085] 实施例2。
[0086] 为了对本发明的效果做进一步验证，采用图2所示的Clock数字体模图像和图5所 示的Shepp-Logan数字体模图像作为本发明的计算机仿真实验对象。
[0087] (1)体模图像大小设定位512X512,模拟CT机的X射线源到旋转中心和探测器的距 离分别为570mm和1040mm，旋转角在[0，2JT]间采样值为1160，每个采样角对应672个探测器 单元，探测器单元的大小为1.407mm。通过CT系统仿真生成大小为1160X672的投影数据q， 对于Clock体模，入射光子数为5.0 X 104;相应的Shepp-Logan体模的入射光子数为2.5 X 1〇5。需要说明的是，在实际的CT数据采集中，投影数据和系统参数即入射光子强度Io可以直 接获取。
[0088] (2)对步骤（1)获取的投影数据q进行Anscombe变换，将服从复合Poission分布的 投影数据q转化为近似服从方差为1的Gaussian分布数据u。
[0089] 对步骤(1)获取的投影数据q进行Anscombe变换，计算式如下：
[0090] ? ·-式（1
[0091] 11=(111，112，'"，_)7表示投影数据经过411% 〇1111^变化后近似服从方差为1的 Gaussian分布数据，其中，T表示转置运算，Ui、U2、···、un是u的分量，N为分量的个数;令f = (负力，"_，况7为待估计投影数据的理想值土上、"_、&是炸勺分量。
[0092] (3)对步骤(2)获取的数据u建立基于全广义变分最小化的理想数据恢复模型，使 用Chambolle-Pock算法求解得到恢复后的投影数据f。
[0093] 对步骤(2)获取的数据u建立全广义变分最小化模型，具体通过如下方法进行： [0094] (3-1)投影数据恢复的一般模型通过式(2)表示：
[0095] …·"式（2);
[0096]其中，第一项为数据保真项，第二项为正则化项，λ>〇是正则化参数。
[0097] (3-2)通过全广义变分正则化代替式(2)中的正则化项R(f)，得到基于全广义变分 最小化的理想数据恢复模型为：
[0098]
[0099]
[0100]
[0115] 其中，τ是步长，τ>〇。
[0116] (4)对步骤(3)得到的恢复后的投影数据f进行Anscombe逆变换，再通过滤波反投 影算法得到CT重建图像。
[0117] 为了验证全广义变分最小化方法在低剂量CT图像重建中的有效性，通过Clock和 Shepp-Logan数值体模进行定量和定性分析。并且与斜波滤波反投影算法、汉宁滤波反投影 算法和非单调全变分最小化算法进行比较。CT成像几何采用弧形探测器的扇形束，其中射 线源到旋转中心和探测器的距离分别为570mm和1040mm，旋转角在[0,2π]间的采样值为 1160，探测器个数为672。体模大小设定为512 X 512，对于Clock体模，入射光子数为5.0 X IO4;相应的Shepp-Logan体模的入射光子数为2 · 5 X IO5。
[0118] 图2是本发明实施例2的数值Clock体模及采用不同方法重建的结果。其中，（a)是 真实的Clock体模图像；（b)表示采用斜波滤波反投影方法(FBP算法，Ramp窗）重建的图像； (c)表示采用汉宁滤波反投影方法(FBP算法，Hanning窗）重建的图像；（d)表示采用非单调 最小化方法(NTVM方法)重建的图像；（e)表示采用本发明的方法重建的图像。通过图2可以 看到，斜波滤波反投影和汉宁滤波反投影算法重建的图像含有大量的噪声和条形伪影，图 像质量严重退化。非单调全变分最小化算法重建的图像虽然在一定程度上减少了 CT图像的 条形伪影和噪声。从视觉效果评价来看，本发明在去除图像的噪声和条形伪影方面具有更 佳的表现。为了进一步展现本方法重建效果，由图3中对应图2的局部放大图可以清晰的看 出，本发明对噪声和伪影的抑制效果。
[0119] 通过计算重建图像的信噪比和均方误差可以定量地分析不同方法的差异。信噪比 计算公式为：
[0120]
[0121]
[0122]
[0123] 其中，Fre3Ji)表示像素点i通过算法得到的重建图像的灰度值，表示重建图像 中所有像素点灰度的平均值;Ft_(i)表示真实图像像素点i处的灰度值。由表1的性噪 比和均方误差分析可以看出，相对于滤波反投影和非单调最小化方法，本发明提出的重建 图像的方法具有最高的信噪比和最小的均方误差。
[0124] 表1图2中重建图像的信噪比和均方误差
[0126]为了进一步验证本文提出方法的优越性，分别给出了四种方法重建结果的线剖面 图，如图4所示，通过图4可以看出本文提出的重建图像的方法与其他三种方法相比，能够较 好的接近真实图像的线剖面。
[0127]图5是本发明实施例2的数值Shepp-Logan体模及采用不同方法重建的结果。其中， (a)是真实的Shepp-Logan体模图像；（b)表示采用斜波滤波反投影方法(FBP算法，Ramp窗） 重建的图像；（c)表示采用汉宁滤波反投影方法(FBP算法，Hanning窗)重建的图像；（d)表示 采用非单调最小化方法(NTVM方法)重建的图像；（e)表示采用本发明的方法重建的图像。通 过图5可以看到，相比滤波反投影和非单调全变分最小化方法，本发明提出的全广义变分最 小化方法在抑制噪声与条形伪影和保持边缘方面都有良好表现，可以更好的保持图像一致 性。为了更清晰地比较四种方法的噪声和伪影的抑制效果，图6给出了低剂量CT重建结果的 一个局部的放大图。从放大图可以明显的看出，本发明提出的全广义变分方法重建的图像 对噪声和伪影有更好的抑制作用。
[0128]由表2的信噪比和均方误差分析可以看出，与其他三种方法相比，本发明提出的全 广义变分方法在提高图像重建的信噪比，降低重建图像均方误差方面均有上佳表现。
[0129]表2图5中重建图像的信噪比和均方误差
[0131] 此外，为了更直观的验证重建图像的质量，画出重建图像的线剖面图，如图7所示， 可以看出本发明重建图像的方法与真实图像的线剖面有更高的吻合度，体现出了全广义变 分方法的优越性。
[0132] 本发明在克服低剂量CT重建的全变分(TV)正则化方法的不足基础上提出了一个 低剂量CT优质重建的全广义变分(TGV)方法性。本发明构造全广义变分正则化项，可以在投 影数据不满足分段常数假设的前提下去除噪声。最后通过Anscombe逆变换和经典的滤波反 投影算法对恢复后的投影数据进行解析重建。数值体模实验结果表明，本发明可以有效地 抑制低剂量CT图像中的噪声和条形伪影，同时可以很好地保持图像的结构信息和空间分辨 率。
[0133] 最后应当说明的是，以上实施例仅用以说明本发明的技术方案而非对本发明保护 范围的限制，尽管参照较佳实施例对本发明作了详细说明，本领域的普通技术人员应当理 解，可以对本发明的技术方案进行修改或者等同替换，而不脱离本发明技术方案的实质和 范围。
【主权项】
1. 一种基于全广义变分正则化的低剂量X射线CT图像重建方法，其特征在于，包括如下 步骤： (1) 获取CT设备的系统参数和低剂量扫描协议下对数变换后的投影数据q; (2) 对步骤（1)获取的投影数据q进行Anscombe变换，将服从复合Poission分布的投影 数据q转化为近似服从方差为1的Gaussian分布数据U; (3) 对步骤（2)获取的数据U建立基于全广义变分最小化的理想数据恢复模型，使用 化ambolle-Pock算法求解得到恢复后的投影数据f; (4) 对步骤(3)得到的恢复后的投影数据f进行Anscombe逆变换，再通过滤波反投影算 法得到CT重建图像。2. 根据权利要求1所述的基于全广义变分正则化的低剂量X射线CT图像重建方法，其特 征在于，所述步骤(1)中获取的CT设备的系统参数包括X射线入射光子强度1〇。3. 根据权利要求2所述的基于全广义变分正则化的低剂量X射线CT图像重建方法，其特 征在于， 所述步骤(2)中对步骤(1)获取的投影数据q进行Anscombe变换，计算式如下：u=(ui，U2，.··，UN)T表示投影数据经过Anscombe变化后近似服从方差为1的Gaussian分 布数据，其中，T表示转置运算，山、化、…、UN是U的分量，N为分量的个数;令 f =(fl，f2，···，fN)T为待估计投影数据的理想值，fl、f2、，，，、fN是f的分量。4. 根据权利要求3所述的基于全广义变分正则化的低剂量X射线CT图像重建方法，其特 征在于， 所述步骤(3)中对步骤(2)获取的数据U建立全广义变分最小化模型，具体通过如下方 法进行： (3-1)投影数据恢复的一般模型通过式(2)表示：其中，第一项为数据保真项，第二项为正则化项，λ>0是正则化参数； (3-2)通过全广义变分正则化代替式(2)中的正则化项R(f)，得到基于全广义变分最小 化的理想数据恢复模型为：其中α〇，αι〉〇是两个正数，V/是f的梯度，ω是对偶变量衣紋 是ω的梯度。5.根据权利要求4所述的基于全广义变分正则化的低剂量X射线CT图像重建方法，其特 征在于， 式(3)的离散形式为：其中，F二护xw，w二R2胃W及微分算子ε，ν由限差分算子近似得到，R胃表示NXN维的实 空间； 根据对偶原理，将式(4)转化为式巧)的鞍点问题：其中，其中，τ是步长，τ〉〇。
【文档编号】G06T11/00GK105844678SQ201610427491
【公开日】2016年8月10日
【申请日】2016年6月15日
【发明人】牛善洲, 李楠, 吴恒, 马建华, 喻高航
【申请人】赣南师范学院