CN103413338B - 一种基于广义变分最小化的少量投影ct图像重建方法 - Google Patents

Abstract

本发明涉及图像处理领域，尤其涉及一种基于广义变分最小化的少量投影CT图像重建方法。本发明针对现有技术中存在的问题，提供基于广义变分最小化的少量投影CT图像重建算法，以减少CT重建时所需的投影数目,从而减少CT对人体的辐射剂量或CT扫描时间,适用于不完全投影数据CT图像重建问题。本发明通过建立CT扫描成像模型的投影矩阵形式，进而建立用于图像重建的优化模型；然后对图像重建的优化模型中的图像进行优化求解，得到重建数据；再建立广义变分最小化的图像修正模型，对重建模型的重建数据进行修正，获得修正后的重建图像；最后判断修正后的重建图像是否满足迭代停止准则等步骤实现，本发明应用于CT图像重建设计领域。

Description

一种基于广义变分最小化的少量投影CT图像重建方法

技术领域

本发明涉及图像处理领域，尤其涉及一种基于广义变分最小化的少量投影CT图像重建方法。

背景技术

计算机层析技术(CT)是一种非常重要的现代无损检测技术，有着十分广泛的应用。如工业上利用CT技术来检测产品内部的缺陷及装配情况；医学上利用CT技术来检测人体头部、腹部等内部器官病变情况；石油、化工等领域中利用CT技术来检测多相流断面，获取管道截面分布的微观信息。图像重建是CT技术中的核心问题，它建立了投影数据与断层结构图像之间的内在关系。

CT技术在实际应用中由于受检测环境、检测时间和检测成本的限制，常常不能获得完备的投影数据，这样的CT图像重建问题称之为少量投影CT图像重建问题。另外，由少量投影数据重建CT图像在医学、工业检测等领域也都有实际的应用背景，如医学CT中需要通过减少投影数据，降低对人体的辐射，提高其安全性；工业CT中需要通过减少投影数据，减少扫描时间，提高其检测效率。因此，研究不完全投影情况下的少量投影CT图像重建方法，具有切实的研究意义和实际的应用价值，并且具有广泛的应用前景。

目前，CT图像重建常用的方法是解析法中的滤波反投影法(FBP)和级数展开法中的代数迭代重建技术“R.Gordon,R.Bender,G.T.Herman,AlgebraicReconstructionTechniques(ART)forthreedimensionalelectronicroscopyandX-rayphotography[J].J.Theor.Biol,1970,29：471-481.”。而由不完全的投影数据重建CT图像的问题在数学上通常被认为是欠定问题。直接用解析重建算法对不完全的投影数据进行重建，重建的效果将严重恶化，得不到较好的重建图像。ART算法将图像重建问题转化为求解线性方程组，当投影数据不完全时，可以把丢失的投影值看做是缺少了的方程，这在某种程度上降低了对数据完备性的要求。但是，当投影数据较少时，ART算法也无法获得较好的重建结果，图像质量依然很差。

由于不完全投影数据的CT图像重建问题具有很强的需求背景和广泛的应用前景，所以一直备受广大学者的关注。Sidky等人“EmilY.Sidky,Chien-MinKaoandXiaochuanPan.Accurateimagereconstructionfromfew-viewsandlimited-angledataindivergent-beamCT[J].JournalofX-RayScienceandTechnology2006,14(2)：119-139.”将图像全变分极小化与代数重建技术(ART)相结合，提出了一种少量投影数据的重建CT图像方法(ART-TV)，得到了较好的重建效果。但是，ART-TV算法采用ART算法来迭代重建图像，采用最速下降算法来最小化TV范数，效率较低，且在进一步减少投影数据时重建效果较差。针对不完投影数据的CT图像重建问题，为了能够进一步减少投影数目，如何建立优化模型、设计求解方法均需进一步研究。

发明内容

本发明所要解决的技术问题是：针对上述存在的问题，提供一种基于广义变分最小化的少量投影CT图像重建算法，保证获得高质量重建图像的前提下，以减少CT重建时所需的投影数目,从而减少CT对人体的辐射剂量或CT扫描时间,适用于不完全投影数据CT图像重建问题。

本发明采用的技术方案如下：

一种基于广义变分最小化的少量投影CT图像重建方法包括：

S1：根据CT扫描成像方式，建立CT扫描的成像模型的投影矩阵；

S2：根据CT扫描的成像模型的投影矩阵，建立用于图像重建的优化模型；

S3：采用无约束最小化的迭代算法对图像重建优化模型中的图像进行优化求解，得到重建模型的重建数据；

S4:建立广义变分最小化的图像修正模型，对S3中重建模型的重建数据进行修正，获得修正后的重建图像；判断修正后的重建图像是否满足迭代停止准则，若满足则停止重建图像，进行S5；否则，返回S3，直至满足迭代停止准则；

S5:所得图像作为结果输出重建结果。

所述S1中CT扫描成像方式包括平行束扫描、扇形束扫描或者锥形束扫描。

所述步骤S2具体步骤包括：

S21：根据公式(1)建立图像重建的正向投影模型，

y＝Au(1)；

其中，u表示待重建的断层图像；CT扫描时，通过射线源S的旋转，可以获得不同角度下的投影数据，A表示一个M×N_image的投影矩阵，其元素A_ij表示断层图像扫描线与像素网格相交的长度，设总的投影数目为M，断层图像的像素数为N_image，y表示由投影数据构成的一个M×1列向量，其元素y_i表示第i个投影数据，1≤i≤M；M表示总的投影数目，其大小由投影角度数目和每个投影角度下的投影数据数，即探测器单元数目决定，即M＝θ_n×D_n，θ_n表示投影角度的数量，D_n表示探测器的单元数目；u表示待重建的断层图像，是一个N_image×1的列向量，其元素u_j表示第j个像素的密度值，1≤j≤N_image，N是断层图像的总像素数，N＝W×H，W和H分别表示断层图像的宽度和高度；

S22：得到重建模型的优化重建的数学模型，即使待重建的断层图像的再投影Au与实际扫描获得的投影数据y的残差最小值minf(u)：

minf(u)＝min‖Au-y‖²(2)。

所述S3中无约束最小化的迭代算法包括牛顿算法，最速下降算法或共轭梯度算法。

所述S3中采用共轭梯度算法具体步骤：

S31：选取初值u₀，计算目标泛函f(u)梯度g₀，g₀＝Au₀-y；设置迭代停止时的允许误差0＜ε＜10^-3；

S32：计算搜索方向p₀，p₀＝-g₀，其中当迭代次数k不为1时，p_k由S37给出；

S33：设置迭代步长α_k，α_k∈4×10^-5～6×10^-5；

S34：更新重建图像u_k+1，u_k+1＝u_k+α_kp_k；

S35：更新目标泛函梯度g_k+1，g_k+1＝Au_k+1-p_k；

S36：计算梯度搜索方向β_k+1，

S37：更新搜索方向p_k+1，p_k+1＝-g_k+1+β_k+1p_k；

S38：如果‖g_k+1‖＜ε，停止，输出u＝u_k+1；否则令k＝0，转到S33；

其中，共轭梯度算法选择Fletcher-Reeves方法计算其搜索方向。

所述S4具体步骤包括：

S41：根据公式(3)计算GV，

$GV=\∫{\left(\sqrt{{u_x}^2+{u_y}^2+\mathrm{\δ}}\right)}^{s}dxdy---\left(3\right)$

并根据公式(3)得到广义变分最小化的图像修正模型

其中，GV表示图像的广义变分范数，u_x和u_x分别表示断层图像在x方向和y方向的梯度值，δ＝10^-6，s为广义调整参数，s的选取范围为0＜s≤1；

S42：通过设置迭代步长α_k∈1×10^-4～2×10^-4的共轭梯度算法求GV最小，对重建数据进行修正，并采用了迭代n步后重新开始的方法，其下降方向是利用上一次搜索方向和本次起始点的负梯度的线性组合生成新的共轭方向，形成新一轮的共轭梯度算法迭代，直至满足收敛条件获得修正后的重建图像。

所述S42具体步骤：

S441：将S3重建模型的重建数据设为初始值u₀，重新开始迭代次数设为n，设置迭代停止时的允许误差10^-8＜ε＜10^-3，n>0；

S442：计算目标函数负梯度，

S443：采用固定迭代步长，α_k∈1×10^-4～2×10^-4；

S444：更新修正值，u_k+1＝u_k+α_kp_k；

S445：判断是否满足迭代停止准则，即是否成立，如成立则停止迭代，输出u＝u_k+1；否则进行S443；

S446：若k＝n-1，令u₀＝u_k+1，k＝0，返回S442；否则进行S447；

S447：利用公式计算F-R搜索方向更新负梯度方向p^k+1＝-g^k+1+β^k+1p_k，若则令u₀＝u_k+1，返回S442；否则k＝k+1，返回S443。

综上所述，由于采用了上述技术方案，本发明的有益效果是：

本发明在很少的投影数据情况下也能重建出高质量的断层图像，适用于不完全投影数据图像重建问题，具有使用投影数据少，重建效率高的特点，可以显著减少人体的辐射剂量或缩短被检测样品的扫描时间，在无损探测成像领域具有广阔的应用前景。

附图说明

本发明将通过例子并参照附图的方式说明，其中：

图1是本发明流程图；

图2是本发明的断层扫描成像模型

图3是本发明数值模拟时所使用的重建对象

具体实施方式

本说明书中公开的所有特征，或公开的所有方法或过程中的步骤，除了互相排斥的特征和/或步骤以外，均可以以任何方式组合。

本说明书(包括任何附加权利要求、摘要和附图)中公开的任一特征，除非特别叙述，均可被其他等效或具有类似目的的替代特征加以替换。即，除非特别叙述，每个特征只是一系列等效或类似特征中的一个例子而已。

本发明相关说明：

1、成像方式：平行束扫描、扇形扫描或锥形扫描。

2、无约束最小化的迭代算法包括牛顿算法，最速下降算法或共轭梯度算法。

3、minf(u)＝min‖Au-y‖²该式主要考虑目标函数应当使重建断层图像的再投影尽可能的接近于实际扫描获得的投影数据，即Au与y的残差最小。

4、步骤6中：在共轭梯度算法中，线性搜索求解迭代步长α_k需要进行大量的运算，从而耗费大量的运算时间，降低了重建速度，在本发明中采用固定的步长因子，从而有效提高了图像重建速度。

如图1所示，本发明的核心思想是将目标泛函优化方法引入到少量投影CT图像重建这个欠定问题中，结合共轭梯度算法进行少量投影CT图像优化重建，提出了一种广义变分最小化的少量投影图像重建新方法。选用著名的Sheep-Logan头部模型(以下简称S-L头模型)作为重建测试对象进行数值模拟试验。该模型有10个位置、大小、方向、密度各异的椭圆组成，象征一个脑断层图像，见图3所示。其中，图3为10个椭圆的分布图，其中英文字母代表10个椭圆的编号(分别是a、b、c、d、e、f、g、h、i、j)，10个椭圆的中心位置、长轴、短轴、旋转角度、及密度分布值由表1给出。

表1S-L头部模型的结构参数

序号	中心坐标	长轴	短轴	旋转角度(度)	密度分布
						a	(0，0)	0.92	0.69	90	1.0
b	(0，-0.0184)	0.874	0.6624	90	0.2
						c	(0.22，0)	0.31	0.11	72	0.0
d	(-0.22，0)	0.41	0.16	108	0.3
						e	(0，0.35)	0.25	0.21	90	0.0
f	(0，0.1)	0.046	0.046	0	0.4
						g	(0，-0.1)	0.046	0.046	0	0.3
h	(-0.08，-0.0605)	0.046	0.023	0	0.3
						i	(0，-0.0605)	0.023	0.023	90	0.3
j	(0.06，-0.0605)	0.046	0.023	90	0.3

下面结合附图，对本发明作进一步描述，如图1所示，该图像重建方法包括以下步骤：

S1：如图2所示，根据CT扫描成像方式，建立CT扫描的成像模型的投影矩阵。图中θ表示扫描角度，S表示射线源，u表示待重建的断层图像，D表示探测器。CT扫描时，通过射线源S的旋转(或被测物的旋转)，可以获得不同角度下的投影数据，设投影数目为M，断层图像的像素数为N_image，则CT扫描的成像模型可以表示为如下的方程组y＝Au。其中，y表示由投影数据构成的一个M×1列向量，其元素y_i表示第i个投影数据，1≤i≤M，M的大小由投影角度数目和每个投影角度下的投影数据数(即探测器单元数目)决定，即M＝θ_n×D_n；θ_n表示投影角度数，此处θ_n＝12；D_n表示每个投影角度下的探测器数目，此处D_n＝256。u表示待重建的断层图像，是一个N_image×1的列向量，其元素u_j表示第j个像素的密度值，1≤j≤N，N_image是断层图像的总像素数，N_image＝W×H＝256×256，W和H分别表示断层图像的宽度和高度(此处W＝H)；A表示一个M×N_image的投影矩阵，其元素A_ij表示断层图像第j个像素对第i个投影数据的贡献，即图中扫描线与像素网格相交的长度。

为了模拟少量投影CT图像重建问题，我们在0-180度内每隔15度均匀选取12个投影数据进行数值模拟试验。当获取投影数据的角度间隔越大时，投影数据越少，其数据越不完备，导致重建质量越差；反之，当获取投影数据的角度间隔越小时，投影数据越多，数据越完备，图像重建质量越好。数值模拟试验中，总的投影数为12×256(M个像元，即为矩阵y的大小)，投影矩阵为256×256×256×12(即为矩阵A的大小)，断层图像大小为256×256(即为矩阵u的大小)。

S2：根据CT扫描的成像模型的投影矩阵，将图像重建问题可转化为最优化模型；建立用于图像重建的优化模型，其优化目标可表示为minf(u)＝min‖Au-y‖²，该式主要考虑目标函数应当使重建断层图像的再投影尽可能的接近于实际扫描获得的投影数据，即Au与y的方差最小。

S3：采用无约束最小化的迭代算法(共轭梯度算法等)进行图像重建模型的优化求解，其迭代步骤为：

S31：选取初值u₀，计算目标泛函f(u)梯度g₀，g₀＝Au₀-y；设置迭代停止时的允许误差0＜ε＜10^-3；

S32：计算搜索方向p₀，p₀＝-g₀，其中当迭代次数为k不为1时，p_k由S37给出；

S33：设置迭代步长α_k，α_k∈4×10^-5～6×10^-5；

S34：更新重建图像u_k+1，u_k+1＝u_k+α_kp_k；

S35：更新目标泛函梯度g_k+1，g_k+1＝Au_k+1-p_k；

S36：计算梯度搜索方向β_k+1，

S37:更新搜索方向p_k+1，p_k+1＝-g_k+1+β_k+1p_k；

S38：如果‖g_k+1‖＜ε，停止，输出u＝u_k+1；否则令k＝0，转到S33。

其中，共轭梯度算法选择Fletcher-Reeves方法计算其搜索方向，在共轭梯度算法中，线性搜索求解迭代步长α_k需要进行大量的运算，从而耗费大量的运算时间，降低了重建速度，在本发明中采用固定的步长因子(因为每次迭代时计算的步长因子都在一个较小的范围内浮动)，从而有效提高了图像重建速度。

S4：建立广义变分最小化的断层图像修正模型对优化重建后的断层图像进行修正，其中广义变分最小化模型表示为其中，GV表示图像的广义变分，计算公式为u_x和u_x分别表示断层图像在x方向和y方向的梯度值，δ＞0，是个很小的数，用以保证GV不为零：s为广义调整参数，本发明中，s的选取范围为0＜s≤1。

S5：采用共轭梯度算法求广义变分最小，共轭梯度算法求GV最小则采用了迭代n步后重新开始的方法，其下降方向是利用上一次搜索方向和本次起始点的负梯度的线性组合生成新的共轭方向，形成新一轮的共轭梯度算法迭代，其具体迭代步骤为：

S51：将上一步骤的优化重建结果设为初始值u₀，重新开始迭代次数设为n，设置迭代停止时的允许误差ε＞0；

S52：计算目标函数负梯度

S53：采用固定迭代步长，α_k∈1×10^-4～2×10^-4；

S54：更新修正值，u_k+1＝u_k+α_kp_k；

S55：判断是否成立，如成立则停止迭代，输出u＝u_k+1否则转S56；

S56：若k＝n-1，令u₀＝u_k+1，k＝0，进行S52；否则进行S57；

S57：利用公式计算F-R搜索方向更新负梯度方向p^k+1＝-g^k+1+β^k+1p_k，若则令u₀＝u_k+1，转S53，否则k＝k+1转S54。

S6：判断修正后的重建图像是否满足迭代停止准则，若满足则停止重建，将S3所得图像作为结果输出重建结果；否则，将步骤5)所得图像作为初值重复S3～S4，直至满足迭代停止准则。

在共轭梯度算法中，线性搜索求解迭代步长α_k需要进行大量的运算，从而耗费大量的运算时间，降低了重建速度，在本发明中采用固定的步长因子，从而有效提高了图像重建速度。

数值模拟试验中，首先利用FBP和ART算法对数值模拟投影数据进行重建，然后采用GV最小化算法进行重建，分别选取广义调整参数s＝1和s＝0.75，最后对这些算法的数值实验结果进行分析和比较。

在投影数据较少的情况，FBP方法重建效果较差，没有重建出图像的细节，密度分辨率不高；ART方法重建结果也不好，图像细节也没有被重建出来，其密度分辨也不高；这说明，在投影数据较少的情况下，采用现有的方法无法得到较好的重建结果；采用广义变分最小化的图像重建算法，可以得到较好的图像重建结果，当广义参数s＝1时，重建算法类似于TV最小化重建算法，其重建效果较好，但是有部分细节没有重建出来；当s进一步减小时，可以获得更好的重建效果，细节分辨能力强，且重建精度高，具有很高密度分辨，这说明采用广义变分最小化的少数图像重建算法是一种解决少数投影数据重建的有效方法，具有非常优越的性能。

采用均方误差(MSE)来客观评价上述几种算法的重建性能，表2列举了上述重建结果对应的均方误差。

表2不同重建方法的重建质量指标比较

由表2可知，广义变分最小化的少量投影CT图像重建方法的重建性能最好，重建质量最高，较之ART算法其均方根误差高6个数量级，较之Sidky等人提出方法(即s＝1时的GV方法)高4个数量级，是一种性能优良的少量投影CT图像重建方法。

本发明并不局限于前述的具体实施方式。本发明扩展到任何在本说明书中披露的新特征或任何新的组合，以及披露的任一新的方法或过程的步骤或任何新的组合。

Claims (5)

1.一种基于广义变分最小化的少量投影CT图像重建方法，其特征在于包括：

S1：根据CT扫描成像方式，建立CT扫描的成像模型的投影矩阵；

S2：根据CT扫描的成像模型的投影矩阵，建立用于图像重建的优化模型；

S3：采用无约束最小化的迭代算法对图像重建优化模型中的图像进行优化求解，得到重建模型的重建数据；

S4:建立广义变分最小化的图像修正模型，对S3中重建模型的重建数据进行修正，获得修正后的重建图像；判断修正后的重建图像是否满足迭代停止准则，若满足则停止重建图像，进行S5；否则，返回S3，直至满足迭代停止准则；S5:所得图像作为结果输出重建结果；

所述S1中CT扫描成像方式包括平行束扫描、扇形束扫描或者锥形束扫描；

所述步骤S2具体步骤包括：

S21：根据公式(1)建立图像重建的正向投影模型，

y＝Au(1)；

其中，u表示待重建的断层图像；CT扫描时，通过射线源S的旋转，可以获得不同角度下的投影数据，A表示一个M×N_image的投影矩阵，其元素A_ij表示断层图像扫描线与像素网格相交的长度，设总的投影数目为M，断层图像的像素数为N_image，y表示由投影数据构成的一个M×1列向量，其元素yi表示第i个投影数据，1≤i≤M；M表示总的投影数目，其大小由投影角度数目和每个投影角度下的投影数据数，即探测器单元数目决定，即M＝θ_n×D_n，θ_n表示投影角度的数量，D_n表示探测器的单元数目；u表示待重建的断层图像，是一个N_image×1的列向量，其元素uj表示第j个像素的密度值，1≤j≤N_image，N是断层图像的总像素数，N＝W×H，W和H分别表示断层图像的宽度和高度；

S22：得到重建模型的优化重建的数学模型，即使待重建的断层图像的再投影Au与实际扫描获得的投影值y的残差最小值minf(u)：

minf(u)＝min||Au-y||²(2)。

2.根据权利要求1所述的一种基于广义变分最小化的少量投影CT图像重建方法，其特征在于所述S3中无约束最小化的迭代算法包括牛顿算法，最速下降算法或共轭梯度算法。

3.据权利要求2所述的一种基于广义变分最小化的少量投影CT图像重建方法，其特征在于所述S3中采用共轭梯度算法具体步骤：

S31：选取初值u₀，计算目标泛函f(u)梯度g₀，g₀＝Au₀-y；设置迭代停止时的允许误差0<ε<10^-3；

S32：计算搜索方向p₀，p0＝-g0，其中当迭代次数k不为1时，pk由S37给出；

S33：设置迭代步长α_k，α_k∈4×10^-5～6×10^-5；

S34：更新重建图像u_k+1，u_k+1＝u_k+α_kp_k；

S35：更新目标泛函梯度g_k+1，g_k+1＝Au_k+1-p_k；

S36：计算梯度搜索方向β_k+1，

S37：更新搜索方向p_k+1，p_k+1＝-g_k+1+β_k+1p_k；

S38：如果||g_k+1||<ε，停止，输出u＝u_k+1；否则令k＝0，转到S33；

其中，共轭梯度算法选择Fletcher-Reeves方法计算其搜索方向。

4.根据权利要求1至3之一所述的一种基于广义变分最小化的少量投影CT图像重建方法，其特征在于所述S4具体步骤包括：

S41：根据公式(3)计算GV，

$GV=\&Integral;{\left(\sqrt{{u_x}^2+{u_y}^2+\mathrm{\&delta;}}\right)}^{s}dxdy---\left(3\right)$

并根据公式(3)得到广义变分最小化的图像修正模型

其中，GV表示图像的广义变分范数，u_x和u_x分别表示断层图像在x方向和y方向的梯度值，δ＝10^-6，s为广义调整参数，s的选取范围为0<s≤1；

S42：通过设置迭代步长α_k∈1×10^-4～2×10^-4的共轭梯度算法求GV最小，对重建数据进行修正，并采用了迭代n步后重新开始的方法，其下降方向是利用上一次搜索方向和本次起始点的负梯度的线性组合生成新的共轭方向，形成新一轮的共轭梯度算法迭代，直至满足收敛条件获得修正后的重建图像。

5.根据权利要求根据权利要求4所述的一种基于广义变分最小化的少量投影CT图像重建方法，其特征在于所述S42具体步骤：

S441：将S3重建模型的重建数据设为初始值u0，重新开始迭代次数设为n，设置迭代停止时的允许误差10^-8<ε<10^-3，n>0；

S442：计算目标函数负梯度，

S443：采用固定迭代步长，α_k∈1×10^-4～2×10^-4；

S444：更新修正值，u_k+1＝u_k+α_kp_k；

S445：判断是否满足迭代停止准则，即是否成立，如成立则停止迭代，输出u＝u_k+1；否则进行S443；

S446：若k＝n-1，令u₀＝u_k+1，k＝0，返回S442；否则进行S447；

S447：利用公式计算F-R搜索方向更新负梯度方向p^k+1＝-g^k+1+β^{k+ 1}p_k，若则令u₀＝u_k+1，返回S442；否则k＝k+1，返回S443。