CN103390285A - 基于边缘引导的锥束ct不完全角度重建方法 - Google Patents

Abstract

本发明涉及一种基于边缘引导的锥束CT不完全角度重建方法，具体含有如下步骤：步骤1：估计初始图像：利用扫描到的投影数据估计初始重建图像；步骤2：图像边缘提取；步骤3：设计加权因子；步骤4：更新优化模型；步骤5：基于稀疏优化锥束CT不完全角度重建；步骤6：判断重建质量是否达到要求？如是，则执行步骤7；如不是，则执行步骤2；步骤7：结束；本发明提供了一种效率高、重建图像质量好的基于边缘引导的锥束CT不完全角度重建方法。

Description

基于边缘引导的锥束CT不完全角度重建方法

（一）、技术领域：本发明涉及一种锥束CT不完全角度重建方法，特别是涉及一种基于边缘引导的锥束CT不完全角度重建方法。

（二）、背景技术:计算机断层成像技术（Computed Tomography,CT）作为一种现代成像技术已经广泛应用于医学、工业等领域。然而，在很多实际应用中，由于受数据采集时间或成像***扫描的几何位置约束，只能在不完全角度范围或在较少的投影角度得到数据，这些都属于不完全角度问题(IncompleteData Problem)。提升在不完全角度扫描下CT图像重建的质量具有重要的理论研究和工程实践意义，如何设计高精度的不完全角度下CT图像重建的方法也是研究的热点和难点问题。

针对锥束CT不完全角度问题，其数据采集不满足精确重建数据完备性条件，解析类重建算法无法获得较高质量的重建图像。迭代类重建算法对数据完备性没有严格的要求，能够取得相对解析类算法较优的重建质量。经典算法为代数迭代技术，代数迭代算法具有一定的抗数据缺失性，通常相同数据量上结果优于解析算法，但是占用计算存储资源，需要较强硬件支持，实际中重建较慢。基于压缩感知理论的CT重建算法通过挖掘待重建物体的先验知识并刻画物体的稀疏特性，能够在压缩采样的情形下比经典迭代算法更优的重建质量。

尽管基于CS理论模型的重建算法能够在较稀疏采样下获得较好的重建结果，但是却没有更深层次地挖掘出能够更加有利于改善重建结果的图像信息。图像的边缘成分是图像的重要要素，已经在图像修复、分割等多种图像应用领域有着重要应用。在图像重建中，有效利用边缘信息能够较大程度改善重建质量。

针对不完全角度CT图像重建问题，经典的处理方法为代数迭代和统计迭代重建算法，例如代数迭代算法（Algebraic Reconstruction Technique,ART）和期望最大（expectation maximization,EM）算法，迭代算法在处理缺失数据较少时通常比解析算法有较好的重建质量，如1984年由Feldkamp等人提出的FDK算法有较好的重建质量。近年来，基于压缩感知理论的指导下，发展起来了一些能够适用于稀疏采样的重建算法，较为典型的有2008年Sidky和Pan等人提出的ASD-POCS(adapt-steepest-descent projection onto convex sets)算法。该算法以TV范数最小化设计优化模型，采用POCS（ART）与TV最速下降交替执行的策略，能够在稀疏采样下重建出较高质量的图像。

在ASD-POCS算法被提出之后，很多学者在模型求解上面又提出了多种优化算法，例如2011年由Vandeghinste等人提出的Split-BregmanTV算法，2012年由Zhang等人提出的ADTVM算法等。其他算法引入另外的先验等等，例如PICCS（prior image constrained compressed sensing,PICCS），2010年王林元等人提出的RRD（reconstruction-reference difference,RRD）等算法均是需求待重建物体在某种先验下的稀疏表达，利用稀疏表达构造优化模型并设计相应的求解算法，从而达到图像重建的目的。另外频域方面的有限角度重建（主要MRI成像）较为强有力的工具为RecPF（reconstruction from partial Fourier space smpling）算法，该算法不涉及大规模矩阵的存储和计算，利用差分算子D_i的良好的性质以及FFT技术实现重建的高效和准确，是CS理论在MRI成功应用的范例。RecPF在CT图像重建领域也有着较强的应用前景，张瀚铭等人基于过采样和稀疏插值技术成功实现较高精度的CT图像重建。

2010年，Yin等人提出利用迭代支集探测（iterative support detection,ISD）的方法在稀疏信号恢复中取得了较优的重建结果。在该方法中，不断对迭代过程中间信号做支撑集探测的处理，根据所探测到的部分支集对优化模型的目标函数进行更新，采取这种迭代方法能在极度欠采样的基础上有效恢复信号。纯粹的基于支集探测的方法并不能直接利用到锥束CT的图像重建中，利用图像稀疏表示集下的支撑集能够为CT图像重建带来重建质量的提升。

（三）、发明内容：

本发明要解决的技术问题是：克服现有技术的缺陷，提供一种效率高、重建图像质量好的基于边缘引导的锥束CT不完全角度重建方法。

本发明的技术方案：

一种基于边缘引导的锥束CT不完全角度重建方法，含有如下步骤：

步骤1：估计初始图像：利用扫描到的投影数据估计初始重建图像；

步骤2：图像边缘提取；

步骤3：设计加权因子；

步骤4：更新优化模型；

步骤5：基于稀疏优化锥束CT不完全角度重建；

步骤6：判断重建质量是否达到要求？如是，则执行步骤7；如不是，则执行步骤2；

步骤7：结束。

步骤1的具体估计方法如下：

步骤1.1：将图像重建问题刻画成如下稀疏模型：

min||x||_TV

s.t.Ax=b

其中，||x||_TV为待重建物体x的总变分（TV）范数，A为***矩阵，向量b为扫描到的投影数据；

步骤1.2：用交替方向法对步骤1.1中的稀疏模型进行求解，迭代次数为N，得到求解公式如下：

$\left\{\begin{array}{c}{x}^{k+1}={(A^{T}A+\underset{j}{\mathrm{\Σ}}{\mathrm{\ρ}}_j{D}_j^T{D}_j)}^{+}(A^{T}b+\underset{j}{\mathrm{\Σ}}{\mathrm{\ρ}}_j{D}_j^T({z_j}^k-{u_j}^k/{\mathrm{\ρ}}_j))\\ {z_j}^{k+1}=\max \left\{\right|D_j{x}^{k+1}+\frac{{u_j}^k}{{\mathrm{\ρ}}_j}|-\frac{\mathrm{\λ}}{{\mathrm{\ρ}}_j},0\}\mathrm{sgn}(D_j{x}^{k+1}+\frac{{u_j}^k}{{\mathrm{\ρ}}_j})\\ {u_j}^{k+1}={u_j}^k+{\mathrm{\ρ}}_j(D_j{x}^{k+1}-{z_j}^{k+1})\end{array}\right.$

其中，(·)⁺为矩阵的Moore-Penrose伪逆，A^T为***矩阵的转置，u_j为j方向的更新因子，u_j ^k为第k次迭代后的j方向的更新因子，u_j ^k+1为第k+1次迭代后的j方向的更新因子，ρ_j为j方向的L2范数惩罚因子，λ为L1范数惩罚因子，D_j为j方向上的梯度算子矩阵，D_j ^T为j方向上的梯度算子矩阵的转置，z_j为j方向梯度图像，z_j ^k为第k次迭代后的j方向梯度图像，z_j ^k+1为第k+1次迭代后的j方向梯度图像，max{·},sgn{·}分别为求最大函数和求符号函数，向量b为扫描到的投影数据，x^k+1为第k+1次迭代后的重建图像，经过N次迭代后的重建图像记作x^(N)，x^(N)为最后估计的初始重建图像。

步骤1.2中的迭代次N为总收敛轮数的0.5～0.2倍。

步骤2的具体方法为：设步骤1.2中某次迭代产生的中间重建图像为x1，利用降噪因子F对中间重建图像x1进行卷积运算：x1_F=x1***F，利用经典边缘算子E求取中间重建图像x1的初步边缘x1_e-mid：x1_e-mid=E[x1]，中间重建图像x1的边缘x1_e为：x1_e=x1_e-mid∩x1_LMI，其中，x1_LMI为中间重建图像x1的局部互信息图像；

步骤3的具体方法为：根据步骤2中的中间重建图像x1的边缘x1_e和L1范数惩罚因子λ确定中间重建图像x1在坐标(i1,j1)处的加权因子w(i1,j1)，该加权因子w(i1,j1)=λx1_e(i1,j1)，其中，x1_e(i1,j1)为中间重建图像x1在坐标(i1,j1)处的边缘；

步骤4的具体方法为：根据步骤3中的加权因子w(i1,j1)计算j方向的加权对角矩阵M_j，M_j＝diag(w(1，1)，w(1，2)，...w(i1，j1)，...，w(N_x，N_y))，那么优化模型更新为如下表达式：

$\min {\left|\right|z_1\left|\right|}_1+{\left|\right|z_2\left|\right|}_1+{\left|\right|z_3\left|\right|}_1$

$s.t.\left\{\begin{array}{c}\mathrm{Ax}=b\\ M_j{D}_{1}x=z_1\\ M_j{D}_{2}x=z_2\\ M_j{D}_{3}x=z_3\end{array}\right.$

其中，N_x为中间重建图像x1在x坐标上的规模，N_y为中间重建图像x1在y坐标上的规模，D₁为1方向上的梯度算子矩阵，D₂为2方向上的梯度算子矩阵，D₃为3方向上的梯度算子矩阵；

步骤5的具体方法为：采用增广拉格朗日乘子法将步骤4中的更新后的带约束的优化模型转为无约束的优化模型，具体公式如下：

$\min \frac{1}{2}{\left|\right|\mathrm{Ax}-b\left|\right|}_2^2+\underset{j}{\mathrm{\Σ}}(u_j^T(M_j{D}_{j}x-z_j)+\frac{{\mathrm{\ρ}}_j}{2}{{\left|\right|M_j{D}_{j}x-z_j\left|\right|}_2^2+\mathrm{\λ}\left|\right|z_j\left|\right|}_1)$

其中，

为j方向的更新因子的转置；

对上述无约束的优化模型采用分离变量，利用交替方向法求最小，迭代公式如下：

$\left\{\begin{array}{c}{x}^{k+1}={(A^{T}A+\underset{j}{\mathrm{\Σ}}{\mathrm{\ρ}}_j{D_j}^T{M}_j{D}_j)}^{+}(A^{T}b+\underset{j}{\mathrm{\Σ}}{\mathrm{\ρ}}_j{D_j}^T{M_j}^T(z_j^k-{u_j}^k/{\mathrm{\ρ}}_j))\\ z_j^{k+1}=\max \left\{\right|M_j{D}_j{x}^{k+1}+\frac{{u_j}^k}{{\mathrm{\ρ}}_j}|-\frac{\mathrm{\λ}}{{\mathrm{\ρ}}_j},0\}\mathrm{sgn}[M_j{D}_j{x}^{k+1}+\frac{u_j^k}{{\mathrm{\ρ}}_j}]\\ u_j^{k+1}=u_j^k+{\mathrm{\ρ}}_j(M_j{D}_j{x}^{k+1}-z_j^{k+1})\end{array}\right.$

x^k+1即为本轮重建图像；

步骤6中的重建质量达到要求是指：本轮重建图像与上一轮重建图像相比无显著变化；首次进行步骤5中的重建时，上一轮重建图像是指步骤1中的初始重建图像。

降噪因子F为高斯核函数或中值滤波器，经典边缘算子E为canny边缘算子、soble边缘算子、prewitt边缘算子、roberts边缘算子、log边缘算子、zeroscross边缘算子中的任一种，L1范数惩罚因子λ的取值在0～1之间。

本锥束CT不完全角度重建方法在基于总变分最小化的基础上，提出将图像的边缘信息与之结合的方法，设计有效的边缘探测算子，在重建的迭代过程中不断利用探测到的边缘信息对重建的优化模型的目标函数进行更新，使得重建算法能够适应更少的采集数据并且提升重建质量。

本发明的有益效果：

本发明利用基于边缘引导的稀疏优化迭代重建技术实现了锥束CT不完角度的高精度图像重建；针对不完全角度或稀疏采样下，传统解析算法重建伴有强烈伪影，严重影响有用信息的辨别的问题；采样迭代重建的技术在迭代的过程中利用中间重建图像的有效的边缘信息，提出了基于边缘引导的锥束CT不完全角度的重建方法；与已有的基于总变分最小化稀疏重建方法不同，本发明结合了图像的边缘与梯度稀疏先验信息，在每轮迭代中选取合适的边缘信息，并设计适当的加权因子，利用加权因子对优化的目标函数进行更新，对更新后的目标函数进行基于交替方向法的迭代策略获取更加高质量的重建图像。本发明利用准确的边缘探测算法，结合高效的优化策略，提升了收敛的效率和重建图像的质量。

（四）、附图说明：

图1为锥束CT扫描模型；

图2为边缘检测示例；

图3为实际三维重建结果；

图4为仿真重建收敛曲线。

（五）、具体实施方式：

基于边缘引导的锥束CT不完全角度重建方法含有如下步骤：

步骤1：估计初始图像：利用扫描到的投影数据估计初始重建图像；

步骤2：图像边缘提取；

步骤3：设计加权因子；

步骤4：更新优化模型；

步骤5：基于稀疏优化锥束CT不完全角度重建；

步骤6：判断重建质量是否达到要求？如是，则执行步骤7；如不是，则执行步骤2；

步骤7：结束。

步骤1的具体估计方法如下：

步骤1.1：将图像重建问题刻画成如下稀疏模型：

min||x||_TV

s.t.Ax=b

其中，||x||_TV为待重建物体x的总变分（TV）范数，A为***矩阵，向量b为扫描到的投影数据；

步骤1.2：用交替方向法对步骤1.1中的稀疏模型进行求解，迭代次数为N，得到求解公式如下：

$\left\{\begin{array}{c}{x}^{k+1}={(A^{T}A+\underset{j}{\mathrm{\Σ}}{\mathrm{\ρ}}_j{D}_j^T{D}_j)}^{+}(A^{T}b+\underset{j}{\mathrm{\Σ}}{\mathrm{\ρ}}_j{D}_j^T({z_j}^k-{u_j}^k/{\mathrm{\ρ}}_j))\\ {z_j}^{k+1}=\max \left\{\right|D_j{x}^{k+1}+\frac{{u_j}^k}{{\mathrm{\ρ}}_j}|-\frac{\mathrm{\λ}}{{\mathrm{\ρ}}_j},0\}\mathrm{sgn}(D_j{x}^{k+1}+\frac{{u_j}^k}{{\mathrm{\ρ}}_j})\\ {u_j}^{k+1}={u_j}^k+{\mathrm{\ρ}}_j(D_j{x}^{k+1}-{z_j}^{k+1})\end{array}\right.$

其中，(·)⁺为矩阵的Moore-Penrose伪逆，A^T为***矩阵的转置，u_j为j方向的更新因子，u_j ^k为第k次迭代后的j方向的更新因子，u_j ^k+1为第k+1次迭代后的j方向的更新因子，ρ_j为j方向的L2范数惩罚因子，λ为L1范数惩罚因子，D_j为j方向上的梯度算子矩阵，D_j ^T为j方向上的梯度算子矩阵的转置，z_j为j方向梯度图像，z_j ^k为第k次迭代后的j方向梯度图像，z_j ^k+1为第k+1次迭代后的j方向梯度图像，max{·},sgn{·｝分别为求最大函数和求符号函数，向量b为扫描到的投影数据，x^k+1为第k+1次迭代后的重建图像，经过N次迭代后的重建图像记作x^(N)，x^(N)为最后估计的初始重建图像。

步骤1.2中的迭代次N为总收敛轮数的0.5～0.2倍。

步骤2的具体方法为：设步骤1.2中某次迭代产生的中间重建图像为x1，利用降噪因子F对中间重建图像x1进行卷积运算：x1_F=x1***F，利用经典边缘算子E求取中间重建图像x1的初步边缘x1_e-mid：x1_e-mid=E[x1]，中间重建图像x1的边缘x1_e为：x1_e=x1_e-mid∩x1_LMI，其中，x1_LMI为中间重建图像x1的局部互信息图像；

步骤3的具体方法为：根据步骤2中的中间重建图像x1的边缘x1_e和L1范数惩罚因子λ确定中间重建图像x1在坐标(i1,j1)处的加权因子w(i1,j1)，该加权因子w(i1,j1)=λx1_e(i1,j1)，其中，x1_e(i1,j1)为中间重建图像x1在坐标(i1,j1)处的边缘；

步骤4的具体方法为：根据步骤3中的加权因子w(i1,j1)计算j方向的加权对角矩阵M_j，M_j＝diag(w(1，1)，w(1，2)，...w(i1，j1)，...，w(N_x，N_y))，那么优化模型更新为如下表达式：

$\min {\left|\right|z_1\left|\right|}_1+{\left|\right|z_2\left|\right|}_1+{\left|\right|z_3\left|\right|}_1$

$s.t.\left\{\begin{array}{c}\mathrm{Ax}=b\\ M_j{D}_{1}x=z_1\\ M_j{D}_{2}x=z_2\\ M_j{D}_{3}x=z_3\end{array}\right.$

其中，N_x为中间重建图像x1在x坐标上的规模，N_y为中间重建图像x1在y坐标上的规模，D₁为1方向上的梯度算子矩阵，D₂为2方向上的梯度算子矩阵，D₃为3方向上的梯度算子矩阵；

步骤5的具体方法为：采用增广拉格朗日乘子法将步骤4中的更新后的带约束的优化模型转为无约束的优化模型，具体公式如下：

$\min \frac{1}{2}{\left|\right|\mathrm{Ax}-b\left|\right|}_2^2+\underset{j}{\mathrm{\Σ}}(u_j^T(M_j{D}_{j}x-z_j)+\frac{{\mathrm{\ρ}}_j}{2}{{\left|\right|M_j{D}_{j}x-z_j\left|\right|}_2^2+\mathrm{\λ}\left|\right|z_j\left|\right|}_1)$

其中，

为j方向的更新因子的转置；

对上述无约束的优化模型采用分离变量，利用交替方向法求最小，迭代公式如下：

$\left\{\begin{array}{c}{x}^{k+1}={(A^{T}A+\underset{j}{\mathrm{\Σ}}{\mathrm{\ρ}}_j{D_j}^T{M}_j{D}_j)}^{+}(A^{T}b+\underset{j}{\mathrm{\Σ}}{\mathrm{\ρ}}_j{D_j}^T{M_j}^T(z_j^k-{u_j}^k/{\mathrm{\ρ}}_j))\\ z_j^{k+1}=\max \left\{\right|M_j{D}_j{x}^{k+1}+\frac{{u_j}^k}{{\mathrm{\ρ}}_j}|-\frac{\mathrm{\λ}}{{\mathrm{\ρ}}_j},0\}\mathrm{sgn}[M_j{D}_j{x}^{k+1}+\frac{u_j^k}{{\mathrm{\ρ}}_j}]\\ u_j^{k+1}=u_j^k+{\mathrm{\ρ}}_j(M_j{D}_j{x}^{k+1}-z_j^{k+1})\end{array}\right.$

x^k+1即为本轮重建图像；

步骤6中的重建质量达到要求是指：本轮重建图像与上一轮重建图像相比无显著变化；首次进行步骤5中的重建时，上一轮重建图像是指步骤1中的初始重建图像。

降噪因子F为高斯核函数或中值滤波器，经典边缘算子E为canny边缘算子、soble边缘算子、prewitt边缘算子、roberts边缘算子、log边缘算子、zeroscross边缘算子中的任一种，L1范数惩罚因子λ的取值在0～1之间。

本锥束CT不完全角度重建方法在基于总变分最小化的基础上，提出将图像的边缘信息与之结合的方法，设计有效的边缘探测算子，在重建的迭代过程中不断利用探测到的边缘信息对重建的优化模型的目标函数进行更新，使得重建算法能够适应更少的采集数据并且提升重建质量。

为了在实际应用中使得不完全角度扫描下的锥束CT能够具有高精度的成像结果，更好地适应实际应用的需求，需对不完全角度扫描下的锥束CT投影数据高精度重建进行研究。本锥束CT不完全角度重建方法基于现有的压缩感知理论和迭代支集探测的研究成果，采用先进稀疏优化算法提出了基于边缘引导的高精度不完全角度扫描锥束CT重建算法。本锥束CT不完全角度重建方法从初始估计图像出发，通过在重建过程中加入边缘探测，并利用所探测边缘对优化的目标函数进行改进，通过不断迭代及循环机制，提高重建精度和收敛速度。

下面通过实例进一步说明本锥束CT不完全角度重建方法：

步骤1：估计初始图像；

锥束CT扫描模型如图1所示，利用扫描到的不完全角度扫描下的锥束CT投影数据，估计初始重建图像，估计方法为：

①不完全角度扫描下的锥束CT图像重建问题：

min||x||_TV

s.t.Ax=b

转为无约束问题：

$\min \frac{1}{2}{\left|\right|\mathrm{Ax}-b\left|\right|}^2+{\mathrm{\λ}\left|\right|x\left|\right|}_{\mathrm{TV}}$

使用ADM方法分离变量可得：

$\min \frac{1}{2}{\left|\right|\mathrm{Ax}-b\left|\right|}^2+\mathrm{\λ}\left(\underset{i}{\mathrm{\Σ}}{\left|\right|z_i\left|\right|}_1\right)$

s.t.D_ix=z_i

②我们在交替方向法框架的基础上，针对以下TV最小化模型给出一种具体的重建算法，即基于增广拉格朗日的交替方向法的总变分最小化算法。

其对应的增广Lagrange函数：

$\min \frac{1}{2}{\left|\right|\mathrm{Ax}-b\left|\right|}^2+\underset{i}{\mathrm{\Σ}}(u_i^T(D_{i}x-z_i)+\frac{{\mathrm{\ρ}}_i}{2}{\left|\right|D_{i}x-z_i\left|\right|}^2{+\mathrm{\λ}\left|\right|z_i\left|\right|}_1)$

令(x^*,z^*)为L_A的最小值点，则乘子更新公式为：

${\tilde{v}}_i=v_i-{\mathrm{\β}}_i(D_i{x}^{*}-z_i^{*})$

$\widetilde{\mathrm{\λ}}=\mathrm{\λ}-\mathrm{\μ}({\mathrm{Ax}}^{*}-b).$

下面使用交替方向法来求解L_A的最小值点，用

表示第k轮迭代的近似最小值点。关于x的部分为：

$\min {\left|\right|\mathrm{Ax}-b\left|\right|}^2+\underset{i}{\mathrm{\Σ}}\left({\mathrm{\ρ}}_i{\left|\right|D_{i}x-z_i+u_i/{\mathrm{\ρ}}_i\left|\right|}^2\right)$

即“x-子问题”，这是一个二次型求最小值的问题，对其进行求导并令d_k(x)=0得到f_k(x)的严格最小值点

$x^{*}={(A^{T}A+\underset{i}{\mathrm{\Σ}}\left({\mathrm{\ρ}}_i{D}_i^T{D}_i\right))}^{+}(A^{T}b+\underset{i}{\mathrm{\Σ}}\left({\mathrm{\ρ}}_i{D}_i^T(z_i-u_i/{\mathrm{\ρ}}_i)\right)),$

其中M⁺表示M的Moore-Penrose伪逆。理论上讲通过求伪逆求解x-子问题是能得到严格的最小值点，然而，在每轮迭代中计算逆或伪逆的数值计算开销是非常巨大的，因此通常使用迭代方法求解。

最速下降方法能够迭代求解x-子问题，然而对于较大规模的问题，也不是一个高效的算法。事实上增广Lagrangian函数是通过交替求解z-子问题和x-子问题达到最小的。因此，在每个步骤精确求解x-子问题是不必要的，只需要使用一步的最速下降结果作为代替。在下降方向的步长选择上可以使用BB（BarzilaiandBorwein）型方法。

求出x^k+1后求解

可通过如下子问题求解：

$\underset{z_i}{\min }{L}_A(x^k,z_i)=\underset{i}{\mathrm{\Σ}}\left(\right|\left|z_i\right||-v_i^T(D_i{x}^k-z_i)+\frac{{\mathrm{\β}}_i}{2}\·{\left|\right|D_i{x}^k-z_i\left|\right|}_2^2)-{\mathrm{\λ}}^T({\mathrm{Ax}}^k-b)+\frac{\mathrm{\μ}}{2}\·{\left|\right|A{x}^k-b\left|\right|}_2^2$ 即“z-子问题”

$\underset{w_i}{\min }\underset{i}{\mathrm{\Σ}}\left(\right|\left|w_i\right||-v_i^T(D_i{\stackrel{\→}{f}}^k-w_i)+\frac{{\mathrm{\β}}_i}{2}\·{\left|\right|D_i{\stackrel{\→}{f}}^k-w_i\left|\right|}_2^2)$

z-子问题可对每个z_i单独求最小，其闭式解：

z_j ^*=shrinkage(D_jx+u_jρ_j,λ/ρ_j)；

ADM的迭代过程为：

$\left\{\begin{array}{c}{x}^{k+1}={(A^{T}A+\underset{i}{\mathrm{\Σ}}{\mathrm{\ρ}}_i{D}_i^T{D}_i)}^{+}(A^{T}b+\underset{i}{\mathrm{\Σ}}{\mathrm{\ρ}}_i{D}_i^T({z_i}^k-{u_i}^k/{\mathrm{\ρ}}_i))\\ {z_j}^{k+1}=\max \left\{\right|D_j{x}^{k+1}+\frac{{u_j}^k}{{\mathrm{\ρ}}_j}|-\frac{\mathrm{\λ}}{{\mathrm{\ρ}}_j},0\}\mathrm{sgn}(D_j{x}^{k+1}+\frac{{u_j}^k}{{\mathrm{\ρ}}_j})\\ {u_j}^{k+1}={u_j}^k+{\mathrm{\ρ}}_j(D_j{x}^{k+1}-{z_j}^{k+1})\end{array}\right.$

其中，||x||_TV为待重建物体x的总变分（TV）范数，(·)⁺为矩阵的Moore-Penrose伪逆，A为***矩阵，b为观测数据，A^T为***矩阵的转置，u_j为更新因子，ρ_j为j方向的L2范数惩罚因子，λ为L1范数惩罚因子，D_j为j方向上的梯度算子矩阵，z_j为j方向梯度图像，max{·},sgn{·}分别为求最大函数和求符号函数。

步骤2：图像边缘提取；

利用降噪因子F对图像x进行二维卷积x_F=x**F，F通常为高斯核函数或者中值滤波器等。利用经典边缘算子E，求取图像初步边缘x_e-mid=E[x]，E通常为canny、soble、prewitt、roberts、log、zeroscross等边缘算子。图像的边缘为：x_e=x_e-mid∩x_LMI，其中，x_LMI为图像局部互信息图像。

例证分析：如图2，所示对于测试体模(a)，未降噪之前图像中已混入了较多高频噪声，利用中值滤波器进行二维滤波后的图像(b)，噪声明显降低能够提升边缘检测的有效性，利用经典边缘检测算子进行检测得到中间过程的边缘图像(c)，中间过程的边缘图像通过局部互信息增强后得到了比较符合所用体模的边缘信息(d)。

步骤3：设计加权因子；

由步骤2所提取边缘x_e，选择比例系数λ，那么加权因子定为w(i,j)=λx_e(i,j)，其中比例系数λ的取值一般为0～1。

步骤4：更新优化模型；

由步骤3所得加权因子w(i,j)=λx_e(i,j)，计算对角矩阵

$M=\mathrm{diag}(d_{\mathrm{1,1}},d_{\mathrm{1,2}},...d_{i,j},...,d_{N_x,N_y}),$ d_i,j=λx_e(i,j)，那么优化模型更新为如下表达式：

$\min {\left|\right|z_1\left|\right|}_1+{\left|\right|z_2\left|\right|}_1$

$s.t.\left\{\begin{array}{c}\mathrm{Ax}=b\\ {\mathrm{MD}}_{i}x=z_i\end{array}\right.$

步骤5：基于稀疏优化的不完全角度扫描下的锥束CT有限角度图像重建；

将步骤4中的更新后的带约束的优化模型采用增广拉格朗日乘子法转为无约束的优化模型：

$\min \frac{1}{2}{\left|\right|\mathrm{Ax}-b\left|\right|}_2^2+\underset{i}{\mathrm{\Σ}}(u_i^T(M_i{D}_{i}x-z_i)+\frac{{\mathrm{\ρ}}_i}{2}{{\left|\right|M_i{D}_{i}x-z_i\left|\right|}_2^2+\mathrm{\λ}\left|\right|z_i\left|\right|}_1)$

对无约束的优化模型的目标函数进行变量分离，分解成两个子问题：

（1）X子问题

$\min {\left|\right|\mathrm{Ax}-b\left|\right|}_2^2+\underset{i}{\mathrm{\Σ}}\left({\mathrm{\ρ}}_i{\left|\right|M_i{D}_{i}x-z_i+u_i/{\mathrm{\ρ}}_i\left|\right|}_2^2\right)$

对上式对x求导开令为零解出x： $x^{*}={(A^{T}A+\underset{i}{\mathrm{\Σ}}\left({\mathrm{\ρ}}_i{M_{i}D}_i^T{D}_i\right))}^{+}(A^{T}b+\underset{i}{\mathrm{\Σ}}\left({\mathrm{\ρ}}_i{M}_i{D}_i^T(z_i-\frac{u_i}{{\mathrm{\ρ}}_i})\right))$

这个解形式上为解析解，实际上计算复杂度还是很高的，在实际实现时采用一步最速下降。

（2）Z子问题

$\min u_j^T(D_{j}x-z_j)+\frac{{\mathrm{\ρ}}_j}{2}{\left|\right|D_{j}x-z_j\left|\right|}^2+\mathrm{\λ}{\left|\right|z_j\left|\right|}_1$

由shrinkage算子得显式解：

$z_j^{*}=\mathrm{shrinkage}\{M_j{D}_{j}x+\frac{u_j}{{\mathrm{\ρ}}_j},\frac{\mathrm{\λ}}{{\mathrm{\ρ}}_j}\}$

解析解形式：

${z_j}^{k+1}=\max \left\{\right|M_j{D}_j{x}^{*}+\frac{u_j}{{\mathrm{\ρ}}_j}|-\frac{\mathrm{\λ}}{{\mathrm{\ρ}}_j},0\}\mathrm{sgn}(M_j{D}_{j}x+\frac{u_j}{{\mathrm{\ρ}}_j})$

综上所述，迭代步骤如下：

$\left\{\begin{array}{c}{x}^{k+1}={(A^{T}A+\underset{i}{\mathrm{\Σ}}{\mathrm{\ρ}}_i{D_i}^T{M}_i{D}_i)}^{+}(A^{T}b+\underset{i}{\mathrm{\Σ}}\left({\mathrm{\ρ}}_i{D_i}^T{M_i}^T(z_i^k-{u_i}^k/{\mathrm{\ρ}}_i)\right))\\ z_j^{k+1}=\max \left\{\right|M_j{D}_j{x}^{k+1}+\frac{{u_j}^k}{{\mathrm{\ρ}}_j}|-\frac{\mathrm{\λ}}{{\mathrm{\ρ}}_j},0\}\mathrm{sgn}[M_j{D}_j{x}^{k+1}+\frac{u_j^k}{{\mathrm{\ρ}}_j}]\\ u_j^{k+1}=u_j^k+{\mathrm{\ρ}}_j(M_j{D}_j{x}^{k+1}-z_j^{k+1})\end{array}\right.$

步骤6：对重建结果进行显示，查看重建效果是否满足要求；

图像质量的考察主要从以下方面进行：1、图像的几何伪影，散射伪影，硬化伪影是否得到有效抑制。2、图像的虚影是否影响或干扰图像细节的辨别。3、重建的分辨率水平是否满足应用需求。

利用该算法对实际数据进行重建，重建结果如图3所示。利用仿真数据测试算法的收敛性能，如图4所示。

Claims (5)

1.一种基于边缘引导的锥束CT不完全角度重建方法，其特征是：含有如下步骤：

步骤1：估计初始图像：利用扫描到的投影数据估计初始重建图像；

步骤2：图像边缘提取；

步骤3：设计加权因子；

步骤4：更新优化模型；

步骤5：基于稀疏优化锥束CT不完全角度重建；

步骤6：判断重建质量是否达到要求？如是，则执行步骤7；如不是，则执行步骤2；

步骤7：结束。

2.根据权利要求1所述的基于边缘引导的锥束CT不完全角度重建方法，其特征是：所述步骤1的具体估计方法如下：

步骤1.1：将图像重建问题刻画成如下稀疏模型：

min||x||_TV

s.t.Ax=b

其中，||x||_TV为待重建物体x的总变分范数，A为***矩阵，向量b为扫描到的投影数据；

步骤1.2：用交替方向法对步骤1.1中的稀疏模型进行求解，迭代次数为N，得到求解公式如下：

$\left\{\begin{array}{c}{x}^{k+1}={(A^{T}A+\underset{j}{\mathrm{\Σ}}{\mathrm{\ρ}}_j{D}_j^T{D}_j)}^{+}(A^{T}b+\underset{j}{\mathrm{\Σ}}{\mathrm{\ρ}}_j{D}_j^T({z_j}^k-{u_j}^k/{\mathrm{\ρ}}_j))\\ {z_j}^{k+1}=\max \left\{\right|D_j{x}^{k+1}+\frac{{u_j}^k}{{\mathrm{\ρ}}_j}|-\frac{\mathrm{\λ}}{{\mathrm{\ρ}}_j},0\}\mathrm{sgn}(D_j{x}^{k+1}+\frac{{u_j}^k}{{\mathrm{\ρ}}_j})\\ {u_j}^{k+1}={u_j}^k+{\mathrm{\ρ}}_j(D_j{x}^{k+1}-{z_j}^{k+1})\end{array}\right.$

其中，(·)⁺为矩阵的Moore-Penrose伪逆，A^T为***矩阵的转置，u_j为j方向的更新因子，u_j ^k为第k次迭代后的j方向的更新因子，u_j ^k+1为第k+1次迭代后的j方向的更新因子，ρ_j为j方向的L2范数惩罚因子，λ为L1范数惩罚因子，D_j为j方向上的梯度算子矩阵，D_j ^T为j方向上的梯度算子矩阵的转置，z_j为j方向梯度图像，z_j ^k为第k次迭代后的j方向梯度图像，z_j ^k+1为第k+1次迭代后的j方向梯度图像，max{·},sgn{·}分别为求最大函数和求符号函数，向量b为扫描到的投影数据，x^k+1为第k+1次迭代后的重建图像，经过N次迭代后的重建图像记作x^(N)，x^(N)为最后估计的初始重建图像。

3.根据权利要求2所述的基于边缘引导的锥束CT不完全角度重建方法，其特征是：所述步骤1.2中的迭代次N为总收敛轮数的0.5～0.2倍。

4.根据权利要求2所述的基于边缘引导的锥束CT不完全角度重建方法，其特征是：所述步骤2的具体方法为：设步骤1.2中某次迭代产生的中间重建图像为x1，利用降噪因子F对中间重建图像x1进行卷积运算：x1_F=x1***F，利用经典边缘算子E求取中间重建图像x1的初步边缘x1_e-mid：x1_e-mid=E[x1]，中间重建图像x1的边缘x1_e为：x1_e=x1_e-mid∩x1_LMI，其中，x1_LMI为中间重建图像x1的局部互信息图像；

步骤3的具体方法为：根据步骤2中的中间重建图像x1的边缘x1_e和L1范数惩罚因子λ确定中间重建图像x1在坐标(i1,j1)处的加权因子w(i1,j1)，该加权因子w(i1,j1)=λx1_e(i1,j1)，其中，x1_e(i1,j1)为中间重建图像x1在坐标(i1,j1)处的边缘；

步骤4的具体方法为：根据步骤3中的加权因子w(i1,j1)计算j方向的加权对角矩阵M_j，M_j=diag(w(1,1),w(1,2),...w(i1,j1),...,w(N_x,N_y))，那么优化模型更新为如下表达式：

$\min {\left|\right|z_1\left|\right|}_1+{\left|\right|z_2\left|\right|}_1+{\left|\right|z_3\left|\right|}_1$

$s.t.\left\{\begin{array}{c}\mathrm{Ax}=b\\ M_j{D}_{1}x=z_1\\ M_j{D}_{2}x=z_2\\ M_j{D}_{3}x=z_3\end{array}\right.$

其中，N_x为中间重建图像x1在x坐标上的规模，N_y为中间重建图像x1在y坐标上的规模，D₁为1方向上的梯度算子矩阵，D₂为2方向上的梯度算子矩阵，D₃为3方向上的梯度算子矩阵；

步骤5的具体方法为：采用增广拉格朗日乘子法将步骤4中的更新后的带约束的优化模型转为无约束的优化模型，具体公式如下：

$\min \frac{1}{2}{\left|\right|\mathrm{Ax}-b\left|\right|}_2^2+\underset{j}{\mathrm{\Σ}}(u_j^T(M_j{D}_{j}x-z_j)+\frac{{\mathrm{\ρ}}_j}{2}{{\left|\right|M_j{D}_{j}x-z_j\left|\right|}_2^2+\mathrm{\λ}\left|\right|z_j\left|\right|}_1)$

其中，

为j方向的更新因子的转置；

对上述无约束的优化模型采用分离变量，利用交替方向法求最小，迭代公式如下：

$\left\{\begin{array}{c}{x}^{k+1}={(A^{T}A+\underset{j}{\mathrm{\Σ}}{\mathrm{\ρ}}_j{D_j}^T{M}_j{D}_j)}^{+}(A^{T}b+\underset{j}{\mathrm{\Σ}}{\mathrm{\ρ}}_j{D_j}^T{M_j}^T(z_j^k-{u_j}^k/{\mathrm{\ρ}}_j))\\ z_j^{k+1}=\max \left\{\right|M_j{D}_j{x}^{k+1}+\frac{{u_j}^k}{{\mathrm{\ρ}}_j}|-\frac{\mathrm{\λ}}{{\mathrm{\ρ}}_j},0\}\mathrm{sgn}[M_j{D}_j{x}^{k+1}+\frac{u_j^k}{{\mathrm{\ρ}}_j}]\\ u_j^{k+1}=u_j^k+{\mathrm{\ρ}}_j(M_j{D}_j{x}^{k+1}-z_j^{k+1})\end{array}\right.$

x^k+1即为本轮重建图像；

步骤6中的重建质量达到要求是指：本轮重建图像与上一轮重建图像相比无显著变化；首次进行步骤5中的重建时，上一轮重建图像是指步骤1中的初始重建图像。

5.根据权利要求4所述的基于边缘引导的锥束CT不完全角度重建方法，其特征是：所述降噪因子F为高斯核函数或中值滤波器，经典边缘算子E为canny边缘算子、soble边缘算子、prewitt边缘算子、roberts边缘算子、log边缘算子、zeroscross边缘算子中的任一种，L1范数惩罚因子λ的取值在0～1之间。

Images