角度受限下稀疏投影数据的改进粒子群优化随机CT重建

高红霞; 罗澜; 骆英浩; 陈展鸿; 马鸽

doi:doi:10.3788/AOS201838.0111003

光学学报, 2018, 38 (1): 0111003, 网络出版: 2018-10-22

角度受限下稀疏投影数据的改进粒子群优化随机CT重建下载： 4158次

Improved Stochastic CT Reconstruction Based on Particle Swarm Optimization for Limited-Angle Sparse Projection Data

论文大纲

高红霞 ¹罗澜 ¹骆英浩 ^1,*陈展鸿 ¹马鸽 ²

作者单位

¹ 华南理工大学自动化科学与工程学院, 广东广州 510640

² 广州大学机械与电气工程学院, 广东广州 510006

AI 词云图 AI一句话精读 AI短摘要

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

摘要

当计算机断层成像(CT)中X射线的采样范围和数量受限时,得到的稀疏投影数据完备性很低,重建算法的搜索空间巨大。基于凸优化思路的迭代求解算法及其改进采用固定搜索路径,难以在有限时间内收敛至全局最优解;粒子群优化具有全局搜索能力,但计算成本和存储代价过高。为解决这类不完备投影数据的重建问题,提出基于粒子群优化的随机稀疏重建算法。首先,通过随机策略生成具有多样性的初始种群,以保证算法的搜索能力;其次,随机选择梯度下降或基于个体历史最优解和全局历史最优解的随机方向进行迭代,以兼顾算法效率和搜索方向的多样性;最后,基于适应度评价,有针对性地重新生成随机初始种群,强制跳离局部最优。针对角度受限下无噪声和含噪声的稀疏投影数据,分别进行重建实验。结果显示,与常见的凸优化迭代和粒子群优化算法相比,本文算法既能保证算法效率,又在重建质量和算法稳健性上具有明显优势。

Abstract

Because of the sampling scope and quantity limitation in the computed tomography (CT), the completeness of sparse projection data is very low, which leads to a huge search space for the reconstruction algorithm. The iterative algorithm based on convex optimization can not converge to the global minima in finite time due to the fixed search path. Particle swarm optimization has global search capability, but costs tremendous computation and memory. To improve the quality of reconstruction from incomplete projection data, a new stochastic sparse reconstruction algorithm based on particle swarm optimization is proposed. Firstly, the initial solutions with diversity are generated by the stochastic strategy to ensure the search capability. Secondly, the proposed algorithm stochastically chooses either gradient descent direction or random direction based on the local best known solution and the global best known solution in the iteration, to ensure the efficiency of this algorithm and the diversity of search directions. Finally, to avoid trapping in local optimum, the random initial populations are generated based on the fitness evaluation, which represents the current situation. The contrast reconstruction experiments are conducted on both noise-free and noisy limited-angle sparse projection data. The experimental results demonstrate that the proposed algorithm is efficient and evidently superior in reconstruction quality and robustness compared to common iterative algorithms based on convex optimization or particle swarm optimization.

1 引言

计机断层成像(CT)是依据从多角度方向检测到的经过被探测物的X射线衰减量,重建出被探测物某一希望层面的轴向X射线图像的透视成像技术^[1]。用于CT重建的投影数据通常需要满足Tuy-Smith数据完备性条件^[2-3],即需要在360°或180°内均匀且密集地采集完整的投影数据,才能获得精确的重建效果。然而受X射线剂量、采集时间、扫描工作范围等的限制,在诸多实际应用中,只能在采集范围受限、采集角度稀疏的投影数据上实现CT重建。例如,脑CT中通过减少投影的采样角度来降低射线对脑的损伤,工业CT中对大型不规则工件进行内部缺陷的无损探伤,都会得到不完备投影数据^[4-5]。

在进行不完备投影数据的重建时,通常要利用待重建图像的先验条件,如全变分(TV)约束等将不完备投影数据下的不适定问题转化为具有约束条件的优化问题,以降低重建对数据完备性的依赖^[6-7]。针对角度受限型投影数据的重建问题,2006年Sidky等^[8]利用CT图像梯度稀疏性先验知识提出基于凸集投影和TV最小化(POCS-TV)的算法,实现了角度受限下90°投影范围内64组投影数据的重建,但在重建轮廓处存在伪影。为了提高算法的重建质量,通过设计自适应松弛因子^[9]或改进TV约束项等方法^[10-11],抑制角度受限下不完备投影数据的重建伪影,但受制于梯度下降等确定型求解算法的局限,上述算法并不能在有限时间内得到高质量的重建结果^[12]。在效率优化方面,研究者针对360°或180°投影范围内角度完备的稀疏采样数据开展了重建算法优化研究,针对梯度下降算法收敛速度慢的缺点,利用二阶梯度近似^[13]、共轭梯度下降^[14]和软阈值滤波^[15]等方法进行改进,提高了角度完备下稀疏投影数据的重建效率。但当面对完备性更低的角度受限下稀疏投影数据的重建问题时,由于数据完备性、对称性更低,重建算法的搜索空间非常巨大,前述这些改进算法无法解决。传统的基于凸优化思路的确定型重建算法,从单一初始解开始沿一定的下降路径迭代求解,无法在有限时间内遍历这类不完备数据带来的巨大搜索空间,找到符合重建要求的解^[16-18]。这使得角度受限下稀疏投影数据的求解成为CT重建的研究难点之一。

随机优化^[19-20]通过随机生成样本保证解的多样性,搜索中用群体的随机采样来代表整个群体^[16],由适应度评价解的品质,最终由种群进化得到全局最优解,更适于不完备数据下的全局搜索。然而,在求解数据完备性很低的不适定问题时,其计算成本过高;因此,虽然在CT图像复原、脑CT图像配准和分割等相关研究上,研究人员已开展了随机优化方法的探索^[21-23],但在不完备投影数据的CT重建上仍未见到采用随机算法的研究报道。

针对角度受限下稀疏投影数据的CT重建问题,本文提出一种基于粒子群优化(PSO)的随机稀疏重建算法。首先,通过随机策略生成具有多样性的初始种群,以保证算法的搜索能力;然后,随机选择梯度下降或基于个体历史最优解和全局历史最优解的随机下降方向进行迭代,以兼顾算法效率和搜索方向的多样性;最后,基于适应度评价,有针对性地重新生成随机初始种群,强制跳离局部最优。角度受限下稀疏投影数据的重建实验表明,与现有CT重建算法相比,本文算法在保证算法效率的同时,明显提升了重建效果和稳健性,尤其在含噪投影数据的重建上,具有更强的适应性和稳定性。

2 不完备投影数据的重建模型

X射线的反投影重建是指X射线源和探测器同时绕待重建对象的中心旋转扫描,采集不同方向的衰减后的X射线投影数据,并据此反推算出反映对象密度分布的灰度图像的过程。离散化的投影过程如图1所示,其中x_j为待重建图像x第j个像素的灰度,p_i为探测器的第i条X射线测量值。

图 1. 离散化的X射线投影过程

Fig. 1. Discrete projection process of X rays

下载图片查看所有图片

X射线的反投影重建可归结为如下的线性方程求解^[1]:

\begin{matrix} p = Ax, (1) \end{matrix}

式中x= $\begin{matrix} (x_{1}, x_{2}, \dots, x_{N})^{T} \end{matrix}$ ∈R^N×¹,为转换成列向量形式的待重建图像,N=W×H,W和H分别为图像的宽和高;p=(p₁,p₂,…,p_M)^T∈R^M×¹,为投影数据,M为投影数据的总个数;A= $\begin{matrix} (a_{ij})_{M \times N} \end{matrix}$ 为投影系数矩阵,a_ij的值定义为第i条X射线与第j个像素相交的长度,由射线源和待重建对象的相对位置确定。因此,CT投影重建可表达为如下所示的目标函数优化问题^[1]:

\begin{matrix} \min_{x} ‖ Ax - p ‖_{2}^{2} 。 (2) \end{matrix}

为了获得精确的重建结果,通常要求在360°或180°内均匀且密集地采集完整的投影数据。然而,如前所述,医学和工业CT中都存在如图2所示的采集范围受限、采集角度稀疏的不完备投影数据重建问题。这时,X射线的投影数据总数M远小于图像像素总数N,投影系数矩阵A严重不满秩,使得不完备投影数据下的重建成为不适定问题^[24-25],目标函数[(2)式]的解不唯一。

图 2. 不完备投影数据示意图。(a)投影角度受限;(b)投影角度稀疏;(c)投影角度受限且稀疏

Fig. 2. Illustration of incomplete projection data. (a) Limited-angle data; (b) sparse view data; (c) limited-angle sparse view data

下载图片查看所有图片

为改善不完备数据带来的求解不适定性,Sidky等^[8]基于CT重建图像的梯度稀疏性先验知识,引入TV正则化项,得到(3)式的不完备投影数据重建模型,可解决投影角度受限、稀疏或完备性更低的角度受限下稀疏投影数据的重建问题。

\begin{matrix} \min_{x} ‖ Ax - p ‖_{2}^{2} + λ ‖ x ‖_{TV}, (3) \end{matrix}

式中‖x‖_TV=∑ $\begin{matrix} \sqrt[]{(x_{m, n} - x_{m - 1, n})^{2} + (x_{m, n} - x_{m, n - 1})^{2}} \end{matrix}$ 为图像TV项,x_m_,_n是待重建图像x中的第[(m-1)×W+n]个元素(W为待重建图像的宽,m、n分别为图像行、列坐标)。

在这一求解框架下,Sidky等^[8]基于凸优化的思路提出了经典的传统的CT重建算法(POCS-TV算法),对图2(a)、(b)所示的角度受限和角度稀疏的投影数据实现了有效重建。该算法基于某一初始值,以目标函数的负梯度为搜索方向,经过反复迭代对目标函数寻优。之后,为提高POCS-TV算法的重建质量,抑制不完备投影数据重建的伪影,研究人员开展了改进研究。Sidky等^[9]利用前一次主迭代和当前主迭代之间的POCS投影和TV梯度下降的结果差值设计步长参数,自适应地修改梯度下降步长,以保证数据一致性与图像梯度稀疏性之间的平衡,提高稀疏投影数据的重建质量。Liu等^[10]用各向异性TV替换各项同性TV,提出自适应权重全变分-凸集投影(AWTV-POCS)算法,加强了TV项对图像边缘保持的作用。

然而,在对图2(c)所示的角度受限下稀疏投影数据进行重建时,上述确定型迭代算法因下降路径单一,难以在有限时间内遍历完备性更低的投影数据带来的巨大解空间^[16],陷入局部最优,得不到高质量的重建结果。因此,设计高效、可行的重建算法对数据完备性更低的角度受限下稀疏投影数据实现有效重建,降低对投影范围、投影间隔的要求,拓展CT重建的应用场合,成为不完备投影数据重建研究中亟待解决的难题。

3 角度受限下稀疏投影数据的改进PSO随机重建算法

3.1 PSO算法

PSO是由Kennedy^[26]提出的一种进化计算技术。其主要思想是利用群体中个体对信息的共享,使整个群体的运动在问题求解空间中从无序到有序演化,从而获得最优解^[26-27]。算法主要步骤^[28]如下:

1) 随机生成一定规模的初始种群,保证样本的多样性。

2) 计算粒子适应度值来评价解的品质,更新个体和全局历史最优解。

3) 粒子依据个体和群体的历史最优解更新自身的运动速度和位置。

4) 重复执行步骤2)~3),直至迭代次数达到设定值或相邻2次迭代的偏差小于给定阈值。

作为一种随机搜索策略,PSO避免了遗传算法中“交叉”和“变异”的复杂规则,种群更新方式简单,实现容易,精度高,收敛快^[29]。但面对不完备数据带来的巨大搜索空间时,种群规模要保持在粒子自身维度的5~10倍以上,才能使粒子分散分布在搜索空间中,提供足够的搜索能力。以CT重建中常见的128 pixel×128 pixel图像为例,粒子群的种群规模要求在10⁵以上,种群中每个粒子的维度在10⁴以上,每次迭代的计算量和存储代价都很大,寻优效率低,并不能完全解决角度受限下稀疏投影数据的重建问题。

3.2 基于PSO的随机稀疏重建算法

对图2(c)所示的角度受限情况下的稀疏投影数据进行重建。基于凸优化思路的确定型迭代从某一初始解出发,以负梯度为搜索方向,保证向目标函数的优化方向迭代,但无法在有限时间内对巨大的搜索空间进行寻优,易陷入局部最优;PSO具有全局搜索能力,但对种群规模的要求,使得计算成本和存储代价过高。本文针对目标函数[(3)式],提出基于PSO的随机稀疏重建算法。该算法融合了PSO和凸优化的梯度下降算法,在初始解生成、算法迭代方向更新和跳离局部最优中引入随机策略,在保证算法效率的同时提高了不完备投影数据的重建质量。该算法的总体思路如下。

3.2.1 基于随机策略的初始解生成

梯度下降算法在对由不完备投影数据带来的巨大搜索空间进行寻优时,算法初始解在解空间中的位置严重影响着最终解的质量。为保证初始解的多样性,本文对目标函数引入随机策略来生成规模为K的初始解种群。首先,对CT投影重建物理模型Ax=p求其伪逆,生成靠近最优解的初始解:

\begin{matrix} x_{0}^{0} = (A^{T} {A)}^{- 1} A^{T} p, (4) \end{matrix}

以保证算法的搜索效率。然后,为增加初始解的多样性,在 $\begin{matrix} x_{0}^{0} \end{matrix}$ 邻域内随机生成初始种群{ $\begin{matrix} x_{k}^{0} \end{matrix}$ },以保证算法的搜索能力,避免算法过早陷入局部最优,公式为

\begin{matrix} x_{k}^{t} = x_{0}^{t} + rand [length (x_{0}^{t})], k = 1,2, \dots, K - 1, (5) \end{matrix}

式中 $\begin{matrix} x_{k}^{t} \end{matrix}$ 表示第t次迭代中第k个种群粒子,此时t=0;rand[length(x)]函数生成维度与自变量x相同的随机向量,其元素为0~1之间的随机数。这里,K值越大,初始解种群的多样性越高,算法全局搜索能力越强,但算法的计算量也会显著增加,总体效率降低。因此,K值的选取主要取决于算法设计时对全局搜索能力或总体效率的偏好。

3.2.2 基于随机策略的算法迭代方向更新

梯度下降算法中,每次的迭代方向即为目标函数的优化方向,但单一下降方向易陷入局部最优。为保证搜索方向的多样性,本文融合粒子群的随机更新策略,在每次迭代中随机更新算法的搜索方向,提高全局搜索能力。算法迭代更新的策略如下。

1) 对解的取值进行约束,保证随机解的有效性:

\begin{matrix} x_{k}^{t} = \{\begin{matrix} 0, & if x_{k}^{t} (j) < 0 \\ 1, & if x_{k}^{t} (j) > 1 \\ x_{k}^{t}, & otherwise \end{matrix} 。 (6) \end{matrix}

2) 依据(3)式得到对应的适应度函数:

\begin{matrix} F (x_{k}^{t}) = {(A x_{k}^{t} - p)}_{2}^{2} + λ {(x_{k}^{t})}_{TV} 。 (7) \end{matrix}

对随机解进行评价,并计算得到当前的个体历史最优解 $\begin{matrix} P_{k}^{t + 1} \end{matrix}$ 和种群全局历史最优解 $\begin{matrix} P_{g}^{k + 1} \end{matrix}$ :

\begin{matrix} \begin{matrix} P_{k}^{t + 1} = \underset{x_{k}^{t}}{arg mi n_{t}} F (x_{k}^{t}), (8) \\ P_{g}^{k + 1} = \underset{x_{k}^{t}}{arg mi n_{k, t}} F (x_{k}^{t}) 。 (9) \end{matrix} \end{matrix}

3) 对样本进行随机迭代更新:

\begin{matrix} \begin{matrix} x_{k}^{t + 1} = x_{k}^{t} + \nabla x_{k}^{t} [rand ()], (10) \\ \begin{matrix} \nabla x_{k}^{t} [rand ()] = \\ \{\begin{matrix} - d \cdot \nabla_{x_{k}^{t}} {‖x_{k}^{t}‖}_{TV}, if rand () \geq V \\ \begin{matrix} \begin{matrix} w \cdot rand [length (x_{k}^{t})] + c_{1} \cdot rand () \cdot \\ (P_{k}^{t + 1} - x_{k}^{t}) + c_{2} \cdot rand () \cdot (P_{g}^{k + 1} - x_{k}^{t}) \end{matrix}, if rand () < V \end{matrix} \end{matrix}, (11) \end{matrix} \end{matrix} \end{matrix}

式中Δ $\begin{matrix} x_{k}^{t} \end{matrix}$ 为粒子当前的迭代方向;rand()生成0到1之间的随机数;d为下降步长参数; $\nabla \begin{matrix} _{x_{k}^{t}} \end{matrix}$ ‖ $\begin{matrix} x_{k}^{t} \end{matrix}$ ‖_TV为 $\begin{matrix} x_{k}^{t} \end{matrix}$ 的TV的梯度;rand[length( $\begin{matrix} x_{k}^{t} \end{matrix}$ )]生成维度与 $\begin{matrix} x_{k}^{t} \end{matrix}$ 相同的随机向量,以保证本次的迭代方向与上一次的无关;w、c₁和c₂为各分量系数,用于调整随机向量、个体历史最优解和全局历史最优解在新的随机方向上的比例,其中,c₁、c₂为常数,w随迭代次数而减小,使算法在迭代后期适当缩小搜索范围;V为阈值。

在这一迭代策略下,当产生的随机数大于等于阈值时,采用梯度下降方向更新样本,降低了大规模空间遍历中的不确定性,使得粒子群规模K可保持在5~10的较低水平,一定程度上保证了算法的重建效率;否则,采用基于个体历史最优和全局历史最优的随机方向进行迭代更新,保证搜索方向的多样性。

阈值越大,粒子群随机方向的选择概率越大,全局搜索能力提高,但算法效率下降;反之,则更多选择确定性的梯度下降方向,以提高算法效率,但却牺牲了全局搜索能力。此外,在基于随机策略的初始解生成中,初始种群规模K越大,粒子群更新的计算量也会相应增加。此时,应减小阈值,增加确定性下降方向的选择概率,以保证算法效率。但如前所述,样本的随机迭代更新中引入了梯度下降方向,降低了对粒子群种群规模和迭代方向随机性的要求;因此,在实际应用中只需选取较小的K和V值,即可在保证算法效率的同时获得好的重建效果。

3.2.3 跳离局部最优的随机策略

对角度受限下的稀疏投影数据进行迭代重建时,常因搜索空间巨大,使解连续多次无法更新,陷入局部最优。实验观察发现,陷入局部最优的重建图像主要表现为平滑区出现噪声、边缘不清晰或整体过平滑3种情况。针对这一问题,本文提出如下的随机策略,来跳离局部最优:

1) 当连续3次迭代没有更新P_g时,执行如下操作

\begin{matrix} \begin{matrix} P_{g_smooth} = g \otimes P_{g}, (12) \\ P_{g_sharpen} = P_{g} + edge (P_{g}), (13) \\ P_{g_random} = P_{g} + rand [length (P_{g})], (14) \end{matrix} \end{matrix}

式中g为3×3的高斯滤波模板,⊗为卷积操作,P_{g_smooth}为滤波后的样本,edge(P_g)为对P_g进行sobel边缘检测的结果,P_{g_sharpen}为边缘增强后的样本,rand[length(P_g)]生成维度与P_g相同的随机向量,其元素为0~1之间的随机数,P_{g_random}为随机扰动后的样本。

2) 用适应度函数F(x)=‖Ax-p $\begin{matrix} ‖_{2}^{2} \end{matrix}$ +λ‖x‖_TV评价P_{g_smooth}、P_{g_sharpen}和P_{g_random},将适应度最优的样本

\begin{matrix} x_{0}^{t + 1} = arg min [F (P_{g_smooth}), F (P_{g_sharpen}), F (P_{g_random})], (15) \end{matrix}

作为新的初始解。

3) 应用(5)式重新生成随机初始种群。

在这一随机策略中,适应度函数是依据算法的目标函数而确定的,兼顾了投影重建的数据保真、噪声滤除和边缘增强。适应度最优保证了后续的种群更新朝向目标函数的下降方向。因此,通过比较图像滤波、边缘增强和随机扰动3种操作的适应度,可有选择性地确定当前局部最优下应采取的图像操作,再重新生成随机初始种群,使得平滑区出现噪声、边缘不清晰或者整体过平滑等情况强制跳离局部最优。

综合以上基于随机策略的初始解生成、算法迭代方向更新和跳离局部最优,本文提出的基于PSO的随机稀疏重建算法具体流程如下。

输入: 投影数据p,投影矩阵A

输出: 重建图像x

1) 初始化。置当前迭代次数t=0, $\begin{matrix} P_{k}^{0} \end{matrix}$ =0, $\begin{matrix} P_{g}^{0} \end{matrix}$ =0;设置T、K、V、w、c₁、c₂;进化停止计数值q=0。

2) 根据(4)式计算初始解。

3) 根据(5)式生成随机初始种群。

4) 根据(6)式对 $\begin{matrix} x_{k}^{t} \end{matrix}$ 进行取值有效性约束。

5) 根据(7)式计算解的适应度。

6) 计算得到当前的个体历史最优解 $\begin{matrix} P_{k}^{t + 1} \end{matrix}$ 和种群全局历史最优解 $\begin{matrix} {P^{t + 1}}_{g} \end{matrix}$ 。

\begin{matrix} \begin{matrix} P_{k}^{t + 1} = \{\begin{matrix} \underset{x_{k}^{t}}{arg mi n_{t}} F (x_{k}^{t}), & if mi n_{t} F (x_{k}^{t}) < F (P_{k}^{t}) \\ P_{k}^{t}, & ifmi n_{t} F (x_{k}^{t}) \geq F (P_{k}^{t}) \end{matrix}, \\ P_{g}^{t + 1} = \{\begin{matrix} \underset{x_{k}^{t}}{arg mi n_{k, t}} F (x_{k}^{t}), & if mi n_{k, t} F (x_{k}^{t}) < F (P_{g}^{t}) \\ P_{g}^{t}, & ifmi n_{k, t} F (x_{k}^{t}) \geq F (P_{g}^{t}) \end{matrix} 。 \end{matrix} \end{matrix}

7) 更新进化停止计数值q。若步骤6中 $\begin{matrix} P_{g}^{t + 1} \end{matrix}$ = $\begin{matrix} P_{g}^{t} \end{matrix}$ ,则q=q+1,否则q=0。

8) 若q<3,执行样本随机迭代更新;否则执行跳离局部最优。

① 根据(11)式对样本 $\begin{matrix} x_{k}^{t} \end{matrix}$ 进行随机迭代更新,得到 $\begin{matrix} x_{k}^{t + 1} \end{matrix}$ 。

② 跳离局部最优。

a 对当前的 $\begin{matrix} P_{g}^{t + 1} \end{matrix}$ 执行(12)、(13)和(14)式,分别得到P_{g_smooth}、P_{g_sharpen}和P_{g_random}。

b 根据(15)式计算上述3种操作后适应度值最优的样本 $\begin{matrix} x_{0}^{t + 1} \end{matrix}$ 。

c 重新初始化种群。

\begin{matrix} x_{k}^{t + 1} = x_{0}^{t + 1} + rand [length (x_{0}^{t + 1})], k = 1,2, \dots, K - 1 。 \end{matrix}

9) 迭代终止条件。若满足t=T或 $\begin{matrix} {(P_{g}^{t + 1} - P_{g}^{t})}_{2} \end{matrix}$ <10^-3,则停止迭代,当前 $\begin{matrix} {P^{t + 1}}_{g} \end{matrix}$ 即为重建结果x;否则,t=t+1,返回步骤4)。

4 角度受限下稀疏投影数据的重建实验结果及分析

以128 pixel×128 pixel的二维Shepp Logan模型、Forbild Head模型为CT重建对象,投影数据生成模型如下:

1) X射线探测器由256个接收器阵列组成,单个接收单元大小为0.5 mm。

2) 投影角度为(0°,90°),投影数据的采集组数为30。

3) 待重建图像的每一个像素的物理尺寸为0.5 mm×0.5 mm。

为了验证本文算法的有效性,在无噪声的X射线投影数据和含噪声的X射线投影数据上进行重建实验,通过重建效果、重建误差和重建效率来对比算法性能。

对比算法包括传统凸优化(L2-TV)算法^[30-32]、基于凸优化的改进算法——随机平均梯度(SAG)算法^[33-34]、POCS-TV算法^[8]和PSO算法^[35]。为更好地比较各算法的重建效果,设置PSO算法的种群规模K=10000,迭代次数小于等于5000。其他各算法的迭代停止条件:相邻2次迭代的重建结果差值小于10^-3,最大迭代次数小于等于1000。

本文算法的参数设置为:T=1000,K=5,V=0.4,c₁=c₂=2,w=w_max-(w_max-w_min)(t/T),式中w_max=0.95,w_min=0.4,t表示当前迭代次数。

本文实验的软硬件环境为:Matlab R2016a,Intel(R)Core(TM)i7-4790K CPU@4.00 GHz,32 GB内存。

4.1 角度受限下无噪声稀疏投影数据的重建实验结果

图3、4分别为各重建算法对投影角度(0°,90°)下的Shepp Logan图像和Forbild Head图像的无噪声稀疏投影数据进行重建的结果。

由图3和图4的重建结果可以看出,因投影角度受限,L2-TV算法及POCS-TV算法重建的图像边缘不够清晰,也存在不同程度的伪影;SAG算法有较好的重建效果,但在重建Shepp Logan图像时边缘恢复能力不足;PSO算法因种群更新的随机性,即使在种群规模远高于本文算法的情况下,经5000次迭代也无法实现有效重建;本文算法在有限的迭代次数内获得了高质量的重建结果,在图像的完整性、边缘恢复和去除伪影能力上都表现出了明显的优势。

图 3. Shepp Logan图像的无噪声稀疏投影数据重建结果。(a) Shepp Logan原图;(b) L2-TV算法;(c) SAG算法;(d) POCS-TV算法;(e) PSO算法;(f)本文算法

Fig. 3. Reconstruction results from noise-free Shepp sparse Logan projection data. (a) Original image of Shepp Logan; (b) L2-TV algorithm; (c) SAG algorithm; (d) POCS-TV algorithm; (e) PSO algorithm; (f) proposed algorithm

下载图片查看所有图片

图 4. Forbild Head图像的无噪声稀疏投影数据重建结果。(a) Forbild Head原图;(b) L2-TV算法;(c) SAG算法;(d) POCS-TV算法;(e) PSO算法;(f)本文算

Fig. 4. Reconstruction results from noise-free Forbild Head sparse projection data. (a) Original image of Forbild Head; (b) L2-TV algorithm; (c) SAG algorithm; (d) POCS-TV algorithm; (e) PSO algorithm; (f) proposed algorithm

下载图片查看所有图片

为更好地展示本文算法的重建效果,原始图像与各重建图的水平剖线及垂直剖线对比结果如图5、6所示。其中,图5为二维Shepp Logan图像的重建剖线对比,图6为二维Forbild Head图像的重建剖线对比。可以看出,本文算法重建的图像水平剖线、垂直剖线几乎与原始图像完全重合,优于其他对比算法。

图 5. 无噪声情况下Shepp Logan图像投影数据的水平剖线(a)及垂直剖线(b)对比

Fig. 5. Comparison of horizontal (a) and vertical (b) profiles with noise-free Shepp Logan projection data

下载图片查看所有图片

图 6. 无噪声情况下Forbild Head图像投影数据的水平剖线(a)及垂直剖线(b)对比

Fig. 6. Comparison of horizontal (a) and vertical (b) profiles with noise-free Forbild Head projection data

下载图片查看所有图片

此外,为客观对比各算法的重建效果,采用重建结果与原始图像的误差值进行算法评价,结果如表1所示。可以看出,本文算法的重建误差明显优于其他算法。

表 1. 无噪声情况下各算法重建误差比较

Table 1. Comparison of reconstruction error with noise-free projection data

Image	L2-TV algorithm	SAG algorithm	POCS-TV algorithm	PSO algorithm	Proposed algorithm
Shepp Logan	24.5674	5.0062	197.4777	3340.6008	0.0439
Forbild Head	26.9956	0.0833	236.7307	3515.6233	0.0379

查看所有表

4.2 角度受限下含噪声稀疏投影数据的重建实验结果

为验证本文算法的稳健性和适应性,对角度受限下含噪声稀疏投影数据的重建结果,从重建图像效果、剖线对比和重建误差3个方面进行算法性能对比。含噪声投影数据由无噪声投影数据加上标准差为1.5的较强高斯噪声得到。

图7和图8分别为各重建算法对(0°,90°)下Shepp Logan图像和Forbild Head图像的含噪声投影数据进行重建的结果。图9和图10分别为Shepp Logan、Forbild Head原始图像与各算法重建图像的水平和垂直剖线对比。表2为各算法重建图像与原始图像的误差值。

从图7~10和表2可以看出,本文算法的重建结果边缘相对清晰,平滑区光滑,较好地抑制了重建伪影,在含噪声投影数据下表现出更好的稳健性和适应性。L2-TV算法、SAG算法和POCS-TV算法都存在明显的伪影和边缘缺失现象,而PSO算法由于其样本更新方式的限制,经过5000次迭代依然无法实现重建。

图 7. Shepp Logan图像的含噪声投影数据重建结果。(a) Shepp Logan原图;(b) L2-TV算法;(c) SAG算法;(d) POCS-TV算法;(e) PSO算法;(f)本文算法

Fig. 7. Reconstruction results from Shepp Logan projection data with Gaussian noise. (a) Original image of Shepp Logan; (b) L2-TV algorithm; (c) SAG algorithm; (d) POCS-TV algorithm; (e) PSO algorithm; (f) proposed algorithm

下载图片查看所有图片

图 8. Forbild Head图像的含噪声投影数据重建结果。(a) Forbild Head原图;(b) L2-TV算法;(c) SAG算法;(d) POCS-TV算法;(e) PSO算法;(f)本文算法

Fig. 8. Reconstruction results from Forbild Head projection data with Gaussian noise. (a) Original image of Forbild Head; (b) L2-TV algorithm; (c) SAG algorithm; (d) POCS-TV algorithm; (e) PSO algorithm; (f) proposed algorithm

下载图片查看所有图片

图 9. 含噪声情况下Shepp Logan图像投影数据的水平剖线(a)及垂直剖线(b)对比

Fig. 9. Comparison of horizontal (a) and vertical (b) profiles with noisy Shepp Logan projection data

下载图片查看所有图片

图 10. 含噪声情况下Forbild Head图像投影数据的水平剖线(a)及垂直剖线(b)对比

Fig. 10. Comparison of horizontal (a) and vertical (b) profiles with noisy Forbild Head projection data

下载图片查看所有图片

表 2. 含噪声情况下各算法重建误差比较

Table 2. Comparison of reconstruction error with noisy projection data

Image	L2-TV algorithm	SAG algorithm	POCS-TV algorithm	PSO algorithm	Proposed algorithm
Shepp Logan	28.2731	38.2937	215.8484	3363.2369	26.1932
Forbild Head	41.3788	45.8357	251.2832	3380.0266	30.1645

查看所有表

4.3 重建算法效率对比

为验证本文算法的效率,对PSO算法以外的各算法,对比其重建误差随迭代次数和时间变化的情况,因为PSO算法无法在有限的迭代次数、计算时间内获得重建结果,所以不参与对比。图11给出了各算法重建误差随迭代次数变化的曲线,图12给出了各算法重建误差随时间变化的曲线。图中(a)和(b)分别为无噪声的Shepp Logan 和Forbild Head投影数据重建误差曲线,图中(c)和(d)分别为含噪声的Shepp Logan和Forbild Head投影数据重建误差曲线。

图 11. 各算法重建误差随迭代次数变化的曲线。(a)无噪声的Shepp Logan投影数据;(b)无噪声的Forbild Head投影数据;(c)含噪声的Shepp Logan投影数据;(d)含噪声的Forbild Head投影数据

Fig. 11. Reconstruction error curves with the change of iteration number of different algorithms. (a) Noise-free Shepp Logan projection data; (b) noise-free Forbild Head projection data; (c) noisy Shepp Logan projection data; (d) noisy Forbild Head projection data

下载图片查看所有图片

图 12. 各算法重建误差随迭代时间变化的曲线。(a)无噪声的Shepp Logan投影数据;(b)无噪声的Forbild Head投影数据;(c)含噪声的Shepp Logan投影数据;(d)含噪声的Forbild Head投影数据

Fig. 12. Reconstruction error curves with the change of iteration time of different algorithms. (a) Noise-free Shepp Logan projection data; (b) noise-free Forbild Head projection data; (c) noisy Shepp Logan projection data; (d) noisy Forbild Head projection data

下载图片查看所有图片

由图11、12可知,本文算法的迭代误差下降快,整体重建效率高。尤其是在重建误差随迭代次数变化的曲线上表现更为明显,同时迭代时间也较满意,说明本文算法单次迭代的计算量与确定型迭代算法基本持平。基于凸优化的L2-TV、SAG和POCS-TV算法也能在有限时间内获得重建结果,尤其是SAG算法在无噪声的Shepp Logan图像重建上时间效率较好,但这类算法采用固定的梯度下降方向,迭代中存在误差振荡,且整体重建误差明显高于本文算法。

5 结论

针对角度受限下稀疏投影数据的重建,将梯度下降法和PSO策略相融合,解决了初始解多样性、算法迭代方向的随机更新和陷入局部最优等问题。其中,梯度下降方向的引入降低了大规模空间遍历中的不确定性,使得算法能够在粒子群的种群规模较小时实现对不完备投影数据的高效、高质量重建。角度受限下无噪声和含噪声稀疏投影数据的对比实验结果显示,本文算法既保证了算法效率,又在重建质量和算法稳健性上有明显的优势。

本文的实验设计对投影的角度范围和稀疏性都提出了非常严苛的要求,得到的投影数据完备性和对称性很低;因此,在含噪声投影数据的重建图像中,与投影方向垂直的边缘上还存在一定的重建误差。后续研究将考虑投影数据在不同方向上的完备性,进一步提高算法的重建能力。算法迭代方向的随机更新中也可探索自适应地调整梯度下降与粒子群随机方向的比重,提高算法的重建效率。

此外,依据目标函数对本文算法中的适应度评价函数做相应修改,可将随机稀疏重建算法推广应用到其他不完备观测数据下的高维信号盲复原问题中,提升算法性能和效率。

参考文献

[1] 庄天戈. CT原理与算法[M]. 上海: 上海交通大学出版社, 1992: 77- 97.

[2] Heang K T. An inversion formula for cone-beam reconstruction[J]. Siam Journal on Applied Mathematics, 1983, 43(3): 546-552.

[3] Bruce D S. Image reconstruction from cone-beam projections: necessary and sufficient conditions and reconstruction methods[J]. IEEE Transactions on Medical Imaging, 1985, 4(1): 14-25.

[4] Hunink M G, Gazelle G S. CT screening: a trade-off of risks, benefits, and costs[J]. Journal of Clinical Investigation, 2003, 111(11): 1612-1619.

[5] 张朝宗. 工业CT技术和原理[M]. 北京: 科学出版社, 2009: 46.

[6] 王林元, 刘宏奎, 李磊, 等. 基于稀疏优化的计算机断层成像图像不完全角度重建综述[J]. 物理学报, 2014, 63(20): 208702.

Wang L Y, Liu H K, Li L, et al. Review of sparse optimization-based computed tomography image reconstruction from few-view projections[J]. Acta Physica Sinica, 2014, 63(20): 208702.

[7] 郭红波, 贺小伟, 侯榆青, 等. 基于非凸稀疏正则的荧光分子断层成像[J]. 光学学报, 2015, 35(7): 0717001.

Guo H B, He X W, Hou Y Q, et al. Fluorescence molecular tomography based on nonconvex sparse regularization[J]. Acta Optica Sinica, 2015, 35(7): 0717001.

[8] 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, 2006, 14(2): 119-139.

[9] Emil Y S, Pan X. Image reconstruction in circular cone-beam computed tomography by constrained, total-variation minimization[J]. Physics in Medicine & Biology, 2008, 53(17): 4777-4807.

[10] Liu Y, Ma J H, Fan Y, et al. Adaptive-weighted total variation minimization for sparse data toward low-dose X-ray computed tomography image reconstruction[J]. Physics in Medicine & Biology, 2012, 57(23): 7923-7956.

[11] 李毅. 有限角度三维CT图像重建算法研究[D]. 太原: 中北大学, 2011.

LiY. The study of limited angle three-dimensional CT image reconstruction algorithm[D]. Taiyuan: North University of China, 2011.

[12] Dhawan A P, Rangayyan R M, Gordon R. Image restoration by Wiener deconvolution in limited-view computed tomography[J]. Applied Optics, 1985, 24(23): 4013-4020.

[13] Park J C, Song B, Jin S K, et al. Fast compressed sensing-based CBCT reconstruction using Barzilai-Borwein formulation for application to on-line IGRT[J]. Medical Physics, 2012, 39(3): 1207-1217.

[14] 邹晶, 孙艳勤, 张朋. 由少量投影数据快速重建图像的迭代算法[J]. 光学学报, 2009, 29(5): 1198-1204.

Zou J, Sun Y Q, Zhang P. A fast iterative image reconstruction algorithm from few-views projection data[J]. Acta Optica Sinica, 2009, 29(5): 1198-1204.

[15] Yu H Y, Wang G. A soft-threshold filtering approach for reconstruction from a limited number of projections[J]. Physics in Medicine & Biology, 2010, 55(13): 3905-3916.

[16] Karp R M. An introduction to randomized algorithms[J]. Discrete Applied Mathematics, 1991, 34(1/2/3): 165-201.

[17] Andersson T, Lathen G, Lenz R, et al. Modified gradient search for level set based image segmentation[J]. IEEE Transactions on Image Processing, 2013, 22(2): 621-630.

[18] BianchiniS. On the Euler-Lagrange equation for a variational problem[M]. [S.l.]: Birkhäuser Basel, 2007: 61- 77.

[19] FreyJ. Introduction to stochastic search and optimization: estimation, simulation, and control[M]. Hoboken: John Wiley & Sons, 2003: 368- 369.

[20] Solis F J. Wets R J B. Minimization by random search techniques[J]. Mathematics of Operations Research, 1981, 6(1): 19-30.

[21] 陈昊, 颜秀铭, 曹福凯, 等. 粒子群优化算法在CT图像复原中的应用研究[ C]. 全国光谱仪器与分析学术研讨会, 2007.

[22] 于泽. 脑CT图像配准中自适应粒子群算法的应用研究[D]. 大连: 大连海事大学, 2010.

YuZ. Research on cerebral CT image registration by adaptive particle swarm optimization[D]. Dalian: Dalian Maritime University, 2010.

[23] Zhao J J, Ji G H, Xia Y, et al. Cavitary nodule segmentation in computed tomography images based on self-generating neural networks and particle swarm optimization[J]. International Journal of Bio-Inspired Computation, 2015, 7(1): 567-580.

[24] IdierJ. Bayesian approach to inverse problems[M]. Hoboken: John Wiley & Sons, 2008.

[25] 卢孝强, 孙恰. 基于乘性正则化的有限角度CT重建算法[J]. 光学学报, 2010, 30(5): 1285-1290.

Lu X Q, Sun Q. Limited angle computed tomography reconstruction algorithm based on multiplicative regularization method[J]. Acta Optica Sinica, 2010, 30(5): 1285-1290.

[26] KennedyJ. Particle swarm optimization[M] ∥ Encyclopedia of Machine Learning. New York: Springer, 2011: 760- 766.

[27] Shi X H, Liang Y C, Lee H P, et al. An improved GA and a novel PSO-GA-based hybrid algorithm[J]. Information Processing Letters, 2005, 93(5): 255-261.

[28] 高璇. 粒子群算法优化及其在图像检索中的应用研究[D]. 西安: 西安电子科技大学, 2013.

GaoX. Research on particle swarm optimization and its application in image retrieval[D]. Xi'an:Xidian University, 2013.

[29] 方峻. 粒子群算法及其应用研究[D]. 成都: 电子科技大学, 2006.

[30] NojimaY, Chen YW, Han XH. Image and video restoration with TV/L2-norm constraint[C]. International Conference on Computer Information Systems and Industrial Applications, 2015: 642- 644.

[31] Liu X W. Efficient algorithms for hybrid regularizers based image denoising and deblurring[J]. Computers & Mathematics with Applications, 2015, 69(7): 675-687.

[32] Liu Y, Liang Z R, Ma J H, et al. Total variation-Stokes strategy for sparse-view X-ray CT image reconstruction[J]. IEEE Transactions on Medical Imaging, 2014, 33(3): 749-763.

[33] Roux N L, Schmidt M, Bach F. A stochastic gradient method with an exponential convergence rate for finite training sets[J]. Advances in Neural Information Processing Systems, 2013, 4: 2663-2671.

[34] SchmidtM, Roux NL, BachF. Minimizing finite sums with the stochastic average gradient[M]. New York: Springer-Verlag, 2017.

[35] Rini D P, Shamsuddin S M, Yuhaniz S S. Particle swarm optimization: technique, system and challenges[J]. International Journal of Computer Applications, 2011, 14(1): 19-26.

1 引言

2 不完备投影数据的重建模型

3 角度受限下稀疏投影数据的改进PSO随机重建算法

3.1 PSO算法

3.2 基于PSO的随机稀疏重建算法

4 角度受限下稀疏投影数据的重建实验结果及分析

4.1 角度受限下无噪声稀疏投影数据的重建实验结果

4.2 角度受限下含噪声稀疏投影数据的重建实验结果

4.3 重建算法效率对比

5 结论

高红霞, 罗澜, 骆英浩, 陈展鸿, 马鸽. 角度受限下稀疏投影数据的改进粒子群优化随机CT重建[J]. 光学学报, 2018, 38(1): 0111003. Hongxia Gao, Lan Luo, Yinghao Luo, Zhanhong Chen, Ge Ma. Improved Stochastic CT Reconstruction Based on Particle Swarm Optimization for Limited-Angle Sparse Projection Data[J]. Acta Optica Sinica, 2018, 38(1): 0111003.

角度受限下稀疏投影数据的改进粒子群优化随机CT重建 下载： 4158次

1 引言

2 不完备投影数据的重建模型

图 1. 离散化的X射线投影过程

Fig. 1. Discrete projection process of X rays

图 2. 不完备投影数据示意图。(a)投影角度受限;(b)投影角度稀疏;(c)投影角度受限且稀疏

Fig. 2. Illustration of incomplete projection data. (a) Limited-angle data; (b) sparse view data; (c) limited-angle sparse view data

3 角度受限下稀疏投影数据的改进PSO随机重建算法

3.1 PSO算法

3.2 基于PSO的随机稀疏重建算法

4 角度受限下稀疏投影数据的重建实验结果及分析

4.1 角度受限下无噪声稀疏投影数据的重建实验结果

图 3. Shepp Logan图像的无噪声稀疏投影数据重建结果。(a) Shepp Logan原图;(b) L2-TV算法;(c) SAG算法;(d) POCS-TV算法;(e) PSO算法;(f)本文算法

Fig. 3. Reconstruction results from noise-free Shepp sparse Logan projection data. (a) Original image of Shepp Logan; (b) L2-TV algorithm; (c) SAG algorithm; (d) POCS-TV algorithm; (e) PSO algorithm; (f) proposed algorithm

图 4. Forbild Head图像的无噪声稀疏投影数据重建结果。(a) Forbild Head原图;(b) L2-TV算法;(c) SAG算法;(d) POCS-TV算法;(e) PSO算法;(f)本文算

Fig. 4. Reconstruction results from noise-free Forbild Head sparse projection data. (a) Original image of Forbild Head; (b) L2-TV algorithm; (c) SAG algorithm; (d) POCS-TV algorithm; (e) PSO algorithm; (f) proposed algorithm

图 5. 无噪声情况下Shepp Logan图像投影数据的水平剖线(a)及垂直剖线(b)对比

Fig. 5. Comparison of horizontal (a) and vertical (b) profiles with noise-free Shepp Logan projection data

图 6. 无噪声情况下Forbild Head图像投影数据的水平剖线(a)及垂直剖线(b)对比

Fig. 6. Comparison of horizontal (a) and vertical (b) profiles with noise-free Forbild Head projection data

表 1. 无噪声情况下各算法重建误差比较

Table 1. Comparison of reconstruction error with noise-free projection data

4.2 角度受限下含噪声稀疏投影数据的重建实验结果

图 7. Shepp Logan图像的含噪声投影数据重建结果。(a) Shepp Logan原图;(b) L2-TV算法;(c) SAG算法;(d) POCS-TV算法;(e) PSO算法;(f)本文算法

Fig. 7. Reconstruction results from Shepp Logan projection data with Gaussian noise. (a) Original image of Shepp Logan; (b) L2-TV algorithm; (c) SAG algorithm; (d) POCS-TV algorithm; (e) PSO algorithm; (f) proposed algorithm

图 8. Forbild Head图像的含噪声投影数据重建结果。(a) Forbild Head原图;(b) L2-TV算法;(c) SAG算法;(d) POCS-TV算法;(e) PSO算法;(f)本文算法

Fig. 8. Reconstruction results from Forbild Head projection data with Gaussian noise. (a) Original image of Forbild Head; (b) L2-TV algorithm; (c) SAG algorithm; (d) POCS-TV algorithm; (e) PSO algorithm; (f) proposed algorithm

图 9. 含噪声情况下Shepp Logan图像投影数据的水平剖线(a)及垂直剖线(b)对比

Fig. 9. Comparison of horizontal (a) and vertical (b) profiles with noisy Shepp Logan projection data

图 10. 含噪声情况下Forbild Head图像投影数据的水平剖线(a)及垂直剖线(b)对比

Fig. 10. Comparison of horizontal (a) and vertical (b) profiles with noisy Forbild Head projection data

表 2. 含噪声情况下各算法重建误差比较

Table 2. Comparison of reconstruction error with noisy projection data

4.3 重建算法效率对比

图 11. 各算法重建误差随迭代次数变化的曲线。(a)无噪声的Shepp Logan投影数据;(b)无噪声的Forbild Head投影数据;(c)含噪声的Shepp Logan投影数据;(d)含噪声的Forbild Head投影数据

Fig. 11. Reconstruction error curves with the change of iteration number of different algorithms. (a) Noise-free Shepp Logan projection data; (b) noise-free Forbild Head projection data; (c) noisy Shepp Logan projection data; (d) noisy Forbild Head projection data

图 12. 各算法重建误差随迭代时间变化的曲线。(a)无噪声的Shepp Logan投影数据;(b)无噪声的Forbild Head投影数据;(c)含噪声的Shepp Logan投影数据;(d)含噪声的Forbild Head投影数据

Fig. 12. Reconstruction error curves with the change of iteration time of different algorithms. (a) Noise-free Shepp Logan projection data; (b) noise-free Forbild Head projection data; (c) noisy Shepp Logan projection data; (d) noisy Forbild Head projection data

5 结论

Article Outline

相关论文

相关资讯

关于本站 Cookie 的使用提示

全站搜索

角度受限下稀疏投影数据的改进粒子群优化随机CT重建下载： 4158次