Image Restoration under Cauchy Noise: A Group Sparse Representation and Multidirectional Total Generalized Variation Approach

Image restoration aims to estimate the original clear image based on an observed image. Cauchy noise [1, 2] is a prevalent impulse noise with the characteristic of a thick tail, which is widely present in radar imaging, astronomical images, biomedical images, and wireless communication systems. Numerous studies have been conducted to remove Cauchy noise, particularly in the wavelet domain. Zhao et al. [3] employed the maximum a posteriori (MAP) estimation of double tree complex wavelet transform for image noise removal. Hill et al. [4] proposed a new method of single-variable shrinkage in the DWT domain, using translation invariant wavelet decomposition for image denoising. Wu et al. [5] utilized a multi-level wavelet convolution neural network to train the denoiser, maintaining more detailed texture during the denoising process and better balancing the size of the receptive field with the computational load.

Recently, variational models based on total variation regularization have been extensively employed for image denoising. Sciacchitano et al. [1] introduced a total variation model [6] for denoising Cauchy noise in images. Mei et al. [2] applied a specialized multiplier alternating direction method [7] to solve the nonconvex total variation minimization problem. However, total variation is influenced by the staircase effect in smooth regions and cannot effectively recover the details and textures of an image, resulting in an oversmoothed image. To address the edge staircase effect, Yang et al. [8] incorporated higher-order total variation to reduce the staircase effect. Parrisoto et al. [9] introduced a high-order anisotropic total variation regularization term that retains and enhances the inherent anisotropy characteristics of an image, recovering the details and textures as much as possible. Shi et al. [10] proposed a nonlinear diffusion equation for denoising to better recover image detail features, including a diffusion coefficient based on gray scale to estimate noise amplitude and a gradient-based diffusion coefficient to control anisotropic diffusion in accordance with the local structure of an image. Liu and Gao [11] proposed a non-blind image deblurring method that combines multi-directional fractional total variation with traditional total variation, addressing the staircase effect and texture loss issues of total variation while effectively preserving texture and edge information.

For images containing not only flat regions but also inclined ones, the total generalized variation (TGV) model demonstrates strong denoising capabilities in single-image denoising. Bredies et al. [12] suggested using total generalized variation as a penalty function for modeling images with edges and smooth changes, exhibiting excellent denoising performance in single-image denoising. Lv [13] proposed a multilook M for total generalized variation-based multiplicative noise removal to eliminate speckle noise from images. Shi et al. [14] introduced a total generalized variation method for reconstructing electrical impedance tomography (EIT), combining the first and second derivative terms as regularizers to address issues related to Tikhonov regularization's excessive smoothness in image reconstruction and the staircase effect in total variation regularization reconstruction. To improve the denoising effect, Zhong et al. [15] presented a high-order TGV regularization variational model based on the automatic selection of spatial adaptive regularization parameters according to local image characteristics. To better restore image texture in a certain direction, Kongskov et al. [16] proposed a texture direction estimation algorithm and a novel direction total generalized variation model (DTGV), which significantly improves texture preservation and noise removal. Li et al. [17] introduced a weighted model of second-order total generalized variation for Gaussian noise removal by incorporating an edge indicator function into the regularizer of the second-order total generalized variation model [12]. To enhance the intensity of the diffusion tensor and improve the image's visual quality, the method employs first and second derivatives. Since the Shearlet transform can sparsely represent an image and generate the best approximation [18], Lv [19] combined the Shearlet transform with a second-order total generalized variation regularization term, proposing a new recovery model for images corrupted with Cauchy noise. Liu [20] developed a hybrid regularization model for image denoising and deblurring by combining the advantages of total generalized variation with the wavelet method. Zhu et al. [21] presented a second-order TGV and wavelet framework-based hybrid regularization method, inheriting the benefits of wavelet framework regularization and second-order TGV. This model not only removes the staircase effect but also preserves sharp edges while maintaining strong sparse estimation capabilities for piecewise smooth functions.

However, since the TGV model processes pixels independently, it disregards the similarity of processed images, resulting in weak robustness against high-amplitude noise. To address this issue, Zhang et al. [22] employed the nonlocal total generalized variation method for image repair and super-resolution reconstruction. To exploit the sparse representation of image block information, Jung and Kang [23] proposed a nonlocal total generalized variation model based on sparse representation, which improves the denoising performance by considering image block similarity. Zhang et al. [24] proposed a nonlocal total generalized variation model based on local similarity, which takes into account the similarity of image patches in the local area and addresses the over-smoothness issue in flat regions. Chen et al. [25] introduced the sparsity based on overlapping group in TGV model, which can improve the denoising effect of TGV by using structural similarity, and the robustness of TGV to strong noise pollution by using the first and second-order gradient information of domain difference.

Compared with the traditional method using local features, the recent variation method using nonlocal information of image has a great improvement of noise removal effect. Ding et al. [26] used the nonlocal self-similarity of natural images to treat a group of similar blocks as an approximate low rank matrix, and expressed the denoising problem of each group as a low rank matrix restoration problem. Laus et al. [27] used a nonlocal, completely unsupervised method to remove the Cauchy noise in the image. Kim et al. [28] used the weighted kernel norm as the regularization term, and used the similar blocks in the image to remove the additive Cauchy noise through nonlocal similarity.

Although the nonlocal method makes use of the similarity among image blocks to improve the denoising effect, the low similarity or dissimilarity among image blocks limits the applicability of the method. In order to overcome these limitations and better consider the structure of the processed image, group or structured sparse representation has been widely studied in the field of image restoration in recent years [29, 30]. The methods based on group or structured sparse representation can capture the internal characteristics of image structure, enhance the inherent local sparsity and nonlocal self-similarity of image, and improve the effect of image deblurring [30-32], image inpainting [33-35], and denoising [36].

Inspired by the fact that group sparse representation could further suppress the staircase effect of the traditional total generalized variational model and multidirectional total generalized variation can better protect the edge structure features of the image and further suppress false edges, we propose a new model to remove noise and/or blur in the image corrupted by Cauchy noise and/or blur. The model uses the prior knowledge of group sparse representation learned from dictionary and higher-order derivative based on nonlocal multi-directional total generalized variation. The prior of group sparse representation can effectively denoise the uniform region and preserve the texture and detail in formation of the image, while the nonlocal multi-directional total generalized variation based high-order derivative prior can denoise the smooth region and preserve the edge in formation of the image. In addition, an effective iterative algorithm is proposed to solve the model.

In this study, a new nonlocal total generalized variation (NLTGV) model for Cauchy noise removal that combines the advantages of total generalized variation and nonlocal means is proposed. The model effectively removes Cauchy noise while preserving image details and textures. The main contributions of this study are as follows:

1. A novel NLTGV model is proposed for Cauchy noise removal. The model combines the strengths of the edge-preserving total generalized variation and the nonlocal means, which considers the similarity of image patches, leading to improved denoising performance.
2. We propose an iterative algorithm for solving the NLTGV model, which effectively converges to the optimal solution while maintaining low computational complexity.
3. Experimental results show that the proposed NLTGV model outperforms several state-of-the-art denoising methods in terms of both visual quality and quantitative evaluation metrics, demonstrating the effectiveness of the proposed model in Cauchy noise removal and image detail preservation.

The rest of this paper is organized as follows: Section 2 describes the proposed NLTGV model and the corresponding optimization algorithm. Section 3 presents experimental results and comparisons with other denoising methods. Finally, Section 4 concludes the paper and discusses future work.

It is assumed that the Cauchy noise $\eta$ in Eq. (1) follows the Cauchy distribution $C(0, \gamma)$ with the location parameter set to $\delta=0$. In order to restore the unknown original image $\mathrm{x}$, the maximum a posteriori (MAP) estimation is commonly employed by maximizing the conditional posterior probability $p(x \mid y)$. Based on the Bayesian principle $p(x \mid y)=\frac{p(y \mid x) p(x)}{p(y)}$ and the application of negative logarithm operation, the MAP estimation $\mathrm{x}$ is obtained through the minimization of $-\log (p(y \mid x))-\log (p(x))$, which $-\log (p(y \mid x))$ serves as a data-fitting term, while $-\log (p(y \mid x))$ also being utilized as a regularizer. It is further assumed that an observed image $y \in \mathbb{R}^{\sqrt{N}\times \sqrt{N}} $ and its corresponding unknown original image $x \in \mathbb{R}^{\sqrt{N} \times \sqrt{N}}$ are transformed into column vectors according to the matrix columns, indexed by $A=\{1,2, \cdots, \sqrt{N}\} \times\{1,2, \cdots, \sqrt{N}\}$. Given that image pixels are considered independent and identically distributed, there are $p(x)=\prod_{j \in A} p\left(x_j\right)$. The relationship between $\delta=0$ and $p(y \mid x)=p(\eta)$ is established, resulting in $p\left(y_j \mid x_j\right)=\frac{\gamma}{\pi\left(\gamma^2+\left((H x)_j-y_j\right)^2\right)}$. By taking logarithms on both sides of this relationship, the following equation is derived:

$\begin{aligned} & -\log (\mathrm{p}(\mathrm{y} \mid \mathrm{x}))=-\sum_{\mathrm{j} \in \mathrm{A}} \log \left(\mathrm{p}\left(\mathrm{y}_{\mathrm{j}} \mid \mathrm{x}_{\mathrm{j}}\right)\right) \\ & =\sum_{\mathrm{j} \in \mathrm{A}} \log \left(\gamma^2+\left((\mathrm{Hx})_{\mathrm{j}}-\mathrm{y}_{\mathrm{j}}\right)^2\right)+\log \pi-\log \gamma .\end{aligned}$      (13)

As discussed in Section 2.2, the nonlocal self-similarity of the image is leveraged, and a more robust geometric structure is preserved through the prior knowledge of group sparse representation. Consequently, the data fidelity term is combined with group sparse representation in the model for Cauchy noise removal. In order to further preserve edge information in multiple directions of an image, the multi-directional total generalized variation method introduces two additional diagonal directions, as depicted in Figure 1, in accordance with the traditional directional total generalized variation model. This allows for the detection of edge features in eight distinct directions during image restoration. In conjunction with overlapping group sparsity [39], the overlapping group sparse multidirectional total generalized variation (MDTGV) regularization term is defined as follows:

$\text{MDTGV}(x)=\underset{w}{\mathop{\min }}\,\left\{ {{\lambda }_{1}}\sum\limits_{j=1}^{4}{\begin{align}  & \phi \left( {{\nabla }_{j}}x-{{w}_{j}} \right) \\ & +{{\lambda }_{0}}\sum\limits_{j=1}^{16}{\phi \left( {{\left[ \varepsilon \left( w \right) \right]}_{j}} \right)} \\\end{align}} \right\},$     (14)

where, $\nabla_1 x$ and $\nabla_2 x$ are described in Section 2.3, $\nabla_3 x$ and $\nabla_4 x$ represent the first-order gradients of $x$ in the diagonal $45^{\circ}$ and $135^{\circ}$ directions, respectively, and $\varepsilon(w)_j=\nabla_j w_j+\sum_{i=1 \text { and } i \neq j}^4 \frac{\nabla_i}{2} w_j(j=\{1,2,3,4\})$ denotes the second-order gradient operator.

1.png

Figure 1. Eight domain space of pixels

Building upon the overlapping group sparse multidirectional total generalized variation model, a method for eliminating Cauchy noise from an image is proposed as follows:

$\begin{align}  & \underset{x,\text{w}}{\mathop{\min }}\,\frac{\lambda }{2}<\log \left( {{\gamma }^{2}}+{{\left( Hx-y \right)}^{2}} \right), \\ & 1>+{{\lambda }_{1}}\sum\limits_{j=1}^{4}{\phi \left( {{\nabla }_{j}}x-{{w}_{j}} \right)+{{\lambda }_{0}}\sum\limits_{j=1}^{16}{\phi \left( {{\left[ \varepsilon \left( w \right) \right]}_{j}} \right)}}. \\\end{align}$         (15)

In summary, two types of priors are combined: group sparse representation priors (GSR) and overlapping group sparse multidirectional total generalized variation (MDTGV). Thus, a new model for denoising images contaminated with Cauchy noise and/or blur is proposed in this study:

$\begin{align}  & \underset{\{{{a}_{{{G}_{i}}}}\},\{{{D}_{{{G}_{i}}}}\},\{{{x}_{{{G}_{i}}}}\},\text{w}}{\mathop{\min }}\,\sum\limits_{i=1}^{n}{\left( \frac{\lambda }{2}<\log \left( {{\gamma }^{2}}+{{\left( H{{x}_{{{G}_{i}}}}-{{y}_{{{G}_{i}}}} \right)}^{2}} \right),\text{ }1>+{{\mu }_{{{G}_{i}}}}{{\left\| {{\alpha }_{{{G}_{i}}}} \right\|}_{0}}+\left\| {{D}_{{{G}_{i}}}}{{\alpha }_{{{G}_{i}}}}-{{x}_{{{G}_{i}}}} \right\|_{2}^{2}+{{\lambda }_{1}}\sum\limits_{j=1}^{4}{\phi \left( {{\nabla }_{j}}{{x}_{{{G}_{i}}}}-{{w}_{j}} \right)+{{\lambda }_{0}}\sum\limits_{j=1}^{16}{\phi \left( {{[\varepsilon \left( w \right)]}_{\text{j}}} \right)}} \right)}, \\ &  \\\end{align}$      (16)

where, $x \in \mathbb{R}^{\sqrt{N} \times \sqrt{N}}, y \in \mathbb{R}^{\sqrt{N} \times \sqrt{N}}, D_{G_i} \in \mathbb{R}^{P \times k}(P<<k), \alpha_{G_i} \in \mathbb{R}^{k \times h}, x_{G_i} \in \mathbb{R}^{P \times h}$, and $\mathrm{n}$ represent the number of image block groups, as described in subsection 2.2. The regularization parameters $\lambda_1$ and $\lambda_0$ are utilized to balance the first-order difference and the second-order difference in multidirectional total generalized variation.

3.1 Decomposition of model (16)

The model (16) proposed in this study is non-convex due to data fidelity terms and the product of $D_{G_i}$ and $\alpha_{G_i}$. Auxiliary variables $p \in \mathbb{R}^N, q_j \in\left(\mathbb{R}^{P h}\right)^2(j=\{1,2,3,4\}), r_j \in\left(\mathbb{R}^{\mathrm{Ph}}\right)^4(j=\{1,2, \cdots, 16\})$ and $z \in \mathbb{R}^N$ are introduced, allowing the unconstrained model (16) to be rewritten as:

$\begin{align}  & \underset{\{{{a}_{{{G}_{i}}}}\},\{{{D}_{{{G}_{i}}}}\},\{{{x}_{{{G}_{i}}}}\},\text{w,p,q,r,z}}{\mathop{\min }}\,\sum\limits_{i=1}^{n}{\left( \frac{\lambda }{2}<\log \left( {{\gamma }^{2}}+{{\left( z-{{y}_{{{G}_{i}}}} \right)}^{2}} \right),\text{ }1>+{{\mu }_{{{G}_{i}}}}{{\left\| {{\alpha }_{{{G}_{i}}}} \right\|}_{0}}+\left\| {{D}_{{{G}_{i}}}}{{\alpha }_{{{G}_{i}}}}-p \right\|_{2}^{2}+{{\lambda }_{1}}\sum\limits_{j=1}^{4}{\phi \left( {{q}_{j}} \right)+{{\lambda }_{0}}\sum\limits_{j=1}^{16}{\phi \left( {{r}_{j}} \right)}} \right)}, \\ & s.t. z=H{{x}_{{{G}_{i}}}},p={{x}_{{{G}_{i}}}},{{q}_{j}}={{\nabla }_{j}}{{x}_{{{G}_{i}}}}-{{w}_{j}},{{r}_{j}}=\varepsilon {{(w)}_{j}}. \\\end{align}$       (17)

By incorporating a penalty term into the constraint condition, the constrained model (17) is relaxed to the following unconstrained model:

$\underset{\{{{a}_{{{G}_{i}}}}\},\{{{D}_{{{G}_{i}}}}\},\{{{x}_{{{G}_{i}}}}\},\text{w,p,q,r,z}}{\mathop{\min }}\,\;\;\sum\limits_{i=1}^{n}{\left\{ \begin{align}  & \frac{\lambda }{2}<\log \left( {{\gamma }^{2}}+{{\left( z-{{y}_{{{G}_{i}}}} \right)}^{2}} \right),\text{ }1>+{{\mu }_{{{G}_{i}}}}{{\left\| {{\alpha }_{{{G}_{i}}}} \right\|}_{0}}+\left\| {{D}_{{{G}_{i}}}}{{\alpha }_{{{G}_{i}}}}-p \right\|_{2}^{2}\text{+}{{\lambda }_{1}}\sum\limits_{j=1}^{4}{\phi \left( {{q}_{j}} \right)+}{{\lambda }_{0}}\sum\limits_{j=1}^{16}{\phi \left( {{r}_{j}} \right)} \\ & +\frac{{{\Upsilon }_{1}}}{2}\sum\limits_{j=1}^{4}{\left\| {{q}_{j}}\text{-(}{{\nabla }_{j}}{{x}_{{{G}_{i}}}}-{{w}_{j}}) \right\|_{2}^{2}\text{+}\frac{{{\Upsilon }_{2}}}{2}\sum\limits_{j=1}^{16}{\left\| {{r}_{j}}-\varepsilon {{(w)}_{j}} \right\|_{2}^{2}+\frac{\beta }{2}\left\| p\text{-}{{x}_{{{G}_{i}}}} \right\|_{2}^{2}+\frac{\xi }{2}\left\| z-H{{x}_{{{G}_{i}}}} \right\|_{2}^{2}}} \\\end{align} \right\}},$          (18)

where, $\beta, \Upsilon_1, \Upsilon_2$ and $\xi$ are parameters greater than 0 . In order to solve the model (18), an alternating minimization algorithm (AMA) [38] is employed, which minimizes one variable while keeping other variables fixed and iterates this process. When AMA is applied to model (18), the following subproblems arise:

$({{a}_{{{G}_{i}}}},{{D}_{{{G}_{i}}}})\in \underset{{{a}_{{{G}_{i}}}},{{D}_{{{G}_{i}}}}}{\mathop{\arg \min }}\,{{\mu }_{{{G}_{i}}}}{{\left\| {{\alpha }_{{{G}_{i}}}} \right\|}_{0}}+\left\| {{D}_{{{G}_{i}}}}{{\alpha }_{{{G}_{i}}}}-p \right\|_{2}^{2},$        (19)

$p\in \underset{p}{\mathop{\arg \min }}\,\left\| {{D}_{{{G}_{i}}}}{{\alpha }_{{{G}_{i}}}}-p \right\|_{2}^{2}+\frac{\beta }{2}\left\| p\text{-}{{x}_{{{G}_{i}}}} \right\|_{2}^{2},$      (20)

$q\in \underset{q}{\mathop{\arg \min }}\,{{\lambda }_{1}}\sum\limits_{j=1}^{4}{\phi \left( {{q}_{j}} \right)}\text{+}\frac{{{\Upsilon }_{1}}}{2}\sum\limits_{j=1}^{4}{\left\| {{q}_{j}}\text{-(}{{\nabla }_{j}}{{x}_{{{G}_{i}}}}-{{w}_{j}}) \right\|_{2}^{2}},$        (21)

$r\in \underset{r}{\mathop{\arg \min }}\,{{\lambda }_{0}}\sum\limits_{j=1}^{16}{\phi \left( {{r}_{j}} \right)}\text{+}\frac{{{\Upsilon }_{2}}}{2}\sum\limits_{j=1}^{16}{\left\| {{r}_{j}}-\varepsilon {{(w)}_{j}} \right\|_{2}^{2}},$      (22)

$\begin{align}  & z\in \underset{z}{\mathop{\arg \min }}\,\frac{\lambda }{2}<\log \left( {{\gamma }^{2}}+{{\left( z-{{y}_{{{G}_{i}}}} \right)}^{2}} \right), \\ & \text{ }1>\text{+}\frac{\xi }{2}\left\| z-H{{x}_{{{G}_{i}}}} \right\|_{2}^{2}, \\\end{align}$        (23)

$\begin{align}  & {{x}_{{{G}_{i}}}}\in \underset{{{x}_{{{G}_{i}}}}}{\mathop{\arg \min }}\,\frac{\beta }{2}\left\| p\text{-}{{x}_{{{G}_{i}}}} \right\|_{2}^{2}\text{+}\frac{{{\Upsilon }_{1}}}{2}\sum\limits_{j=1}^{4}{\left\| {{q}_{j}}\text{-(}{{\nabla }_{j}}{{x}_{{{G}_{i}}}}-{{w}_{j}}) \right\|_{2}^{2}} \\ & \text{+}\frac{\xi }{2}\left\| z-H{{x}_{{{G}_{i}}}} \right\|_{2}^{2}, \\\end{align}$          (24)

$\begin{align}  & w\in \underset{w}{\mathop{\arg \min }}\,\frac{{{\Upsilon }_{1}}}{2}\sum\limits_{j=1}^{4}{\left\| {{q}_{j}}\text{-(}{{\nabla }_{j}}{{x}_{{{G}_{i}}}}-{{w}_{j}}) \right\|_{2}^{2}} \\ & \text{+}\frac{{{\Upsilon }_{2}}}{2}\sum\limits_{j=1}^{16}{\left\| {{r}_{j}}-\varepsilon {{(w)}_{j}} \right\|_{2}^{2}}. \\\end{align}$        (25)

3.2 Solving $\left(a_{G_i}, D_{G_i}\right)$​  sub-problems

With $\mathrm{p}$ fixed, the subproblem (19) for $\left(a_{G_i}, D_{G_i}\right)$ is resolved. By fixing $D_{G_i}$ and $p$, the minimization problem of $a_{G_i}$ becomes:

$\underset{{{a}_{{{G}_{i}}}}}{\mathop{\min }}\,{{\mu }_{{{G}_{i}}}}{{\left\| {{\alpha }_{{{G}_{i}}}} \right\|}_{0}}+\left\| {{\text{D}}_{{{G}_{i}}}}{{\alpha }_{{{G}_{i}}}}-p \right\|_{2}^{2}.$        (26)

The Orthogonal Matching Pursuit (OMP) algorithm [40] can be utilized to solve the subproblem (26). OMP, an iterative greedy algorithm, selects the most relevant column to the current residuals in $D_{G_i}$ at each step. The OMP algorithm stops when the error $\left\|D_{G_i} \alpha_{G_i}-p\right\|_2^2$ falls below $\theta^2$ ( $\theta$ is a very small number). In order to obtain the dictionary $D_{G_i}$, the KSVD algorithm [41] is employed to solve the subproblem (19). In fact, the sparse coding and dictionary are repeatedly updated $\mathrm{J}$ times within the KSVD algorithm.

3.3 Solving p and z subproblems

The subproblem of solving $\mathrm{p}$ in (20) is a least square problem. Assuming $L=\min _p\left\|D_{G_i} \alpha_{G_i}-p\right\|_2^2+$ $\frac{\beta}{2}\left\|p-x_{G_i}\right\|_2^2$, the partial derivatives of $\mathrm{p}$ for $\mathrm{L}$ is $\frac{\partial L}{\partial p}=2\left(p-D_{G_i} a_{G_i}\right)+\beta\left(p-x_{G_i}\right)$, and letting $\frac{\partial L}{\partial p}=0$, the closed-form solution of the $p$ subproblem in (20) is as follows:

$p\text{=}{{\left( 2I+\beta I \right)}^{-1}}\left( 2{{D}_{{{G}_{i}}}}{{a}_{{{G}_{i}}}}+\beta {{x}_{{{G}_{i}}}} \right),$        (27)

where, I represents the identity matrix, and I maintains the same meaning in subsequent formulas unless stated otherwise.

When $8 \gamma \geq \frac{\lambda}{\xi}$, it can be demonstrated that the minimization problem $(23)$ is strictly convex, and its minimum can be obtained by solving its Euler-Lagrange equation. Assuming $L=\min _z \frac{\lambda}{2}<\log \left(\gamma^2+\left(z-y_{G_i}\right)^2\right), 1>$ $+\frac{\xi}{2}\left\|z-H x_{G_i}\right\|_2^2$, the partial derivatives of $z$ for $\mathrm{L}$ is $\frac{\partial L}{\partial z}=\lambda \frac{z-y_{G_i}}{\gamma^2+\left(z-y_{G_i}\right)^2}+\xi\left(z-H x_{G_i}\right)$, and letting $G(z)=$ $\frac{\partial L}{\partial z}=0$, the following equation can be deduced:

$\text{G(}z)=\lambda \frac{z-{{y}_{{{G}_{i}}}}}{{{\gamma }^{2}}+{{(z-{{y}_{{{G}_{i}}}})}^{2}}}+\xi (z-H{{x}_{{{G}_{i}}}})\text{=}0.$         (28)

In order to solve Eq. (28), Newton's method is employed to obtain by iterating the following Eq. (29), which converges after several iterations:

${{z}^{t+1}}={{z}^{t}}-\frac{G({{z}^{t}})}{{{G}^{'}}({{z}^{t}})}.$           (29)

3.4 Solving q and r sub-problems

A comprehensive analysis of this type of problem is provided in study [42]. The Majorization-Minimization (MM) algorithm [43, 44] aims to iteratively solve Eq. (30) in order to address the minimization problem of $q_j=\underset{q_j}{\operatorname{argmin}} \lambda_1 \phi\left(q_j\right)+\frac{\Upsilon_1}{2}\left\|q_j-\left(\nabla_j x_{G_i}-w_j\right)\right\|_2^2$:

$\mathrm{q}_{\mathrm{j}}^{\mathrm{t}+1}=\underset{\mathrm{q}_{\mathrm{j}}}{\arg \min } \frac{\lambda_1}{2}\left\|\Lambda\left(\mathrm{q}_{\mathrm{j}}^{\mathrm{t}}\right) \mathrm{q}_{\mathrm{j}}\right\|_2^2$$+\frac{\Upsilon_1}{2}\left\|\mathrm{q}_{\mathrm{j}}-\left(\nabla_{\mathrm{j}} \mathrm{x}_{\mathrm{G}_{\mathrm{i}}}-\mathrm{w}_{\mathrm{j}}\right)\right\|_2^2$, $\mathrm{j}=1,2,3,4, \mathrm{t}=0,1, \cdots$,          (30)

where, $\Lambda\left(q_j^t\right)$ is the diagonal matrix, $\left[\Lambda\left(q_j^t\right)\right]_{l, l}=\sqrt{\left.\sum_{i_1, i_2=-m_1}^{m_2}\left[\sum_{k_1, k_2=-n_1}^{n_2} \mid q_{t_1-i_1+k_1, t_2-i_2+k_2}^{(j)}\right]^2\right]^{-\frac{1}{2}}}$ are diagonal elements, with $l=\left(t_2-1\right) P+t_1, m_1=\left[\frac{P-1}{2}\right], m_2=\left[\frac{P}{2}\right], n_1=\left[\frac{h-1}{2}\right], n_2=\left[\frac{h}{2}\right], t_1 \in\{1,2, \cdots, P\}$, $t_2 \in\{1,2, \cdots, h\}, \mathrm{P}$ and $\mathrm{h}$ as explained in Section 2.2. The symbol $q_{t_1-i_1+k_1, t_2-i_2+k_2}^{(j)}$ denotes the element values of the $t_1-i_1+k_1$ th row and the $t_2-i_2+k_2$ th column in the matrix $q_j$. The problem (30) is a least square problem. Setting the partial derivative of $q_j$ to 0 , the solution of $q_j$ is obtained:

$q_{j}^{t+1}={{\left( {{\Upsilon }_{1}}I+{{\lambda }_{1}}\Lambda {{\left( q_{j}^{t} \right)}^{T}}\Lambda \left( q_{j}^{t} \right) \right)}^{-1}}{{\Upsilon }_{1}}\left( {{\nabla }_{j}}{{x}_{{{G}_{i}}}}-{{w}_{j}} \right),$          (31)

where, $q_j^0=\nabla_j x_{G_i}-w_j$ and $j=\{1,2,3,4\}$.

Likewise, the solution of r_j can be derived as follows:

$r_{j}^{t+1}={{\left( {{\Upsilon }_{2}}I+{{\lambda }_{0}}\Lambda {{\left( r_{j}^{t} \right)}^{T}}\Lambda \left( r_{j}^{t} \right) \right)}^{-1}}{{\Upsilon }_{2}}\varepsilon {{(w)}_{j}},$           (32)

where, $r_j^0=\varepsilon(w)_j, j=\{1,2, \cdots, 16\}, \Lambda\left(r_j^t\right)$ and $\left[\Lambda\left(r_j^t\right)\right]_{l, l}=\sqrt{\sum_{i_1, i_2=-m_1}^{m_2}\left[\sum_{k_1, k_2=-n_1}^{n_2}\left|r_{t_1-i_1+k_1, t_2-i_2+k_2}^{(j)}\right|^2\right]^{-\frac{1}{2}}}$ are diagonal matrices and diagonal elements, respectively. The terms $1, m_1, m_2, n_1, n_2, t_1, t_2$ have the same meaning as in equation (30) above, and $r_{t_1-i_1+k_1, t_2-i_2+k_2}^{(j)}$ represents the element value of the $t_1-i_1+k_1$ th row and the $t_2-i_2+k_2$ th column in the matrix $r_j$.

3.5 Solving $x_{G_i}$ sub-problem

$x_{G_i}$ sub-problem in (24) is a least square problem. Suppose $L=\min _{x_{G_i}} \frac{\beta}{2}\left\|p-x_{G_i}\right\|_2^2+\frac{r_1}{2} \sum_{j=1}^4 \| q_j-\left(\nabla_j x_{G_i}-\right.$ $\left.w_j\right)\left\|_2^2+\frac{\xi}{2}\right\| z-H x_{G_i} \|_2^2$, the partial derivative $x_{G_i}$ for $L$ is $\frac{\partial L}{\partial x_{G_i}}=\beta\left(x_{G_i}-p\right)+\Upsilon_1 \sum_{j=1}^4 \nabla_j^T\left(\nabla_j x_{G_i}-w_j-q_j\right)+$ $\xi H^T H x_{G_i}-\xi H^T z$. Let $\frac{\partial L}{\partial x}=0$, the closed form solution of the $x_{G_i}$ sub-problem in (24) is as follows:

${{x}_{{{G}_{i}}}}\text{=}{{\left( \beta I+{{\Upsilon }_{1}}\sum\limits_{j=1}^{4}{\nabla _{j}^{T}{{\nabla }_{j}}}+\xi {{H}^{T}}H \right)}^{-1}}\left( \beta p+{{\Upsilon }_{1}}\sum\limits_{j=1}^{4}{\nabla _{j}^{T}({{w}_{j}}+{{q}_{j}})}+\xi {{H}^{T}}z \right).$        (33)

3.6 Solving w₁~w₄ sub-problem

Suppose , the partial derivative w₁ for L is . Let , we can deduce The solution of w₁ is as follows:

$\operatorname{Suppose} L=\min _w \frac{\Upsilon_1}{2} \sum_{j=1}^4\left\|q_j-\left(\nabla_j x_{G_i}-w_j\right)\right\|_2^2+\frac{\Upsilon_2}{2} \sum_{j=1}^{16}\left\|r_j-\varepsilon(w)_j\right\|_2^2$, the partial derivative $w_1$ for L is $\frac{\partial L}{\partial w_1}=$$\Upsilon_1\left(w_1-\nabla_1 x_{G_i}+q_1\right)+\Upsilon_2 \varepsilon^T\left(\varepsilon(w)_1-r_1\right)$. Let $\frac{\partial L}{\partial w_1}=0$, we can deduce $\left(\Upsilon_1+\Upsilon_2 \nabla_1^T \nabla_1+\frac{\Upsilon_2}{2} \nabla_2^T \nabla_2+\frac{\Upsilon_2}{2} \nabla_3^T \nabla_3+\right.$ $\left.\frac{\Upsilon_2}{2} \nabla_4^T \nabla_4\right) w_1=\Upsilon_1\left(\nabla_1 x_{G_i}-q_1\right)+\Upsilon_2 r_1-\frac{\Upsilon_2}{2} \nabla_2^T\left(\nabla_1 w_2-2 r_5\right)-\frac{\Upsilon_2}{2} \nabla_3^T\left(\nabla_1 w_3-2 r_9\right)-\frac{\Upsilon_2}{2} \nabla_4^T\left(\nabla_1 w_4-2 r_{13}\right)$. The solution of $w_1$ is as follows:

$\begin{align}  & {{w}_{1}}={{({{\Upsilon }_{1}}I+{{\Upsilon }_{2}}\nabla _{1}^{T}{{\nabla }_{1}}+\frac{{{\Upsilon }_{2}}}{2}\nabla _{2}^{T}{{\nabla }_{2}}+\frac{{{\Upsilon }_{2}}}{2}\nabla _{3}^{T}{{\nabla }_{3}}+\frac{{{\Upsilon }_{2}}}{2}\nabla _{4}^{T}{{\nabla }_{4}})}^{-1}} \\ & [{{\Upsilon }_{1}}({{\nabla }_{1}}{{x}_{{{G}_{i}}}}-{{q}_{1}})+{{\Upsilon }_{2}}{{r}_{1}}-\frac{{{\Upsilon }_{2}}}{2}\nabla _{2}^{T}({{\nabla }_{1}}{{w}_{2}}-2{{r}_{5}})-\frac{{{\Upsilon }_{2}}}{2}\nabla _{3}^{T}({{\nabla }_{1}}{{w}_{3}}-2{{r}_{9}})-\frac{{{\Upsilon }_{2}}}{2}\nabla _{4}^{T}({{\nabla }_{1}}{{w}_{4}}-2{{r}_{13}})]. \\\end{align}$           (34)

Through the above method, w₂~w₄ can be solved similarly as follows:

$\begin{align}  & {{w}_{2}}={{({{\Upsilon }_{1}}I+\frac{{{\Upsilon }_{2}}}{2}\nabla _{1}^{T}{{\nabla }_{1}}+{{\Upsilon }_{2}}\nabla _{2}^{T}{{\nabla }_{2}}+\frac{{{\Upsilon }_{2}}}{2}\nabla _{3}^{T}{{\nabla }_{3}}+\frac{{{\Upsilon }_{2}}}{2}\nabla _{4}^{T}{{\nabla }_{4}})}^{-1}} \\ & [{{\Upsilon }_{1}}({{\nabla }_{2}}{{x}_{{{G}_{i}}}}-{{q}_{2}})+{{\Upsilon }_{2}}{{r}_{2}}-\frac{{{\Upsilon }_{2}}}{2}\nabla _{1}^{T}({{\nabla }_{2}}{{w}_{1}}-2{{r}_{6}})-\frac{{{\Upsilon }_{2}}}{2}\nabla _{3}^{T}({{\nabla }_{2}}{{w}_{3}}-2{{r}_{10}})-\frac{{{\Upsilon }_{2}}}{2}\nabla _{4}^{T}({{\nabla }_{2}}{{w}_{4}}-2{{r}_{14}})]. \\\end{align}$           (35)

$\begin{align}  & {{w}_{3}}={{({{\Upsilon }_{1}}I+\frac{{{\Upsilon }_{2}}}{2}\nabla _{1}^{T}{{\nabla }_{1}}+\frac{{{\Upsilon }_{2}}}{2}\nabla _{2}^{T}{{\nabla }_{2}}+{{\Upsilon }_{2}}\nabla _{3}^{T}{{\nabla }_{3}}+\frac{{{\Upsilon }_{2}}}{2}\nabla _{4}^{T}{{\nabla }_{4}})}^{-1}} \\ & [{{\Upsilon }_{1}}({{\nabla }_{3}}{{x}_{{{G}_{i}}}}-{{q}_{3}})+{{\Upsilon }_{2}}{{r}_{3}}-\frac{{{\Upsilon }_{2}}}{2}\nabla _{1}^{T}({{\nabla }_{3}}{{w}_{1}}-2{{r}_{7}})-\frac{{{\Upsilon }_{2}}}{2}\nabla _{2}^{T}({{\nabla }_{3}}{{w}_{2}}-2{{r}_{11}})-\frac{{{\Upsilon }_{2}}}{2}\nabla _{4}^{T}({{\nabla }_{3}}{{w}_{4}}-2{{r}_{15}})]. \\\end{align}$         (36)

$\begin{align}  & {{w}_{4}}={{({{\Upsilon }_{1}}I+\frac{{{\Upsilon }_{2}}}{2}\nabla _{1}^{T}{{\nabla }_{1}}+\frac{{{\Upsilon }_{2}}}{2}\nabla _{2}^{T}{{\nabla }_{2}}+\frac{{{\Upsilon }_{2}}}{2}\nabla _{3}^{T}{{\nabla }_{3}}+{{\Upsilon }_{2}}\nabla _{4}^{T}{{\nabla }_{4}})}^{-1}} \\ & [{{\Upsilon }_{1}}({{\nabla }_{4}}{{x}_{{{G}_{i}}}}-{{q}_{4}})+{{\Upsilon }_{2}}{{r}_{4}}-\frac{{{\Upsilon }_{2}}}{2}\nabla _{1}^{T}({{\nabla }_{4}}{{w}_{1}}-2{{r}_{8}})-\frac{{{\Upsilon }_{2}}}{2}\nabla _{2}^{T}({{\nabla }_{4}}{{w}_{2}}-2{{r}_{12}})-\frac{{{\Upsilon }_{2}}}{2}\nabla _{3}^{T}({{\nabla }_{4}}{{w}_{3}}-2{{r}_{16}})]. \\\end{align}$          (37)

Ultimately, to ensure that the solution of Eq. (18) converges to the solution of Eq. (17), and that the solution of Eq. (17) is consistent with the solution of Eq. (16), the value of the parameter $\beta, \Upsilon_1, \Upsilon_2, \xi$ should be set to a very large value. However, setting the parameter values to very large values initially may lead to numerical stability issues (as discussed in Chapter 17 of study [45]). Therefore, following the idea of the NAMA method [46], the parameter values are set to small values at the beginning and are gradually increased during the iteration process, enabling the solution of the model to converge to Eq. (17). Algorithm 1 provides a summary of the entire algorithm for solving model (16).

In this study, a model fusing the prior knowledge of GSR and MDTGV is proposed to restore images contaminated with Cauchy noise and/or blur. The model is solved using a penalty method, variable splitting strategy, and alternating minimization scheme. The prior knowledge of GSR leverages the nonlocal self-similarity of images by considering the non-zero coefficients appearing in the form of clustering in sparse representation signals. This approach takes into account the sparsity of the group structure, preserves more geometric structures, and requires less computation time than global sparse representation. MDTGV regularization can describe the edge information in 8 directions of the image and reconstruct clearer detailed features, resulting in ideal visual results and better visual quality than TGV. Experimental results demonstrate that the proposed method outperforms comparison methods in terms of PSNR and SSIM values in most cases. It is important to note that in these methods, the Cauchy noise level and blur kernel are known; however, in real corrupted images, the Cauchy noise level and blur kernel may be unknown.

To address this limitation, future work will focus on improving the GSR + MDTGV method to adapt to nonparametric blind super-resolution, where the noise level and blur kernel are unknown. Initially, methods will be employed to estimate the noise level and blur kernel size. Subsequently, given the estimated noise level and blur kernel, the SR method based on GSR+MDTGV will be used to super-resolve the final high-resolution image. This extension of the proposed method has the potential to enhance its applicability in various real-world scenarios where the noise level and blur kernel are not readily available.