基于改进广义全变分的稀疏图像重建算法

班晓征; 李志华; 李贝贝; 徐敏达

doi:doi:10.3788/LOP55.111103

激光与光电子学进展, 2018, 55 (11): 111103, 网络出版: 2019-08-14

基于改进广义全变分的稀疏图像重建算法下载： 1460次

Sparse Image Reconstruction Based on Improved Total Generalized Variation

论文大纲

班晓征李志华 ^*李贝贝徐敏达

作者单位

江南大学物联网工程学院, 江苏无锡 214122

AI 词云图 AI语音精读 AI语音超短摘要

注：本部分内容由 AI 自动生成，请您知悉。

摘要

为了提升稀疏采样环境下的图像重建质量,针对广义全变分模型重建图像时不能充分利用图像本身结构自相似性信息的不足,建立了一个非局部约束下的改进广义全变分图像重建模型。该模型引入了变化域非局部自相似性作为图像重建的先验信息,同时在八邻域空间计算多方向的广义全变分正则化约束,从而更好地保护了图像的结构特征,进一步地,使用增广拉格朗日理论对模型进行去约束化、求解,提出了基于改进广义全变分的图像重建算法。仿真实验结果表明,所提出的重建模型和图像重建算法可以有效地去除图像中的伪影和噪声,满足稀疏采样情形下对图像重建质量的要求。与其他重建算法进行比较可知,本文算法所重建的图像不论是主观视觉效果,还是各个客观评价指标均有不同程度的改善和提高。

Abstract

Regarding to the total generalized variational model cannot fully utilize the self-similarity information of the image structure when reconstructing images, an improved generalized variational image reconstruction model under non-local constraints is established to improve the quality of image reconstruction in the sparse sampling situation. This model introduces a non-local self-similarity of the transform domain as a priori information for image reconstruction. And the multi-directional total generalized variational regularization constraint is calculated in the eight-neighborhood space to protect the structural characteristics of the image. Further, the augmented Lagrangian theory is used to remove the constraint and solve the model, and an image reconstruction algorithm based on the improved total generalized variation is proposed. Simulation experimental results show that the proposed reconstruction model and image reconstruction algorithm can effectively remove the artifacts and noise in the image and meet the requirements of image reconstruction quality under sparse sampling condition. Compared with the famous reconstruction algorithms, the images reconstructed by proposed algorithm has significant improvement in both subjective visual effects and all objective evaluation indicators.

1 引言

图像重建算法主要是通过类似反演计算等手段将探测器捕获到的投影信息恢复成原始图像中的像素^[1]。稀疏采样是指在采集投影数据时受到空间或时间限制,所获取的投影数据完备性不足。稀疏采样下图像重建主要有解析类图像重建算法^[2-3]和统计迭代类图像重建算法^[4-9]。改进解析类图像重建算法通过提升投影数据的质量,进而提升重建图像的质量,但是此类算法对投影数据的完备性要求高,稀疏采样情况下极易丢失图像细节。统计迭代类重建算法向传统迭代重建模型中添加优化约束限制,利用投影数据的极大似然函数特征对求解目标进行最大近似处理,文献[ 4-5]引入纹理惩罚模型作为先验知识,并提出通过贝叶斯理论推理出最大后验图像作为输出的重建图像,使重建图像边缘效果得到显著的增益。统计迭代方法充分考虑了图像中的噪声情况并建立了数学优化模型,所以得到了许多研究的跟踪结果^[6-10]。由于在稀疏采样环境下进行图像重建可以等同看作是欠采样数据恢复问题,文献[ 6]根据压缩感知理论,提出了最小化图像总变差迭代(TV-ART)算法,在每次代数重建后对图像的梯度图像进行了最小化约束,抑制了图像中高频噪声。在TV-ART算法基础上,文献[ 7]提出了一种结合先验图像的压缩感知(PICCS)图像重建算法,PICCS算法充分利用了图像自身的特征,将自身作为下次迭代的先验知识,使得重建的图像具有良好的结构性。但是ART算法迭代求解速度慢,通过迭代并不能获得最优解。为了更加快速地对问题进行求解,文献[ 8]基于TV最小化模型,使用Split-Bregman算法缩短了重建时间。

另一方面,针对TV模型容易产生阶梯效应的问题,文献[ 9]在TV重建模型中引入了非局部平均滤波修正思想,每次迭代时使用改进的非局部平均滤波修正当前重建的图像。该算法由于使用了非局部约束的思想,重建图像的边缘结构保持的更好,并且信噪比也更高。文献[ 10]将广义全变分(TGV)正则化应用于求解三维图像重建问题,有效避免了TV重建模型中的阶梯效应,保持了重建图像边缘的完整性,但没有很好地区分图像的不同结构,在图像光滑部分出现噪点。

针对目前稀疏采样情况下的图像重建过程中容易出现系统噪声和虚假边缘等问题,本文提出了一种结合图像非局部信息的多方向广义全变分图像重建算法。该算法首先通过非局部自相似性原理计算先验图像,并将其作为图像重建问题的先验知识,再使用增广拉格朗日方法对模型进行快速的迭代求解。实验结果表明,本文算法不仅可以保证重建图像具有良好的结构性质,还可以有效去除噪声。

2 广义全变分重建模型

稀疏采样情况下的图像重建问题可以建立模型,表示为

\begin{matrix} \min_{X} ({(AX - b)}_{2}^{2} + {(ψX)}_{0}), (1) \end{matrix}

式中第一项为求解目标项,第二项为添加的惩罚约束项,ψ表示稀疏变换。通过进行稀疏变换可以捕获图像中的特征,并将其作为图像重建的先验信息。

TGV模型由Bredies等^[11]提出,借助高阶多项式进行了最小化逼近。对于任意的X∈ $\begin{matrix} L_{loc}^{1} \end{matrix}$ (Ω),其k阶TGV定义为

\begin{matrix} TG V_{α}^{k} (X) = \sup \{\begin{matrix} \int X di v^{k} wdx | w \in C_{c}^{k} [Ω, S_{Sym}^{k} (R^{d})], \\ {(di v^{l} w)}_{\infty} \leq α_{l}, l = 0, \dots, k - 1 \end{matrix}\}, (2) \end{matrix}

式中Ω∈R^d并且是一个非空的开区间, $\begin{matrix} S_{Sym}^{k} \end{matrix}$ (R^d)是k阶对称张量空间,k≥1且α=(α₀,…,α_k-₁)是加权系数, $\begin{matrix} C_{c}^{k} \end{matrix} \begin{matrix} [Ω, {S^{k}}_{Sym} (R^{d})] \end{matrix}$ 是紧支撑对称向量空间,div表示散度算子。

当k=2时,根据Lengendre-Fenchel对偶定理和离散变换,可以得到离散二阶TGV模型,即:

\begin{matrix} TGV (X) = \arg \min_{p} \{α_{1} \overset{2}{\sum_{j = 1}} {(\nabla X_{j} - p_{j})}_{1} + α_{0} \overset{4}{\sum_{j = 1}} {({[ε (p)]}_{j})}_{1}\}, (3) \end{matrix}

式中ε(p)=Ñp+Ñp^T/2为对称梯度算子,α₁与α₀为梯度权重,用于平衡一阶差分算子与二阶差分算子。

二阶TGV模型包含的低阶差分算子有助于去除图像中的噪声,同时高阶差分算子有助于保持图像的边缘信息。为此,将(3)式所示的二阶TGV范数作为惩罚项代入到(1)式中,可以得到广义全变分图像重建模型,即:

\begin{matrix} \underset{X}{argmin} [{(AX - b)}_{2}^{2} + TGV (X)] 。 (4) \end{matrix}

3 改进的广义全变分图像重建模型

3.1 多方向TGV模型

传统TGV模型中仅在水平和垂直两个方向上计算梯度信息,但是由于图像中的边缘信息是多方向的,因此本文在图1所示的像素的八邻域空间中计算一阶梯度算子,不仅可以检测到水平和垂直方向上的边缘,还可以考虑到如图1所示的对角线上的45°、135°等8个方向的边缘特征,有利于尽可能多地保留图像中的细节。

图 1. 像素点的八邻域空间

Fig. 1. Eight neighborhoods of pixel points

下载图片查看所有图片

对于图像X,其在八邻域空间中计算得到的多方向广义全变分约束(MTGV)定义为

\begin{matrix} MTGV (X) = \min_{p} \{α_{1} \overset{4}{\sum_{j = 1}} {(\nabla X_{j} - p_{j})}_{1} + α_{0} \overset{16}{\sum_{j = 1}} {({[ε (p)]}_{j})}_{1}\}, (5) \end{matrix}

式中ÑX₁与ÑX₂分别表示水平与垂直方向的一阶梯度,ÑX₃与ÑX₄分别表示45°和135°对角线方向的一阶梯度信息。

3.2 非局部先验约束

虽然MTGV模型可以更多地探测到图像中边缘信息,但是其仍然仅考虑了图像的局部信息,忽视了图像中的非局部结构信息。文献[ 12]利用图像的非局部自相似性性质提出了著名的三维块匹配(BM3D)算法,在变换域中对图像进行非局部约束,有效地去除了图像中的高频噪声。考虑到BM3D算法中第二次协同滤波处理对图像信噪比的提升是有限的,因此对BM3D算法进行改进,仅对图像做一次协同滤波处理,提出了利用图像块相似性的非局部滤波算法,称之块组相似性(PGS)方法,PGS方法仅对图像进行一次块匹配和协同硬阈值滤波处理,这样做的好处是显著缩短了非局部先验的计算时间。PGS方法基本步骤概括如下:

1) 将大小为N×N的图像分为j个等大的n×n图像块X_j(j=1,2,…,N/n)。

2) 对于每个图像块X_j,在其周围大小为Ns×Ns的搜索框中,采用欧氏距离作为相似性度量,聚合与其最为相似的前C个图像块,并存储为三维的图像相似块矩阵S_Cj。

3) 对步骤2)获得的相似块矩阵S_Cj进行三维正交变换,随后在变换域上使用协同的硬阈值方式进行滤波处理。

4) 对处理过的相似块矩阵S_Cj进行三维反变换,并通过权重平均的方式更新每个图像块在原图像的像素值。

PGS方法保护了图像结构和纹理细节,因此使用PGS方法对先验图像进行非局部约束,同时通过减法操作对优化目标进行稀疏化处理,得到非局部先验约束为

\begin{matrix} ϕ_{PGS} (X) = {(X - PGS (X_{up}))}_{1}, (6) \end{matrix}

式中X为待求解目标,X_up代表先验图像,在本文中为上一步迭代出的重建图像。

3.3 非局部信息约束下的多方向广义全变分模型

进一步地,使用3.1节提出的MTGV约束代替(4)式中的TGV约束项,并加入3.2节所提出PGS非局部先验约束方法,得到非局部信息约束下的广义全变分图像重建模型(PGS-MTGV),表示为

\begin{matrix} \begin{matrix} X = \underset{X}{argmin} \{\frac{β}{2} {(AX - b)}_{2}^{2} + λ_{PGS} {(X - PGS (X_{up}))}_{1} + \\ λ_{MTGV} \{α_{1} \overset{4}{\sum_{j = 1}} {(\nabla X_{j} - p_{j})}_{1} + α_{0} \overset{16}{\sum_{j = 1}} {({[ε (p)]}_{j})}_{1}\}\}, (7) \end{matrix} \end{matrix}

式中第一项为求解的目标项,第二项是非局部结构先验项,第三项是多方向广义全变分约束项。λ_PGS代表非局部先验项的权重,λ_MTGV代表多方向广义全变分约束项的权重,并且两个权重满足λ_PGS+λ_MTGV=1,通过调节权重可以平衡两个不同的约束对图像重建效果的影响。

4 基于增广拉格朗日算法的模型求解

通过引入辅助变量y、r_j(j=1,…,4)、s_j(j=1,…,16)分别表示X-PGS(X_up)、 $\begin{matrix} \nabla_{j} \end{matrix}$ X-p_j、[ε(p)]_j,公式(7)可以等价转换成

\begin{matrix} \begin{matrix} X = \underset{X}{argmin} [\frac{β}{2} {(AX - b)}_{2}^{2} + λ_{PGS} {(y)}_{1} + λ_{MTGV} (α_{1} \overset{4}{\sum_{j = 1}} {(r_{j})}_{1} + α_{0} \overset{16}{\sum_{j = 1}} {(s_{j})}_{1})] \\ s . t . y = X - PGS (X_{up}), r_{j} = \nabla_{j} X - p_{j}, s_{j} = [ε {(p)]}_{j} 。 (8) \end{matrix} \end{matrix}

使用中间变量V代替非局部先验图像PGS(X_up),并且添加增广拉格朗日乘子yy、rr以及ss后,根据增广拉格朗日定理,(8)式可以进一步转化成无约束优化问题,即

\begin{matrix} \begin{matrix} X = \underset{y, rj, sj, X, pj}{argmin} \{\frac{β}{2} {(AX - b)}_{2}^{2} + λ_{PGS} {(y)}_{1} + \\ λ_{MTGV} (α_{1} \overset{4}{\sum_{j = 1}} {(r_{j})}_{1} + α_{0} \overset{16}{\sum_{j = 1}} {(s_{j})}_{1}) + \frac{λ_{PGS} μ_{1}}{2} {(y - (X - V) - yy)}_{2}^{2} + \\ \frac{λ_{MTGV} α_{1} μ_{2}}{2} \overset{4}{\sum_{j = 1}} {(r_{j} - (\nabla_{j} X - p_{j}) - r r_{j})}_{2}^{2} + \frac{λ_{MTGV} α_{0} μ_{3}}{2} \overset{16}{\sum_{j = 1}} {(s_{j} - [ε {(p)]}_{j} - s s_{j})}_{2}^{2}\} 。 (9) \end{matrix} \end{matrix}

(9)式可以看作由三个极小化子问题组成,问题1是对辅助变量y、r_j、s_j的更新,问题2是求解待重建图像X,问题3是求解MTGV约束中的梯度信息p_j。每次迭代时,首先需要对三个问题依次进行更新,然后对拉格朗日乘子进行更新,更新方法为

\begin{matrix} \{\begin{matrix} yy = yy + γ (X - V - y) \\ rr = rr + γ (\nabla X - p - r) \\ ss = ss + γ [ε (p) - s] \end{matrix} 。 (10) \end{matrix}

4.1 辅助变量的更新

从(9)式中提取辅助变量y、r_j、s_j的相关表达式,得到

\begin{matrix} \{\begin{matrix} y = \underset{y}{argmin} [{(y)}_{1} + \frac{μ_{1}}{2} {(y - (X - V) - yy)}_{2}^{2}] \\ r_{j} = \underset{r_{j}}{argmin} [{(r_{j})}_{1} + \frac{μ_{2}}{2} {(r_{j} - (\nabla_{j} X - p_{j}) - r r_{j})}_{2}^{2}], j = 1,2, 3,4 \\ s_{j} = \underset{s_{j}}{argmin} [{(s_{j})}_{1} + \frac{μ_{3}}{2} {(s_{j} - [ε {(p)]}_{j} - s s_{j})}_{2}^{2}], j = 1, \dots, 16 \end{matrix} 。 (11) \end{matrix}

(11)式可以看作是对辅助变量y、r_j、s_j的L1范数优化问题,具有shrinkage解析解,所以可以直接采用软阈值收缩函数来更新,具体为

\begin{matrix} \{\begin{matrix} y = sgn [(X - V) - yy] \max \{|(X - V) - yy| - \frac{1}{μ_{1}}, 0\} \\ r_{j} = sgn [(\nabla_{j} X - p_{j}) - r r_{j}] \max \{|(\nabla_{j} X - p_{j}) - r r_{j}| - \frac{1}{μ_{2}}, 0\} \\ s_{j} = sgn \{\{{[ε (p)]}_{j} - s s_{j}\} - yy\} \max \{|[ε {(p)]}_{j} - s s_{j}| - \frac{1}{μ_{3}}, 0\} \end{matrix} 。 (12) \end{matrix}

4.2 X方向子问题

从(9)式中提取变量X的有关项,并且使用循环差分算子D_j近似替代梯度算子 $\begin{matrix} \nabla_{j} \end{matrix}$ ,继而X方向的子问题可以表示为

\begin{matrix} \begin{matrix} X = \underset{X}{argmin} [\frac{β}{2} {(AX - b)}_{2}^{2} + \frac{λ_{PGS} μ_{1}}{2} {(y - (X - V) - yy)}_{2}^{2} + \\ \frac{λ_{MTGV} α_{1} μ_{2}}{2} \overset{4}{\sum_{j = 1}} {(r_{j} - (D_{j} X - p_{j}) - r r_{j})}_{2}^{2}] 。 (13) \end{matrix} \end{matrix}

易见,(13)式对变量X是一个带有罚函数的最小二乘问题,可以通过求偏导得到X的闭合解。因为投影矩阵A是一个大型稀疏矩阵,对矩阵A计算会耗费大量的内存和时间,为此利用非精确近似思想^[13],使用ρ=A^T(AX_k-b)表示X_k时的梯度,τ满足条件τ>0,得到 $\begin{matrix} {(AX - b)}_{2}^{2} \end{matrix}$ 的近似替代式ϕ_similar= $\begin{matrix} {(A X_{k} - b)}_{2}^{2} \end{matrix}$ +2ρ^T $\begin{matrix} (X - X_{k}) \end{matrix}$ +1/τ $\begin{matrix} {(X - X_{k})}_{2}^{2} \end{matrix}$ 。将ϕ_similar代入到(13)式中,再对(13)式中的X求偏导,从而可以得到

\begin{matrix} \begin{matrix} X = {(λ_{PGS} μ_{1} + λ_{MTGV} α_{1} μ_{2} \overset{4}{\sum_{j = 1}} {D^{T}}_{j} D + \frac{β}{τ} I)}^{+} \cdot \\ [λ_{PGS} μ_{1} (y + V - yy) + λ_{MTGV} α_{1} μ_{2} \overset{4}{\sum_{j = 1}} {D^{T}}_{j} (r_{j} + p_{j} - r r_{j}) + \frac{β}{τ} X - βρ] 。 (14) \end{matrix} \end{matrix}

通过快速傅里叶变换简化(14)式的求逆运算,得到求解X的傅里叶解析表达式

\begin{matrix} X = F^{- 1} \{\frac{F [λ_{PGS} μ_{1} (y + V - yy) + λ_{MTGV} α_{1} μ_{2} \overset{4}{\sum_{j = 1}} {D^{T}}_{j} (r_{j} + p_{j} - r r_{j}) + \frac{β}{τ} X - βρ]}{diag [F (λ_{PGS} μ_{1} + λ_{MTGV} α_{1} μ_{2} \overset{4}{\sum_{j = 1}} {D^{T}}_{j} D + \frac{β}{τ} I) F^{- 1}]}\}, (15) \end{matrix}

式中F代表傅里叶变换,diag表示取对角线元素,F^-1代表傅里叶逆变换。当获得重建图像X后,使用非负约束(POS)方法^[14]保证重建图像的非负性,具体为

\begin{matrix} POS (X_{i, j}) = \{\begin{matrix} X_{i, j}, & X_{i, j} > 0 \\ 0, & X_{i, j} < 0 \end{matrix} 。 (16) \end{matrix}

4.3 p1~p4方向子问题

从(9)式中提取p_j的相关项,并使用循环差分算子D_j近似替代梯度算子 $\begin{matrix} \nabla_{j} \end{matrix}$ 后,得到求解p_j的最优化表达式为

\begin{matrix} p = \underset{p_{j}}{argmin} \{\frac{λ_{MTGV} α_{1} μ_{2}}{2} \overset{4}{\sum_{j = 1}} {(r_{j} - (D_{j} X - p_{j}) - r r_{j})}_{2}^{2} + \frac{λ_{MTGV} α_{0} μ_{3}}{2} \overset{16}{\sum_{j = 1}} {(s_{j} - [ε {(p)]}_{j} - s s_{j})}_{2}^{2}\} 。 (17) \end{matrix}

对(17)式求p₁方向的偏导,得到p₁的求解式为

\begin{matrix} \begin{matrix} p_{1} = {(α_{1} μ_{2} \cdot I + α_{0} μ_{3} {D^{T}}_{1} D_{1} + \frac{α_{0} μ_{3}}{2} {D^{T}}_{2} D_{2} + \frac{α_{0} μ_{3}}{2} {D^{T}}_{3} D_{3} + \frac{α_{0} μ_{3}}{2} {D^{T}}_{4} D_{4})}^{+} \\ [α_{1} μ_{2} (D_{1} X - r_{1} + r r_{1}) + α_{0} μ_{3} (s_{1} - s s_{1}) - \frac{α_{0} μ_{3}}{2} {D^{T}}_{2} (D_{1} p_{2} - 2 s_{5} + 2 s s_{5}) - \\ \frac{α_{0} μ_{3}}{2} {D^{T}}_{3} (D_{1} p_{3} - 2 s_{8} + 2 s s_{8}) - \frac{α_{0} μ_{3}}{2} {D^{T}}_{4} (D_{1} p_{4} - 2 s_{10} + 2 s s_{10})] 。 (18) \end{matrix} \end{matrix}

由于(18)式中没有复杂矩阵,可以直接使用快速傅里叶变换将(18)式等价变化为

\begin{matrix} p_{1} = F^{- 1} [\frac{F (R_{1})}{α_{1} μ_{2} + α_{0} μ_{3} diag (F {D^{T}}_{1} D_{1} F^{- 1} + F \frac{{D^{T}}_{2} D_{2}}{2} F^{- 1} + F \frac{{D^{T}}_{3} D_{3}}{2} F^{- 1} + F \frac{{D^{T}}_{4} D_{4}}{2} F^{- 1})}], (19) \end{matrix}

式中R₁=α₁μ₂ $\begin{matrix} (D_{1} X - r_{1} + r r_{1}) \end{matrix}$ +α₀μ₃ $\begin{matrix} (s_{1} - s s_{1}) \end{matrix}$ -α₀μ₃/2 $\begin{matrix} D_{2}^{T} \end{matrix} \begin{matrix} (D_{1} p_{2} - 2 s_{5} + 2 s s_{5}) \end{matrix}$ -α₀μ₃/2 $\begin{matrix} D_{3}^{T} \end{matrix} \begin{matrix} (D_{1} p_{3} - 2 s_{8} + 2 s s_{8}) \end{matrix}$ -α₀μ₃/2 $\begin{matrix} D_{4}^{T} \end{matrix} \begin{matrix} (D_{1} p_{4} - 2 s_{10} + 2 s s_{10}) \end{matrix}$ ,通过上述方法,可以类似地对p₂~p₄进行求解。

4.4 基于改进广义全变分模型的图像重建算法

概括上述求解、推导过程,本文提出基于改进广义全变分的稀疏图像重建(PGS-MTGV)算法。

Name:PGS-MTGV algorithm;
Input:projection matrix A、projection data b;
Output:reconstructed image X。

1. Initialize reconstruction parameters:λ_DTGV,λ_MTGV,α₀,α₁,μ₁,μ₂,μ₃,γ;
2. Initialize X_up=0 as the iterated initial prior image;
3. While (not convergence & i<I_iter)do
4. X_up=X;
5. Calculate non-local constraints on prior image by PGS method;
6. Update y^k+¹、

\begin{matrix} r_{j}^{k + 1} \end{matrix}

、

\begin{matrix} s_{j}^{k + 1} \end{matrix}

by formula (12);
7. Update X^k⁺¹ by formula (15) and formula (16);
8. Update

\begin{matrix} p_{j}^{k + 1} \end{matrix}

by formula (17);
9. Update yy^k+¹、r

\begin{matrix} r_{j}^{k + 1} \end{matrix}

、s

\begin{matrix} s_{j}^{k + 1} \end{matrix}

by formula (10);
10. i=i+1;
11. End;
12. Return X;

5 实验结果和分析

5.1 实验环境及参数设定

本文实验投影光线模拟为平行光光束,并且为了对所提出的稀疏图像重建算法重建效果的高效性和有效性进行验证,在有限投影角度和稀疏投影角度两种采样环境下,使用图2所示的Sheep-Logan、Forbild-Head、Forbild-Abdomen三个大小为256 mm×256 mm的图像进行实验。文献[ 15]通过研究发现,当采样不足时,实际投影过程会产生近似服从非平稳高斯分布的噪声,为了更好地模拟实际投影重建情况,实验中向投影数据添加了10%高斯随机噪声。

图 2. 测试用的图像数据。(a) Sheep-Logan 图像;(b) Forbild-Head 图像;(c) Forbild-Abdomen 图像

Fig. 2. Image data used in the test. (a) Sheep-Logan image; (b) Forbild-Head image; (c) Forbild-Abdomen image

下载图片查看所有图片

选取迭代次数为500,将本文算法与L1-Bregman算法[16]、交替方向全变分最小化算法(TV-ADM)TV-ADM算法[17]、基于拉格朗日的全变分最小化算法(TVAL3)TVAL3算法[18]以及交替方向广义全变分最小化算法(TGV-ADM)TGV-ADM算法[10]进行比较,初始化实验参数设置如表1所示。实验配置为64位的windows10,8 GB内存,双核CPU,500 GB硬盘,matlabR2016a作为仿真开发工具。

表 1. 初始化参数设置

Table 1. Initialization parameters

Parameter	β	μ_i(i=1,2,3)	γ	α₀	α₁
Value	64	32	1	0.2	2

查看所有表

由于随着迭代次数的增加,非局部约束对于图像重建效果的影响逐渐减小,所以适当的降低非局部约束的权重,增加MTGV正则项所占的权重。设置λ_MTGV=1/1+exp[-0.001(i+1)],其中i代表迭代的次数。

5.2 评价指标

为了对图像重建效果进行公平地评估,本文采用的定量评价指标为结构相似性(SSIM)、峰值信噪比(PSNR)和均方根误差(RMSE)。SSIM指标可以识别出图像中的结构信息,考虑到人眼的视觉特性对两幅图像在亮度和对比度上进行对比,SSIM越接近1表明对比的图像在视觉上越相似。PSNR和RMSE的定义分别为

\begin{matrix} \begin{matrix} R_{psnr} (X, X_{u}) = 10 \lg \frac{255^{2} \times M \times N}{{(X - X_{u})}_{2}^{2}}, (20) \\ R_{psnr} (X, X_{u}) = \sqrt[]{\frac{\overset{M \times N}{\sum_{j = 0}} (X_{j} - X_{uj})}{M \times N}}, (21) \end{matrix} \end{matrix}

式中X_u表示原始图像,X表示重建后的图像,M和N是图像的长和宽,M×N为图像中像素单元的总数,X_j和X_uj表示图像中的第j个像素单元。PSNR反映了对比的图像噪声水平的差异,并且值越大表明图像重建的质量越好,而RMSE的值越接近0表明对比的图像在像素级别越相似。

另外,为了测试图像的细节区域,分别选取了两个感兴趣区(ROI),如图3所示。通过对感兴趣区进行放大可以更清晰地观察到图像的细节,从视觉上对各算法重建图像的优劣进行比较,有助于准确分析各重建算法的有效性。

图 3. 测试图像设置的ROI区。(a) Sheep-Logan图像;(b) Forbild-Head图像;(c) Forbild-Abdomen图像

Fig. 3. ROI regions of the test images. (a) Sheep-Logan image; (b) Forbild-Head image; (c) Forbild-Abdomen image

下载图片查看所有图片

5.3 算法有效性分析

首先设置采样空间为135°,一共进行15次采样,模拟有限投影角度采样情形进行图像重建实验。图4~6是在有限角度采样情形下各算法重建的图像及图像中两个ROI细节的放大图像。

图 4. 有限角度各算法重建图像及ROI对比(Sheep-Logan)。(a)原始图像;(b) L1-Bregman 算法;(c) TV-ADM 算法;(d) TVAL3 算法;(e) TGV-ADM 算法;(f)本文算法

Fig. 4. Reconstructed image and ROI comparison of different algorithms in finite angle (Sheep-Logan). (a) Original image; (b) L1-Bregman algorithm; (c) TV-ADM algorithm; (d) TVAL3 algorithm; (e) TGV-ADM algorithm; (f) proposed algorithm

下载图片查看所有图片

图 5. 有限角度各算法重建图像及ROI对比(Forbild-Head)。 (a)原始图像;(b) L1-Bregman 算法;(c) TV-ADM 算法;(d) TVAL3 算法;(e) TGV-ADM 算法;(f)本文算法

Fig. 5. Reconstructed image and ROI comparison of different algorithms in finite angle (Forbild-Head). (a) Original image; (b) L1-Bregman algorithm; (c) TV-ADM algorithm; (d) TVAL3 algorithm; (e) TGV-ADM algorithm; (f) proposed algorithm

下载图片查看所有图片

图 6. 有限角度各算法重建图像及ROI对比(Forbild-Abdomen)。 (a)原始图像;(b) L1-Bregman 算法;(c) TV-ADM 算法;(d) TVAL3 算法;(e) TGV-ADM 算法;(f)本文算法

Fig. 6. Reconstructed image and ROI comparison of different algorithms in finite angle (Forbild-Abdomen). (a) Original image; (b) L1-Bregman algorithm; (c) TV-ADM algorithm; (d) TVAL3 algorithm; (e) TGV-ADM algorithm; (f) proposed algorithm

下载图片查看所有图片

观察图4~6可以发现,由于缺乏部分角度的投影信息,各算法重建的图像出现了不同程度的信息缺失效果。其中,L1-Bregman算法重建的图像存在严重的颗粒感,这是由于L1范数空间解的不连续性造成的,其余几种算法通过对图像进行稀疏化处理获得了更多的先验约束,防止了过拟合情况的出现。采用TV正则化模型的两种算法对噪声有很好的平滑效果,其中TV-ADM算法重建的图像虽然没有噪声但是在图像中的细节部分出现了明显的阶梯效应;TVAL3算法重建出来的图像质量相较TV-ADM算法有所提升,这是由于TVAL3算法改进了迭代算法中参数选取的策略,重建出的图像轮廓更加完整,但是TVAL3算法仍然不能很好地区分图像的光滑区域和边缘区域。TGV-ADM算法恢复出了较多的图像细节,但从视觉上看本文提出的算法更能区分图像中的光滑部分和细节部分,比如图6(f)中伪影较少,且边缘轮廓基本逼近了原图。

为了进一步地分析本文算法的可用性,下面模拟稀疏角度采样环境,设置总的采样空间为180°,同样地设置15个采样点。图7~9显示了本文算法与其他几种对比算法的重建图像和ROI区域差异。

图 7. 稀疏角度各算法重建图像及ROI对比(Sheep-Logan)。 (a)原始图像;(b) L1-Bregman 算法;(c) TV-ADM 算法;(d) TVAL3 算法;(e) TGV-ADM算法;(f)本文算法

Fig. 7. Reconstructed image and ROI comparison of different algorithms in sparse angle (Sheep-Logan). (a) Original image; (b) L1-Bregman algorithm; (c) TV-ADM algorithm; (d) TVAL3 algorithm; (e) TGV-ADM algorithm; (f) proposed algorithm

下载图片查看所有图片

图 8. 稀疏角度各算法重建图像及ROI对比(Forbild-Head)。 (a)原始图像;(b) L1-Bregman 算法;(c) TV-ADM 算法;(d) TVAL3 算法;(e) TGV-ADM算法;(f)本文算法

Fig. 8. Reconstructed image and ROI comparison of different algorithms in sparse angle (Forbild-Head). (a) Original image; (b) L1-Bregman algorithm; (c) TV-ADM algorithm; (d) TVAL3 algorithm; (e) TGV-ADM algorithm; (f) proposed algorithm

下载图片查看所有图片

图 9. 稀疏角度各算法重建图像及ROI对比(Forbild-Abdomen)。 (a)原始图像;(b) L1-Bregman 算法;(c) TV-ADM 算法;(d) TVAL3 算法;(e) TGV-ADM算法;(f)本文算法

Fig. 9. Reconstructed image and ROI comparison of different algorithms in sparse angle (Forbild-Abdomen). (a) Original image; (b) L1-Bregman algorithm; (c) TV-ADM algorithm; (d) TVAL3 algorithm; (e) TGV-ADM algorithm; (f) proposed algorithm

下载图片查看所有图片

观察图7~9可以看出,从主观视觉上看,在稀疏角度采样环境中这些图像重建算法都可以完整地重建出原图像的总体轮廓。这是由于稀疏角度采样空间大于有限角度,可以比较完整地获得图像中每个像素的信息,并且像素越靠近图像中间区域,重建信息获取的越充分。L1-Bregman算法的重建效果要优于有限角度下,但是图像中仍然存在很明显的颗粒感,模糊掉了图像中的细节。虽然TVAL3算法相比TV-ADM算法重建出了更多图像中的细节,但两种算法重建的图像边缘出现了很明显的梯度纹路。相比而言,TGV模型重建细节的能力比TV模型更加出色,例如对于Forbild-Abdomen图像,TGV-ADM算法重建出了最多的点状圆形单元,如图9(e)所示。对比ROI放大图可以看出,这些圆形单元的边缘并不清晰,仅图9(f)更接近原图像。综合对比各算法重建图像和ROI放大图可知,本文算法重建图像的边缘阶梯效应更弱,重建效果更好。

总的来说,不管是在有限角度采样情形还是稀疏角度采样情形,本文算法重建的图像光滑度更好,细节上与原图更加接近。

5.4 算法高效性分析

为了表现本文算法重建图像的高效性,利用5.2节各评价指标对几种算法重建的图像进行客观评价。表2和表3分别为在有限投影角度下和稀疏投影角度下分别采用L1-Bregman算法、TV-ADM算法、TVAL3算法、TGV-ADM算法以及所提出的PGS-MTGV算法重建出的图像相对原始图像的SSIM、PSNR、RMSE的客观质量参数和执行时间对比结果。本文算法的约束项比对比算法要多,尤其非局部先验约束需要较长的计算时间,导致算法在执行时间上并不占优势,但各个客观质量参数具有不同程度的优势。

表 2. 有限投影角度下各算法重建图像的客观评价参数对比

Table 2. Comparison of objective evaluation parameters of reconstructed images of various algorithms under finite projection angle

Algorithm	Image	SSIM	PSNR /dB	RMSE	Time /s
L1-Bregman	Sheep-Logan	0.549	15.748	0.163	10.969
	Forbild-Head	0.430	13.175	0.219	11.063
	Forbild-Abdomen	0.535	15.672	0.165	11.328
TV-ADM	Sheep-Logan	0.891	22.203	0.078	16.467
	Forbild-Head	0.893	21.916	0.080	16.453
	Forbild-Abdomen	0.967	33.024	0.022	15.781
TVAL3	Sheep-Logan	0.894	22.807	0.073	19.625
	Forbild-Head	0.899	23.072	0.070	20.109
	Forbild-Abdomen	0.968	33.643	0.021	22.234
TGV-ADM	Sheep-Logan	0.899	25.039	0.056	70.016
	Forbild-Head	0.935	27.420	0.043	69.594
	Forbild-Abdomen	0.967	35.312	0.018	68.609
Proposed	Sheep-Logan	0.903	26.091	0.054	150.682
	Forbild-Head	0.954	29.229	0.035	149.531
	Forbild-Abdomen	0.977	39.231	0.011	147.000

查看所有表

表2中本文算法重建图像的SSIM值相对于其他算法更接近1,并且RMSE值更接近0,这也反映了在有限角度采样数据情形下,本文算法重建出的图像与原始图像更加接近。另外,本文算法的PSNR值也更大,这说明本文算法对图像重建过程中的噪声具有很好的抑制效果。

从表3中可知,稀疏采样角度下各算法重建图像的SSIM、PSNR、RMSE各评价指标均优于有限角度,这是由于稀疏采样的采样范围更广,投影数据更完备。在稀疏角度采样情况下,本文算法的SSIM、PSNR与RMSE指标仍有不同程度的优势,这反映了本文算法的重建图更接近原图像。

表 3. 稀疏投影角度下各算法重建图像的客观评价参数对比

Table 3. Comparison of objective evaluation parameters of reconstructed images of various algorithms under sparse projection angle

Algorithm	Image	SSIM	PSNR /dB	RMSE	Time /s
L1-Bregman	Sheep-Logan	0.619	17.449	0.134	11.625
	Forbild-Head	0.185	14.441	0.190	11.078
	Forbild-Abdomen	0.543	16.340	0.152	12.953
TV-ADM	Sheep-Logan	0.995	45.429	0.005	16.797
	Forbild-Head	0.993	39.136	0.011	16.531
	Forbild-Abdomen	0.981	36.339	0.015	17.250
TVAL3	Sheep-Logan	0.995	46.075	0.005	20.953
	Forbild-Head	0.995	43.074	0.007	20.234
	Forbild-Abdomen	0.985	38.522	0.012	22.000
TGV-ADM	Sheep-Logan	0.991	44.582	0.006	69.031
	Forbild-Head	0.995	43.810	0.006	67.406
	Forbild-Abdomen	0.978	39.608	0.011	71.922
Proposed	Sheep-Logan	0.999	50.018	0.003	144.250
	Forbild-Head	0.999	52.172	0.003	146.359
	Forbild-Abdomen	0.991	39.930	0.010	149.141

查看所有表

为了更进一步地分析稀疏采样角度下各算法的重建效果,图10选取各算法重建图像与原图像的第135列的纵向剖面图进行对比。从图中可以看出,本文算法重建图像剖面线与原图像有些许偏差,但像素值的波动很小,这说明本文算法具有比较好的细节保留能力与抗噪能力。

图 10. 各算法重建图像第135列纵向剖面图对比。(a) Sheep-Logan 图像;(b) Forbild-Head 图像;(c) Forbild-Abdomen 图像

Fig. 10. Comparison of vertical cross-sectional views of 135th column of each reconstruction image. (a) Sheep-Logan image; (b) Forbild-Head image; (c) Forbild-Abdomen image

下载图片查看所有图片

6 结论

提出了基于多方向广义全变分和非局部先验的PGS-MTGV图像重建算法,该算法充分利用了图像非局部自相似性的特点,从而有效地提升重建图像的结构特性。与现有的4种经典算法进行实验对比,实验结果表明,在有限角度采样情况下,本文算法重建图像的SSIM、PSNR和RMSE指标都较对比算法有优势;在稀疏角度采样情况下,本文算法的SSIM指标更接近1,重建图像有较好的结构特性。本文算法可以在稀疏采样环境重建出高质量的图像,并有较高的通用性,但是本文使用的非局部约束方法过于繁琐,因此在未来工作中考虑进一步改进非局部先验方法,提高算法在稀疏采样坏境下的性能。

参考文献

[1] Herman GT. Fundamentals of computerized tomography: image reconstruction from projections[M]. Berlin: Springer Science & Business Media, 2009.

[2] Kachelrieß M, Watzke Oliver, Kalender W A. Generalized multi-dimensional adaptive filtering for conventional and spiral single-slice, multi-slice, and cone-beam CT[J]. Medical Physics, 2001, 28(4): 475-490.

[3] 马晨欣, 胡君杰, 闫镔. CT扇形束滤波反投影图像重建算法优化[J]. 激光与光电子学进展, 2012, 49(9): 091103.

Ma C X, Hu J J, Yan B. Optimization of fan-beam CT filtered backprojection reconstruction algorithm[J]. Laser & Optoelectronics Progress, 2012, 49(9): 091103.

[4] Zhang H, Han H, Wang J, et al. Deriving adaptive MRF coefficients from previous normal-dose CT scan for low-dose image reconstruction via penalized weighted least-squares minimization[J]. Medical Physics, 2014, 41(4): 041916.

[5] Zhang H, Han H. LiangJ, et al. Extracting information from previous full-dose CT scan for knowledge-based Bayesian reconstruction of current low-dose CT images[J]. IEEE Transactions on Medical Imaging, 2015, 35(3): 860-870.

[6] 练秋生, 郝鹏鹏. 基于压缩传感和代数重建法的CT图像重建[J]. 光学技术, 2009, 35(3): 422-425.

Lian Q S, Hao P P. Image reconstruction for CT based on compressed sensing and ART[J]. Optical Technique, 2009, 35(3): 422-425.

[7] Lauzier P T, Chen G H. Characterization of statistical prior image constrained compressed sensing (PICCS): II. Application to dose reduction[J]. Medical Physics, 2013, 40(2): 021902.

[8] 张煜林, 孔慧华, 潘晋孝, 等. 基于Split-Bregman算法的能谱计算机层析图像重建[J]. 光学学报, 2017, 37(4): 0411003.

Zhang Y L, Kong H H, Pan J X, et al. Spectral computed tomographic image reconstruction based on split-bregman algorithm[J]. Acta Optica Sinica, 2017, 37(4): 0411003.

[9] ErtasM, AkanA, YildirimI, et al. An iterative reconstruction for tomosynthesis imaging using non-local means[C]∥Proceedings of IEEE International Conference on Imaging Systems and Techniques, 2014: 325- 328.

[10] Chen J L, Wang L Y, Yan B, et al. Efficient and robust 3D CT image reconstruction based on total generalized variation regularization using the alternating direction method[J]. Journal of X-Ray Science and Technology, 2015, 23(6): 683-699.

[11] Bredies K, Kunisch K, Pock T. Total generalized variation[J]. SIAM Journal on Imaging Sciences, 2010, 3(3): 492-526.

[12] Dabov K, Foi A, Egiazarian K. Image denoising with block-matching and 3D filtering[J]. Proceedings of SPIE, 2006, 6064: 354-365.

[13] 王林元, 张瀚铭, 蔡爱龙, 等. 非精确交替方向总变分最小化重建算法[J]. 物理学报, 2013, 62(19): 198701.

Wang L Y, Zhang H M, Cai A L, et al. Image reconstruction algorithm based on inexact alternating direction total-variation minimization[J]. Acta Physica Sinica, 2013, 62(19): 198701.

[14] Sidky E Y, Kao C M, Pan X. Accurate image reconstruction from few-views and limited-angle data in divergent-beam CT[J]. Journal of X-ray Science and Technology, 2009, 14(2): 119-139.

[15] Li T F, Li X, Wang J, et al. Nonlinear sinogram smoothing for low-dose X-ray CT[J]. IEEE Transactions on Nuclear Science, 2004, 51(5): 2505-2513.

[16] Wen Z W, Yin W T, Goldfarb D, et al. A fast algorithm for sparse reconstruction based on shrinkage, subspace optimization, and continuation[J]. SIAM Journal on Scientific Computing, 2010, 32(4): 1832-1857.

[17] Xiao Y H, Song H N. Aninexact alternating directions algorithm for constrained total variation regularized compressive sensing problems[J]. Journal of Mathematical Imaging and Vision, 2012, 44(2): 114-127.

[18] Li C B, Yin W T, Jiang H, et al. An efficient augmented Lagrangian method with applications to total variation minimization[J]. Computational Optimization and Applications, 2013, 56(3): 507-530.

3.3 非局部信息约束下的多方向广义全变分模型

4.4 基于改进广义全变分模型的图像重建算法

班晓征, 李志华, 李贝贝, 徐敏达. 基于改进广义全变分的稀疏图像重建算法[J]. 激光与光电子学进展, 2018, 55(11): 111103. Xiaozheng Ban, Zhihua Li, Beibei Li, Minda Xu. Sparse Image Reconstruction Based on Improved Total Generalized Variation[J]. Laser & Optoelectronics Progress, 2018, 55(11): 111103.

基于改进广义全变分的稀疏图像重建算法 下载： 1460次

1 引言

2 广义全变分重建模型

3 改进的广义全变分图像重建模型

3.1 多方向TGV模型

图 1. 像素点的八邻域空间

Fig. 1. Eight neighborhoods of pixel points

3.2 非局部先验约束

3.3 非局部信息约束下的多方向广义全变分模型

4 基于增广拉格朗日算法的模型求解

4.1 辅助变量的更新

4.2 X方向子问题

4.3 p1~p4方向子问题

4.4 基于改进广义全变分模型的图像重建算法

5 实验结果和分析

5.1 实验环境及参数设定

图 2. 测试用的图像数据。(a) Sheep-Logan 图像;(b) Forbild-Head 图像;(c) Forbild-Abdomen 图像

Fig. 2. Image data used in the test. (a) Sheep-Logan image; (b) Forbild-Head image; (c) Forbild-Abdomen image

表 1. 初始化参数设置

Table 1. Initialization parameters

5.2 评价指标

图 3. 测试图像设置的ROI区。(a) Sheep-Logan图像;(b) Forbild-Head图像;(c) Forbild-Abdomen图像

Fig. 3. ROI regions of the test images. (a) Sheep-Logan image; (b) Forbild-Head image; (c) Forbild-Abdomen image

5.3 算法有效性分析

图 4. 有限角度各算法重建图像及ROI对比(Sheep-Logan)。(a)原始图像;(b) L1-Bregman 算法;(c) TV-ADM 算法;(d) TVAL3 算法;(e) TGV-ADM 算法;(f)本文算法

Fig. 4. Reconstructed image and ROI comparison of different algorithms in finite angle (Sheep-Logan). (a) Original image; (b) L1-Bregman algorithm; (c) TV-ADM algorithm; (d) TVAL3 algorithm; (e) TGV-ADM algorithm; (f) proposed algorithm

图 5. 有限角度各算法重建图像及ROI对比(Forbild-Head)。 (a)原始图像;(b) L1-Bregman 算法;(c) TV-ADM 算法;(d) TVAL3 算法;(e) TGV-ADM 算法;(f)本文算法

Fig. 5. Reconstructed image and ROI comparison of different algorithms in finite angle (Forbild-Head). (a) Original image; (b) L1-Bregman algorithm; (c) TV-ADM algorithm; (d) TVAL3 algorithm; (e) TGV-ADM algorithm; (f) proposed algorithm

图 6. 有限角度各算法重建图像及ROI对比(Forbild-Abdomen)。 (a)原始图像;(b) L1-Bregman 算法;(c) TV-ADM 算法;(d) TVAL3 算法;(e) TGV-ADM 算法;(f)本文算法

Fig. 6. Reconstructed image and ROI comparison of different algorithms in finite angle (Forbild-Abdomen). (a) Original image; (b) L1-Bregman algorithm; (c) TV-ADM algorithm; (d) TVAL3 algorithm; (e) TGV-ADM algorithm; (f) proposed algorithm

图 7. 稀疏角度各算法重建图像及ROI对比(Sheep-Logan)。 (a)原始图像;(b) L1-Bregman 算法;(c) TV-ADM 算法;(d) TVAL3 算法;(e) TGV-ADM算法;(f)本文算法

Fig. 7. Reconstructed image and ROI comparison of different algorithms in sparse angle (Sheep-Logan). (a) Original image; (b) L1-Bregman algorithm; (c) TV-ADM algorithm; (d) TVAL3 algorithm; (e) TGV-ADM algorithm; (f) proposed algorithm

图 8. 稀疏角度各算法重建图像及ROI对比(Forbild-Head)。 (a)原始图像;(b) L1-Bregman 算法;(c) TV-ADM 算法;(d) TVAL3 算法;(e) TGV-ADM算法;(f)本文算法

Fig. 8. Reconstructed image and ROI comparison of different algorithms in sparse angle (Forbild-Head). (a) Original image; (b) L1-Bregman algorithm; (c) TV-ADM algorithm; (d) TVAL3 algorithm; (e) TGV-ADM algorithm; (f) proposed algorithm

图 9. 稀疏角度各算法重建图像及ROI对比(Forbild-Abdomen)。 (a)原始图像;(b) L1-Bregman 算法;(c) TV-ADM 算法;(d) TVAL3 算法;(e) TGV-ADM算法;(f)本文算法

Fig. 9. Reconstructed image and ROI comparison of different algorithms in sparse angle (Forbild-Abdomen). (a) Original image; (b) L1-Bregman algorithm; (c) TV-ADM algorithm; (d) TVAL3 algorithm; (e) TGV-ADM algorithm; (f) proposed algorithm

5.4 算法高效性分析

表 2. 有限投影角度下各算法重建图像的客观评价参数对比

Table 2. Comparison of objective evaluation parameters of reconstructed images of various algorithms under finite projection angle

表 3. 稀疏投影角度下各算法重建图像的客观评价参数对比

Table 3. Comparison of objective evaluation parameters of reconstructed images of various algorithms under sparse projection angle

图 10. 各算法重建图像第135列纵向剖面图对比。(a) Sheep-Logan 图像;(b) Forbild-Head 图像;(c) Forbild-Abdomen 图像

Fig. 10. Comparison of vertical cross-sectional views of 135th column of each reconstruction image. (a) Sheep-Logan image; (b) Forbild-Head image; (c) Forbild-Abdomen image

6 结论

Article Outline

相关论文

相关资讯

关于本站 Cookie 的使用提示

全站搜索

基于改进广义全变分的稀疏图像重建算法下载： 1460次