CN115880387A - 基于离散nlhtv正则项的稀疏角度采样ct迭代重建方法 - Google Patents

Abstract

本发明提出基于离散NLHTV正则项的稀疏角度采样CT迭代重建方法，属于图像处理领域的计算机断层射线成像(CT)技术。本发明同时考虑CT图像的非局部一阶和非局部二阶变分的稀疏性，以解决稀疏角度CT重建问题。该方法的优点是不仅可以利用非局部混合阶数导数信息保护图像结构特征，而且可以有效利用非局部自相似性恢复重建图像的细节。由于所提出的重建模型存在非光滑凸问题，因此引入交替方向乘法器(ADMM)优化方法，通过拆分为多个子问题来有效求解优化模型，构建重建技术方案。在稀疏采样条件下，该发明能有效提高重建图像的边缘与细节质量，从而提高CT重建图像的整体质量，实用价值显著。

Description

基于离散NLHTV正则项的稀疏角度采样CT迭代重建方法

技术领域

本发明属于图像处理领域的计算机断层射线成像(CT)技术，涉及一种基于离散NLHTV正则项的稀疏角度采样CT重建迭代重建方法。

背景技术

作为一种先进的成像技术，计算机断层扫描(CT)在不损坏物体情况下可以通过重建切片图像来可视化物体的内部结构，它已广泛应用于医学影像诊断、工业无损检测、公安检查和逆向工程等领域。在CT成像领域，采样足够的投影数据是传统重建算法获得高质量断层图像的前提条件。但是，X射线对人体有害是众所周知的。因此，降低CT的辐射剂量日益成为人们普遍关注的问题。通过将传统的密集角度投影采集更改为稀疏角度采样模式，可实现有效降低辐射剂量。但是，稀疏角度采样模式可以轻易导致重建图像产生条带伪影，使得CT图像重建质量明显退化，非常不利于诊断。如何保证稀疏采样投影数据情况下获取高质量CT重建图像具有重要意义和实用价值。

为了提高重建质量，有必要引入先验信息进行正则化。压缩传感理论表明，稀疏先验在稀疏角度采样CT重建中发挥重要作用。目前，变分正则化是一种广泛用于稀疏角度CT重建，因为它们成功促进稀疏性，其中总变分(TV)是应用最为广泛。与传统的ART和EM相比，TV方法的重建性能由于其边缘保留能力而得到显着提高。虽然TV方法已经取得良好的重建性能，但仍有一些局限性。TV重建图像可能会受到不希望的斑块伪影的影响，这对保护一些平稳过渡结构并不友好。为了解决TV的弱点，有必要在正则化项中引入重建图像的高阶导数信息，如基于总广义变异(TGV)的重建算法，Hessian矩阵正则化，混合阶变分导数正则化等。除了利用随着不同阶数变化的稀疏性，但现有的方法仍然没有考虑高阶变分的非局部自相似特性。

本发明提出一种由非局部一阶变分和二阶变分组成的新的非局部混合变分正则项(NLHTV)，能有效提高变分CT图像重建模型的性能。相比现有的变分正则化项，NLHTV通过建模重建图像的非局部“分段常数”和“分段线性”特征，能有效用图像的非局部自相似性来恢复图像细节。该方法能改善重建图像的边缘质量，从而提高CT重建图像的整体质量。此外，提出了一种ADMM算法的高效求最优解方案，CT图像重建算法如图1所示。

发明内容

为了解决现有技术存在的问题，本发明提供了一种基于离散NLHTV正则项的稀疏角度采样CT迭代重建方法。

本发明的技术方案：

基于离散NLHTV正则项的稀疏角度采样CT迭代重建方法，包括以下步骤：

步骤1CT成像模型建立；

CT成像的基本原理是通过测量X射线从不同方向穿过物体的衰减量，从而重构物体内部的横切面。假设X射线是单色的，X射线沿传播路径L穿过物体的衰减，由在数学上由比尔定律描述，公式为

其中，I₀是X射线源发射的初始强度。I_i是X射线穿过物体后的X射线强度，μ(l)是沿着透射路径的X射线衰减系数，如图2所示。

为了重构断面图像，被重建区域需要离散为N₁×N₂大小的矩形栅格。被重构断面图像表达为u＝[u₁,u₂,···,u_j,···,u_N]^T，u_j表示第j个像素值，N＝N₁×N₂是像素总数。经过离散化后，公式(1)中的线积分方程表达为迭代方程。针对第i个射线，获得如下表达式

其中，a_i,j表达为对投影数值p_i做出贡献的像素u_j的权系数，如图2所示。考虑X射线穿过物体的所有路径，可表达为以下线性方程组

Au＝p+e(3)

其中，A∈R^M×N是投影矩阵，a_i,j≥0是矩阵A的第i行，第j列的元素。矢量p＝[p₁,p₂,···,p_i,···,p_M]^T代表所有投影数据，M是被测量投影数据的总数。e∈R^M×1表达为测量误差，包括探测器噪声和X射线的光子噪声。

步骤2定义NLHTV正则项；

在全广义变差正则化定义基础上，引入加性稀疏正则化项，提出一种非局部混合变分正则化(NLHTV)项J_NLHTV(u)，如下所示：

其中，Ω为定义域，符号u(x)和u(y)分别代表空间位置x和y处的像素值，g(x)和g(y)分别表示x和y处的非局部一阶变分，g^m(x)代表g(x)的第m个分量，同理g^m(y)代表g(y)的第m个分量。w₁,w₂,w₃表示支撑权值，定义如下：

w₂(x,y)＝α·w₁(x,y)

w₃(x)＝β·∫w₁(x,y)dy

其中，C是归一化因子，h_c,h_l,h_g均为大于0的常数，权值w₁,w₂,w₃衰减快慢由h_c,h_l,h_g决定，α，β均为大于0的常数。

步骤3离散化NLHTV

在NLHTV正则项的离散定义基础上，基于离散NLHTV正则项的稀疏角度CT图像重建方法能使得损失函数最小化，表达为下式

/>

其中，参数μ＞0用来平衡NLHTV正则项和数据拟合项；i，j为像素索引；d_ij＝|i-j|；u_i,u_j代表索引为i，j的像素值，即u(x)在iΔl与u(y)在jΔl处的像素值(Δl为变量x或y的离散化的间隔)；g_i，g_j代表一阶索引为i，j处像素的非局部一阶变分，即g(x)在iΔl处的值与g(y)在jΔl处的值；w₁(i,j),w₂(i,j)表示支撑权值w₁(x,y)，w₂(x,y)的离散形式，即在w₁(x,y)，w₂(x,y)在(iΔl,jΔl)处的值；w₃(n)代表∑_mw₃(nΔl,mΔl)，m与n为像素索引，w₃(nΔl,mΔl)表示w₃(x,y)在(nΔl,mΔl)处的值。

所提重构模型可以同时利用非局部一阶、非局部二阶变分正则化项，可有效地利用稀疏度和非局部自相似冗余先验信息，提高噪声抑制能力，结构细节恢复能力，并保持重建图像边缘。

步骤4求解NLHTV重建模型

为了便于公式推导，以矩阵表示形式重写离散NLHTV正则项，表示为

其中，u是重建图像，F是非局部的一阶差分矩阵，

表示由整个图像的非局部一阶导数构成的向量。矩阵D代表与非局部差分相对应的距离矩阵，其元素由d_i,j构成。类似地，W₁,W₂,W₃分别代表与非局部差分相对应的权重矩阵，其元素分别由{w₁}，{w₂}，{w₃}构成，⊙表示两个矩阵的元素乘法。

步骤5求解子问题；

为求解问题(6)，设计了基于ADMM的优化求解方案，将整个优化问题分解为几个相对简单的子问题，从而迭代求解。通过引入辅助变量X,Y,Z，将上述问题表述为

优化问题(7)对应的拉格朗日函数为

其中，ρ＞0为二次项惩罚系数，λ₁、λ₂、λ₃为拉格朗日乘子。

ADMM方法通过迭代优化增广拉格朗日函数中的X,Y,Z变量求解原始问题，对应的迭代公式表达如下

其中，上标k表示迭代次数。

步骤6子问题求解

下面给出每个子问题求解的细节。

1)求解X,Y,Z子问题

关于X,Y,Z的问题相当于以下形式：

上述问题解析解表达可表示为

其中，sign()表示符号函数。换言之，给定x≥0,sign(x)＝1，否则sign(x)＝0。max(x)函数表示取x中的最大值，ρ的意义同公式(8)中。

2)求解子问题u

固定X^k+1,Y^k+1,Z^k+1，子问题u是最小平方问题，可令

式中关于u的梯度等于0，即/>

(符号/>

代表求梯度)，可以进一步写为下式

将公式(12)表达为矩阵方程形式，得

式(13)为求解线性方程组问题，涉及求解大规模矩阵求逆。由于直接获得逆矩阵非常困难，因此采用迭代求解方法。为提高求解子问题u的效率，以加快收敛速度快，提数值稳定性高，采用共轭梯度法求解方程(13)。

3)更新乘子

在迭代过程中，设置误差容忍度ε＞0，判断如果||u^k+1-u^k||≥ε，更新迭代次数k＝k+1；如果||u^k+1-u^k||＜ε，则输出最终重建结果u。

本发明所提出基于NLHTV正则化的稀疏角度CT重建方法，可以同时利用非局部一阶、非局部二阶变分正则化项，对图像特征具有很强的细节描述能力，因而擅长保留重建图像中的丰富细节，使得CT图像细节质量明显改善。同时，本发明的重建方法引入非局部一阶、非局部二阶变分的非局部自相似性的机制，可有效利用重建图像不同区域的医学组织自相似冗余先验信息，来压制重建图像中的噪声和因欠采样导致的伪影，保护重建图像边缘。该发明能有效解决现有基于变分正则化的CT图像重建方法不能利用非局部高阶变分信息来重建高质量的图像细节问题。在相同的稀疏角度采样模式下，该发明能相对于现有的基于变分正则化的CT图像重建方法，重建图像质量提升明显。因而，在保证CT图像的前提下，相对于现有的变分CT重建方法，所发明能进一步有效缩短病人受X射线照射的时间，降低X射线对病人的辐射危害，实用价值显著。

附图说明

图1为本发明的CT图像重建流程的示意图。

图2为X射线源发射X射线沿着路径穿过物体到达探测器，其中I₀表示X射线的最初密度，I_i表示X射线穿过物体后的密度。

图3为被重建区域需要离散为N₁×N₂大小的矩形栅格，被重构断面图像表达为u＝[u₁,u₂,···,u_j,···,u_N]^T，u_j表示第j个像素值，N＝N₁×N₂是像素总数。经过离散化后，公式(1)中的线积分方程表达为迭代方程。针对第i条射线，获得如下表达式

其中，a_i,j表达为对投影数值p_i做出贡献的像素u_j的权系数。

图4(a)是参考图像，图4(b)是所提方法重建结果，该方法去除掉重建结果中的阶梯效应和块状伪影，并有效保留边缘和平滑过渡部分。为了观察重建结果的细节差异，感兴趣区域(ROI)为放大并显示在每张重建图像的右上角。

具体实施方式

以下根据实施例和附图对本发明的技术方案进行进一步说明。

实施例：

基于离散NLHTV正则项的稀疏角度采样CT迭代重建方法，步骤如下：

(1)射线源到探测器距离是1400mm。射线源到旋转轴距离是1000mm。探测器像元个数为1024。重建图像尺寸为256×256，且每个像素大小为0.5mm×0.5mm；

(2)X射线管电压和电流分别设置为90kvp、0.1mA，探测器曝光时间设置为500ms，在[0°,360°)范围内以等角度间隔采样90个投影。真实参考图像是利用360个投影数据重建获得，此时X射线管电压是90kvp，电流是0.5mA。探测器曝光时间是500ms；

(3)实验环境是matlab2014a，计算机操作系统是window10，CPU模型是Intel i7-7700，内存大小为32G。每个像素的搜索窗设置为9×9，算法最大迭代次数设为500；

(4)初始估计图像u_init；

(5)用TGV正则项方法重建初始估计图像u_init，从而重构权矩阵W₁,W₂,W₃

w₂(x,y)＝α·w₁(x,y)

w₃(x)＝β·∫w₁(x,y)dy

(6)设置初始化参数N_max，容忍误差ε＞0，k＝0,u⁰＝0；

(7)更新参数X^k+1,Y^k+1,Z^k+1；

(8)通过下式更新中间重建结果u；

(9)更新拉格朗日乘子

(10)如果||u^k+1-u^k||＜ε中断；

(11)更新迭代次数k＝k+1；

(12)输出最终重建结果u，如图4(a)和图4(b)所示。

Claims (1)

1.基于离散NLHTV正则项的稀疏角度采样CT迭代重建方法，其特征在于，包括以下步骤：

步骤1CT成像模型建立；

设X射线是单色的，X射线沿传播路径L穿过物体的衰减，由在数学上由比尔定律描述，公式为

其中，I₀是X射线源发射的初始强度；I_i是X射线穿过物体后的X射线强度，μ(l)是沿着透射路径的X射线衰减系数；

被重建区域需要离散为N₁×N₂大小的矩形栅格；被重构断面图像表达为u＝[u₁,u₂,···,u_j,···,u_N]^T，u_j表示第j个像素值，N＝N₁×N₂是像素总数；经过离散化后，公式(1)中的线积分方程表达为线性方程；针对第i个射线，获得如下表达式

其中，a_i,j表达为对投影数值p_i做出贡献的像素u_j的权系数，考虑X射线穿过物体的所有路径，可表达为以下线性方程

Au＝p+e(3)

其中，A∈R^M×N是投影矩阵，a_i,j≥0是矩阵A的第i行，第j列的元素；矢量p＝[p₁,p₂,···,p_i,···,p_M]^T代表所有投影数据，M是被测量投影数据的总数；e∈R^M×1表达为测量误差，包括探测器噪声和X射线的光子噪声；

步骤2定义NLHTV正则项；

在全广义变差正则化定义基础上，引入加性稀疏正则化项，提出一种非局部混合变分正则化NLHTV，如下所示

其中，Ω为定义域，符号u(x)和u(y)分别代表空间位置x和y处的像素值，g(x)和g(y)分别表示x和y处的非局部一阶变分，g^m(x)代表g(x)的第m个分量，同理g^m(y)代表g(y)的第m个分量；w₁,w₂,w₃表示支撑权值，定义如下：

w₂(x,y)＝α·w₁(x,y)

w₃(x)＝β·∫w₁(x,y)dy

其中，C是归一化因子，h_c,h_l,h_g均为大于0的常数，权值w₁,w₂,w₃衰减快慢由h_c,h_l,h_g决定，α，β为大于0的常数；

步骤3离散化NLHTV

在NLHTV正则项的离散定义基础上，基于离散NLHTV正则项的稀疏角度CT图像重建方法能使得损失函数最小化，表达为下式

其中，参数μ＞0用来平衡NLHTV正则项和数据拟合项；i，j为像素索引；d_ij＝|i-j|；u_i,u_j代表索引为i，j的像素值，即u(x)在iΔl与u(y)在jΔl处的像素值，Δl为变量x或y的离散化的间隔；g_i，g_j代表一阶索引为i，j处像素的非局部一阶变分，即g(x)在iΔl处的值与g(y)在jΔl处的值；w₁(i,j),w₂(i,j)表示支撑权值w₁(x,y)，w₂(x,y)的离散形式，即在w₁(x,y)，w₂(x,y)在(iΔl,jΔl)处的值；w₃(n)代表∑_mw₃(nΔl,mΔl)，m与n为像素索引，w₃(nΔl,mΔl)表示w₃(x,y)在(nΔl,mΔl)处的值；

上述CT图像重建模型的优点在于可同时利用非局部一阶、非局部二阶变分正则化项；

步骤4求解NLHTV重建模型

为了便于公式推导，以矩阵表示形式重写离散NLHTV正则项，表示为

其中，u是重建图像，F是非局部的一阶差分矩阵，

表示由整个图像的非局部一阶导数构成的向量；矩阵D代表与非局部差分相对应的距离矩阵，其元素由d_i,j构成；W₁,W₂,W₃分别代表与非局部差分相对应的权重矩阵，其元素分别由{w₁}，{w₂}，{w₃}构成，⊙表示两个矩阵的元素乘法；

步骤5求解子问题；

为求解问题(6)，设计了基于ADMM的优化求解方案，将整个优化问题分解为几个相对简单子问题，从而迭代求解；通过引入辅助变量X,Y,Z，将上述问题表述为

优化问题(7)对应的拉格朗日函数为

其中，ρ＞0为二次项惩罚系数，λ₁、λ₂、λ₃为拉格朗日乘子；

ADMM方法通过迭代优化增广拉格朗日函数中的X,Y,Z变量求解原始问题，对应的迭代公式表达如下

其中，上标k表示迭代次数；

步骤6子问题求解

下面给出每个子问题求解的细节；

1)求解X,Y,Z子问题

关于X,Y,Z的问题相当于以下形式：

上述问题解析解表达可表示为

其中，其中，sign()表示符号函数；换言之，给定x≥0,sign(x)＝1，否则sign(x)＝0；max(x)函数表示取x中的最大值，ρ的意义同公式(8)中；

2)求解子问题u

固定X^k+1,Y^k+1,Z^k+1，子问题u是最小平方问题，可令

式中关于u的梯度等于0，即/>

符号/>

代表求梯度，可以进一步写为下式：

将公式(12)表达为矩阵方程形式，得

式(13)为求解线性方程组问题，涉及求解大规模矩阵求逆；由于直接获得逆矩阵非常困难，因此采用迭代求解方法；为提高求解子问题u的效率，以加快收敛速度快，提数值稳定性高，采用共轭梯度法求解方程(13)；

3)更新乘子

在迭代过程中，设置误差容忍度ε＞0，判断如果||u^k+1-u^k||≥ε，更新迭代次数k＝k+1；如果||u^k+1-u^k||＜ε，则输出最终重建结果u。

Images