CN110717959B - 基于曲率约束的x射线有限角ct图像重建方法和装置 - Google Patents

Abstract

本发明公开了一种基于曲率约束的X射线有限角CT图像重建方法和装置，包括：步骤1，输入已知数据集；步骤2，初始化；步骤3，迭代处理估计图像u^(k)，利用有限角CT扫描数据集p，通过与扫描几何参数集G相关的图像重建算子R_G来更新估计图像u^(k)，得到u^(k+1/3)＝R_G(u^(k)，p)，R_G为扫描几何参数集G相关的图像重建算子；步骤4，利用梯度的稀疏性约束图像u^(k+1/3)，得到的u^(k+2/3)＝P₁(u^(k+1/3))，P₁为约束图像梯度稀疏性的算子；步骤5，利用边界的光滑性约束图像u^(k+2/3)，得到的图像u^(k+1)＝P₂(u^(k+2/3))，完成一次迭代，P₂为约束图像边界光滑性的算子；步骤6，判断相邻两次迭代图像间的差别是否小于给定阈值，若不满足，转至步骤3，若满足，则终止迭代。本发明能够克服现有的有限角CT图像重建算法中边界模糊或存在阶梯效应的问题。

Description

基于曲率约束的X射线有限角CT图像重建方法和装置

技术领域

本发明涉及CT图像重建技术领域，特别是关于一种基于曲率约束的X射线有限角CT图像重建方法和装置。

背景技术

工业无损检测中，常遇到电路板、机翼等的检测问题，由于待测物品在某个方向具有较长的边界，可能出现以下情况：1、为了能够重建高分辨率的图像，要求射线源到物品之间的距离较小，为了避免射线源与物品之间的碰撞，扫描过程中转台无法旋转整圈；2、虽然转台能够旋转整圈，但部分X射线由于穿过较长的物体，衰减严重，造成对应角度的探测数据信噪比差而无法使用。这些情况均导致仅能在一定的角度范围内采集得到有效的投影数据。

使用一定角度范围内采集的数据重建得到图像即有限角CT图像重建，隶属于不完全数据图像重建问题。采用传统的图像重建算法(如FDK、SART等)得到的图像，存在部分边界的严重模糊，影响图像质量，进而影响判断。

因此，希望有一种技术方案来克服或至少减轻现有技术的上述缺陷中的至少一个。

发明内容

本发明的目的在于提供一种基于曲率约束的X射线有限角CT图像重建方法来克服或至少减轻现有技术的上述缺陷中的至少一个。

为实现上述目的，本发明提供一种基于曲率约束的X射线有限角CT图像重建方法，该方法包括：

步骤1，输入已知数据集：有限角CT扫描获得的数据集p，CT扫描几何参数集G；

步骤2，初始化：初始估计图像u⁽⁰⁾，确定迭代终止阈值ε以及迭代次数上限N；

步骤3，迭代处理估计图像u^(k)，k＝0，1，…，N，利用有限角CT扫描数据集p，通过与扫描几何参数集G相关的图像重建算子R_G来更新估计图像u^(k)，得到u^(k+1/3)＝R_G(u^(k)，p)，其中，R_G表示扫描几何参数集G相关的图像重建算子；

步骤4，利用梯度的稀疏性约束图像u^(k+1/3)，得到的u^(k+2/3)＝P₁(u^(k+1/3))，其中，P₁表示约束图像梯度稀疏性的算子；

步骤5，利用边界的光滑性约束图像u^(k+2/3)，得到的图像u^(k+1)＝P₂(u^(k+2/3))，完成一次迭代，其中，P₂表示约束图像边界光滑性的算子；

步骤6，判断相邻两次迭代图像间的差别是否小于给定阈值，即||u^(k+1)-u^(k)||≤ε或是否达到迭代次数上限N，若不满足，转至步骤3，若满足，则终止迭代。

进一步地，所述步骤4中，P₁由如下式(1)的最优化问题所定义：

式(1)中，

为正则化项，

为数据保真项，υ表示输入图像，

表示输入图像水平方向的梯度，

表示输入图像竖直方向的梯度，||·||₀表示l₀范数，||·||₁表示l₁范数，||·||₂表示l₂范数，T表示转置，λ₁、ω₁、ω₂为调节参数。

进一步地，所述步骤4中，P₁由如下式(2)的最优化问题所定义：

式(2)中，

为正则化项，

为数据保真项，υ表示输入图像，

表示输入图像水平方向的梯度，

表示输入图像竖直方向的梯度，||·||₀表示l₀范数，||·||₁表示l₁范数，||·||₂表示l₂范数，T表示转置，λ₂、ω₃、ω₄为调节参数。

进一步地，所述步骤5中，P₂由如下式(3)的最优化问题所定义：

式(3)中，argmin_υλ₃||Φ(C(υ))||₁为正则化项，

为数据保真项，Φ(·)表示惩罚函数作用于图像中的每一点处的平均曲率，||·||₁表示l₁范数，||·||₂表示l₂范数，λ₃为调节参数，C(·)表示计算图像每一点的平均曲率，对于第j点，其平均曲率为式(4)：

该优化问题的求解通过将υ视作三维曲面υ(x，y)、将u^(k+2/3)视作三维曲面u^(k+2/3)(x，y)进行采样，通过求解如下偏微分方程得到：

其中，x和y均为空间变量，t表示时间变量，υ(0；x，y)表示t＝0时的υ的取值，υ_t(t，x，y)表示v(t，x，y)关于t的导函数，

表示梯度，

·表示散度，λ₄为调节参数，Φ′表示Φ的导函数。

进一步地，所述步骤2中，R_G选择迭代类重建算法或解析类重建算法获得。

本发明还提供一种基于曲率约束的X射线有限角CT图像重建装置，其包括：

数据输入模块，其用于输入已知数据集，分别为：有限角CT扫描获得的数据集p以及CT扫描几何参数集G；

初始化模块，其用于初始估计图像u⁽⁰⁾，确定迭代终止阈值ε以及迭代次数上限N；

迭代处理估计模块，其用于迭代处理估计图像u^(k)，k＝0，1，…，N，利用有限角CT扫描数据集p，通过与扫描几何参数集G相关的图像重建算子R_G来更新估计图像u^(k)，得到u^{(k +1/3)}＝R_G(u^(k)，p)，其中，R_G表示扫描几何参数集G相关的图像重建算子；

梯度稀疏性约束模块，其用于利用梯度的稀疏性约束图像u^(k+1/3)，得到的u^(k+2/3)＝P₁(u^(k+1/3))，其中，P₁表示约束图像梯度稀疏性的算子；

边界光滑性约束模块，其用于利用边界的光滑性约束图像u^(k+2/3)，得到的图像u^{(k +1)}＝P₂(u^(k+2/3))，完成一次迭代，其中，P₂表示约束图像边界光滑性的算子；

判断模块，其用于判断相邻两次迭代图像间的差别是否小于给定阈值，即||u^(k+1)-u^(k)||≤ε或是否达到迭代次数上限N，若不满足，转至步骤3，若满足，则终止迭代。

进一步地，所述梯度稀疏性约束模块中的P₁由如下式(1)的最优化问题所定义：

式(1)中，

为正则化项，

为数据保真项，υ表示输入图像，

表示输入图像水平方向的梯度，

表示输入图像竖直方向的梯度，||·||₀表示l₀范数，||·||₁表示l₁范数，||·||₂表示l₂范数，T表示转置，λ₁、ω₁、ω₂为调节参数。

进一步地，所述梯度稀疏性约束模块中的P₁由如下式(2)的最优化问题所定义：

式(2)中，

为正则化项，

为数据保真项，υ表示输入图像，

表示输入图像水平方向的梯度，

表示输入图像竖直方向的梯度，||·||₀表示l₀范数，||·||₁表示l₁范数，||·||₂表示l₂范数，T表示转置，λ₂、ω₃、ω₄为调节参数。

进一步地，所述边界光滑性约束模块中的P₂由如下式(3)的最优化问题所定义：

式(3)中，argmin_υλ₃||Φ(C(υ))||₁为正则化项，

为数据保真项，Φ(·)表示惩罚函数作用于图像中的每一点处的平均曲率，||·||₁表示l₁范数，||·||₂表示l₂范数，λ₃为调节参数，C(·)表示计算图像每一点的平均曲率，对于第j点，其平均曲率为式(4)：

该优化问题的求解通过将υ视作三维曲面υ(x，y)、将u^(k+2/3)视作三维曲面u^(k+2/3)(x，y)进行采样，通过求解如下偏微分方程得到：

其中，x和y均为空间变量，t表示时间变量，υ(0；x，y)表示t＝0时的υ的取值，υt(t，x，y)表示v(t，x，y)关于t的导函数，

表示梯度，

·表示散度，λ₄为调节参数，Φ′表示Φ的导函数。

本发明提供的方法利用图像梯度的稀疏性以及边界的光滑性，通过迭代的方式，逐步恢复图像的边界信息，能够有效抑制图像的边界模糊以及阶梯效应，提高重建质量，进而克服了现有的有限角CT图像重建算法中边界模糊或存在阶梯效应的问题。

附图说明

图1为本发明实施例的流程图；

图2为本发明实施例中扫描的样品照片，其中示出了多个不同形状与材质的模体；

图3为本发明实施例中SART的重建图像；

图4为本发明实施例中AEDS的重建图像；

图5为利用本发明实施例获得的重建图像。

具体实施方式

下面结合附图和实施例对本发明进行详细的描述。

本发明实施例提供的基于曲率约束的X射线有限角CT图像重建方法包括如下步骤：

步骤1，输入已知数据集：有限角CT扫描获得的数据集p，CT扫描几何参数集G。其中，CT扫描几何参数集G包括射线源到探测器中心的距离(SDD)、射线源到转台中心的距离(SOD)、探测器单元个数、探测器单元尺寸、扫描角度数和角度采样间隔。

步骤2，初始化：初始估计图像u⁽⁰⁾，一般将其所有像素值设置为0，确定迭代终止阈值ε以及迭代次数上限N。

步骤3，迭代处理估计图像u^(k)，k＝0，1，…，N，利用有限角CT扫描数据集p，通过与扫描几何参数集G相关的图像重建算子R_G来更新估计图像u^(k)，得到u^(k+1/3)＝R_G(u^(k)，p)，其中，R_G表示扫描几何参数集G相关的图像重建算子，u^(k+1/3)表示输入图像。

步骤4，利用梯度的稀疏性约束图像u^(k+1/3)，得到的图像u^(k+2/3)＝P₁(u^(k+1/3))，其中，P₁表示约束图像梯度稀疏性的算子，u^(k+2/3)表示进行约束图像梯度稀疏性处理后得到的图像。其中，图像通过求梯度得到梯度域图像(其包括灰度图像的X方向的梯度，Y方向的梯度和图像像素值)，梯度的稀疏性指的则是梯度域图像中非零元素的多少。通常，梯度域图像中非零元素越少，梯度的稀疏性越好。本实施例利用图像梯度的稀疏性处理图像，一方面可以简化模型，避免过拟合，另一方面可以特征选择，选择参数非零的特征即可，而且使整个模型获得更好的可解释性，只有几个非零参数对结果有影响，更好解释。

步骤5，利用边界的光滑性约束图像u^(k+2/3)，得到的图像u^(k+1)＝P₂(u^(k+2/3))，完成一次迭代，其中，P₂表示约束图像边界光滑性的算子，u^(k+1)表示进行约束图像边界的光滑性处理后得到的图像。

步骤6，判断相邻两次迭代图像间的差别是否小于给定阈值，即||u^(k+1)-u^(k)||≤ε或是否达到迭代次数上限N，若不满足，转至步骤3，若满足，则终止迭代，基于曲率约束的X射线有限角CT图像重建方法结束。

本实施例提供的方法利用图像梯度的稀疏性以及边界的光滑性，通过迭代的方式，逐步恢复图像的边界信息，能够有效抑制图像的边界模糊以及阶梯效应，提高重建质量，进而克服了现有的有限角CT图像重建算法中边界模糊或存在阶梯效应的问题。

上述实施例中，所述步骤4中，P₁由如下式(1)的最优化问题所定义：

式(1)中，

为正则化项，

为数据保真项，υ表示输入图像，

表示输入图像水平方向的梯度，

表示输入图像竖直方向的梯度，||·||₀表示l₀范数，||·||₁表示l₁范数，||·||₂表示l₂范数，T表示转置。λ₁、ω₁、ω₂为调节参数，其中：λ₁用于控制正则化项和数据保真项权重，具体取值也是通过实验得到。原则上，要求ω₁大于十倍的ω₂，具体数值通过实验确定，如λ₁＝0.1,ω₁＝1，ω₂＝0。

在一个实施例中，所述步骤4中，P₁由如下式(2)的最优化问题所定义：

式(2)中，

为正则化项，

为数据保真项，υ表示输入图像，

表示输入图像水平方向的梯度，

表示输入图像竖直方向的梯度，||·||₀表示l₀范数，||·||₁表示l₁范数，||·||₂表示l₂范数，T表示转置。λ₂、ω₃、ω₄为调节参数，其中：λ₂用于控制正则化项和数据保真项权重，具体取值也是通过实验得到。原则上，要求ω₃大于十倍的ω₄，具体数值通过实验确定，如λ₂＝0.1、ω₃＝1、ω₄＝0。

在一个实施例中，所述步骤5中，P2由如下式(3)的最优化问题所定义：

式(3)中，argmin_υλ₃||Φ(C(υ))||₁为正则化项，

为数据保真项，Φ(·)表示惩罚函数作用于图像中的每一点处的平均曲率，||·||₁表示l₁范数，||·||₂表示l₂范数，λ₃为调节参数，λ₃用于控制正则化项和数据保真项权重，具体取值也是通过实验得到，如λ₃＝1E-4。C(·)表示计算图像每一点的平均曲率，对于第j点，其平均曲率为式(4)：

通过将υ，u^(k+2/3)视作三维曲面υ(x，y)，u^(k+2/3)(x，y)的采样，进而求解如下偏微分方程得到：

其中，x和y均为空间变量，t表示时间变量，υ(0；x，y)表示t＝0时的υ的取值，υ_t表示v(t，x，y)关于t的导函数，

表示梯度，

·表示散度，λ₄为调节参数，λ₄用于控制正则化项和数据保真项权重，具体取值也是通过实验得到,如λ₄＝1E-3。Φ′表示Φ的导函数。该偏微分方程采用现有技术中的经典数值方法求解。

所述步骤2中，R_G选择迭代类重建算法或解析类重建算法获得。其中，所述迭代类重建算法为ART、SART、EM；所述解析类重建算法为FBP、BPF。

本发明还提供一种基于曲率约束的X射线有限角CT图像重建装置，其包括：数据输入模块、初始化模块、迭代处理估计模块、梯度稀疏性约束模块、边界光滑性约束模块和判断模块，其中：

数据输入模块用于输入已知数据集，分别为：有限角CT扫描获得的数据集p以及CT扫描几何参数集G。其中，CT扫描几何参数集G包括射线源到探测器中心的距离(SDD)、射线源到转台中心的距离(SOD)、探测器单元个数、探测器单元尺寸、扫描角度数和角度采样间隔。

初始化模块用于初始估计图像u⁽⁰⁾，确定迭代终止阈值ε以及迭代次数上限N。一般将其所有像素值设置为0。

迭代处理估计模块用于迭代处理估计图像u^(k)，k＝0，1，…，N，利用有限角CT扫描数据集p，通过与扫描几何参数集G相关的图像重建算子R_G来更新估计图像u^(k)，得到u^(k+1/3)＝R_G(u^(k)，p)，其中，表示扫描几何参数集G相关的图像重建算子，u^(k+1/3)表示输入图像。R_G可以选择现有的迭代类重建算法或解析类重建算法获得。其中，所述迭代类重建算法为ART、SART、EM；所述解析类重建算法为FBP、BPF。

梯度稀疏性约束模块用于利用梯度的稀疏性约束图像u^(k+1/3)，得到的图像u^(k+2/3)＝P₁(u^(k+1/3))，其中，P₁表示约束图像梯度稀疏性的算子，u^(k+2/3)表示进行约束图像梯度稀疏性处理后得到的图像。其中，图像通过求梯度得到梯度域图像(其包括灰度图像的X方向的梯度，Y方向的梯度和图像像素值)，梯度的稀疏性指的则是梯度域图像中非零元素的多少。通常，梯度域图像中非零元素越少，梯度的稀疏性越好。本实施例利用图像梯度的稀疏性处理图像，一方面可以简化模型，避免过拟合，另一方面可以特征选择，选择参数非零的特征即可，而且使整个模型获得更好的可解释性，只有几个非零参数对结果有影响，更好解释。

边界光滑性约束模块用于利用边界的光滑性约束图像u^(k+2/3)，得到的图像u^(k+1)＝P₂(u^(k+2/3))，完成一次迭代，其中，P₂表示约束图像边界光滑性的算子，u^(k+1)表示进行约束图像边界的光滑性处理后得到的图像。

判断模块用于判断相邻两次迭代图像间的差别是否小于给定阈值，即||u^(k+1)-u^(k)||≤ε或是否达到迭代次数上限N，若不满足，转至步骤3，若满足，则终止迭代，基于曲率约束的X射线有限角CT图像重建方法结束。

本实施例提供的装置利用图像梯度的稀疏性以及边界的光滑性，通过迭代的方式，逐步恢复图像的边界信息，能够有效抑制图像的边界模糊以及阶梯效应，提高重建质量，进而克服了现有的有限角CT图像重建算法中边界模糊或存在阶梯效应的问题。

在一个实施例中，所述梯度稀疏性约束模块中的P₁由如下式(1)的最优化问题所定义：

式(1)中，

为正则化项，

为数据保真项，υ表示输入图像，

表示输入图像水平方向的梯度，

表示输入图像竖直方向的梯度，||·||₀表示l₀范数，||·||₁表示l₁范数，||·||₂表示l₂范数，T表示转置。λ₁、ω₁、ω₂为调节参数，其中：λ₁用于控制正则化项和数据保真项权重，具体取值也是通过实验得到。原则上，要求ω₁大于十倍的ω₂，具体数值通过实验确定，如λ₁＝0.1,ω₁＝1，ω₂＝0。

在一个实施例中，所述梯度稀疏性约束模块中的P₁由如下式(2)的最优化问题所定义：

式(2)中，

为正则化项，

为数据保真项，υ表示输入图像，

表示输入图像水平方向的梯度，

表示输入图像竖直方向的梯度，||·||₀表示l₀范数，||·||₁表示l₁范数，||·||₂表示l₂范数，T表示转置。λ₂、ω₃、ω₄为调节参数，其中：λ₂用于控制正则化项和数据保真项权重，具体取值也是通过实验得到。原则上，要求ω₃大于十倍的ω₄，具体数值通过实验确定，如λ₂＝0.1、ω₃＝1、ω₄＝0。

在一个实施例中，所述边界光滑性约束模块中的P₂由如下式(3)的最优化问题所定义：

式(3)中，argmin_υλ₃||Φ(C(υ))||₁为正则化项，

为数据保真项，Φ(·)表示惩罚函数作用于图像中的每一点处的平均曲率，||·||₁表示l₁范数，||·||₂表示l₂范数，λ₃为调节参数，λ₃用于控制正则化项和数据保真项权重，具体取值也是通过实验得到，如λ₃＝1E-4。C(·)表示计算图像每一点的平均曲率，对于第j点，其平均曲率为式(4)：

通过将υ，u^(k+2/3)视作三维曲面υ(x，y)，u^(k+2/3)(x，y)的采样，进而求解如下偏微分方程得到：

其中，x和y均为空间变量，t表示时间变量，υ(0；x，y)表示t＝0时的υ的取值，υ_t表示v(t，x，y)关于t的导函数，

表示梯度，

·表示散度，λ₄为调节参数，λ₄用于控制正则化项和数据保真项权重，具体取值也是通过实验得到,如λ₄＝1E-3。Φ′表示Φ的导函数。该偏微分方程采用现有技术中的经典数值方法求解。

上述各实施例中，l₀范数表示参数中非零元素的个数。l₁范数表示各参数绝对值的和，可以防止过拟合，l₁范数不会使得参数等于0而是接近于0，有参数平滑的作用，优化后的参数向量往往比较稀疏。l₂范数表示各个参数的平方的和的开方值，也就是欧氏距离，l₂范数越小，可以使得参数的每个元素都很小，接近于0，正则化项处处可导。

下面为了更好地体现本发明提供的有限角CT图像重建方法在重建效果方面的优势，下面结合一具体实施例将本发明所述的算法与已存在的典型算法SART和AEDS进行比较。

本实施例数据采集自工业CT系统，扫描样品为多个不同形状与材质的模体，如图2所示。实验电压与电流分别为140kV以及100uA，在120度范围内每隔0.5度采集一次投影数据。

分别采用SART、AEDS和本发明方法对扫描数据重建图像，结果分别如图3、图4、图5所示。其中：图3是SART算法的重建结果；图4是AEDS算法的重建结果；图5是本发明方法的重建结果。可以看出，SART算法的结果中，多个方向的边界存在模糊；AEDS算法的结果中，正方形与水平方向夹角较小的边界附近存在较严重的阶梯效应；而本发明的结果中，各个方向的边界均得到了清晰的重建。

最后需要指出的是：以上实施例仅用以说明本发明的技术方案，而非对其限制。本领域的普通技术人员应当理解：可以对前述各实施例所记载的技术方案进行修改，或者对其中部分技术特征进行等同替换；这些修改或者替换，并不使相应技术方案的本质脱离本发明各实施例技术方案的精神和范围。

Claims (7)

1.一种基于曲率约束的X射线有限角CT图像重建方法，其特征在于，包括：

步骤1，输入已知数据集：有限角CT扫描获得的数据集p，CT扫描几何参数集G；

步骤2，初始化：初始估计图像u⁽⁰⁾，确定迭代终止阈值ε以及迭代次数上限N；

步骤3，迭代处理估计图像u^(k)，k＝0，1，…，N，利用有限角CT扫描数据集p，通过与扫描几何参数集G相关的图像重建算子R_G来更新估计图像u^(k)，得到u^(k+1/3)＝R_G(u^(k)，p)，其中，R_G表示扫描几何参数集G相关的图像重建算子；

步骤4，利用梯度的稀疏性约束图像u^(k+1/3)，得到的u^(k+2/3)＝P₁(u^(k+1/3))，其中，P₁表示约束图像梯度稀疏性的算子；

步骤5，利用边界的光滑性约束图像u^(k+2/3)，得到的图像u^(k+1)＝P₂(u^(k+2/3))，完成一次迭代，其中，P₂表示约束图像边界光滑性的算子；

步骤6，判断相邻两次迭代图像间的差别是否小于给定阈值，即||u^(k+1)-u^(k)||≤ε或是否达到迭代次数上限N，若不满足，转至步骤3，若满足，则终止迭代；

所述步骤5中，P₂由如下式(3)的最优化问题所定义：

式(3)中，argmin_υλ₃||Φ(C(υ))||₁为正则化项，

为数据保真项，Φ(·)表示惩罚函数作用于图像中的每一点处的平均曲率，‖·‖₁表示l₁范数，‖·‖₂表示l₂范数，λ₃为调节参数，υ表示输入图像，C(·)表示计算图像每一点的平均曲率，对于第j点，其平均曲率为式(4)：

该优化问题的求解通过将υ视作三维曲面υ(x，y)、将u^(k+2/3)视作三维曲面u^(k+2/3)(x，y)进行采样，通过求解如下偏微分方程得到：

其中，x和y均为空间变量，t表示时间变量，υ(0；x，y)表示t＝0时的υ的取值，υ_t(t，x，y)表示v(t,x,y)关于t的导函数，

表示梯度，

表示散度，λ₄为调节参数，Φ表示Φ的导函数。

2.如权利要求1所述的基于曲率约束的X射线有限角CT图像重建方法，其特征在于，所述步骤4中，P₁由如下式(1)的最优化问题所定义：

式(1)中，

为正则化项，

为数据保真项，υ表示输入图像，

表示输入图像水平方向的梯度，

表示输入图像竖直方向的梯度，‖·‖₀表示l₀范数，‖·‖₁表示l₁范数，‖·‖₂表示l₂范数，T表示转置，λ₁、ω₁、ω₂为调节参数。

3.如权利要求1所述的基于曲率约束的X射线有限角CT图像重建方法，其特征在于，所述步骤4中，P₁由如下式(2)的最优化问题所定义：

式(2)中，

为正则化项，

为数据保真项，υ表示输入图像，

表示输入图像水平方向的梯度，

表示输入图像竖直方向的梯度，‖·‖₀表示l₀范数，‖·‖₁表示l₁范数，‖·‖₂表示l₂范数，T表示转置，λ₂、ω₃、ω₄为调节参数。

4.如权利要求1所述的基于曲率约束的X射线有限角CT图像重建方法，其特征在于，所述步骤2中，R_G选择迭代类重建算法或解析类重建算法获得。

5.一种基于曲率约束的X射线有限角CT图像重建装置，其特征在于，包括：

数据输入模块，其用于输入已知数据集，分别为：有限角CT扫描获得的数据集p以及CT扫描几何参数集G；

初始化模块，其用于初始估计图像u⁽⁰⁾，确定迭代终止阈值ε以及迭代次数上限N；

迭代处理估计模块，其用于迭代处理估计图像u^(k)，k＝0，1，…，N，利用有限角CT扫描数据集p，通过与扫描几何参数集G相关的图像重建算子R_G来更新估计图像u^(k)，得到u^(k+1/3)＝R_G(u^(k)，p)，其中，R_G表示扫描几何参数集G相关的图像重建算子；

梯度稀疏性约束模块，其用于利用梯度的稀疏性约束图像u^(k+1/3)，得到的u^(k+2/3)＝P₁(u^(k+1/3))，其中，P₁表示约束图像梯度稀疏性的算子；

边界光滑性约束模块，其用于利用边界的光滑性约束图像u^(k+2/3)，得到的图像u^(k+1)＝P₂(u^(k+2/3))，完成一次迭代，其中，P₂表示约束图像边界光滑性的算子；

判断模块，其用于判断相邻两次迭代图像间的差别是否小于给定阈值，即||u^(k+1)-u^(k)||≤ε或是否达到迭代次数上限N，若不满足，转至步骤3，若满足，则终止迭代；

所述边界光滑性约束模块中的P₂由如下式(3)的最优化问题所定义：

式(3)中，argmin_υλ₃||Φ(C(υ))||₁为正则化项，

为数据保真项，Φ(·)表示惩罚函数作用于图像中的每一点处的平均曲率，‖·‖₁表示l₁范数，‖·‖₂表示l₂范数，λ₃为调节参数，υ表示输入图像，C(·)表示计算图像每一点的平均曲率，对于第j点，其平均曲率为式(4)：

该优化问题的求解通过将υ视作三维曲面υ(x，u)、将u^(k+2/3)视作三维曲面u^(k+2/3)(x，y)进行采样，通过求解如下偏微分方程得到：

其中，x和y均为空间变量，t表示时间变量，υ(0；x，y)表示t＝0时的υ的取值，υ_t(t，x，y)表示v(t,x,y)关于t的导函数，

表示梯度，

表示散度，λ₄为调节参数，Φ′表示Φ的一阶导数。

6.如权利要求5所述的基于曲率约束的X射线有限角CT图像重建装置，其特征在于，所述梯度稀疏性约束模块中的P₁由如下式(1)的最优化问题所定义：

式(1)中，

为正则化项，

为数据保真项，υ表示输入图像，

表示输入图像水平方向的梯度，

表示输入图像竖直方向的梯度，‖·‖₀表示l₀范数，‖·‖₁表示l₁范数，‖·‖₂表示l₂范数，T表示转置，λ₁、ω₁、ω₂为调节参数。

7.如权利要求5所述的基于曲率约束的X射线有限角CT图像重建装置，其特征在于，所述梯度稀疏性约束模块中的P₁由如下式(2)的最优化问题所定义：

式(2)中，

为正则化项，

为数据保真项，υ表示输入图像，

表示输入图像水平方向的梯度，

表示输入图像竖直方向的梯度，‖·‖₀表示l₀范数，‖·‖₁表示l₁范数，‖·‖₂表示l₂范数，T表示转置，λ₂、ω₃、ω₄为调节参数。

Images