Magnetic resonance image restoration via least absolute deviations measure with isotropic total variation constraint

1.   Introduction

Magnetic resonance (MR) imaging is a very useful medical imaging method because of the ability to render high anatomical resolution of soft tissues ^[1]. However, due to the inherent limitations of imaging equipment, the MR image will be disturbed by noise during the process of acquisition, resulting in the degradation of image quality. The effect of noise is mainly observed as the blurred regions and the false edges, which degrades the performance of computer-aided diagnosis. Therefore, it is necessary to remove noise while preserving the edges without introducing artifacts to further improve the accuracy of clinical diagnosis.

Restoration methods for MR image can be broadly classified into four types: filtering methods ^[2,3], partial differential equation-based methods ^[4,5], deep learning-based methods ^[6,7] and total variation (TV)-based methods ^[8,9]. For more information about different MR image restoration techniques, please refer to the comprehensive review ^[10]. Here, we are mainly concerned with the TV-based methods. Let $\Omega$ be the domain of image definition; then, the TV of $f$ can be defined as $\int_\Omega \|\nabla f\|_2$, where $\nabla$ is a gradient operator. The famous TV-based deblurring and denoising model is the following minimization problem:

where $A$ is a blurring operator and $y$ denotes the observed image. In Eq (1.1), $\|Af-y\|_2^2$ is a data fidelity term that controls the difference between the degraded image and the observed image. Additionally, the regularization parameter $\lambda$ is used to balance the contributions of the fidelity term and the TV term. Minimizing Eq (1.1) can preserve sharp discontinuities while removing noise and other unwanted fine-scale detail ^[11]. We assume that the MR images are the matrices of size $n \times n$ and denote by $X$ the Euclidean space $R^{n \times n}$; if $f\in X$, then the discrete TV of $f$ can be expressed as

Equation (1.2) is the so-called anisotropic TV using Neumann (symmetric) boundary conditions ^[12], but it could easily yield metrication artifacts. Therefore, one usually uses the isotropic TV ^[12]:

Numerous methods for solving Eq (1.1) have been proposed. In ^[13], Goldfarb and Yin formulated TV minimization as a second-order cone program which was then solved by interior point algorithm. In ^[14,15], Chambolle proposed several gradient descent methods, which have some superior features, including easy implementation and quick convergence. However, for some very ill-conditioned problems, Chambolle's method approaches the minimum too slowly. Yu et al. ^[16] extended Chambolle's method by using an adaptive stepsize with a nonmonotone line search scheme instead of the original constant stepsize. In ^[17], Dahl et al. presented public software for TV denoising, inpainting and deblurring by using some first-order optimization algorithms. In ^[18], Zhu et al. proposed the application of gradient projection algorithms with different step length selection and line search strategies for the dual formulation of a TV model. In ^[19], Bonettini and Ruggiero established the convergence of a general primal-dual method for nonsmooth convex optimization problems in the TV image restoration framework. Based on the TV model and the alternating direction method of multipliers, the authors of ^[20] proposed a fast algorithm that is able to simultaneously estimate the regularization parameter and restore the degraded image. In ^[21], based on the primal-dual TV model, Yu et al. proposed a nonmonotone adaptive projected gradient method, which implements image restoration by deblurring and denoising alternatively. In ^[22], Prasath et al. proposed a spatially adaptive multiscale variable exponent-based anisotropic variational partial differential equation method that overcomes over smoothing and staircasing artifacts, while still retaining and enhancing edge structures across scale.

According to the literature survey, it was found that most of the existing methods mainly focus on the design of solving algorithm for Eq (1.1), and the effect of noise is ignored. In order to further improve the restoration performance of TV-based model, we propose a new MR image deblurring and denoising model based on an isotropic TV regularized least absolute deviations measure. Due to the inherent nonsmoothness of the least absolute deviations measure and the TV term, optimizing the proposed model is challenging. To tackle this problem efficiently, we further introduce an alternating solving algorithm. The proposed model is not sensitive to noise and possesses better restoration ability, which is consistent with the experimental results.

The rest of this paper is organized as follows. We introduce the proposed restoration model and the solving algorithm in Section 2. Comparative experiments on real data are presented and analyzed in Section 3. Finally, we conclude this paper in Section 4.

2.   Model and algorithm

In this section, we present the proposed MR image restoration model, as well as its solving algorithm.

2.1.   Proposed model

Based on the least absolute deviations measure, we propose the following restoration model:

which is more robust than the traditional one given by Eq (1.1), because it can better handle outliers (noise) in the image. For example, if the error $(Af-y)$ is greater than $1$, then the $l_2$-norm enlarges the error, resulting in a larger error than the least absolute deviations measure. In other words, Eq (1.1) is more sensitive to noise than Eq (2.1).

Note that the two terms in Eq (2.1) are both nonsmooth. To transfer $f$ out of the nonsmooth term $\int_\Omega \|\nabla f\|_2$, an auxiliary variable $x \in R^{n^2}$ is introduced, and the difference between $x$ and $f$ is penalized, yielding the following approximation model to Eq (2.1):

where $\alpha > 0$ and $\beta > 0$ are the model parameters.

2.2.   Solving algorithm

By introducing an artificial variable $z$, Eq (2.2) can be rewritten as

Equation (2.3) is an equality-constrained minimization problem, and to make it easier to solve, we transform it into the following unconstrained minimization problem by introducing a Lagrange multiplier $\xi$:

where $\sigma > 0$ is a penalty parameter. The objective function $L$ in Eq (2.4) is known as an augmented Lagrangian function, and now, this minimization problem can be efficiently solved by using an alternating optimization algorithm.

● First, given $x_k$, $z_k$, and $\xi_k$, we get $f_{k+1}$ by solving the following minimization problem:

Taking the derivative of Eq (2.5) with respect to $f$ and forcing the result to zero, we obtain

Thus, we have

● Second, given $f_{k+1}$, we get $x_{k+1}$ by solving the following minimization problem:

Note that the TV of $x$ in Eq (2.7) has the dual form ^[16], i.e., $\int_\Omega \|\nabla x \|_2 = \max_{w\in C_c^1(\Omega), |w|\le1} \int_\Omega \nabla x \cdot w = \max_{|w|\le1} \int_\Omega -x \nabla \cdot w$, where $w: \Omega\rightarrow R^2$ is the dual variable and $\nabla \cdot$ is a divergence operator. Based on this dual model, optimizing Eq (2.7) is equivalent to solving the following optimization problem:

where $\mu = \alpha / \beta$. Here, for convenience, we introduce the nonmonotone adaptive projected gradient method proposed in ^[21] to solve Eq (2.8), which is described as follows. Further, for the details of this method, we refer the interested readers to ^[21].

The full description of the nonmonotone adaptive projected gradient method for solving Eq (2.8).

Lemma 1 (Theorem 3.1 ^[21]). Suppose that the sequence $\{w_k\}$ is generated by the nonmonotone adaptive projected gradient method; then, any accumulation point of $\{w_k\}$ is a constrained stationary point.

When the minimum $w^*$ of Eq (2.8) is determined by the nonmonotone adaptive projected gradient method, then $x_{k+1}$ can be generated by

● Third, given $f_{k+1}$ and $\xi_k$, we get $z_{k+1}$ by solving the following minimization problem:

which is equivalent to the following minimization problem:

By using the soft-thresholding operator $S$, we have

where $S$ is defined as $S_b(a) = \max\{a-b, 0 \} - \max\{-a-b, 0 \}$.

● Finally, we update the Lagrange multiplier $\xi$ by

The full steps of the proposed method LADTV (short for isotropic total variation regularized least absolute deviations measure) are presented in Algorithm 1. Further, for this algorithm, we have some remarks.

Table Algorithm 1.  LADTV method.

1 Initialization: $A$, $\alpha$, $\beta$, $x_0$, $z_0$, and $\xi_0$.
2 while not converged do 7 end

---

 | Show Table

DownLoad: CSV

Remark 1. If $A$ is set to be an identity matrix, then Eq (2.1) is reduced to a denoising problem.

Remark 2. In Algorithm 1, we use this order, $f_{k+1} \rightarrow x_{k+1} \rightarrow z_{k+1} \rightarrow \xi_{k+1}$, to update the model parameters. Different orders may yield different results and also show different restoration performance.

Remark 3. For convenience of optimization, we directly introduce the traditional soft-thresholding operator to solve Eq (2.11). Some other advanced algorithms can be considered to further improve the solving accuracy, such as the primal dual hybrid gradient method ^[23], proximal gradient method ^[24] and block coordinate descent method ^[25].

3.   Experiments

In this section, we conduct comparative experiments to show the restoration performance of our approach from blurred and noisy observation.

We test two kinds of point spread function (psf) provided by using the Matlab function $fspecial$. One is the Gaussian blurring, the other is the motion blurring. Gaussian blurring is a rotationally symmetric Gaussian low-pass filter. Motion blurring shifts image pixels linearly in a certain direction with a fixed length. For simplicity, $G(hsize, sigma)$ denotes the Gaussian blurring kernel with a blurring size ($hsize$) and a standard deviation ($sigma$), and $M(len, theta)$ is the motion blurring kernel with a motion length ($len$) and an angle ($theta$). Then, each observed image is generated by adding noise to the blurred image by using the Matlab function $randn$ with the variance parameter $var$.

The quality of restoration performance is measured by the signal-to-noise ratio (SNR for short), which is defined as ${\rm SNR} = 10\log_{10} \|f_0-M(f_0)\|_2^2/\|f_0-f\|_2^2$. Here, $f_0$ denotes the original image, $f$ denotes the restored image and $M(f_0)$ is the mean intensity value of $f_0$.

3.1.   Experimental setup

3.1.1.   Data description

We use clinical data including eight parts to show the restoration performance of the proposed method. The original test images are shown in Figure 1. All test images are from the Department of Radiology, Maanshan People's Hospital, Maanshan, China, and the resolution of each image is $256 \times 256$.

Figure 1.  Original images: (a) abdomen, (b) ankle, (c) breast, (d) cervical vertebra, (e) head, (f) liver, (g) pelvis, (h) sacroiliac joint.

DownLoad: Full-Size Img PowerPoint

3.1.2.   Baseline methods

In order to highlight the restoration performance of LADTV, we compare it with other three restoration methods including waveletFISTA ^[26], tvFISTA ^[27] and NAPG ^[21]. The waveletFISTA is a fast iterative shrinkage-thresholding algorithm (FISTA) for solving wavelet-based image deblurring problem, and the tvFISTA is a fast algorithm for the constrained TV-based image deburring problem. The third method NAPG is a nonmonotone adaptive projected gradient algorithm based on a primal-dual TV model for solving image restoration problem from noisy and blurred observation. In LADTV, we set $\alpha = 0.05$, $\beta = 2$ and $\sigma = 1$. The parameters of other methods are set to default values.

3.2.   Results and analysis

First, we set a uniform $3 \times 3$ blur and let psf = (1/9)$\ast$ones(3, 3); then, each observed image is generated by adding noise to the blurred image with $var = 0.001$ and $var = 0.01$, respectively. Tables 1 and 2 report the obtained SNR results.

Table 1.  The SNR (dB) results obtained by four methods when $var = 0.001$.

Test image number	waveletFISTA	tvFISTA	NAPG	LADTV
(a)	27.17	28.04	26.20	31.79
(b)	17.76	18.06	24.47	28.13
(c)	26.04	27.83	22.40	26.96
(d)	20.71	21.15	27.05	31.50
(e)	30.33	32.04	26.56	31.90
(f)	28.75	29.59	26.41	31.75
(g)	19.91	20.68	22.04	26.51
(h)	19.96	20.72	20.15	24.33
mean value	23.83	24.76	24.41	29.11

---

 | Show Table

DownLoad: CSV

Table 2.  The SNR (dB) results obtained by four methods when $var = 0.01$.

Test image number	waveletFISTA	tvFISTA	NAPG	LADTV
(a)	25.91	26.97	26.14	29.53
(b)	17.63	17.98	24.43	27.29
(c)	22.63	25.39	22.35	25.36
(d)	20.51	20.00	27.01	30.01
(e)	28.49	30.41	26.49	30.10
(f)	26.32	27.98	26.32	29.37
(g)	19.20	20.23	22.00	25.08
(h)	19.36	20.22	20.11	23.16
mean value	22.51	23.65	24.36	27.49

---

 | Show Table

DownLoad: CSV

In each table, the red font indicates the maximum value, and the second best value is marked with blue. It is observed that our method gets the suboptimal results on images (c) and (e), and it obtains the highest values on all the other test images. However, from the perspective of mean value, our method LADTV has the best restoration performance.

Further, we use the Gaussian blurring and the motion blurring to blur the test images, and then each observed image is generated by adding noise to the blurred image with $var = 0.001$. Figures 2 and 3 show the results when the original images are degraded by $G(3, 3)$. Figures 4 and 5 show the results when the original images are degraded by $M(3, 3)$. From the restoration results displayed in Figures 2–5, we see that LADTV has very obvious advantages compared to other methods.

Figure 2.  Restoration results with $G(3, 3)$.

DownLoad: Full-Size Img PowerPoint

Figure 3.  Restoration results with $G(3, 3)$ (con't).

DownLoad: Full-Size Img PowerPoint

Figure 4.  Restoration results with $M(3, 3)$.

DownLoad: Full-Size Img PowerPoint

Figure 5.  Restoration results with $M(3, 3)$ (con't).

DownLoad: Full-Size Img PowerPoint

Second, to learn more about the restoration performance of the least absolute deviations measure (i.e., LADTV) and the traditional Euclidean distance measure (i.e., NAPG), we add more serious blurred interference. As examples, images "abdomen" and "head" are taken to show the restoration performance. Figures 6 and 7 show the trend of SNR values as the number of iterations increases. We can see that on the whole the SNR values obtained by LADTV are better than that obtained by NAPG.

Figure 6.  SNR vs. Iterations on "abdomen".

DownLoad: Full-Size Img PowerPoint

Figure 7.  SNR vs. Iterations on "head".

DownLoad: Full-Size Img PowerPoint

Tables 3 and 4 shows the computing time of LADTV and NAPG when the two images are degraded by Gaussian blurring and Motion blurring, respectively. Because our method LADTV involves more variable update calculations, it takes more time to complete the restoration task. Figures 8 and 9 display the restoration results corresponding to Tables 3 and 4. Figure 8 is the results with $G(10, 10)$ and $G(15, 15)$, and Figure 9 is the results with $M(30, 30)$ and $M(50, 50)$. It is observed that in the case of severely interference (i.e., very poor visual effects), both LADTV and NAPG show good ability in deblurring and denoising simultaneously. However, overall, the proposed method LADTV can also get higher quality than NAPG. Given all that, there has to be a balance in the restoration performance and computational efficiency.

Table 3.  Computing time (s) of LADTV and NAPG on "abdomen".

Method	abdomen, G(10, 10)	abdomen, G(15, 15)	head, G(10, 10)	head, G(15, 15)
LADTV	16.61	17.67	16.89	17.73
NAPG	15.95	13.20	13.84	14.45

---

 | Show Table

DownLoad: CSV

Table 4.  Computing time (s) of LADTV and NAPG on "head".

Method	abdomen, M(30, 30)	abdomen, M(50, 50)	head, M(30, 30)	head, M(50, 50)
LADTV	18.03	16.84	19.80	18.98
NAPG	16.30	14.14	13.63	16.36

---

 | Show Table

DownLoad: CSV

Figure 8.  The first line is the restoration results for "abdomen" with $G(10, 10)$, and the second line is the restoration results for "abdomen" with $G(15, 15)$. The third line is the restoration results for "head" with $G(10, 10)$, and the fourth line is the restoration results for "head" with $G(15, 15)$.

DownLoad: Full-Size Img PowerPoint

Figure 9.  The first line is the restoration results for "abdomen" with $M(30, 30)$, and the second line is the restoration results for "abdomen" with $M(50, 50)$. The third line is the restoration results for "head" with $M(30, 30)$, and the fourth line is the restoration results for "head" with $M(50, 50)$.

DownLoad: Full-Size Img PowerPoint

Finally, we carry out some experiments on images "cervical vertebra" and "sacroiliac joint" to show the sensitivity of the proposed model to its parameters $\alpha$, $\beta$, and $\sigma$. For simplicity, we set $\alpha = \{0.01, 0.05, 1\}$, $\beta = \{1.5, 2, 2.5\}$ and $\sigma = \{0.5, 1, 1.5\}$. Preliminary experimental results are presented in Figures 10 and 11. From the perspective of SNR, the best SNR values is obtained when $\alpha = 0.05$, $\beta = 2$ and $\sigma = 1$. In this situation, the ratio of the best values is 41.67% (i.e., 5/12). The ratio of the best values in the situation that $\alpha = 0.05$, $\beta = 2$ and $\sigma = 0.5$ is 33.33% (i.e., 4/12). The ratio of the best values in the situation that $\alpha = 0.05$, $\beta = 2.5$ and $\sigma = 1$ is 25% (i.e., 3/12). However, overall, the fluctuation of these SNR values is not very large. Besides, it seems that the computing time changes without certain regulations, but the gap is not obvious. These results suggest that the proposed method has better stability and reliability.

Figure 10.  Sensitivity analysis on "cervical vertebra".

DownLoad: Full-Size Img PowerPoint

Figure 11.  Sensitivity analysis on "sacroiliac joint".

DownLoad: Full-Size Img PowerPoint

4.   Conclusions

We propose in this paper a new MR image restoration model for recovering high quality image from a blurred and noisy image. We show that combining least absolute deviations measure and isotropic TV is a feasible solution to the MR image restoration problem. This combination brings the ability to suppress noise and preserve image smoothness together for efficiently recoverying degraded image. Both qualitative and quantitative analytical approaches are introduced in our experiments to show the feasibility and effectiveness of our method.

The proposed method is a little computationally expensive compared with other methods. The major reason is that the program takes some time to calculate the dual variable, which is implemented using a projected gradient decent method with nonmonotone line search technique. In future work, we will use C++ to speed up the current program. Besides, we would like to clarify that the main purpose of this paper is to demonstrate that the TV-based MR image restoration can be improved by combining least absolute deviations measure. Once more an efficient solver is developed and deployed, the computational efficiency of the proposed method could also be improved.

The proposed method in this paper is mainly developed for two-dimensional magnetic resonance image restoration, in the future, we intend to extend the proposed method to MR image reconstruction, registration, three-dimensional image restoration and so on. In addition, Moreau Envelope ^[28,29] will also be considered to maintain the details of image more effectively, such as texture and edges.

Acknowledgments

This work was supported in part by the Anhui Provincial Natural Science Foundation (Grant No. 2208085MF168).

Conflict of interest

The authors declare that there is no conflict of interests regarding the publication of this article.