CN112070704B - 一种基于紧小波框架的双正则化有限角ct图像重建方法 - Google Patents

Abstract

本发明涉及一种基于紧小波框架的双正则化有限角CT图像重建方法，属于图像处理领域。该方法包括以下步骤：S1：采集投影数据；S2：建立双正则化有限角CT重建模型；S3：有限角CT迭代重建；S4：输出重建图像；当步骤S3中的迭代重建算法收敛时，输出重建图像。本发明公开的基于紧小波框架的双正则化有限角CT图像重建方法包含对高频部分进行硬阈值处理，和对低频部分进行图像TV最小化光滑处理，经过这种双正则化机制处理后，不仅CT图像的有限角伪影和噪声能够有效地抑制，而且能够有效保护边界，从而提高CT重建图像的质量。

Description

一种基于紧小波框架的双正则化有限角CT图像重建方法

技术领域

本发明属于图像处理领域，涉及一种基于紧小波框架的双正则化有限角CT图像重建方法。

背景技术

受扫描环境、被扫描目标自身结构等因素限制，有时候只能够在一定的角度范围内对待检目标进行扫描，此时采集投影数据是不完备的。例如：在役管道成像、大尺寸物体检测、C型臂CT(Computed Tomography)等。这种扫描情形称为有限角CT扫描。针对有限角CT扫描的重建问题，采用传统的滤波反投影重建算法将使得重建图像出现许多伪影，导致有些重要结构信息丢失或被掩盖，严重地影响无损检测的精度或医生对于病症的诊断。因此针对有限角CT图像重建问题，如何重建出符合无损检测标准或医生诊断要求的高质量CT图像具有较大的实际意义。

现有技术，传统的代数重建算法ART(Algebraic Reconstruction Technique)、SART(Simultaneous Algebraic Reconstruction Technique)和SIRT(SimultaneousIterative Reconstruction Technique)等算法，使得重建出来的图像出现明显的伪影。基于TV(Total Variation)的图像重建算法处理稀疏角CT重建有较好的效果，但是对较小的有限角的扫描，不能有效抑制伪影和较好地保护边界结构。余维引进图像梯度L0正则化，一定程度上有效地抑制伪影并避免过度光滑。王成祥等人将图像小波变换的L0拟范数作为正则化，利用小波变换的多尺度、多分辨率特点来处理有限角CT重建问题，该方法能够一定程度地抑制伪影并保护边界。

公开号为CN 107978005A专利申请公开“一种基于保边界扩散和平滑的有限角CT图像重建算法”。该方法首先利用投影数据进行重建，然后对重建出来的图像x，y方向分别进行梯度L0保边界扩散修正，最后重复迭代修正直到达到满足停机标准。虽然上述专利申请所述方法能够保边界平滑、消除扩散可能引入的线状伪影。但是仍然存在如下的缺陷：上述专利申请所述方法仅考虑图像x，y轴进行梯度L0保边界扩散修正，然而没有考虑到有限角CT图像的伪影具有多方向性的特点。图像的梯度变换只有图像的高频部分高频信息，带有限角伪影的图像不仅高频部分有伪影，低频部分有出现失真，上述专利申请所述方法没有考虑到图像多尺度和多分辨率的特性。

公开号为CN 110717959A专利申请公开“基于曲率约束的x射线有限角CT图像重建方法和装置”。该方法首先利用投影数据进行重建，其次对重建出的图像进行图像梯度L0正则化稀疏约束，然后对稀疏约束的结果进行曲率约束，最后重复迭代直到满足停机标准。虽然上述专利申请所述方法能够克服现有的有限角CT重建算法中边界模糊或存在阶梯效应问题。但是仍然存在如下的缺陷：(1)上述专利申请所述方法仅考虑图像的梯度变换，只有图像的高频部分高频信息；(2)上述专利申请所述方法考虑曲率约束，只考虑到高频轮廓的变化。带有限角伪影的图像不仅高频部分有伪影，低频部分有出现失真，上述专利申请所述方法没有考虑到图像多尺度和多分辨率的特性。

公开号为CN 109697691A专利申请公开“基一种基于L0范数和奇异值阈值分解的双正则项优化的有限角投影重建方法”。该方法首先利用SART算法和误差进行修正进行重建，其次对重建的图像进行图像梯度L0正则化稀疏约束，然后对稀疏约束后的图像进行奇异值分解和软阈值进一步处理。最后重复迭代直到达到满足停机标准。虽然上述专利申请所述方法能够恢复CT图像轮廓，减少有限角伪影。但是仍然存在如下的缺陷：上述专利申请所述方法仅考虑x，y轴方向伪影，然而没有考虑到有限角CT图像的伪影具有多方向性的特点。图像的梯度变换只有图像的高频部分高频信息，带有限角伪影的图像不仅高频部分有伪影，低频部分有出现失真，上述专利申请所述方法没有考虑到图像多尺度和多分辨率的特性。

目前存在的大多数有限角CT优化重建方法要么只考虑对图像高频部分进行约束，要么没有考虑到有限角CT重建图像的伪影在多尺度、多分辨率分解下，不仅存在高频部分，还存在低频部分，要么没有考虑对低频部分的伪影进行有效的处理。本发明考虑图像在紧小波框架变换下的稀疏性，利用B样条紧小波框架对图像进行多尺度、多分辨率分解，采用双正则化约束机制，对高频部分进行L0稀疏约束，对低频部分采用TV正则化进行光滑和稀疏约束，从而提高CT重建图像的质量。

发明内容

有鉴于此，本发明的目的在于提供一种基于紧小波框架的双正则化有限角CT图像重建方法，能够有效地抑制重建图像的伪影并保护图像边界，从而提高CT重建图像的质量。

为达到上述目的，本发明提供如下技术方案：

一种基于紧小波框架的双正则化有限角CT图像重建方法，该方法包括以下步骤：

S1：采集投影数据；

S2：建立双正则化有限角CT重建模型；

S3：有限角CT迭代重建；

S4：输出重建图像；当步骤S3中的迭代重建算法收敛时，输出重建图像。

可选的，所述S1具体为：

在控制与图像处理***的控制下，首先将射线源和线阵探测器绕待检的中心沿着轨道旋转有限的角度来获得不完备的投影数据，然后传送到控制与图像处理***中存储。

可选的，所述S2具体为：

采用离散模型进行重建时，首先需要将所有(x,y)对应的重建像素f(x,y)按照y的维度将其转变成一个长长的列向量f，列向量f的维数为N×1，其中N＝n₁×n₂，n₁为f(x,y)在x方向的维数，n₂为f(x,y)在y方向的维数；

然后，将所有投影视角指标s对应线阵探测器位置a的投影数据g^δ(a,s)按照s的维度将其转变成一个长长的列向量g^δ，列向量g^δ的维数为M×1，其中M＝m₁×m₂，m₁为g^δ(a,s)在a方向的维数，m₂为g^δ(a,s)在s方向的维数，即总的投影视角数；

采用B样条紧小波框架变换将重建图像分解成低频部分和高频部分，使得图像进行多尺度、多分辨率分解；

为抑制高频中的噪声和伪影，通过对B样条紧小波框架变换高频部分进行L0稀疏正则化约束，为了使得重建图像变得光滑和抑制低频部分的伪影，对B样条紧小波框架变换的低频部分进行TV正则化约束；

建立的模型如下：

其中A∈R^M×N是有限角CT***矩阵，f∈R^N×1是待重建图像，g^δ∈R^M×1是有限角CT投影数据，Ω是凸集(Ω:＝{f|f≥0})，||x||_D＝<Dx,x>；D是一个对角矩阵，其对角元素为且对所有i′＝1,2,...,M，/>λ_i是正则化参数，W是B样条紧小波框架；Ω₁是高频子带的指标集，Ω₂是低频子带的指标集，Ω₁∪Ω₂表示所有B样条紧小波框架变换子带的指标集；||β||₀是统计β的非0元素个数，

在B样条紧小波框架变换时，将f中(x,y)的y维度排成1个2维矩阵f(x,y)，然后对该2维矩阵做B样条紧小波框架变换，当做完B样条紧小波框架变换反变换后，重新将f(x,y)按照y的维度将其转变成一个长长的列向量f。

可选的，所述S3具体为：

根据建立的模型(1)，采用交替方向迭代ADMM的方法求解模型(1)；

其中步骤S3有限角CT迭代重建具体过程为：

首先，将模型(1)通过交替方向迭代的方法转换如下迭代格式：

迭代格式(2)中k表示迭代次数，ρ是松弛参数，α是辅助变量，v是对偶变量，t是ADMM***变量时引入的参数；

其次，为避免关于第一个变量f的子问题中求***矩阵A的逆或者采用迭代法求解子问题的不足，将迭代格式(2)转换为临近交替线性化的ADMM迭代格式，具体如下：

迭代格式(3)中μ是临近线性化引入的松弛参数；(3)中为ART迭代更新格式；通过临近交替线性化巧妙地将经典的ART迭代算法融入其中；

然后，为求解子问题α，将α分解成高频和低频分别求解，迭代格式变成如下等价形式：

最后，求出迭代格式(4)的子问题最优解，迭代格式如下所示：

其中W^T表示B样条紧小波框架反变换，表示图像TV最小化光滑算法/>图像TV最小化光滑算法迭代格式如下：

其中τ为SplitBregman算法引入的参数，Gⁿ的分量形式为：

n表示TV最小化光滑算法的迭代次数。

可选的，所述迭代格式中，包括以下步骤：

S31：ART迭代重建和B样条紧小波框架反变换线性组合并确保组合后图像非负，得到初步重建结果，即公式(5)的第1个方程；

S32：在B样条紧小波框架变换域对高频部分进行硬阈值处理，抑制高频部分的噪声和方向性伪影，即公式(5)的第2个方程；

S33：在B样条紧小波框架变换域对低频部分进行TV最小化光滑处理，抑制低频部分的伪影和保持图像的平滑性，即公式(5)的第3个方程；

S34：对偶变量v更新，即公式(5)的第4个方程；当达到一定的迭代次数时停止迭代，否则重复步骤S31-S34。

本发明的有益效果在于：本发明公开了基于紧小波框架的双正则化有限角CT图像重建方法，涉及到一种有限角CT重建技术。本发明考虑到有限角CT重建图像的伪影在多尺度、多分辨率分解下，不仅存在高频部分，还存在低频部分，采用了B样条紧小波框架变换将重建图像分解成低频部分和高频部分，使得伪影被多尺度分解。为了抑制高频中的噪声和方向性伪影，通过对B样条紧小波框架变换的高频部分进行L0稀疏正则化约束，为了使得重建图像变得光滑和抑制低频部分的伪影，对B样条紧小波框架变换的低频部分进行图像TV正则化约束。本发明公开的基于紧小波框架的双正则化有限角CT图像重建方法包含对高频部分进行硬阈值处理，和对低频部分进行图像TV最小化光滑处理，经过这种双正则化机制处理后，不仅CT图像的有限角伪影和噪声能够有效地抑制，而且能够有效保护边界，从而提高CT重建图像的质量。

本发明的其他优点、目标和特征在某种程度上将在随后的说明书中进行阐述，并且在某种程度上，基于对下文的考察研究对本领域技术人员而言将是显而易见的，或者可以从本发明的实践中得到教导。本发明的目标和其他优点可以通过下面的说明书来实现和获得。

附图说明

为了使本发明的目的、技术方案和优点更加清楚，下面将结合附图对本发明作优选的详细描述，其中：

图1为本发明的待检目标的扫描结构示意图；

图2为本发明的有限角CT重建算法几何结构示意图；

图3为基于紧小波框架的双正则化有限角CT图像重建的流程图。

附图标记：1-射线源，2-线阵探测器，3-待检目标，4-轨道，5-图像处理***。

具体实施方式

以下通过特定的具体实例说明本发明的实施方式，本领域技术人员可由本说明书所揭露的内容轻易地了解本发明的其他优点与功效。本发明还可以通过另外不同的具体实施方式加以实施或应用，本说明书中的各项细节也可以基于不同观点与应用，在没有背离本发明的精神下进行各种修饰或改变。需要说明的是，以下实施例中所提供的图示仅以示意方式说明本发明的基本构想，在不冲突的情况下，以下实施例及实施例中的特征可以相互组合。

其中，附图仅用于示例性说明，表示的仅是示意图，而非实物图，不能理解为对本发明的限制；为了更好地说明本发明的实施例，附图某些部件会有省略、放大或缩小，并不代表实际产品的尺寸；对本领域技术人员来说，附图中某些公知结构及其说明可能省略是可以理解的。

本发明实施例的附图中相同或相似的标号对应相同或相似的部件；在本发明的描述中，需要理解的是，若有术语“上”、“下”、“左”、“右”、“前”、“后”等指示的方位或位置关系为基于附图所示的方位或位置关系，仅是为了便于描述本发明和简化描述，而不是指示或暗示所指的装置或元件必须具有特定的方位、以特定的方位构造和操作，因此附图中描述位置关系的用语仅用于示例性说明，不能理解为对本发明的限制，对于本领域的普通技术人员而言，可以根据具体情况理解上述术语的具体含义。

图1为本发明的待检目标的扫描结构示意图：采集数据前，将射线源1和线阵探测器2分别放置于待检目标3的两侧，使得射线源1产生的扇束射线束能够覆盖待检目标3。采集数据时，射线源1沿着轨道4旋转，与此同时线阵探测器2同步旋转，由于扫描场景限制，导致射线源1和线阵探测器2只能在有限角度内旋转采集数据。采集好数据传输至控制及图像处理***5进行重建。

图2为本发明的有限角CT重建算法几何结构示意图：以射线源1到待检目标3中心轴的垂足为坐标原点O建立空间右手直角坐标系O-xy，y轴为原点与射线源1的连线，且从原点指向射线源1为正方向，x轴为垂直于y轴的轴，且向右为正方向。以坐标原点O为旋转中心，(x,y)表示被重建像素J的坐标，射线源1位于S，SL表示扇束射线的中心射线，SF表示经过被重建点J的一条射线。a轴表示线阵探测器中的探测单元位置且正方向与空间坐标系的x轴一致，(a,s)表示投影视角指标s对应线阵探测器位置a。

图3为基于紧小波框架的双正则化有限角CT图像重建的流程图：基于紧小波框架的双正则化有限角CT图像重建方法包括以下步骤：

S1.采集投影数据：在控制与图像处理***5的控制下，首先将射线源1和线阵探测器2绕待检的中心沿着轨道4旋转有限的角度来获得不完备的投影数据，然后传送到控制与图像处理***5中存储；

S2.建立双正则化有限角CT重建模型：采用离散模型进行重建时，首先需要将如图2所示的所有(x,y)对应的重建像素f(x,y)按照y的维度将其转变成一个长长的列向量f，列向量f的维数为N×1，其中N＝n₁×n₂，n₁为f(x,y)在x方向的维数，n₂为f(x,y)在y方向的维数。然后，将如图2所示的所有投影视角指标s对应线阵探测器位置a的投影数据g^δ(a,s)按照s的维度将其转变成一个长长的列向量g^δ，列向量g^δ的维数为M×1，其中M＝m₁×m₂，m₁为g^δ(a,s)在a方向的维数，m₂为g^δ(a,s)在s方向的维数，也就是总的投影视角数。按照本发明考虑采用B样条紧小波框架变换将重建图像分解成低频部分和高频部分，使得图像进行多尺度、多分辨率分解。为了抑制高频中的噪声和伪影，通过对B样条紧小波框架变换高频部分进行L0稀疏正则化约束，为了使得重建图像变得光滑和抑制低频部分的伪影，对B样条紧小波框架变换的低频部分进行TV正则化约束。本发明建立的模型如下：

其中A∈R^M×N是有限角CT***矩阵，f∈R^N×1是待重建图像，g^δ∈R^M×1是有限角CT投影数据，Ω是凸集(Ω:＝{f|f≥0})，||x||_D＝<Dx,x>。D是一个对角矩阵，其对角元素为且对所有i′＝1,2,...,M，/>λ_i是正则化参数，W是B样条紧小波框架。Ω₁是高频子带的指标集，Ω₂是低频子带的指标集，Ω₁∪Ω₂表示所有B样条紧小波框架变换子带的指标集。||β||₀是统计β的非0元素个数, 当我们做B样条紧小波框架变换时，我们将f按照图2中(x,y)的y维度排成1个2维矩阵，然后对该2维矩阵做B样条紧小波框架变换，当做完B样条紧小波框架变换反变换后，我们重新将f(x,y)按照y的维度将其转变成一个长长的列向量f。

S3.有限角CT迭代重建：根据建立的模型(1)，采用交替方向迭代(ADMM,alternating direction method of multipliers)的方法求解模型(1)；

其中步骤S3有限角CT迭代重建具体过程为：

首先,将模型(1)通过交替方向迭代的方法转换如下迭代格式：

迭代格式(2)中k表示迭代次数，ρ是松弛参数，α是辅助变量，v是对偶变量，t是ADMM***变量时引入的参数。

其次，为了避免关于第一个变量f的子问题中求***矩阵A的逆或者采用迭代法求解子问题的不足，本发明嵌入临近交替线性化的思想，将迭代格式(2)转换为临近交替线性化的ADMM迭代格式，具体如下：

迭代格式(3)中μ是临近线性化引入的松弛参数。(3)中为ART迭代更新格式。通过临近交替线性化巧妙地将经典的ART迭代算法融入其中。

然后，为了求解子问题α，我们将α分解成高频和低频分别求解，迭代格式变成如下等价形式：

最后，求出迭代格式(4)的子问题最优解，迭代格式如下所示：

其中(W^T表示B样条紧小波框架反变换)，表示图像TV最小化光滑算法/>图像TV最小化光滑算法迭代格式如下：

其中τ为SplitBregman算法引入的参数，Gⁿ的分量形式为：

n表示TV最小化光滑算法的迭代次数。

根据公式(5),步骤S4有限角CT迭代重建包含如下4个子步骤：S31.ART迭代重建和B样条紧小波框架反变换线性组合并确保组合后图像非负，得到初步重建结果(公式(5)的第1个方程)；S32.在B样条紧小波框架变换域对高频部分进行硬阈值处理，抑制高频部分的噪声和方向性伪影(公式(5)的第2个方程)；S33.在B样条紧小波框架变换域对低频部分进行TV最小化光滑处理，抑制低频部分的伪影和保持图像的平滑性(公式(5)的第3个方程)；S34.对偶变量v更新(公式(5)的第4个方程)。当达到一定的迭代次数时停止迭代，否则重复步骤S31-S34。

S4.输出重建图像。当步骤S3中的迭代重建算法收敛时，输出重建图像。

最后说明的是，以上实施例仅用以说明本发明的技术方案而非限制，尽管参照较佳实施例对本发明进行了详细说明，本领域的普通技术人员应当理解，可以对本发明的技术方案进行修改或者等同替换，而不脱离本技术方案的宗旨和范围，其均应涵盖在本发明的权利要求范围当中。

Claims (2)

1.一种基于紧小波框架的双正则化有限角CT图像重建方法，其特征在于：该方法包括以下步骤：

S1：采集投影数据；

S2：建立双正则化有限角CT重建模型；

S3：有限角CT迭代重建；

S4：输出重建图像；当步骤S3中的迭代重建算法收敛时，输出重建图像；

所述S1具体为：

在控制与图像处理***的控制下，首先将射线源和线阵探测器绕待检的中心沿着轨道旋转有限的角度来获得不完备的投影数据，然后传送到控制与图像处理***中存储；

所述S2具体为：

采用离散模型进行重建时，首先需要将所有(x,y)对应的重建像素f(x,y)按照y的维度将其转变成一个列向量f，列向量f的维数为N×1，其中N＝n₁×n₂，n₁为f(x,y)在x方向的维数，n₂为f(x,y)在y方向的维数；

然后，将所有投影视角指标s对应线阵探测器位置a的投影数据g^δ(a,s)按照s的维度将其转变成一个列向量g^δ，列向量g^δ的维数为M×1，其中M＝m₁×m₂，m₁为g^δ(a,s)在a方向的维数，m₂为g^δ(a,s)在s方向的维数，即总的投影视角数；

采用B样条紧小波框架变换将重建图像分解成低频部分和高频部分，使得图像进行多尺度、多分辨率分解；

为抑制高频中的噪声和伪影，通过对B样条紧小波框架变换高频部分进行L0稀疏正则化约束，为了使得重建图像变得光滑和抑制低频部分的伪影，对B样条紧小波框架变换的低频部分进行TV正则化约束；

建立的模型如下：

其中A∈R^M×N是有限角CT***矩阵，f∈R^N×1是待重建图像，g^δ∈R^M×1是有限角CT投影数据，Ω是凸集，Ω:＝{f|f≥0}，||x||_D＝<Dx,x>；D是一个对角矩阵，其对角元素为且对所有i′＝1,2,...,M，/>λ_i是正则化参数，W是B样条紧小波框架；Ω₁是高频子带的指标集，Ω₂是低频子带的指标集，Ω₁∪Ω₂表示所有B样条紧小波框架变换子带的指标集；||β||₀是统计β的非0元素个数，(▽_xβ)_i′,j′＝β_i′,j′-β_i′-1,j′,(▽_yβ)_i′，j′＝β_i′,j′-β_i′,j′-1；

在B样条紧小波框架变换时，将f中(x,y)的y维度排成1个2维矩阵f(x,y)，然后对该2维矩阵做B样条紧小波框架变换，当做完B样条紧小波框架变换反变换后，重新将f(x,y)按照y的维度将其转变成一个列向量f；

所述S3具体为：

根据建立的模型(1)，采用交替方向迭代ADMM的方法求解模型(1)；

其中步骤S3有限角CT迭代重建具体过程为：

首先，将模型(1)通过交替方向迭代的方法转换如下迭代格式：

迭代格式(2)中k表示迭代次数，α是辅助变量，v是对偶变量，t是ADMM***变量时引入的参数；

其次，为避免关于第一个变量f的子问题中求***矩阵A的逆或者采用迭代法求解子问题的不足，将迭代格式(2)转换为临近交替线性化的ADMM迭代格式，具体如下：

迭代格式(3)中γ，μ是临近线性化引入的松弛参数；(3)中/>为ART迭代更新格式；通过临近交替线性化巧妙地将经典的ART迭代算法融入其中；

然后，为求解子问题α，将α分解成高频和低频分别求解，迭代格式变成如下等价形式：

最后，求出迭代格式(4)的子问题最优解，迭代格式如下所示：

其中W^T表示B样条紧小波框架反变换， 表示图像TV最小化光滑算法/>图像TV最小化光滑算法迭代格式如下：

其中τ为Split Bregman算法引入的参数，Gⁿ的分量形式为：

n表示TV最小化光滑算法的迭代次数。

2.根据权利要求1所述的一种基于紧小波框架的双正则化有限角CT图像重建方法，其特征在于：所述迭代格式中，包括以下步骤：

S31：ART迭代重建和B样条紧小波框架反变换线性组合并确保组合后图像非负，得到初步重建结果，即公式(5)的第1个方程；

S32：在B样条紧小波框架变换域对高频部分进行硬阈值处理，抑制高频部分的噪声和方向性伪影，即公式(5)的第2个方程；

S33：在B样条紧小波框架变换域对低频部分进行TV最小化光滑处理，抑制低频部分的伪影和保持图像的平滑性，即公式(5)的第3个方程；

S34：对偶变量v更新，即公式(5)的第4个方程；当达到一定的迭代次数时停止迭代，否则重复步骤S31-S34。