ISSN 
下载:
-
图像重建属于图像处理中的一个重要分支，其原理是根据现实场景的投影数据来获取场景中的物质分布信息。目前已广泛应用于物体内部结构图像的检测和观察中，尤其在医学影像方面有广泛应用[1]。通过CT成像，可利用射线从不同的方向穿过物体从而得到该物体横截面信息[2]。
图像重建算法大体上可以分为解析类算法和迭代类算法。解析类算法计算量小，对完全投影数据能快速获得较好质量的重建，但要求投影数据具有很高的完整性[3]。迭代类算法的最大优势是在投影数据信噪比低、数据严重缺损的情况下获得高质量的重建图像[4]。常用的代数重建法采用级数迭代的思路，可用于不完全投影数据的图像重建[5]。首先对一幅图像进行离散化，把整个图像区域看作由一个正方形网格组成，一个网格对应的是一个像素，每个像素内部的灰度值视为常数[2]。而后将每个重建区域的投影视为一个方程，通过解线性方程组求出整个重建图像的灰度分布。
由于在不完全投影的情况下，方程组未知数的个数要大于方程个数，是典型的不适定问题。同时投影矩阵一般是大规模的稀疏矩阵，其秩是未知的，因此用传统的矩阵计算方法求解十分困难。一般可采用基于全变差的正则化方法来求解上述线性方程组，优势为在图像重建的过程中能很好地保持图像的边缘信息[6-7]。
用优化的方法解决稀疏角度CT图像重建中的问题当前已受到广泛关注[8-9]，其基本方法都是把CT图像重建问题归结为一个加权最小二乘(penalized weighted least squares, PWLS)问题[10]，然后用代数迭代算法求解。文献[11]把Chambolle-Pock的原始对偶算法[12]应用在CT图像的重建问题中。文献[13]提出在有限角度和稀疏角度下的CT图像重建问题中引入全变差最小化目标函数，通过将上一步迭代重建的图像作为反馈加入到新的迭代中。此类算法大多采用正则化方法进行代数/统计迭代重建来解决稀疏角度问题，可较好地提高重建图像的质量，减少伪影，但若正则化设计不当，存在不收敛的可能性，重建速度较慢。自适应临近点方法在基于全变差的图像恢复问题中具有出色的数值表现[14]。
为平衡稀疏角度CT图像重建过程中的重建质量与重建效率，本文把CT图像重建中原始的全变差模型转变成一个基于全变差的鞍点问题，然后再通过自适应临近点算法求解该鞍点问题。本文算法可在临近点迭代的每一步中自适应地选取临近参数矩阵，而且该临近参数矩阵在每一步迭代中是变动的、正定的，但不要求对称。同时在算法中有一步校正步以保证其收缩性，而且预测步等价于一个加速的原始对偶算法[12]。在收缩算法的框架下, 该算法具有全局收敛性。
-
在不完全投影矩阵是大规模稀疏矩阵的情况下，通过基于全变差的正则化方法来进行图像重建，可描述成如下形式的优化问题[7]：
$$ \mathop {\min }\limits_x \frac{1}{2}\left\| {\mathit{\boldsymbol{Kx}} - \mathit{\boldsymbol{f}}} \right\|_2^2 + \lambda {\left\| {\left| {\nabla \mathit{\boldsymbol{x}}} \right|} \right\|_1} $$ (1) 式中，x为需重建的目标图像向量；K为投影矩阵；f为射线投影值向量；λ为正则化系数；$\nabla $是一个线性算子。
定义 1 假设X，Y都为多维向量空间，存在X到Y的线性转换K，规定G和F为将相应的X和Y空间映射为非负实数字的凸函数[11]。
令$\mathit{\boldsymbol{A}}{\rm{ = }}(_\nabla ^\mathit{\boldsymbol{K}})$，则原始的全变差模型可以写成[11]:
$$ \mathop {\min }\limits_x \{ F(\mathit{\boldsymbol{Ax}}) + G(\mathit{\boldsymbol{x}})\} $$ (2) 式中，$F(\mathit{\boldsymbol{Ax}}) = \frac{1}{2}\left\| {\mathit{\boldsymbol{Kx}} - \mathit{\boldsymbol{f}}} \right\|_2^2 + \lambda {\left\| {\left| {\nabla \mathit{\boldsymbol{x}}} \right|} \right\|_1}$；$G(\mathit{\boldsymbol{x}}) = 0$。
定义 2 指示函数是定义在某集合W上的函数，定义在集合${\rm{Box}}(\omega )$上的指示函数定义为：
$$ {\delta _{{\rm{Box}}(\lambda )}}(|\omega |) = \left\{ {{\rm{ }}_{\infty {\rm{ }}\left\| \omega \right\| > \lambda }^{0{\rm{ }}\left\| \omega \right\| \le \lambda }} \right. $$ (3) 定义 3 设R是赋范线性空间，W是R中的凸集，$\varphi $是W上的凸泛函，规定集合：
$$ {W^*} = \left\{ {{r^*} \in {R^*}|{{\sup }_{r \in W}}( < r, {r^*} > - \varphi (r)) < \infty } \right\} $$ (4) 称其为W的共轭集合。在W*上规定泛函：
$$ {\varphi ^*}({r^*}) = {\sup _{r \in W}}\left\{ { < r, {r^*} > - \varphi (r)} \right\} $$ (5) 称其为$\varphi $的共轭泛函。
设$F(\mathit{\boldsymbol{Ax}}) = {F_1}(\mathit{\boldsymbol{Kx}}) + {F_2}(\nabla \mathit{\boldsymbol{x}}) = {F_1}(\mathit{\boldsymbol{p)}} + {F_2}(\mathit{\boldsymbol{q}})$，则通过定义2和定义3可给出${F_1}$和${F_2}$的共轭函数分别为$F_1^*(\mathit{\boldsymbol{p}}) = \frac{1}{2}\left\| \mathit{\boldsymbol{p}} \right\|_2^2 + {\mathit{\boldsymbol{p}}^{\rm{T}}}\mathit{\boldsymbol{f}}$和$F_2^*(q) = {\delta _{{\rm{Box}}(\lambda )}}(\left| q \right|)$，则${F^*}(\mathit{\boldsymbol{p}}, q) = {F^*}(\mathit{\boldsymbol{Ax}}) = F_1^*(\mathit{\boldsymbol{p}}) + F_2^*(q)$。
可把式(2)转变成如下的鞍点问题：
$$ \begin{array}{c} \mathop {\min }\limits_x \mathop {\max }\limits_{p, q} \Phi (\mathit{\boldsymbol{x}}, \mathit{\boldsymbol{p}}, q): = \\ \left\{ { < \mathit{\boldsymbol{Kx}}, \mathit{\boldsymbol{p}} > - \frac{1}{2}\left\| \mathit{\boldsymbol{p}} \right\|_2^2 - {\mathit{\boldsymbol{p}}^{\rm{T}}}\mathit{\boldsymbol{f}} + < \nabla \mathit{\boldsymbol{x}}, q > - {\delta _{{\rm{Box}}(\lambda )}}(|q|)} \right\} \end{array} $$ (6) 式(6)可以转换成一个变分不等式问题${\rm{VI(}}\mathit{\Omega } {\rm{, }}F{\rm{)}}$，需找到${\mathit{\boldsymbol{u}}^*}$满足如下不等式：
$$ {\rm{VI(}}\mathit{\Omega } {\rm{, }}F{\rm{):}}{(\mathit{\boldsymbol{u}} - {\mathit{\boldsymbol{u}}^*})^{\rm{T}}}F({\mathit{\boldsymbol{u}}^*}) \ge 0 $$ (7) 式中，$\mathit{\boldsymbol{u}}{\rm{ = }}\left( {_\mathit{\boldsymbol{y}}^\mathit{\boldsymbol{x}}} \right){\rm{; }}\mathit{\boldsymbol{y}}{\rm{ = }}\left( {_q^p} \right){\rm{; }}F{\rm{(}}\mathit{\boldsymbol{u}}{\rm{) = }}\left( {_{\nabla {\Phi _y}}^{\nabla {\Phi _x}}} \right);{\rm{ }}\mathit{\Omega } : = \mathit{\boldsymbol{x}} \times \mathit{\boldsymbol{y}}$。
显然，${\rm{ }}F{\rm{(}}\mathit{\boldsymbol{u}}{\rm{)}}$是单调的，因此，${\rm{VI(}}\mathit{\Omega } {\rm{, }}F{\rm{)}}$为单调变分不等式，从而${\rm{VI(}}\mathit{\Omega }{\rm{, }}F{\rm{)}}$的解集${\mathit{\Omega }^ * }$非空。
定义 4 假设$\psi $是一个实值凸泛函，$\psi $的近似点算子定义为：
$$ {\rm{pro}}{{\rm{x}}_\sigma }[\psi ](z) = \arg {\min _v}\left\{ {\psi (v) + \frac{{\left\| {v - z} \right\|_2^2}}{{2\sigma }}} \right\} $$ (8) 本文提出的自适应临近点算法可求解式(7)。具体的算法实施步骤如下。
1) 初始化：$\sigma \leftarrow 1/\left\| \mathit{\boldsymbol{A}} \right\|$，$\tau \leftarrow 1/\left\| \mathit{\boldsymbol{A}} \right\|$，$\gamma = 1.2$，$\lambda = 5 \times {10^{ - 4}}$，$k = 1$，${t_1} = 1$。
2) 根据定义4给出预测步：
$$ \left\{ \begin{array}{l} {{\mathit{\boldsymbol{\tilde x}}}^k} = {\rm{pro}}{x_\sigma }[G]({\mathit{\boldsymbol{x}}^k} - \sigma {\mathit{\boldsymbol{A}}^{\rm{T}}}{\mathit{\boldsymbol{y}}^k})\\ {{\mathit{\boldsymbol{\bar x}}}^k} = {{\mathit{\boldsymbol{\tilde x}}}^k} + {\theta _k}({{\mathit{\boldsymbol{\tilde x}}}^k} - {\mathit{\boldsymbol{x}}^k})\\ ({{\mathit{\boldsymbol{\tilde p}}}^k}, {{\mathit{\boldsymbol{\tilde q}}}^k}) = {\rm{pro}}{x_\tau }[{F^*}]({\mathit{\boldsymbol{p}}^k} + \tau \mathit{\boldsymbol{H}}{{\mathit{\boldsymbol{\bar x}}}^k}, {\mathit{\boldsymbol{q}}^k} + \tau \nabla {{\mathit{\boldsymbol{\bar x}}}^k}) \end{array} \right. $$ 式中，${\theta _k} = \frac{{{t_{k - 1}} - 1}}{{{t_k}}}, {\rm{ }}{t_k} = \frac{1}{2}(1 + \sqrt {1 + 4t_{k - 1}^2} )$。
3) 校正步：对预测点进行校正产生下一个迭代点：
$$ {\mathit{\boldsymbol{u}}^{k + 1}} = {\mathit{\boldsymbol{u}}^k} - \gamma \alpha _k^*{\mathit{\boldsymbol{H}}^{ - 1}}{\mathit{\boldsymbol{M}}_k}({\mathit{\boldsymbol{u}}^k} - {\mathit{\boldsymbol{\tilde u}}^k}) $$ 式中，${\mathit{\boldsymbol{M}}_k} = \left( {_{ - \mathit{\boldsymbol{A}}{\rm{ }}\frac{{\rm{1}}}{\sigma }\mathit{\boldsymbol{E}}}^{\frac{1}{\tau }\mathit{\boldsymbol{E}}{\rm{ }} - {\theta _k}{\mathit{\boldsymbol{A}}^{\rm{T}}}}} \right)$；$\mathit{\boldsymbol{H}} = \left( {_{{\rm{ }}\frac{{\rm{1}}}{\sigma }\mathit{\boldsymbol{E}}}^{\frac{1}{\tau }\mathit{\boldsymbol{E}}}} \right)$，$\mathit{\boldsymbol{E}}$为单位矩阵。
预测步的一阶最优性条件是一个临近点算法(PPA)的结构，即若预测点${\mathit{\boldsymbol{\tilde u}}^k} = ({\mathit{\boldsymbol{\tilde x}}^k}, {\mathit{\boldsymbol{\tilde y}}^k})$是预测步产生的，则${\mathit{\boldsymbol{\tilde u}}^k} \in \mathit{\Omega }$且满足：
$$ {(\mathit{\boldsymbol{u'}} - {\mathit{\boldsymbol{\tilde u}}^k})^{\rm{T}}}(F({\mathit{\boldsymbol{\tilde u}}^k}) + {\mathit{\boldsymbol{M}}_k}\mathit{\boldsymbol{(}}{\mathit{\boldsymbol{\tilde u}}^k} - {\mathit{\boldsymbol{u}}^k})) \ge 0{\rm{ }}\forall \mathit{\boldsymbol{u'}} \in \mathit{\Omega } $$ (9) 文献[15]中给出原始对偶过程的步长τ和σ的选取满足条件为：
$$ \tau \sigma \frac{{{{(1 + {\theta _k})}^2}}}{4}||{\mathit{\boldsymbol{A}}^{\rm{T}}}\mathit{\boldsymbol{A}}|| < 1{\rm{ }}{\theta _k} \in ( - 1, 1) $$ (10) 在文献[11]的基础上提出简化算法，方便、快捷地求出$||\mathit{\boldsymbol{A}}{\rm{|| = }}||\mathit{\boldsymbol{H}}{\rm{, }}\nabla ||$，设定算法运行N步结束，随着N的增大$s \to ||\mathit{\boldsymbol{A}}{\rm{||}}$，具体的算法执行过程如下。
1) 初始化：输入一张非零的图像${\mathit{\boldsymbol{x}}_0}, 0 \leftarrow n$；
2) 循环：
$$ \left\{ \begin{array}{l} {\mathit{\boldsymbol{x}}_{n + 1}} = {\mathit{\boldsymbol{K}}^{\rm{T}}}\mathit{\boldsymbol{K}}{\mathit{\boldsymbol{x}}_n} - {\rm{div}}\nabla {\mathit{\boldsymbol{x}}_n}\\ {\mathit{\boldsymbol{x}}_{n + 1}} = {{{\mathit{\boldsymbol{x}}_{n + 1}}} {\left/ { { {||{\mathit{\boldsymbol{x}}_{n + 1}}|{|_2}}}} \right. } }\\ s = \sqrt {||\mathit{\boldsymbol{K}}{\mathit{\boldsymbol{x}}_{n + 1}}||_2^2 + ||\nabla {\mathit{\boldsymbol{x}}_{n + 1}}||_2^2} \end{array} \right. $$ 直到$n \ge N$。
-
在收缩算法的框架下，给出参数的选取满足原始对偶过程条件下的自适应临近点算法的收敛性的分析。
引理 1 若自适应临近点算法中的参数$\sigma > 0, {\rm{ }}\tau > 0, {\rm{ }}{\theta _k} \in ( - 1, 1)$满足式(10)中的原始对偶过程条件，则${\mathit{\boldsymbol{M}}_k}$是一致正定的，即：
$$ {(\mathit{\boldsymbol{u}} - \mathit{\boldsymbol{\tilde u}})^{\rm{T}}}{\mathit{\boldsymbol{M}}_k}(\mathit{\boldsymbol{u}} - \mathit{\boldsymbol{\tilde u}}) \ge \frac{{{\delta _k}}}{{1 + {\delta _k}}}||\mathit{\boldsymbol{u}} - \mathit{\boldsymbol{\tilde u}}||_H^2 > 0{\rm{ }}\forall \mathit{\boldsymbol{u}} \ne \mathit{\boldsymbol{\tilde u}} $$ 式中，${\rm{ }}{\delta _k} = \frac{2}{{1 + {\theta _k}}}\sqrt {\frac{1}{{\tau \sigma ||{\mathit{\boldsymbol{A}}^{\rm{T}}}\mathit{\boldsymbol{A}}||}}} - 1 > 0$。
证明：由${\rm{ }}{\delta _k} = \frac{2}{{1 + {\theta _k}}}\sqrt {\frac{1}{{\tau \sigma ||{\mathit{\boldsymbol{A}}^{\rm{T}}}\mathit{\boldsymbol{A}}||}}} - 1 > 0$可得：
$$ \tau (1 + {\delta _k})||{\mathit{\boldsymbol{A}}^{\rm{T}}}\mathit{\boldsymbol{A}}||\frac{{(1 - {\theta _k})}}{4} = \frac{1}{{\sigma (1 + {\delta _k})}} $$ 因此，对任意的$\mathit{\boldsymbol{u}} \ne \mathit{\boldsymbol{\tilde u}}$有：
$$ {(\mathit{\boldsymbol{u}} - \mathit{\boldsymbol{\tilde u}})^{\rm{T}}}{\mathit{\boldsymbol{M}}_k}(\mathit{\boldsymbol{u}} - \mathit{\boldsymbol{\tilde u}}) = ||\mathit{\boldsymbol{u}} - \mathit{\boldsymbol{\tilde u}}||_H^2 - $$ $$ (1 + {\theta _k}){(\mathit{\boldsymbol{y}} - \mathit{\boldsymbol{\tilde y}})^{\rm{T}}}A(\mathit{\boldsymbol{x}} - \mathit{\boldsymbol{\tilde x}}) = $$ $$ ||\mathit{\boldsymbol{u}} - \mathit{\boldsymbol{\tilde u}}||_H^2 - \frac{1}{{\sigma (1 + {\delta _k})}}||\mathit{\boldsymbol{y}} - \mathit{\boldsymbol{\tilde y}}|{|^2} - $$ $$ \begin{array}{c} \frac{1}{{\tau (1 + {\delta _k})||{\mathit{\boldsymbol{A}}^{\rm{T}}}\mathit{\boldsymbol{A}}||}}||A(\mathit{\boldsymbol{x}} - \mathit{\boldsymbol{\tilde x}})|{|^2} \ge ||\mathit{\boldsymbol{u}} - \mathit{\boldsymbol{\tilde u}}||_H^2 - \\ \frac{1}{{1 + {\delta _k}}}\left( {\frac{1}{\tau }||\mathit{\boldsymbol{x}} - \mathit{\boldsymbol{\tilde x}}|{|^2}{\rm{ + }}\frac{1}{\sigma }{\rm{||}}\mathit{\boldsymbol{y}} - \mathit{\boldsymbol{\tilde y}}|{|^2}} \right){\rm{ = }}\frac{{{\delta _k}}}{{1 + {\delta _k}}}||\mathit{\boldsymbol{u}} - \mathit{\boldsymbol{\tilde u}}||_H^2 \end{array} $$ 由于$||\mathit{\boldsymbol{u}} - \mathit{\boldsymbol{\tilde u}}||_H^2 > 0$及$\delta > 0$，因此${\mathit{\boldsymbol{M}}_k}$是一致正定的。
引理 2 ${\mathit{\boldsymbol{M}}_k}$和H的定义如算法1校正步中所示，则${\rm{||}}\mathit{\boldsymbol{M}}_k^{\rm{T}}{\mathit{\boldsymbol{H}}^{ - 1}}{\mathit{\boldsymbol{M}}_k}||$一致有界，即存在$m > 0$，满足$|{\rm{|}}\mathit{\boldsymbol{M}}_k^{\rm{T}}{\mathit{\boldsymbol{H}}^{ - 1}}{\mathit{\boldsymbol{M}}_k}|| \le m$。
证明：
$$ \mathit{\boldsymbol{M}}_k^{\rm{T}}{\mathit{\boldsymbol{H}}^{ - 1}}{\mathit{\boldsymbol{M}}_k} = \left( {\begin{array}{*{20}{c}} {\frac{1}{\tau }\mathit{\boldsymbol{E}} + \sigma {\mathit{\boldsymbol{A}}^{\rm{T}}}\mathit{\boldsymbol{A}}}&{ - ({\theta _k} + 1){\mathit{\boldsymbol{A}}^{\rm{T}}}}\\ { - ({\theta _k} + 1)\mathit{\boldsymbol{A}}}&{\frac{1}{\sigma }\mathit{\boldsymbol{E}} + \tau \theta _k^2\mathit{\boldsymbol{A}}{\mathit{\boldsymbol{A}}^{\rm{T}}}} \end{array}} \right) $$ 又由于$\mathop {\lim }\limits_{k \to \infty } {\theta _k} = 1$，易得：
$$ \mathop {\lim }\limits_{k \to \infty } \mathit{\boldsymbol{M}}_k^{\rm{T}}{\mathit{\boldsymbol{H}}^{ - 1}}{\mathit{\boldsymbol{M}}_k} = \left( {\begin{array}{*{20}{c}} {\frac{1}{\tau }\mathit{\boldsymbol{E}} + \sigma {\mathit{\boldsymbol{A}}^{\rm{T}}}\mathit{\boldsymbol{A}}}&{ - 2{\mathit{\boldsymbol{A}}^{\rm{T}}}}\\ { - 2\mathit{\boldsymbol{A}}}&{\frac{1}{\sigma }\mathit{\boldsymbol{E}} + \tau \mathit{\boldsymbol{A}}{\mathit{\boldsymbol{A}}^{\rm{T}}}} \end{array}} \right) $$ 记，由范数的连续性可得：
$$ \mathop {\lim }\limits_{k \to \infty } ||\mathit{\boldsymbol{M}}_k^{\rm{T}}{\mathit{\boldsymbol{H}}^{ - 1}}{\mathit{\boldsymbol{M}}_k}|| = ||\mathit{\boldsymbol{M}}|| $$ 则存在$N > 0$，当$k > N$有：
$$ \mathop {\lim }\limits_{k \to \infty } \left\| {\mathit{\boldsymbol{M}}_k^{\rm{T}}{\mathit{\boldsymbol{H}}^{ - 1}}{\mathit{\boldsymbol{M}}_k}} \right\| \le \left\| \mathit{\boldsymbol{M}} \right\| + 1 $$ 令$m = \max \{ ||\mathit{\boldsymbol{M}}_1^{\rm{T}}{\mathit{\boldsymbol{H}}^{ - 1}}{\mathit{\boldsymbol{M}}_1}||, \cdots ||\mathit{\boldsymbol{M}}_N^{\rm{T}}{\mathit{\boldsymbol{H}}^{ - 1}}{\mathit{\boldsymbol{M}}_N}||, ||\mathit{\boldsymbol{M}} + 1||\} $，则引理得证。
定理 1 对任意$\{ {\mathit{\boldsymbol{u}}^*}\} \in {\mathit{\Omega }^*}$，由自适应临近点生成的序列$\{ {\mathit{\boldsymbol{u}}^k}\} $满足不等式：
$$ \left\| {{\mathit{\boldsymbol{u}}^{k + 1}} - {\mathit{\boldsymbol{u}}^*}} \right\|_H^2 \le \left\| {{\mathit{\boldsymbol{u}}^k} - {\mathit{\boldsymbol{u}}^*}} \right\|_H^2 - C{\left\| {{\mathit{\boldsymbol{u}}^k} - {{\mathit{\boldsymbol{\tilde u}}}^k}} \right\|^2} $$ 式中，$C > 0$是个正常数。
证明：由式(6)和式(9)得：
$$ \begin{array}{c} ||{\mathit{\boldsymbol{u}}^{k + 1}} - {\mathit{\boldsymbol{u}}^*}||_H^2{\rm{ }} = \\ {\rm{ }}||({\mathit{\boldsymbol{u}}^k} - {\mathit{\boldsymbol{u}}^*}) - {\alpha _k}{\mathit{\boldsymbol{H}}^{ - 1}}{\mathit{\boldsymbol{M}}_k}({\mathit{\boldsymbol{u}}^k} - {{\mathit{\boldsymbol{\tilde u}}}^k})||_H^2{\rm{ = }}\\ ||{\mathit{\boldsymbol{u}}^k} - {\mathit{\boldsymbol{u}}^*}||_H^2 + 2{\alpha _k}{({\mathit{\boldsymbol{u}}^k} - {\mathit{\boldsymbol{u}}^*})^{\rm{T}}}{\mathit{\boldsymbol{M}}_k}({\mathit{\boldsymbol{u}}^k} - {{\mathit{\boldsymbol{\tilde u}}}^k}) + \\ \alpha _k^2||{\mathit{\boldsymbol{H}}^{ - 1}}{\mathit{\boldsymbol{M}}_k}({\mathit{\boldsymbol{u}}^k} - {{\mathit{\boldsymbol{\tilde u}}}^k})||_H^2 \le \\ ||\mathit{\boldsymbol{u}} - {\mathit{\boldsymbol{u}}^*}||_H^2 - 2{\alpha _k}{({\mathit{\boldsymbol{u}}^k} - {{\mathit{\boldsymbol{\tilde u}}}^k})^{\rm{T}}}{\mathit{\boldsymbol{M}}_k}({\mathit{\boldsymbol{u}}^k} - {{\mathit{\boldsymbol{\tilde u}}}^k}) + \\ \alpha _k^2||{\mathit{\boldsymbol{H}}^{ - 1}}{\mathit{\boldsymbol{M}}_k}({\mathit{\boldsymbol{u}}^k} - {{\mathit{\boldsymbol{\tilde u}}}^k})||_H^2 = \\ ||\mathit{\boldsymbol{u}} - {\mathit{\boldsymbol{u}}^*}||_H^2 - \gamma (2 - \gamma )\alpha _k^*{({\mathit{\boldsymbol{u}}^k} - {{\mathit{\boldsymbol{\tilde u}}}^k})^{\rm{T}}}\mathit{\boldsymbol{M}}({\mathit{\boldsymbol{u}}^k} - {{\mathit{\boldsymbol{\tilde u}}}^k}) \end{array} $$ 由引理1和引理2可知分别存在正常数J和L满足：
$$ \begin{array}{c} {({\mathit{\boldsymbol{u}}^k} - {{\mathit{\boldsymbol{\tilde u}}}^k})^{\rm{T}}}{\mathit{\boldsymbol{M}}_k}({\mathit{\boldsymbol{u}}^k} - {{\mathit{\boldsymbol{\tilde u}}}^k}) \ge J||{\mathit{\boldsymbol{u}}^k} - {{\mathit{\boldsymbol{\tilde u}}}^k}|{|^2}\\ ||\mathit{\boldsymbol{M}}_k^{\rm{T}}{\mathit{\boldsymbol{H}}^{ - 1}}{\mathit{\boldsymbol{M}}_k}|| \le L \end{array} $$ 因此有：
$$ \begin{array}{c} \alpha _k^*{({\mathit{\boldsymbol{u}}^k} - {{\mathit{\boldsymbol{\tilde u}}}^k})^{\rm{T}}}{\mathit{\boldsymbol{M}}_k}({\mathit{\boldsymbol{u}}^k} - {{\mathit{\boldsymbol{\tilde u}}}^k}) = \\ \frac{{{{({\mathit{\boldsymbol{u}}^k} - {{\mathit{\boldsymbol{\tilde u}}}^k})}^{\rm{T}}}{\mathit{\boldsymbol{M}}_k}({\mathit{\boldsymbol{u}}^k} - {{\mathit{\boldsymbol{\tilde u}}}^k})}}{{||{\mathit{\boldsymbol{H}}^{ - 1}}{\mathit{\boldsymbol{M}}_k}({\mathit{\boldsymbol{u}}^k} - {{\mathit{\boldsymbol{\tilde u}}}^k})||_H^2}}{({\mathit{\boldsymbol{u}}^k} - {{\mathit{\boldsymbol{\tilde u}}}^k})^{\rm{T}}}{\mathit{\boldsymbol{M}}_k}({\mathit{\boldsymbol{u}}^k} - {{\mathit{\boldsymbol{\tilde u}}}^k}) \ge \\ \frac{{{J^{\rm{2}}}||{\mathit{\boldsymbol{u}}^k} - {{\mathit{\boldsymbol{\tilde u}}}^k}|{|^4}}}{{||\mathit{\boldsymbol{M}}_k^{\rm{T}}{\mathit{\boldsymbol{H}}^{ - 1}}{\mathit{\boldsymbol{M}}_k}||||{\mathit{\boldsymbol{u}}^k} - {{\mathit{\boldsymbol{\tilde u}}}^k}|{|^2}}} \ge \frac{{{J^{\rm{2}}}||{\mathit{\boldsymbol{u}}^k} - {{\mathit{\boldsymbol{\tilde u}}}^k}|{|^2}}}{L} \end{array} $$ 令$C = \frac{{\gamma (2 - \gamma ){J^{\rm{2}}}}}{L}$则定理得证。
上面的定理说明由自适应临近点算法生成的序列$\{ {\mathit{\boldsymbol{u}}^k}\} $是Fejèr单调的。因此本文提出的算法在$\mathit{\boldsymbol{H}}$模下是一个收缩算法，下面给出算法的主要收敛性定理。
定理2由自适应临近点算法产生的序列$\{ {\mathit{\boldsymbol{u}}^k}\} $收敛到变分不等式${\rm{VI(}}\mathit{\Omega }, F{\rm{)}}$的解${\mathit{\boldsymbol{u}}^ * }$。
证明：由定理1可知$\{ {\mathit{\boldsymbol{u}}^k}\} $有界，则易得$\mathop {\lim }\limits_{k \to \infty } ||{\mathit{\boldsymbol{u}}^k} - {\mathit{\boldsymbol{\tilde u}}^k}|| = 0$
由此可知$\{ {\mathit{\boldsymbol{\tilde u}}^k}\} $有界。令${\mathit{\boldsymbol{u}}^*}$是$\{ {\mathit{\boldsymbol{\tilde u}}^k}\} $的一个聚点且存在一个子列$\{ \mathit{\boldsymbol{\tilde u}}_j^k\} $收敛到${\mathit{\boldsymbol{u}}^*}$。则对子列$\mathit{\boldsymbol{\tilde u}}_j^k$可得：
$$ {(\mathit{\boldsymbol{u}} - \mathit{\boldsymbol{\tilde u}}_j^k)^{\rm{T}}}\{ F(\mathit{\boldsymbol{\tilde u}}_j^k) + {\mathit{\boldsymbol{M}}_k}(\mathit{\boldsymbol{u}}_j^k - \mathit{\boldsymbol{\tilde u}}_j^k)\} \ge 0 $$ 式中，${\mathit{\boldsymbol{\tilde u}}^k} \in \mathit{\Omega }$，$\forall \mathit{\boldsymbol{u}} \in \mathit{\Omega }$，又因为$\{ \mathit{\boldsymbol{\tilde u}}_j^k\} \to {\mathit{\boldsymbol{u}}^*}$及$\mathop {\lim }\limits_{_{j \to \infty }} ||\mathit{\boldsymbol{u}}_j^k - \mathit{\boldsymbol{\tilde u}}_j^k|| = 0$，则有：
$$ {\rm{ (}}\mathit{\boldsymbol{u}} - {\mathit{\boldsymbol{u}}^*}{{\rm{)}}^{\rm{T}}}F({\mathit{\boldsymbol{u}}^*}) \ge 0 $$ 式中，${\mathit{\boldsymbol{u}}^*} \in \mathit{\Omega }$，$\forall \mathit{\boldsymbol{u}} \in \mathit{\Omega }$，故${\mathit{\boldsymbol{u}}^ * }$是变分不等式${\rm{VI(}}\mathit{\Omega }, F{\rm{)}}$的解，由于变分不等式${\rm{VI(}}\mathit{\Omega }, F{\rm{)}}$的所有解都满足定理1，因此：
$$ ||{\mathit{\boldsymbol{u}}^{k + 1}} - {\mathit{\boldsymbol{u}}^*}|{|^2} \le ||{\mathit{\boldsymbol{u}}^k} - {\mathit{\boldsymbol{u}}^*}|{|^2}{\rm{ }}\forall k > 0 $$ 故定理得证。
-
为验证自适应临近点算法在CT图像恢复中的有效性，选取Shepp-Logan数值体模[16]进行定量和定性分析。用均方误差和信噪比来量化重建的效果，从而衡量一个算法的优劣。所有算法都是由Matlab软件实现并且设定100步终止。若信噪比${\rm{SN}}{{\rm{R}}_x} = 20{\lg _{10}}\frac{{\left\| {{\mathit{\boldsymbol{x}}_{{\rm{true}}}}} \right\|}}{{\left\| {\mathit{\boldsymbol{\hat x}} - {\mathit{\boldsymbol{x}}_{{\rm{true}}}}} \right\|}}$越大，同时均方误差${\rm{MS}}{{\rm{E}}_x} = \frac{{||\mathit{\boldsymbol{\hat x}} - {\mathit{\boldsymbol{x}}_{{\rm{true}}}}||}}{{||{\mathit{\boldsymbol{x}}_{{\rm{true}}}}||}}$越小表明重建的图像越好，其中$\mathit{\boldsymbol{\hat x}}$是由算法重建出来的图像，${\mathit{\boldsymbol{x}}_{{\rm{true}}}}$是原始的图像。
图 1展示了由滤波反投影和自适应临近点算法分别从60个角度重建出来的图像，体模大小设定为256×256。滤波反投影(filtered back projection, FBP)类算法属于经典的解析重建算法之一，被广泛应用于CT重建中[3]。其基本思想是先对投影数据进行滤波，再对滤波后的投影数据进行反投影运算从而得到重建图像。FBP算法和APPA算法的MSE分别为0.202 1和0.177 9，而SNR分别为13.889 0和14.994 7。与FBP算法相比，APPA方法信噪比提高了7.96%，而相应的均方误差下降了11.9%。可看出自适应临近点算法重建出来的图像信息恢复较为完整，细节清晰，去掉了较多的伪影，具备实用性和有效性。
图 1 两种算法从60个角度重建出来的图像
图 2展示了自适应临近点算法与基于投影数据全广义变分最小化的低剂量CT重建(TGV)算法[16]从180个角度重建出来的图像，体模大小设定为512×512。与APPA算法类似，TGV算法中也根据对偶原理将系统模型转化为鞍点问题。在图中，APPA算法和TGV算法的SNR分别为24.612 7和24.035 2，APPA算法的信噪比提高了2.40%。可见本文提出的方法在提高图像重建的信噪比方面有较好表现。
图 2 两种算法从180个角度重建出来的图像
-
本文主要讨论了稀疏角度CT图像重建问题中的自适应临近点算法，在原始对偶过程条件下建立了其全局收敛性结果。数值实验结果表明，本文的自适应临近点算法对去除稀疏角度CT图像重建中图像的噪声和条形伪影方面具有较佳的表现。
An Adaptive Proximal Point Algorithm for Sparse-View CT Image Reconstruction
-
摘要: 针对全变差正则化的稀疏角度CT图像重建问题，提出了一种自适应临近点算法。该算法在临近点迭代的每一步中自适应地选取临近参数矩阵，且该矩阵是变动且不对称的。在收缩算法的框架和一定条件下可建立该算法的全局收敛性。对广泛使用的滤波反投影算法和自适应临近点算法进行了二维的仿真数据实验，实验结果表明自适应临近点算法在图像重建中是有效实用的。Abstract: This paper presents an adaptive proximal point algorithm (APPA) for sparse-view computed tomography (CT) image reconstruction based on total variation regularization. The proposed algorithm chooses an adaptive proximal parameter matrix which is neither necessary symmetric nor constant in each iteration. By using the framework of contraction method, the global convergence result could be established for the proposed algorithm under suitable conditions. Numerical results for 2-D CT reconstruction from simulated digital Shepp-Logan phantom data show that APPA method is effective and practical.
-
[1] 高河伟, 张丽, 陈志强, 等.有限角度CT图像重建算法综述[J]. CT理论与应用研究, 2006, 15(1):46-50. doi: 10.3969/j.issn.1004-4140.2006.01.011 GAO He-wei, ZHANG Li, CHEN Zhi-qiang, et al. Reviews of image reconstruction from limited-angle[J]. Computerized Tomography Theory and Applications, 2006, 15(1):46-50. doi: 10.3969/j.issn.1004-4140.2006.01.011 [2] HSIEH J, 张朝宗.计算机断层成像技术:原理, 设计, 伪像和进展[M].北京:科学出版社, 2006. HSIEH J, ZHANG Chao-zong. Computer tomography:Principle, design, artifacts and recent advances[M]. Beijing:Science Press, 2006. [3] YU Han, BIN Yan, LEI Li, et al. Rebinned filtered back-projection reconstruction from truncated data for half-covered helical cone-beam computed tomography[J]. IEEE Transactions on Nuclear Science, 2014, 61(5):2753-2763. doi: 10.1109/TNS.2014.2351817 [4] NASSI M, BRODY W R, MEDOFF B P, et al. Iterative reconstruction-reprojection:an algorithm for limited data cardiac-computed tomography[J]. IEEE Transactions on Biomedical Engineering, 1982, 29(5):333-341. doi: 10.1109-TBME.1982.324900/ [5] HERMAN G T, MEYER L B. Algebraic reconstruction techniques can be made computationally efficient[positron emission tomography application] [J]. IEEE Transactions on Medical Imaging, 1993, 12(3):600-609. doi: 10.1109/42.241889 [6] FANG Li, JUAN F P, ABASCAL J. et al. Total variation regularization with split Bregman-based method in magnetic induction tomography using experimental data[J]. IEEE Sensors Journal, 2017, 17(4):976-985. doi: 10.1109/JSEN.2016.2637411 [7] PERSSON M, BONE D, ELMQVIST H. Total variation norm for three-dimensional iterative reconstruction in limited view angle tomography[J]. Physics in Medicine and Biology, 2001, 46(3):853-866. doi: 10.1088/0031-9155/46/3/318 [8] SIDKY E Y, PAN X. Image reconstruction in circular cone-beam computed tomography by constrained, total-variation minimization[J]. Physics in Medicine and Biology, 2008, 53(17):4777-4807. doi: 10.1088/0031-9155/53/17/021 [9] SIDKY E Y, CHARTRAND R, BOONE J M, et al. Constrained minimization for enhanced exploitation of gradient sparsity:Application to CT image reconstruction[J]. IEEE Journal of Translational Engineering in Health and Medicine, 2014, 2:1-18. http://ieeexplore.ieee.org/xpl/articleDetails.jsp?arnumber=6714374& [10] FESSLER J A. Penalized weighted least-squares image reconstruction for positron emission tomography[J]. IEEE Transactions on Medical Imaging, 1994, 13(2):290-300. doi: 10.1109/42.293921 [11] SIDKY E Y, JØRGENSEN J H, PAN X. Convex optimization problem prototyping for image reconstruction in computed tomography with the Chambolle-Pock algorithm[J]. Physics in Medicine and Biology, 2012, 57(10):3065-3091. doi: 10.1088/0031-9155/57/10/3065 [12] CHAMBOLLE A, POCK T. A first-order primal-dual algorithm for convex problems with applications to imaging[J]. Journal of Mathematical Imaging and Vision, 2011, 40(1):120-145. doi: 10.1007/s10851-010-0251-1 [13] 蔺鲁萍, 王永革.不完全角度CT图像重建的模型与算法[J].北京航空航天大学学报, 2017, 43(4):823-830. http://d.old.wanfangdata.com.cn/Periodical/bjhkhtdxxb201704023 LIN Lu-ping, WANG Yong-ge. CT image reconstruction model and algorithm from few views[J]. Journal of Beijing University of Aeronautics and Astronautics, 2017, 43(4):823-830. http://d.old.wanfangdata.com.cn/Periodical/bjhkhtdxxb201704023 [14] CHEN Ying, WU Jian, YU Gao-hang. Adaptive proximal point algorithms for total variation image restoration[J]. Statistics, Optimization & Information Computing, 2015, 3(1):15-29. http://www.iapress.org/index.php/soic/article/download/20150302/180 [15] HE B, YUAN X. Convergence analysis of primal-dual algorithms for a saddle-point problem:From contraction perspective[J]. SIAM Journal on Imaging Sciences, 2012, 5(1):119-149. doi: 10.1137/100814494 [16] 牛善洲, 吴恒, 喻泽峰, 等.基于投影数据全广义变分最小化的低剂量CT重建[J].南方医科大学学报, 2017, 37(12):1585-1591. doi: 10.3969/j.issn.1673-4254.2017.12.04 NIU Shan-zhou, WU Heng, YU Ze-feng, et al. Total generalized variation minimization based on projection data for low-dose CT reconstruction[J]. Journal of Southern Medical University, 2017, 37(12):1585-1591. doi: 10.3969/j.issn.1673-4254.2017.12.04 -
点击查看大图
图(2)
计量
- 文章访问数: 4070
- HTML全文浏览量: 1288
- PDF下载量: 91
- 被引次数: 0
