CN102737392A - 一种低剂量x线ct图像的非局部正则化先验重建方法 - Google Patents

Abstract

本发明公开了一种低剂量X线CT图像的非局部正则化先验重建方法，包括(1)利用X线CT成像设备采集病人先前扫描的标准剂量图像；(2)利用X线CT成像设备在Low-mAs扫描协议下采集病人的CT投影数据，并同时获取相应的校正参数及系统矩阵；(3)根据步骤(2)中获得的投影数据满足的统计分布，构建用于图像重建的数学模型；(4)利用步骤(1)中的标准剂量图像构建先前标准剂量图像导引的非局部正则化先验，并采用最大后验估计方法进行模型转化，结合步骤(3)得到的数学模型构建用于重建图像的目标函数；(5)对步骤(4)中构建的用于CT图像重建目标函数采用迭代算法进行求解，完成图像重建。本发明可以实现Low-mAs扫描协议下低剂量CT图像的重建。

Description

一种低剂量X线CT图像的非局部正则化先验重建方法

技术领域

本发明涉及一种医学影像的图像重建方法，具体来说是指一种低剂量X线CT图像的非局部正则化先验重建方法。

背景技术

X线CT图像的质量与辐射剂量密切相关，辐射剂量越高，图像质量就越好，然而，过高剂量的X射线照射又会诱发癌症、白血病或其他遗传性疾病。以最少的剂量获得最佳的CT影像诊断效果是CT成像领域研究的重要目标之一。

在其他扫描参数不变的情况下，直接降低X射线低毫安秒（Low-mAs）就可以实现低剂量CT成像，然而此时由于探测数据中将引入大量噪声，直接导致重建图像的严重退化，继而影响临床诊断。大量实验证明，对于低剂量CT图像重建，基于统计建模的重建方法优于传统滤波反投影方法。重要的是，基于统计建模的重建方法可以引入目标图像的先验信息，以保证重建图像的质量。但是，传统的统计建模方法中仅引入图像自身的局部邻域的约束作为先验信息，未能考虑同一病人先前扫描图像所能提供的先验信息，同时传统先验依赖的图像局部邻域内的像素灰度信息，无法有效的区分边缘信息和噪声，导致重建图像出现过平滑或阶梯状效应。

发明内容

本发明的目的在于提供一种低剂量X线CT图像的非局部正则化先验重建方法，新构建的非局部均值先验为各类基于统计建模的CT图像重建模型提供正则化项，可以实现Low-mAs扫描协议下低剂量CT图像的优质重建。

本发明的上述目的可以通过以下的技术措施来实现：一种低剂量X线CT图像的非局部正则化先验重建方法，包括以下步骤，其中本申请文件中相同的符号代表相同的含义：

(1)获取成像对象在标准剂量射线下的CT医学图像；

(2)获取同一成像对象在低剂量射线下的CT医学图像的投影数据，并同时获取相应的校正参数及系统矩阵，其中，低剂量射线的射线剂量为标准剂量射线的射线剂量的1/7至1/20，线积分投影数据是指对数变换后的投影数据p，获取的校正参数是指X线入射光子数N₀；

(3)根据步骤(2)中的线积分投影数据满足的统计特性，构建用于CT医学图像重建的数学模型；

(4)利用步骤(1)中的标准剂量射线下的医学图像构建先前标准剂量图像导引的非局部正则化先验，并采用最大后验估计方法进行模型转化，结合步骤(3)得到的数学模型构建用于重建图像的带先验约束项的目标函数；其中，所述非局部正则化先验构建的具体过程为：

(41)在标准剂量射线下的CT医学图像内选择一个包含图像几何信息的大方形邻域N_i，在大方形邻域内选择一个小方形邻域，该小方形邻域的中心位于像素点m处；在待重建图像中选择一个与上述小方形邻域大小相同的小方形邻域，该小方形邻域的中心位于像素点i，其中，像素点i和平扫图像中大方形邻域中心点几何位置相对应；

(42)通过相似性测度计算两个小方形邻域的距离，相似性测度的方法采用两像素点m、i邻域内所有像素点灰度值的加权欧几里德距离的反比例函数；利用两像素点m、i间相似性，正则化先验项的权值w_m，i定义为：

$w_{m,i}=\exp (-{\left|\right|\mathrm{\μ}\left(n_i\right)-{\mathrm{\μ}}^{\mathrm{nd}}\left(n_m\right)\left|\right|}_2^2/h^2)/{\mathrm{\Σ}}_{m\∈N_i}\exp (-{\left|\right|\mathrm{\μ}\left(n_i\right)-{\mathrm{\μ}}^{\mathrm{nd}}\left|\right|}_2^2/h^2)$

μ(n_i)＝{μ_c:c∈n_i}, ${\mathrm{\μ}}^{\mathrm{nd}}\left(n_m\right)=\{{\mathrm{\μ}}_c^{\mathrm{nd}}:c\∈n_m\}$

其中，μ(n_i)代表待重建图像μ中像素点i处相似邻域的灰度数组，μ^nd(n_m)代表标准剂量射线下的医学图像μ^nd中像素点m处相似邻域的灰度数组，n_i和n_m表示以像素点i和像素点m为中心的邻域，c为像素点表示的索引，h为计算像素点间权值量的指数函数衰减的控制系数；

(43)由上述(42)得到的权值矩，构建如下形式的非局部先验

N_i表示图像域内像素点i处的方形邻域，w_m，i表示正则化先验项的权值，v(·)表示正则化先验项的势函数；

(5)对步骤(4)中构建的用于CT图像重建的目标函数采用迭代算法进行求解，完成病人低剂量射线下的CT医学图像的重建。

本发明中，所述步骤(1)中标准剂量射线下的医学图像为临床标准mAs（毫安秒）扫描协议下通过滤波反投影方法重建的医学图像。

本发明中，所述步骤（3）中的统计分布特性为：

(a)对数变换后的线积分投影数据满足近似高斯分布，对应的数学重建模型中数据项为：(y-Hμ)′∑^-1(y-Hμ)；其中，μ表示待重建图像，H为系统矩阵；y表示本发明步骤(2)中获得的对数变换后的投影数据；

∑表示对角线元素为

的对角矩阵，

表示相应于i探测单元的方差；

(b)对数变换前投影数据满足近似泊松分布，对应的数学重建模型中数据项为：

其中，y表示所述步骤(2)中获得的对数变换前的投影数据，b_i表示成像过程X线球管中的曝光光子数，r_i表示相应的背景噪声。

本发明中，所述步骤(5)中的迭代算法可为梯度下降法或共轭梯度下降法或超松弛迭代算法。

所述梯度下降法的具体实现过程为：第一步，首先基于上一步迭代获得的重建图像及标准剂量射线下的医学图像获得待重建图像的逐像素点正则化先验项的权值量w_m，i，以作为下一步迭代之用；第二步，在第一步获取的权值量基础上，利用梯度下降算法进行迭代重建；第三步，交替进行第一、二步迭代直至收敛，获得最终重建的低剂量射线下的CT医学图像。

与现有技术相比，本发明公开的重建方法突破常规单纯利用待重建图像本身信息，引入由先前标准剂量图像导引的非局部先验，提高了X线CT低剂量成像条件下重建图像的质量，且本发明公开算法在图像噪声抑制、边缘保持等方面均有上佳表现

附图说明

下面结合附图和具体实施方式对本发明做进一步详细说明。

图1是本发明重建方法的流程图；

图2a为修订的Shepp-Logan体模图像；

图2b为仿真低剂量数据滤波反投影方法重建的图像；

图2c为仿真低剂量数据采用基于传统非局部先验迭代方法重建的图像；

图2d为仿真低剂量数据采用本发明重建方法重建的图像；

图3a为重建图像部分水平中线剖面图；

图3b为重建图像部分垂直中线剖面图；

图4a为原始的Shepp-Logan体模图像；

图4b为图4a仿射变换后的图像；

图4c为低剂量数据本发明公开的重建方法重建的图像，其中图4a图像作为标准剂量图像；

图4d为低剂量数据本发明公开的重建方法重建的图像，其中图4b图像作为标准剂量图像；

图5a为真实物理体模标准剂量数据采用滤波反投影方法重建的图像；

图5b为真实物理体模低剂量数据采用滤波反投影方法重建的图像；

图5c为真实物理体模低剂量数据采用基于传统非局部先验先验的迭代方法重建的图像；

图5d为真实物理体模低剂量数据采用本发明公开的重建方法重建的图像。

具体实施方式

本发明一种低剂量X线CT图像的非局部正则化先验重建方法，其具体实施步骤如图1所示，包括以下步骤：

1、数据采集：利用X线CT成像设备采集病人先前扫描的标准剂量射线下的CT医学图像，该标准剂量射线下的医学图像为临床标准毫安秒扫描协议下通过滤波反投影方法重建的医学图像；利用X线CT成像设备在低毫安秒扫描协议下采集同一病人低剂量射线的CT医学图像的线积分投影数据，并同时获取相应的校正参数及系统矩阵，其中，低剂量射线的射线剂量为标准剂量射线的射线剂量的1/7至1/20，可以由用户灵活设定，其中线积分投影数据是指对数变换后的投影数据y，获取的校正参数是指X线入射光子数N₀；

2、数据校正：利用获取的系统参数进行探测数据较正；

3、构建图像重建模型：对步骤2得到的低剂量线积分投影数据满足高斯分布的统计特性进行数学建模，并采用最大后验方法，完成标准剂量图像导引的先验项的设计，构建出用于CT增强扫描图像重建的带约束的目标函数Φ(μ)， $\mathrm{\Φ}\left(\mathrm{\μ}\right)={(y-\mathrm{H\μ})}^{\′}{\mathrm{\Σ}}^{-1}(y-\mathrm{H\μ})+\mathrm{\β}\underset{i}{\mathrm{\Σ}}{\mathrm{\Σ}}_{m\∈N_i}{w}_{m,i}v({\mathrm{\μ}}_i-{\mathrm{\μ}}_m^{\mathrm{nd}}),$ 其中μ表示待重建图像，H为系统矩阵；y表示线积分投影数据；∑表示对角线元素为

的对角矩阵，

表示相应于i探测单元的方差；R(μ)表示正则化先验项，其形式为

N_i表示图像域内像素点i处的方形邻域，w_m,i表示正则化先验项的权值，v(·)表示正则化先验项的势函数，β为正则化参数，可采用模拟退火算法确定，其中，上述的非局部正则化先验构建的具体过程为：

(31)在标准剂量射线下的CT医学图像内选择一个包含图像几何信息的大方形邻域N_i，在大方形邻域内选择一个小方形邻域，该小方形邻域的中心位于像素点m处；在待重建图像中选择一个与上述小方形邻域大小相同的小方形邻域，该小方形邻域的中心位于像素点i，其中，像素点i和平扫图像中大方形邻域中心点几何位置相对应；

(32)通过相似性测度计算两个小方形邻域的距离，相似性测度的方法采用两像素点m、i邻域内所有像素点灰度值的加权欧几里德距离的反比例函数；利用两像素点m、i间相似性，正则化先验项的权值w_m，i定义为：

$w_{m,i}=\exp (-{\left|\right|\mathrm{\μ}\left(n_i\right)-{\mathrm{\μ}}^{\mathrm{nd}}\left(n_m\right)\left|\right|}_2^2/h^2)/{\mathrm{\Σ}}_{m\∈N_i}\exp (-{\left|\right|\mathrm{\μ}\left(n_i\right)-{\mathrm{\μ}}^{\mathrm{nd}}\left|\right|}_2^2/h^2)$

μ(n_i)＝{μ_c:c∈n_i}, ${\mathrm{\μ}}^{\mathrm{nd}}\left(n_m\right)=\{{\mathrm{\μ}}_c^{\mathrm{nd}}:c\∈n_m\}$

其中，μ(n_i)代表待重建图像μ中像素点i处相似邻域的灰度数组，μ^nd(n_m)代表标准剂量射线下的医学图像μ^nd中像素点m处相似邻域的灰度数组，n_i和n_m表示以像素点i和像素点m为中心的邻域，c为像素点表示的索引，h为计算像素点间权值量的指数函数衰减的控制系数；

(33)由上述(32)得到的权值矩，构建如下形式的非局部先验

N_i表示图像域内像素点i处的方形邻域，w_m，i表示正则化先验项的权值，v(·)表示正则化先验项的势函数；

4、完成重建：在步骤3构建的相关模型基础上，采用梯度下降法的迭代算法进行图像重建。具体实现过程为：第一步，首先基于上一步迭代获得的重建图像及标准剂量图像获得待重建图像的逐像素点正则化先验项的权值量w_m，i，以作为下一步迭代之用；第二步，在第一步获取的权值量基础上，利用梯度下降算法进行迭代重建；第三步，交替进行第一、二步直至收敛，获得最终重建的病人低剂量射线下的CT医学图像。

本发明中步骤4中的迭代算法也可以采用共轭梯度下降法或超松弛迭代算法。

本发明技术与现有技术的实验比较及结果如下：

利用修定的Shepp-Logan体模模拟生成低剂量CT投影数据进行本发明算法实验评估。实验中，射源到旋转中心和探测器的距离分别为：570.00mm和1040.00mm，探测元的个数为672，大小为1.407mm；旋转一周的探测角向采样个数为1160。Shepp-Logan体模图像大小为512×512。图2中a为理想Shepp-Logan体模图像，该理想条件下的图像结构清晰不含任何伪影和噪声；b为低剂量数据采用滤波反投影算法直接重建图像，可以看到由于剂量的降低导致重建图像出现严重的统计噪声；c为低剂量数据采用基于传统非局部先验最大后验重建图像，易见该方法重建图像较滤波反投影重建图像有较大改善，然而均质区域仍包含明显的伪影；d为本专利方法重建图像，由图2d重建图像可看出，利用本发明方法重建获得的结果在抑制噪声和伪影方面明显优于其它算法重建结果。

图3中绘出对应于图2中重建图像的水平中线及垂直中线剖面图，鉴于整个剖面图中含有512个像素点，全部显示则难以区分各重建方法，故仅显示时仅截取其中一段，其区间为[200，300]。由图3可以看出，无论背景区域还是在目标区域，本发明公开方法的重建值更接近于理想值。

图4给出了配准和非配准条件下本发明算法重建图像。其中图4a为原始的Shepp-Logan体模图像；图4b为图4a经过仿射变换后的图像，作为非配准条件下标准剂量图像；图4c为低剂量数据本发明公开算法利用图4a图像作为导引图像的重建结果；图4d为低剂量数据本发明公开算法利用图4b图像作为导引图像的重建结果。从图中我们可以看出图4c和图4d视觉上无本质区别，在抑制噪声和伪影方面无明显区别。也就是说，本发明算法对标准剂量图像与待重建低剂量图像间的配准精度无严格要求。

图5给出真实物理体模数据的实验重建结果。标准剂量扫描参数为：管电流125mA，管电压：120kVp；低剂量扫描参数为：管电流17mA，管电压：120kVp。图5a为标准剂量数据滤波反投影算法重建图像；图5b为低剂量数据滤波反投影算法重建图像，可以看到相对于图5a图像中存在明显的噪声及条状伪影；图5c为低剂量数据采用基于传统非局部先验最大后验重建图像，易见部分的噪声和伪影被抑制，同时图像的分辨率也有所下降；图5d为低剂量数据采用本发明公开算法重建所得图像，可以看到噪声和伪影得到了很好的抑制，且分辨率得到较好保持。

本发明的实施方式不限于此，可以根据实际需要进行修改，以适应不同的实际需求，因此，在本发明上述基本技术思想前提下，按照本领域的普通技术知识和惯用手段对本发明内容所做出其它多种形式的修改、替换或变更，均落在本发明权利保护范围之内。

Claims (5)

1.一种低剂量X线CT图像的非局部正则化先验重建方法，包括以下步骤：

(1)获取成像对象在标准剂量射线下的CT医学图像；

(2)获取同一成像对象在低剂量射线下的CT医学图像的投影数据，并同时获取相应的校正参数及系统矩阵，其中，低剂量射线的射线剂量为标准剂量射线的射线剂量的1/7至1/20；

(3)根据步骤(2)中的投影数据满足的统计特性，构建用于CT医学图像重建的数学模型；

(4)利用步骤(1)中的标准剂量射线下的医学图像构建先前标准剂量图像导引的非局部正则化先验，并采用最大后验估计方法进行模型转化，结合步骤(3)得到的数学模型构建用于重建图像的带先验约束项的目标函数；其中，所述非局部正则化先验构建的具体过程为：

(41)在标准剂量射线下的CT医学图像内选择一个包含图像几何信息的大方形邻域N_i，在大方形邻域内选择一个小方形邻域，该小方形邻域的中心位于像素点m处；在待重建图像中选择一个与上述小方形邻域大小相同的小方形邻域，该小方形邻域的中心位于像素点i，其中，像素点i和平扫图像中大方形邻域中心点几何位置相对应；

(42)通过相似性测度计算两个小方形邻域的距离，相似性测度的方法采用两像素点m、i邻域内所有像素点灰度值的加权欧几里德距离的反比例函数；利用两像素点m、i间相似性，正则化先验项的权值w_m，i定义为：

$w_{m,i}=\exp (-{\left|\right|\mathrm{\μ}\left(n_i\right)-{\mathrm{\μ}}^{\mathrm{nd}}\left(n_m\right)\left|\right|}_2^2/h^2)/{\mathrm{\Σ}}_{m\∈N_i}\exp (-{\left|\right|\mathrm{\μ}\left(n_i\right)-{\mathrm{\μ}}^{\mathrm{nd}}\left|\right|}_2^2/h^2)$

μ(n_i)＝{μ_c:c∈n_i}, ${\mathrm{\μ}}^{\mathrm{nd}}\left(n_m\right)=\{{\mathrm{\μ}}_c^{\mathrm{nd}}:c\∈n_m\}$

其中，μ(n_i)代表待重建图像μ中像素点i处相似邻域的灰度数组，μ^nd(n_m)代表标准剂量射线下的医学图像μ^nd中像素点m处相似邻域的灰度数组，，n_i和n_m表示以像素点i和像素点m为中心的邻域，c为像素点表示的索引，h为计算像素点间权值量的指数函数衰减的控制系数；

(43)由上述(42)得到的权值矩，构建如下形式的非局部先验N_i表示图像域内像素点i处的方形邻域，w_m，i表示正则化先验项的权值，v(·)表示正则化先验项的势函数；

(5)对步骤(4)中构建的用于CT图像重建的目标函数采用迭代算法进行求解，完成低剂量射线下的CT医学图像的重建。

2.根据权利要求1所述的低剂量X线CT图像的非局部正则化先验重建方法，其特征在于：所述步骤(1)中标准剂量射线下的医学图像为临床标准毫安秒扫描协议下通过滤波反投影方法重建的医学图像。

3.根据权利要求1所述的低剂量X线CT图像的非局部正则化先验重建方法，其特征在于：所述步骤（3）中的统计分布特性为：

(a)对数变换后的线积分投影数据满足近似高斯分布对数变换前投影数据满足近似泊松分布，对应的数学重建模型中数据项为：(y-Hμ)′∑^-1(y-Hμ)；其中，μ表示待重建图像，H为系统矩阵；y表示本发明步骤(2)中获得的对数变换后的投影数据；

∑表示对角线元素为的对角矩阵，

表示相应于i探测单元的方差；

(b)对数变换前投影数据满足近似泊松分布，对应的数学重建模型中数据项为：

其中，y表示所述步骤(2)中获得的对数变换前的投影数据，b_i表示成像过程X线球管中的曝光光子数，r_i表示相应的背景噪声。

4.根据权利要求1所述的低剂量X线CT图像的非局部正则化先验重建方法，其特征在于：所述步骤(5)中的迭代算法为梯度下降法或共轭梯度下降法或超松弛迭代算法。

5.根据权利要求4所述的低剂量X线CT图像的非局部正则化先验重建方法，其特征在于：所述梯度下降法的具体实现过程为：第一步，首先基于上一步迭代获得的重建图像及标准剂量射线下的医学图像获得待重建图像的逐像素点正则化先验项的权值量w_m，i，以作为下一步迭代之用；第二步，在第一步获取的权值量基础上，利用梯度下降算法进行迭代重建；第三步，交替进行第一、二步迭代直至收敛，获得最终重建的低剂量射线下的CT医学图像。

Images