Compressive total variation for image reconstruction and restoration

Abstract

In this paper, we make use of the fact that the matrix $\mathit{u}$ is (approximately) low-rank in image inpainting, and the corresponding gradient transform matrices ${\mathcal{D}}_x\mathit{u},{\mathcal{D}}_y\mathit{u}$ are sparse in image reconstruction and restoration. Therefore we consider that these gradient matrices ${\mathcal{D}}_x\mathit{u},{\mathcal{D}}_y\mathit{u}$ also are (approximately) low-rank, and also verify it by numerical test and theoretical analysis. We propose a model called compressive total variation (CTV) to characterize the sparsity and low-rank prior knowledge of an image. In order to solve the proposed model, we design a concrete algorithm with provably convergence, which is based on inertial proximal ADMM. The performance of the proposed model is tested for magnetic resonance imaging (MRI) reconstruction, image denoising and image deblurring. The proposed method not only recovers edges of the image but also preserves fine details of the image. And our model is much better than the existing regularization models based on the TGV, Shearlet-TGV, ${\ell }_1-\alpha {\ell }_2$TV and BM3D in test for images with piecewise constant regions. And it visibly improves the performances of TV, ${\ell }_1-\alpha {\ell }_2$TV and TGV, and is comparable to Shearlet-TGV in test for natural images.

Introduction

Many image processing problems can be formulated as an inverse problem, in which the data $\mathit{b}$ is assumed to be obtained approximately by applying a linear operator $\mathit{A}$ on an image $\mathit{u}$ with the additive noise $\mathit{e}$. In most of the cases, solving $\mathit{u}$ can be via

$$\mathit{b}=\mathit{Au}\text{with perturbations }\mathit{e}.$$

However, this problem is ill-posed in the sense that directly inverting $\mathit{A}$ would lead to bad and possibly multiple solutions. It is necessary and even desirable to constrain the solutions through a regularization, which provides the prior knowledge of images that one wants to reconstruct or recover. The regularization is usually used to avoid non-uniqueness of solutions and improve the quality of solutions. A general model for such an inverse problem is

$$\hat{\mathit{u}}≔arg\underset{\mathit{u}}{min}J\left(\mathit{u}\right)+\rho \Phi \left(\mathit{u}\right),$$

where $\Phi \left(\mathit{u}\right)$ is the data fitting term, and $J\left(\mathit{u}\right)$ is the regularization term (or penalty term), and $\rho >0$ is a positive parameter to balance two terms, and $\hat{\mathit{u}}$ is an optimal solution of the model.

Usually, the data fitting term $\Phi \left(\mathit{u}\right)$ is taken as ${\Vert \mathit{A}\mathit{u}-\mathit{b}\Vert }_2^2∕2$. And in this subsection, we only give a simple review of the regularization term $J\left(\mathit{u}\right)$.

A classical regularization is the total variation (TV) as follows

$$TV\left(\mathit{u}\right)={\Vert \mathcal{D}\mathit{u}\Vert }_{2,1}=\Vert {\sqrt{{\left|{\mathcal{D}}_x\mathit{u}\right|}^2+{\left|{\mathcal{D}}_y\mathit{u}\right|}^2}\Vert }_1,$$

where ${\mathcal{D}}_x,{\mathcal{D}}_y$ denote the horizontal and vertical partial derivative operators (with certain boundary conditions assumed), respectively, and $\mathcal{D}=[{\mathcal{D}}_x;{\mathcal{D}}_y]$ is the gradient operator $\nabla$ in the discrete setting. Here we give the mathematical definition only in the discrete setting. It is originated by Rudin-Osher-Fatemi [1], which is referred to as the ROF model. It is widely used in image processing applications, such as denoising [1], deconvolution [2], MRI reconstruction [3], recognition [4], inpainting [5], and super-resolution [6] and so on. This TV model is isotropic, and later an anisotropic formulation was addressed in the literature [7] (see also [8] and other references). For any two-dimensional (2D) image, the anisotropic TV is defined by

$$TV\left(\mathit{u}\right)={\Vert \mathcal{D}\mathit{u}\Vert }_1={\Vert {\mathcal{D}}_x\mathit{u}\Vert }_1+{\Vert {\mathcal{D}}_y\mathit{u}\Vert }_1.$$

There exist some different penalties as alternatives to TV. A few notable examples of nonconvex replacements of TV are ${\ell }_p$ quasi-norm of $\mathcal{D}\mathit{u}$ for $0<p<1$ [9], [10], ${\ell }_1-\alpha {\ell }_2$ norm of $\mathcal{D}\mathit{u}$ for $0<\alpha \le 1$ [11], [12], [13]. And there also exists higher order TV-total generalized variation (TGV) (see Section 2.2) [14], [15], which is more precise in describing intensity variations in smooth regions, and thus reduces oil painting artifacts while still being able to preserve sharp edges like TV does. Except these TV and its variants, there also exists some combination of TV (or TGV) and different waves, for example TV+wavelet [16], TGV+shearlet [17], [18] and so on. The connection between TV (or even TGV) and wavelet frames also has been analyzed in [19].

Although TV regularization can preserve edge information, it does not take full advantage of the similarities in images. More regularizations, for example the sparse regularization based dictionary learning and the low-rank regularization based on similar image patches, have been introduced in [20], [21], [22] and  [23], [24], respectively. In here, we do not recall them in detail.

In this subsection, we delineate the motivation. First, we recall the interpretation of TV from the perspective of compressive sensing [25], [26]. The compressive sensing aims to reconstruct a signal or an image from an underdetermined system of linear equations $\mathit{b}=\mathit{A}\mathit{u}$, provided that the signal or image is sufficiently sparse or sparse in a transform domain. For instance, an image is mostly sparse after taking the gradient transform. Mathematically, this problem is equal to minimize the ${\ell }_0$ norm of the image gradient, i.e., $J\left(\mathit{u}\right)={\Vert \mathcal{D}\mathit{u}\Vert }_0$. However, such method is NP-hard and thus computationally infeasible in high dimensional settings. Candès et al. [25] replace ${\ell }_0$ by a convex relaxation ${\ell }_1$, i.e., the ${\ell }_1$ norm of $\mathcal{D}\mathit{u}$, which is the TV. They showed that images could be reconstructed via ${\ell }_1$ norm through numerical tests.

Now, we show the motivation.

• Low-rank Matrix Completion: Matrix completion has been a valuable but difficult task in image inpainting and recommender systems. Given an image with missing parts or corrupted regions, the goal is to find the missing pixel values while ensuring the completed results visually reasonable. For an image, it is probably of (approximately) low rank structure [27]. Candès et al. [28] first solved the matrix completion problem by approximating the rank with nuclear norm as follows

  $$\underset{\mathit{X}}{min}{\Vert \mathit{X}\Vert }_{\ast }$$$$\text{s.t.}{\mathbb{P}}_{\Omega }\left(\mathit{X}\right)={\mathbb{P}}_{\Omega }\left(\mathit{M}\right),$$

  where ${\mathbb{P}}_{\Omega }$ is a projection operator on index set $\Omega \subset {\mathbb{R}}^{m\times n}$, i.e., ${\mathbb{P}}_{\Omega }{\left(\mathit{X}\right)}_{i,j}=X_{i,j}$ for $(i,j)\in \Omega$, and ${\mathbb{P}}_{\Omega }{\left(\mathit{X}\right)}_{i,j}=0$ for $(i,j)\in {\mathbb{R}}^{m\times n}\setminus \Omega$. More works about low-rank matrix completion , readers can refer to an overview [29].

• Total Nuclear Variation (TNV): For a separable image $\mathit{u}\in {\mathbb{R}}^{m\times n\times l}$, all channels at any point should share a common gradient direction. For example, color images $\mathit{u}\in {\mathbb{R}}^{m\times n\times 3}$ using three channels (red, green, blue), share a common gradient direction. Furthermore, having shared directions is equivalent to having a rank-one Jacobian. Therefore Holt [30] thus proposed using nuclear norm in the Jacobian framework, which was called total nuclear variation.

• Low-rank Prior: The similar image patches have similar underlying structures. Thus the matrix constructed from stacking the similar patches together has low rank. Those successful works in [23], [24] show that the methods based on low rank minimization can capture the underlying multi-scale structure and provide good representation for images.

• Robust Principal Component Analysis (RPCA): Some 2D data, has the superposition of a low-rank component and a sparse component, for example video surveillance, face recognition, latent semantic indexing and ranking and collaborative filtering. In [31], [32], the authors proposed a very convenient convex program called robust principal component analysis, to recover both the low-rank component and the sparse component as follows

  $$\underset{\mathit{X},\mathit{Y}}{min}{\Vert \mathit{X}\Vert }_{\ast }+\lambda {\Vert \mathit{Y}\Vert }_1$$$$\text{s.t.}\mathit{X}+\mathit{Y}=\mathit{C},$$

  where $\lambda >0$ is the balancing parameter and $\mathit{C}$ is the observations. More works about RPCA, readers can refer to [33], [34].

• Compressive Phase Retrieval (CPR): In order to solve sparse phase retrieval problem, Ohlsson et al. [35] proposed a model via the lifting technique, which involves a symmetric positive semi-definite, sparse and low-rank matrix $\mathit{X}$, as follows

  $$\underset{\mathit{X}⪰\mathbf{O}}{min}\text{tr}\left(\mathit{X}\right)+\lambda {\Vert \mathit{X}\Vert }_1$$$$\text{s.t.}\mathcal{A}\left(\mathit{X}\right)=\mathit{b},$$

  where $\mathit{b}$ is the observations, $\mathcal{A}:{\mathbb{H}}^{n\times n}\to {\mathbb{H}}^m$ with $\mathbb{H}=\mathbb{R}$ or $\mathbb{H}=ℂ$, and $\mathcal{A}{\left(\mathbf{X}\right)}_j=〈{\mathbf{a}}_j{\mathbf{a}}_j^H,\mathbf{X}〉$. This model is called compressive phase retrieval via the lifting (CPRL) in [35]. Since then, compressive Phase Retrieval has been studied by many scholars [36], [37], [38]. And this model has been extended to compressive affine phase retrieval in [39].

Note that gradient matrices ${\mathcal{D}}_x\mathit{u},{\mathcal{D}}_y\mathit{u}$ of an image $\mathit{u}$ are sparse. And the image $\mathit{u}$ also has (approximately) low rank structure. Motivated by robust principal component analysis and compressive phase retrieval, we conjecture that these gradient matrices also are (approximately) low-rank.

The main contributions of this paper are three folds.

• We consider that the gradient matrices ${\mathcal{D}}_x\mathit{u},{\mathcal{D}}_y\mathit{u}$ of an image $\mathit{u}$ are not only sparse but also (approximately) low-rank, and verify it by numerical tests (Section 2.1) and theoretical analysis (Section 2.2). In order to characterize the low-rank prior information, we introduce the Nuclear norm Total (Generalized) Variation (NT(G)V) (Section 2.3). And we establish a model-compressive total variation (CTV), which reflects these prior informations of the image (Section 2.4).

• Based on inertial proximal ADMM, we design a concrete algorithm to solve our model (Section 3).

• To demonstrate the effectiveness of the proposed method, we test many numerical examples, including MRI reconstruction from incomplete data (Section 4.1), denoising from noisy data (Section 4.2), and deblurring from blurring data (Section 4.3). The tests show that proposed CTV method is better than that of classic TV, TGV, Shearlet-TGV, ${\ell }_1-\alpha {\ell }_2$TV and BM3D methods in testing for piecewise constant images. Our method can not only hold sharper edges but also preserve various image features. And proposed method is comparable to TGV, Shearlet-TGV, and ${\ell }_1-\alpha {\ell }_2$TV methods in testing for natural images.

Our notation is standard, as used above in this section. The standard inner product is denoted by $〈\cdot ,\cdot 〉$. Let ${\Vert \cdot \Vert }_p(0<p\le \infty )$ denote the ${\ell }_p$ norm (or quasi-norm) of matrix or vector, i.e., ${\Vert \mathbf{X}\Vert }_p={\left({\sum }_{i,j}{\left|X_{i,j}\right|}^p\right)}^{1∕p}$ or ${\Vert \mathbf{x}\Vert }_p={\left({\sum }_j{\left|x_j\right|}^p\right)}^{1∕p}$. For any matrix $\mathit{X}\in {\mathbb{R}}^{m\times n}$ with rank $r$, let $\mathit{X}=\mathit{U}\mathit{\Sigma }{\mathit{V}}^T=\mathit{U}\text{diag}\left(\mathit{\sigma }\right){\mathit{V}}^T$ be the singular value decomposition (SVD) of $\mathit{X}$, where $\mathit{\sigma }≔\mathit{\sigma }\left(\mathit{X}\right)={({\sigma }_1,\dots ,{\sigma }_r)}^T$ is the singular vector with ${\sigma }_1\ge {\sigma }_2\ge \cdots \ge {\sigma }_r\ge 0$. We denote ${\Vert \mathit{X}\Vert }_{\ast }≔{\Vert \mathit{\sigma }\Vert }_1$ as the nuclear norm of the matrix $\mathit{X}$. The superscript ${}^T$ denotes the matrix/vector transpose operator. And we denote ${\mathbf{I}}_n$ by $n\times n$ identity matrix, $\mathit{O}$ by zero matrix and $S^n$ denotes the set of $n\times n$ symmetric matrices. In all of this paper, we use boldfaced letter to denote matrix or vector.

The rest of the paper is organized as follows. In Section 2, we test the rank of gradient matrices of images and find that gradient transform matrices are not only sparse but also low-rank or approximately/relatively low-rank, and establish a model-compressive total variation (CTV), which reflects these priors of the image. In Section 3, we design a concrete algorithm based inertial proximal ADMM with provable convergence to solve our model. In Section 4, we compare our CTV method with some other existing methods in MRI reconstruction, image denoising and image deblurring. Conclusions and discussions are given in Section 5.