Undersampled CS image reconstruction using nonconvex nonsmooth mixed constraints

Abstract

Compressed sensing magnetic resonance imaging (CS-MRI) has attracted considerable attention due to its great potential in reducing scanning time and guaranteeing high-quality reconstruction. In conventional CS-MRI framework, the total variation (TV) penalty and L1-norm constraint on wavelet coefficients are commonly combined to reduce the reconstruction error. However, TV sometimes tends to cause staircase-like artifacts due to its nature in favoring piecewise constant solution. To overcome the model-dependent deficiency, a hybrid TV (TV^1,2) regularizer is introduced in this paper by combining TV with its second-order version (TV²). It is well known that the wavelet coefficients of MR images are not only approximately sparse, but also have the property of tree-structured hierarchical sparsity. Therefore, a L0-regularized tree-structured sparsity constraint is proposed to better represent the measure of sparseness in wavelet domain. In what follows, we present our new CS-MRI framework by combining the TV^1,2 regularizer and L0-regularized tree-structured sparsity constraint. However, the combination makes CS-MRI problem difficult to handle due to the nonconvex and nonsmooth natures of mixed constraints. To achieve solution stability, the resulting composite minimization problem is decomposed into several simpler subproblems. Each of these subproblems has a closed-form solution or could be efficiently solved using existing numerical method. The results from simulation and in vivo experiments have demonstrated the good performance of our proposed method compared with several conventional MRI reconstruction methods.

Appendix

In this Appendix, we describe how the proximal function (20) associated with TV and TV² denoising can be efficiently solved using a fast projection method.

1.1 A.1 Fast projection method for tv denoising

The subproblem $\mathbf {f}_{\text {TV}} = \text {prox}_{\pi } \left (2 \alpha \rho \left \| \mathbf {f} \right \|_{\text {TV}} \right ) \left (\bar {\mathbf {f}} \right )$ is equivalent to solving the following unconstrained optimization problem

$$ \mathbf{f}_{\text{TV}} = \min_{\mathbf{f}} \left\{ \left\| \partial_{h}\mathbf{f} \right\|_{1} + \left\| \partial_{v}\mathbf{f} \right\|_{1} + \frac{\pi}{2} \left\| {\mathbf{f}} - \bar{\mathbf{f}} \right\|_{2}^{2} \right\}, $$

(21)

where $\pi = \frac {1}{2 \alpha \pi \rho }$. The TV-based optimization problem (21) remains difficult to solve due to its non-smooth and non-differential natures. The split Bregman technique [21] has been widely used to solve this unconstrained optimization problem. To further speed up the algorithm, Jia and Zhao [33] proposed a fast algorithm based on split Bregman method for the solution of (21). This fast algorithm referred to as fast projection method is shown in Algorithm 3.

The projection method could achieve satisfactory denoising performance with convergence rate $\mathcal {O} (k^{-1} )$. For the sake of illustration, we use the two-dimensional version of image f with size m × n in this subsection. The first order discrete differential operator ∂_hf is given by $\left (\partial _{h} \mathbf {f} \right )_{i,1} = 0$ for i = 1,⋯ ,m and $\left (\partial _{h} \mathbf {f} \right )_{i,j} = \mathbf {f}_{i,j} - \mathbf {f}_{i,j-1}$ for i = 1,⋯ ,m and j = 2,⋯ ,n. Analogous to the discrete version of ∂_hf, ∂_vf is given by $\left (\partial _{v} \mathbf {f} \right )_{1,j} = 0$ for j = 1,⋯ ,n and $\left (\partial _{v} \mathbf {f} \right )_{i,j} = \mathbf {f}_{i,j} - \mathbf {f}_{i-1,j}$ for i = 2,⋯ ,m and j = 1,⋯ ,n. For $\mathbf {w} \in \mathcal {R}^{m \times n}$, the conjugate operators of ∂_hw and ∂_vw are defined respectively as

$$\left( {\partial_{h}^{T}} \mathbf{w} \right)_{i,j} = \left\{\begin{array}{lll} -\mathbf{w}_{i,2} & \text{if}~ j = 1,\\ \mathbf{w}_{i,j} - \mathbf{w}_{i,j + 1} & \text{if}~1 < j < n, \\ \mathbf{w}_{i,j} & \text{if}~ j = n. \end{array}\right. $$

and

$$\left( {\partial_{v}^{T}} \mathbf{w} \right)_{i,j} = \left\{\begin{array}{llll} -\mathbf{w}_{2,j} & \text{if}~ i = 1,\\ \mathbf{w}_{i,j} - \mathbf{w}_{i + 1,j} & \text{if}~ 1 < i < m, \\ \mathbf{w}_{i,j} & \text{if}~ i = m. \end{array}\right. $$

For τ > 0 and $\mathbf {c} \in \mathcal {R}^{m \times n}$, the $\text {Cut} \left (\cdot , \cdot \right )$ operator is defined as

$$ \text{Cut} \left( \mathbf{c}, \frac{1}{\tau} \right) = \mathbf{c} - \frac{\mathbf{c}}{\left| \mathbf{c} \right|} \circ \max \left\{ \left| \mathbf{c} \right| - \frac{1}{\tau}, 0 \right\}. $$

(22)

According to the $\text {Shrinkage} \left (\cdot , \cdot \right )$ operator defined in (19), it is clear that $\text {Cut} \left (\mathbf {c}, \frac {1}{\tau } \right ) + \text {Shrinkage} \left (\mathbf {c}, \frac {1}{\tau } \right ) = \mathbf {c}$.

1.2 A.2 Fast projection method for TV² denoising

The subproblem $\mathbf {f}_{\text {TV}^{2}} = \text {prox}_{\pi } \left (2 \alpha \left (1 - \rho \right ) \left \| \mathbf {f} \right \|_{\text {TV}^{2}} \right ) \left (\bar {\mathbf {f}} \right )$ is equivalent to solving the following unconstrained optimization problem

$$ \mathbf{f}_{\text{TV}^{2}} = \min_{\mathbf{f}} \left\{ \sum\nolimits_{{\Lambda} \in {\Theta}} \left\| \partial_{{\Lambda}} \mathbf{f} \right\|_{1} + \frac{\pi}{2} \left\| {\mathbf{f}} - \bar{\mathbf{f}} \right\|_{2}^{2} \right\}, $$

(23)

where ${\Theta } = \left \{ hh, hv, vh, vv \right \}$ and $\pi = \frac {1}{2 \alpha \pi (1 - \rho )}$. Analogous to the fast projection method for TV denoising, the pseudocode of fast projection method for TV² denoising is shown in Algorithm 4.

We now define the second order discrete differential operators (∂_hhf, ∂_hvf, ∂_vhf and ∂_vvf) as follows

$$\left( \partial_{hh} \mathbf{f} \right)_{i,j} = \left\{\begin{array}{llll} \mathbf{f}_{i,n} - 2\mathbf{f}_{i,j} + \mathbf{f}_{i,j + 1} & \text{if}~ 1 \leq i \leq m, j = 1,\\ \mathbf{f}_{i,j-1} - 2\mathbf{f}_{i,j} + \mathbf{f}_{i,j + 1} & \text{if}~ 1 \leq i \leq m, 1 < j < n,\\ \mathbf{f}_{i,j-1} - 2\mathbf{f}_{i,j} + \mathbf{f}_{i,1} & \text{if}~ 1 \leq i \leq m, j = n. \end{array}\right. $$

$$\left( \partial_{vv} \mathbf{f} \right)_{i,j} = \left\{\begin{array}{lll} \mathbf{f}_{m,j} - 2\mathbf{f}_{i,j} + \mathbf{f}_{i + 1,j} & \text{if}~ i = 1, 1 \leq j \leq n,\\ \mathbf{f}_{i-1,j} - 2\mathbf{f}_{i,j} + \mathbf{f}_{i + 1,j} & \text{if}~ 1 < i < m, 1 \leq j \leq n,\\ \mathbf{f}_{i-1,j} - 2\mathbf{f}_{i,j} + \mathbf{f}_{1,j} & \text{if}~ i = m, 1 \leq j \leq n. \end{array}\right. $$

and

$$\left( \partial_{hv} \mathbf{f} \right)_{i,j} = \left\{\begin{array}{lll} \mathbf{f}_{i,j} - \mathbf{f}_{i + 1,j} - \mathbf{f}_{i,j + 1} + \mathbf{f}_{i + 1,j + 1} & \text{if}~ 1 \leq i < m, 1 \leq j < n,\\ \mathbf{f}_{i,j} - \mathbf{f}_{1,j} - \mathbf{f}_{i,j + 1} + \mathbf{f}_{1,j + 1} & \text{if}~ i = m, 1 \leq j < n,\\ \mathbf{f}_{i,j} - \mathbf{f}_{i + 1,j} - \mathbf{f}_{i,1} + \mathbf{f}_{i + 1,1} & \text{if}~ 1 \leq i < m, j = n,\\ \mathbf{f}_{i,j} - \mathbf{f}_{1,j} - \mathbf{f}_{i,1} + \mathbf{f}_{1,1} & \text{if}~ i = m, j = n. \end{array}\right. $$

where ∂_hvf = ∂_vhf. Let $\partial _{hh}^{T}$, $\partial _{vv}^{T}$, $\partial _{hv}^{T}$ and $\partial _{vh}^{T}$ be conjugate operators of ∂_hh, ∂_vv, ∂_hv and ∂_vh, we then have $\partial _{hh}^{T} = \partial _{hh}$, $\partial _{vv}^{T} = \partial _{vv}$ and $\partial _{hv}^{T} = \partial _{vh}^{T}$ with

$$\left( \partial_{hv}^{T} \mathbf{w} \right)_{i,j} = \left\{\begin{array}{llll} \mathbf{w}_{i,j} - \mathbf{w}_{i,n} - \mathbf{w}_{m,j} + \mathbf{w}_{m,n} &\text{if}~ i = 1, j = 1,\\ \mathbf{w}_{i,j} - \mathbf{w}_{i,j-1} - \mathbf{w}_{m,j} + \mathbf{w}_{m,j-1} &\text{if}~ i = 1, 1 <j \leq n,\\ \mathbf{w}_{i,j} - \mathbf{w}_{i-1,j} - \mathbf{w}_{i,n} + \mathbf{w}_{i-1,n} &\text{if}~ 1 < i \leq m, j = 1,\\ \mathbf{w}_{i,j} - \mathbf{w}_{1,j-1} - \mathbf{w}_{i-1,j} + \mathbf{w}_{i-1,j-1} &\text{if}~ 1 < i \leq m, 1 < j \leq n. \end{array}\right. $$