CN104103086B - 一种稀疏采样角度下基于变分不等式的ct图像重建方法 - Google Patents

Abstract

本发明公开了一种稀疏采样角度下基于变分不等式的CT图像重建方法，包括以下步骤：(1)针对扇束CT，在0到180度角度范围内进行等角度间隔投影扫描，投影方向角度间隔的大小在2到6度之间，获取稀疏投影数据y；(2)通过X射线源，探测器以及待重建物体的位置信息计算投影矩阵A；(3)根据步骤(1)中得到的投影数据y以及步骤(2)中得到的投影矩阵A，同时引入图像梯度的稀疏性以及非负性作为先验知识得到稀疏采样角度下图像重建问题的重建模型；(4)将步骤(3)中的重建模型转化为变分不等式形式；(5)求解步骤(4)中的变分不等式，得到重建后的图像。本发明在保证图像重建质量的前提下，可以提升重建的收敛速度并降低单次迭代的时间。

Description

一种稀疏采样角度下基于变分不等式的CT图像重建方法

技术领域

本发明涉及精密电子封装过程中元器件内部图像重建以及医学CT图像重建领域，特别是涉及稀疏采样角度下基于变分不等式的CT图像重建方法。

背景技术

计算机断层成像技术(Computed Tomography,CT)是X射线照相技术与复杂的计算机信号处理方法相结合的产物，在不同的方向对物体进行X射线投影，通过测量的投影数据获得物体的横截面信息，可以准确而又直观地重构出物体内部的结构。目前CT技术已被广泛应用于安全检查，工业无损探伤以及医学诊断等领域。

在CT重建的过程中，我们往往只能获得不完全的，稀疏采样角度下的投影数据。如考虑到集成电路大规模生产的需求，要求短时间内完成封装元器件的内部图像的重建并进行质量分析，因此，投影数据采集以及重建的速度对于提高集成电路检测的效率至关重要，相比于每隔1度进行一次投影，增大投影角度的间隔，采集稀疏角度的投影数据无疑会大幅度提升数据采集的速度，有助于缺陷检测过程的快速进行；在医学领域，考虑到X光辐射对人体的危害性，高剂量的X光辐射会增加人体患疾病的风险，在满足CT重建质量要求的情况下，希望能够尽可能减少投影的角度，降低X光射线穿过的次数来降低X光的辐射剂量。对于这类稀疏采样角度下图像重建问题，通常使用迭代算法进行求解，但是使用传统的迭代算法，如最大期望(Expectation Maximization，EM)算法，有序子集最大期望(OrderedSubset Expectation Maximization，OSEM)算法，代数重建算法(AlgebraicReconstruction Technique，ART)、联合迭代重建算法(Simultaneous IterativeReconstruction Technique，SIRT)并不能得到高质量的重建图像。为了提升重建的质量，在重建的过程中考虑正则化的概念，以约束、罚项等形式引入图像的先验知识，如图像梯度的稀疏性以及非负性质。然而这些考虑了正则化的迭代重建方法如凸集投影-全变分(Projection On Convex Sets-Total Variation,POCS-TV)算法由于投影数据的不完全性，对重建的速度产生了重大的影响，往往需要经过很多次的迭代才能获得高质量的重建图像，收敛速度慢，单次迭代时间长。因此，对于稀疏采样角度下的图像重建问题，必须发明新的算法在保证重建质量的前提下，提高重建的收敛速度，降低单次迭代的时间。

发明内容

本发明的目的在于克服现有方法的缺点与不足，提供了稀疏采样角度下基于变分不等式的CT图像重建方法，能够在保证重建质量的前提下，有效提高收敛的速度，并降低单次迭代的时间。

为了达到上述目的，本发明采用以下技术方案：

一种稀疏采样角度下基于变分不等式的CT图像重建方法，包括下述步骤：

(1)针对扇束CT，在0到180度角度范围内进行等角度间隔投影扫描，投影方向角度间隔的大小在2到6度之间，获取稀疏投影数据y；

(2)通过X射线源，探测器以及待重建物体的位置信息计算投影矩阵A；

(3)根据步骤(1)中得到的投影数据y以及步骤(2)中得到的投影矩阵A，同时引入图像梯度的稀疏性以及非负性作为先验知识进行约束得到稀疏采样角度下图像重建问题的重建模型：

其中u为待重建图像，被离散化为一个n×n的图像矩阵；表示图像梯度；矩阵A中的每一个元素A_i,j表示为第i条射线穿过第j个像素的长度；N为待重建图像像素的个数；M为投影数据的个数；

(4)将步骤(3)中的稀疏采样角度下图像重建问题的重建模型转化为变分不等式形式：

<f(u^*),u-u^*>≥0 (2)

其中u^*,u∈C，u^*为变分不等式的解，F(u)＝||u||_TV，f(u)＝F'(u)，

C₁＝{u|||Au-y||₂≤ε}，C₂＝{u|u≥0}，C＝C₁∩C₂；

(5)求解步骤(4)中的变分不等式，得到重建图像，即找到一个待重建的图像u^*，使得对于任意的u∈C，<f(u^*),u-u^*>≥0成立。

优选的，步骤(3)中，||u||_TV的梯度是一个图像，该图像的每个像素值都是||u||_TV相对应像素的偏微分，并用如下表示：

优选的，在重建的过程中，待重建图像u以及||u||_TV的梯度都被转化为列向量，图像中第s行，第t列的像素点对应于列向量中第(s-1)·n+t个像素点。

优选的，步骤(4)中稀疏采样角度下图像重建模型与变分不等式之间的转化，转化过程包括以下步骤：

变分不等式的定义：设G:Rⁿ→Rⁿ为连续映射，为非空闭凸集，存在x^*∈S，使得对于<G(x^*),x-x^*>≥0成立；

(4-1)令C₁＝{u|||Au-y||₂≤ε}，C₂＝{u|u≥0}，C＝C₁∩C₂，F(u)＝||u||_TV，f(u)＝F'(u)

公式(1)变为

argminF(u)s.t.u∈C (3)

其中F(u)是一个凸函数，因为C₁，C₂皆为凸集，所以C＝C₁∩C₂也为凸集

(4-2)令u^*为公式(3)的最优解，那么对于F(u)≥F(u^*)；

(4-3)对于任何的γ∈[0,1]，当u，u^*属于凸集C时，u^*+γ(u-u^*)∈C，于是根据步骤(4-2)可得:

F(u^*+γ(u-u^*))-F(u^*)≥0 (4)

因此对于公式(2)成立；

这表明，u^*也是变分不等式，即公式(2)的解，公式(1)等价于公式(2)。

优选的，步骤(5)所述的变分不等式求解方法，求解过程具体包含以下步骤：

(5-1)初始化：u⁰＝0，d⁰＝1，ζ＝0.2，设置总的迭代次数N_Total，当前迭代次数k＝0；

(5-2)对当前的迭代解u^k应用投影算法，得到中间解

(5-3)计算下一次迭代中的梯度下降法步长d^k+1；

(5-4)将列向量转化为二维图像矩阵二维图像矩阵中第s行，第t列的像素点与列向量中第(s-1)·n+t个像素点对应；

(5-5)对进行平滑操作，得到平滑过后的结果

(5-6)将转化为列向量得到u^k+1，图像中第s行，第t列的像素点对应于列向量中第(s-1)·n+t个像素点；

(5-7)k＝k+1，判断k是否小于设置的总迭代次数N_Total；若k小于N_Total，跳转至步骤(5-2)继续迭代；否则重建结束，将重建结果转化为二维图像矩阵。

优选的，步骤(5-2)所述的投影算法，具体为：

(5-2-1)p^k＝u^k-τf(u^k)；

(5-2-2)z^k＝P_C(p^k)，其中为投影算子，计算p^k在可行集C上的投影。

优选的，步骤(5-2-1)所述p^k的计算，p^k＝u^k-τf(u^k)是梯度下降法，使用如下步骤来实现梯度下降法的过程，其中t为当前迭代的次数，设置固定迭代次数为N_GD，d^k为下降的步长，具体为：

(5-2-1-1)初始化：s⁰＝u^k,t＝0；

(5-2-1-2)计算梯度：

(5-2-1-3)梯度下降法：s^t+1＝s^t-d^k·v/|v|；

(5-2-1-4)t＝t+1，判断t是否小于设置的迭代次数N_GD；若t小于N_GD，跳转至步骤(5-2-1-2)继续迭代；否则迭代停止，

优选的，步骤(5-2-2)所述对p^k在C上投影的计算，C＝C₁∩C₂，由两部分组成，因此求p^k在可行集C上的投影也分为两部分：先求p^k在可行集C₁上的投影z^k，λ为正则化参数，并设置迭代次数t为当前迭代的次数；再限制z^k为非负，得到具体过程为：

(5-2-2-1)初始化:x⁰＝p^k,t＝0；

(5-2-2-2)计算初始梯度G⁰＝2(x⁰-p^k)+λA^T(Ax⁰-y)；

(5-2-2-3)计算初始步长

(5-2-2-4)x^t+1＝x^t-a·G^t；

(5-2-2-5)更新梯度G^t+1＝2(x^t+1-p^k)+λA^T(Ax^t+1-y)；

(5-2-2-6)更新迭代的步长

(5-2-2-7)t＝t+1，判断t是否小于设置的迭代次数若t小于跳转至步骤(5-2-2-4)继续迭代；否则迭代停止，并跳转至步骤(5-2-2-8)；

(5-2-2-8)正定性限制，限制其为非负：

优选的，步骤(5-3)所述梯度下降法步长d^k+1的计算：

优选的，步骤(5-5)所述的对进行平滑操作，得到平滑过后的结果具体为：

(5-5-1)针对中的每一个像素如计算与其八邻域像素之间的绝对值差e_i-1,j-1，e_i-1,j，e_i-1,j+1，e_i,j-1，e_i,j+1，e_i+1,j-1，e_i+1,j，e_i+1,j+1

(5-5-2)根据绝对值差计算相应的权重w_i-1,j-1，w_i-1,j，w_i-1,j+1，w_i,j-1，w_izj+1，w_i+1,j-1，w_i+1,j，w_i+1,j+1，其中η为尺度因子，可调节像素的权重值；

(5-5-3)根据步骤(5-5-1)中得到的绝对值差判断其八邻域像素点与当前像素点是否属于同一个集合，若不属于同一集合，则设置相应的权重为0，其中thr为判断的阈值：

(5-5-4)进行平滑操作

本发明与现有技术相比，具有如下优点和有益效果：

1、本发明将稀疏采样角度下图像重建问题转化为变分不等式形式，并结合投影算法以及图像平滑操作对变分不等式进行求解。另外，本发明中投影算法计算量小，适用于解决大规模问题；图像平滑操作可以有效降低投影数据噪声对重建的影响，保证图像的分段光滑性。

2、本发明在稀疏采样角度的投影数据下，能够在保证图像重建质量的前提下，提高重建的收敛速度，同时降低单次迭代的时间，实现重建过程的快速性。

附图说明

图1为本发明所述的稀疏采样角度下基于变分不等式的CT图像重建方法的流程示意图；

图2使用的Shepp-Logan模型原图；

图3(a)为POCS-TV算法得到的重建结果图；

图3(b)为本发明方法得到的重建结果图；

图4(a)为图3(a)中POCS-TV算法重建结果部分区域放大图；

图4(b)为图3(b)中通过本发明方法重建结果部分区域放大图；

图5为模拟投影数据下每次迭代得到的重建图像与原始图像的误差曲线图。

具体实施方式

下面结合实施例及附图对本发明作进一步详细的描述，但本发明的实施方式不限于此。

如图1所示，稀疏采样角度下基于变分不等式的CT图像重建方法，包括以下步骤：

(1)针对扇束CT，在0到180度角度范围内进行等角度间隔投影扫描，投影方向角度间隔的大小在2到6度之间，获取稀疏投影数据y；

(2)通过X射线源，探测器以及待重建物体的位置信息计算投影矩阵A；

(3)根据步骤(1)中得到的投影数据y以及步骤(2)中得到的投影矩阵A，同时引入图像梯度的稀疏性以及非负性作为先验知识进行约束得到稀疏采样角度下图像重建问题的重建模型：

其中u为待重建图像，被离散化为一个n×n的图像矩阵；表示图像梯度；矩阵A中的每一个元素A_i,j表示为第i条射线穿过第j个像素的长度；N为待重建图像像素的个数；M为投影数据的个数；||u||_TV的表达式中引入的ξ为一个较小的正数，使得对任意一个像素，||u||_TV的偏微分都是存在的，选择ξ＝10^-8；||u||_TV的梯度是一个图像，该图像的每个像素值都是||u||_TV相对应像素的偏微分，并用如下表示：

在重建的过程中，待重建图像u以及||u||_TV的梯度都被转化为列向量，图像中第s行，第t列的像素点对应于列向量中第(s-1)·n+t个像素点；

(4)将步骤(3)中的稀疏采样角度下图像重建问题的重建模型转化为变分不等式形式：

<f(u^*),u-u^*>≥0 (2)

变分不等式的定义：设G:Rⁿ→Rⁿ为连续映射，为非空闭凸集，存在x^*∈S，使得对于<G(x^*),x-x^*>≥0成立。

(4-1)令C₁＝{u|||Au-y||₂≤ε}，C₂＝{u|u≥0}，C＝C₁∩C₂，F(u)＝||u||_TV，f(u)＝F'(u)

公式(1)变为：

argminF(u)s.t.u∈C (3)

其中F(u)是一个凸函数，因为C₁，C₂皆为凸集，所以C＝C₁∩C₂也为凸集；

(4-2)令u^*为公式(3)的最优解，那么对于F(u)≥F(u^*)；

(4-3)对于任何的γ∈[0,1]，当u，u^*属于凸集C时，u^*+γ(u-u^*)∈C，于是根据步骤(4-2)可得:

F(u^*+γ(u-u^*))-F(u^*)≥0 (4)

因此对于公式(2)成立。

这表明，u^*也是变分不等式，即公式(2)的解，公式(1)等价于公式(2)。

(5)求解步骤(4)中的变分不等式，得到重建图像，即找到一个u^*(待重建的图像)，使得对于任意的u∈C，<f(u^*),u-u^*>≥0成立。

(5-1)初始化：u⁰＝0，d⁰＝1，ζ＝0.2，设置总的迭代次数N_Total，当前迭代次数k＝0；

(5-2)对当前的迭代解u^k应用投影算法，得到中间解

(5-2-1)p^k＝u^k-τf(u^k)；

p^k＝u^k-τf(u^k)是梯度下降法，使用如下步骤来实现梯度下降法的过程，其中t为当前迭代的次数，设置梯度下降法实现过程中的固定迭代次数N_GD＝20，d^k为下降的步长：

(5-2-1-1)初始化：s⁰＝u^k,t＝0；

(5-2-1-2)计算梯度：

(5-2-1-3)梯度下降法：s^t+1＝s^t-d^k·v/|v|；

(5-2-1-4)t＝t+1，判断t是否小于设置的迭代次数N_GD。若t小于N_GD，跳转至步骤(5-2-1-2)继续迭代；否则迭代停止，(5-2-2)z^k＝P_C(p^k)，其中为投影算子，计算p^k在可行集C上的投影。C＝C₁∩C₂，由两部分组成，因此求p^k在可行集C上的投影也分为两部分：先求p^k在可行集C₁上的投影z^k，正则化参数λ＝100，并设置其迭代次数为10，t为当前迭代的次数；再限制z^k为非负，得到具体过程为：

(5-2-2-1)初始化:x⁰＝p^k,t＝0；

(5-2-2-2)计算初始梯度G⁰＝2(x⁰-p^k)+λA^T(Ax⁰-y)；

(5-2-2-3)计算初始步长

(5-2-2-4)x^t+1＝x^t-a·G^t；

(5-2-2-5)更新梯度G^t+1＝2(x^t+1-p^k)+λA^T(Ax^t+1-y)；

(5-2-2-6)更新迭代的步长

(5-2-2-7)t＝t+1，判断t是否小于设置的迭代次数若t小于跳转至步骤(5-2-2-4)继续迭代；否则迭代停止，并跳转至步骤(5-2-2-8)；

(5-2-2-8)正定性限制，限制其为非负：

(5-3)计算下一次迭代中的梯度下降法步长d^k+1：

(5-4)将列向量转化为二维图像矩阵二维图像矩阵中第s行，第t列的像素点与列向量中第(s-1)·n+t个像素点对应；

(5-5)对进行平滑操作，得到平滑过后的结果具体包括以下步骤：

(5-5-1)针对中的每一个像素如计算与其八邻域像素之间的绝对值差e_i-1,j-1，e_i-1,j，e_i-1,j+1，e_i,j-1，e_i,j+1，e_i+1,j-1，e_i+1,j，e_i+1,j+1

(5-5-2)根据绝对值差计算相应的权重w_i-1,j-1，w_i-1,j，w_i-1,j+1，w_i,j-1，w_i,j+1，w_i+1,j-1，w_i+1,j，w_i+1,j+1，其中尺度因子η设置为0.1，可调节像素的权重值；

(5-5-3)根据步骤(5-5-1)中得到的绝对值差判断其八邻域像素点与当前像素点是否属于同一个集合，若不属于同一集合，则设置相应的权重为0，其中thr为判断的阈值，设置thr＝0.035：

(5-5-4)进行平滑操作

将转化为列向量得到u^k+1，图像中第s行，第t列的像素点对应于列向量中第(s-1)·n+t个像素点；

(5-7)k＝k+1，判断k是否小于设置的总迭代次数N_Total。若k小于N_Total，跳转至步骤(5-2)继续迭代；否则重建结束，将重建结果转化为二维图像矩阵。

为了说明本发明的重建效果及优势，采用模拟投影数据对Shepp-Logan模型进行重建，并在相同的迭代次数下与POCS-TV算法进行比较。二维Shepp-Logan图像大小为128×128。假设模拟投影数据获取及重建过程的参数如下：

(1)X-Ray射线源与待重建物体的中心距离为256mm；

(2)线探测器的中心线经过待重建物体的中心，且距离为256mm；

(3)探测器的个数为256个，每个探测器单元的宽度为0.5mm；

(4)投影的角度范围为(0，180°),投影方向的间隔角度为6°；

(5)待重建物体的断层图像的每个像素的大小为0.5mm×0.5mm；

(6)设置总的迭代次数N_Total＝100；

(7)扫描方式为扇束扫描。

重建结果如图3(a)、图3(b)所示，对各自图中方框内区域进行放大得到图4(a)和图4(b)。从图3(a)、图3(b)以及图4(a)和图4(b)可知，本发明方法得到的重建结果在边缘位置更加的清晰，模糊程度远低于POCS-TV算法，重建质量好于POCS-TV算法。

计算Error＝||u^k-u_true||₂来表明本发明方法和POCS-TV算法在模拟投影数据下每次迭代得到的重建图像与原始图像的误差，如图5所示，可知本发明方法的收敛速度明显快于POCS-TV算法。且本发明方法的单次迭代时间只需要0.2389s，远低于POCS-TV算法的2.5456s，保证了重建过程的快速性。

上述实施例为本发明较佳的实施方式，但本发明的实施方式并不受上述实施例的限制，其他的任何未背离本发明的精神实质与原理下所作的改变、修饰、替代、组合、简化，均应为等效的置换方式，都包含在本发明的保护范围之内。

Claims (10)

1.一种稀疏采样角度下基于变分不等式的CT图像重建方法，其特征在于，包括下述步骤：

(1)针对扇束CT，在0到180度角度范围内进行等角度间隔投影扫描，投影方向角度间隔的大小在2到6度之间，获取稀疏投影数据y；

(2)通过X射线源，探测器以及待重建物体的位置信息计算投影矩阵A；

(3)根据步骤(1)中得到的投影数据y以及步骤(2)中得到的投影矩阵A，同时引入图像梯度的稀疏性以及非负性作为先验知识进行约束得到稀疏采样角度下图像重建问题的重建模型：

$\begin{array}{cc}\underset{u}{\mathrm{argmin}}\left|\right|u|{|}_{TV}& s.t.\left|\right|Au-y|{|}_2\&le;\mathrm{\&epsiv;},u\&GreaterEqual;0\end{array}---\left(1\right)$

其中u为待重建图像，被离散化为一个n×n的图像矩阵；表示图像梯度；矩阵A中的每一个元素A_i,j表示为第i条射线穿过第j个像素的长度；N为待重建图像像素的个数；M为投影数据的个数；

(4)将步骤(3)中的稀疏采样角度下图像重建问题的重建模型转化为变分不等式形式：

<f(u^*),u-u^*>≥0 (2)

其中u^*,u∈C，u^*为变分不等式的解，F(u)＝||u||_TV，f(u)＝F'(u)，

C₁＝{u|||Au-y||₂≤ε}，C₂＝{u|u≥0}，C＝C₁∩C₂；

(5)求解步骤(4)中的变分不等式，得到重建图像，即找到一个待重建的图像u^*，使得对于任意的u∈C，<f(u^*),u-u^*>≥0成立。

2.根据权利要求1所述的稀疏采样角度下基于变分不等式的CT图像重建方法，其特征在于，步骤(3)中，||u||_TV的梯度是一个图像，该图像的每个像素值都是||u||_TV相对应像素的偏微分，并用如下表示：

$\begin{array}{c}\frac{\&part;\left|\right|u|{|}_{TV}}{\&part;u_{s,t}}=\frac{(u_{s,t}-u_{s-1,t})+(u_{s,t}-u_{s,t-1})}{\sqrt{\mathrm{\&xi;}+{(u_{s,t}-u_{s-1,t})}^2+{(u_{s,t}-u_{s,t-1})}^2}}\\ -\frac{(u_{s+1,t}-u_{s,t})}{\sqrt{\mathrm{\&xi;}+{(u_{s+1,t}-u_{s,t})}^2+{(u_{s+1,t}-u_{s+1,t-1})}^2}}\\ -\frac{(u_{s,t+1}-u_{s,t})}{\sqrt{\mathrm{\&xi;}+{(u_{s,t+1}-u_{s,t})}^2+{(u_{s,t+1}-u_{s-1,t+1})}^2}}.\end{array}$

3.根据权利要求1所述的稀疏采样角度下基于变分不等式的CT图像重建方法，其特征在于，在重建的过程中，待重建图像u以及||u||_TV的梯度都被转化为列向量，图像中第s行，第t列的像素点对应于列向量中第(s-1)·n+t个像素点。

4.根据权利要求1所述的稀疏采样角度下基于变分不等式的CT图像重建方法，其特征在于，步骤(4)中稀疏采样角度下图像重建模型与变分不等式之间的转化，转化过程包括以下步骤：

变分不等式的定义：设G:Rⁿ→Rⁿ为连续映射，为非空闭凸集，存在x^*∈S，使得对于<G(x^*),x-x^*>≥0成立；

(4-1)令C₁＝{u|||Au-y||₂≤ε}，C₂＝{u|u≥0}，C＝C₁∩C₂，F(u)＝||u||_TV，f(u)＝F'(u)公式(1)变为

argminF(u)s.t.u∈C (3)

其中F(u)是一个凸函数，因为C₁，C₂皆为凸集，所以C＝C₁∩C₂也为凸集

(4-2)令u^*为公式(3)的最优解，那么对于F(u)≥F(u^*)；

(4-3)对于任何的γ∈[0,1]，当u，u^*属于凸集C时，u^*+γ(u-u^*)∈C，于是根据步骤(4-2)可得:

F(u^*+γ(u-u^*))-F(u^*)≥0 (4)

$\begin{array}{c}\underset{\mathrm{\&gamma;}\&RightArrow;0^{+}}{\lim }\frac{F(u^{*}+\mathrm{\&gamma;}(u-u^{*}\left)\right)-F\left(u^{*}\right)}{\mathrm{\&gamma;}}\\ =\underset{\mathrm{\&gamma;}\&RightArrow;0^{+}}{\lim }\frac{F(u^{*}+\mathrm{\&gamma;}(u-u^{*}\left)\right)-F\left(u^{*}\right)}{\mathrm{\&gamma;}(u-u^{*})}(u-u^{*})\\ =F^{\&prime;}\left(u^{*}\right)(u-u^{*})=f\left(u^{*}\right)(u-u^{*})\&GreaterEqual;0\end{array}---\left(5\right)$

因此对于公式(2)成立；

这表明，u^*也是变分不等式，即公式(2)的解，公式(1)等价于公式(2)。

5.根据权利要求1所述的稀疏采样角度下基于变分不等式的CT图像重建方法，其特征在于，步骤(5)所述的变分不等式求解方法，求解过程具体包含以下步骤：

(5-1)初始化：u⁰＝0，d⁰＝1，ζ＝0.2，设置总的迭代次数N_Total，当前迭代次数k＝0；

(5-2)对当前的迭代解u^k应用投影算法，得到中间解

(5-3)计算下一次迭代中的梯度下降法步长d^k+1；

(5-4)将列向量转化为二维图像矩阵二维图像矩阵中第s行，第t列的像素点与列向量中第(s-1)·n+t个像素点对应；

(5-5)对进行平滑操作，得到平滑过后的结果

(5-6)将转化为列向量得到u^k+1，图像中第s行，第t列的像素点对应于列向量中第(s-1)·n+t个像素点；

(5-7)k＝k+1，判断k是否小于设置的总迭代次数N_Total；若k小于N_Total，跳转至步骤(5-2)继续迭代；否则重建结束，将重建结果转化为二维图像矩阵。

6.根据权利要求5所述的稀疏采样角度下基于变分不等式的CT图像重建方法，其特征在于，步骤(5-2)所述的投影算法，具体为：

(5-2-1)p^k＝u^k-τf(u^k)；

(5-2-2)z^k＝P_C(p^k)，其中为投影算子，计算p^k在可行集C上的投影。

7.根据权利要求6所述的稀疏采样角度下基于变分不等式的CT图像重建方法，其特征在于，步骤(5-2-1)所述p^k的计算，p^k＝u^k-τf(u^k)是梯度下降法，使用如下步骤来实现梯度下降法的过程，其中t为当前迭代的次数，设置固定迭代次数为N_GD，d^k为下降的步长，具体为：

(5-2-1-1)初始化：s⁰＝u^k,t＝0；

(5-2-1-2)计算梯度：

(5-2-1-3)梯度下降法：s^t+1＝s^t-d^k·v/|v|；

(5-2-1-4)t＝t+1，判断t是否小于设置的迭代次数N_GD；若t小于N_GD，跳转至步骤(5-2-1-2)继续迭代；否则迭代停止，

8.根据权利要求6所述的稀疏采样角度下基于变分不等式的CT图像重建方法，其特征在于，步骤(5-2-2)所述对p^k在C上投影的计算，C＝C₁∩C₂，由两部分组成，因此求p^k在可行集C上的投影也分为两部分：先求p^k在可行集C₁上的投影z^k，λ为正则化参数，并设置迭代次数t为当前迭代的次数；再限制z^k为非负，得到具体过程为：

(5-2-2-1)初始化:x⁰＝p^k,t＝0；

(5-2-2-2)计算初始梯度G⁰＝2(x⁰-p^k)+λA^T(Ax⁰-y)；

(5-2-2-3)计算初始步长

(5-2-2-4)x^t+1＝x^t-a·G^t；

(5-2-2-5)更新梯度G^t+1＝2(x^t+1-p^k)+λA^T(Ax^t+1-y)；

(5-2-2-6)更新迭代的步长

(5-2-2-7)t＝t+1，判断t是否小于设置的迭代次数若t小于跳转至步骤(5-2-2-4)继续迭代；否则迭代停止，并跳转至步骤(5-2-2-8)；

(5-2-2-8)正定性限制，限制其为非负：

9.根据权利要求5所述的稀疏采样角度下基于变分不等式的CT图像重建方法，其特征在于，步骤(5-3)所述梯度下降法步长d^k+1的计算：

$d^{k+1}=\mathrm{\&zeta;}\left|\right|z_{pos}^k-z^k\left|\right|.$

10.根据权利要求5所述的稀疏采样角度下基于变分不等式的CT图像重建方法，其特征在于，步骤(5-5)所述的对进行平滑操作，得到平滑过后的结果具体为：

(5-5-1)针对中的每一个像素如计算与其八邻域像素之间的绝对值差e_i-1,j-1，e_i-1,j，e_i-1,j+1，e_i,j-1，e_i,j+1，e_i+1,j-1，e_i+1,j，e_i+1,j+1；

(5-5-2)根据绝对值差计算相应的权重w_i-1,j-1，w_i-1,j，w_i-1,j+1，w_i,j-1，w_i,j+1，w_i+1,j-1，w_i+1,j，w_i+1,j+1，其中η为尺度因子，可调节像素的权重值；

$w=\exp (-{\left(\frac{e}{\mathrm{\&eta;}}\right)}^2)$

(5-5-3)根据步骤(5-5-1)中得到的绝对值差判断其八邻域像素点与当前像素点是否属于同一个集合，若不属于同一集合，则设置相应的权重为0，其中thr为判断的阈值：

$w=\left\{\begin{array}{cc}w& e<thr\\ 0& e>thr\end{array}\right.$

(5-5-4)进行平滑操作

$z_{Smooth}^k(i,j)=Smooth\left(z_{Matrix}^k\right(i,j\left)\right)=\frac{\begin{array}{c}{z}_{Matrix}^k(i-1,j-1)\&CenterDot;w_{i-1,j-1}+z_{Matrix}^k(i-1,j)\&CenterDot;w_{i-1,j}\\ +z_{Matrix}^k(i-1,j+1)\&CenterDot;w_{i-1,j+1}+z_{Matrix}^k(i,j-1)\&CenterDot;w_{i,j-1}\\ +z_{Matrix}^k(i,j+1)\&CenterDot;w_{i,j+1}+z_{Matrix}^k(i+1,j-1)\&CenterDot;w_{i+1,j-1}\\ +z_{Matrix}^k(i+1,j)\&CenterDot;w_{i+1,j}+z_{Matrix}^k(i+1,j+1)\&CenterDot;w_{i+1,j+1}\end{array}}{\begin{array}{c}{w}_{i-1,j-1}+w_{i-1,j}+w_{i-1,j+1}+w_{i,j-1}\\ +w_{i,j+1}+w_{i+1,j-1}+w_{i+1,j}+w_{i+1,j+1}\end{array}}.$