© Wuhan University 2024

0 Introduction

Image deblurring is a classical problem in low-level vision. Image deblurring aims to estimate the latent image from its degraded observation. A blurred image can be considered as a convolution between the latent image and the point spread function (PSF) of the image capture device. The blurring process can be modeled as

$\mathit{g}=\mathit{H}\mathit{f}+\mathit{n}$(1)

where $\mathit{f}$ is an original image without any form of degradation, $\mathit{H}$ represents the PSF, $\mathit{g}$ is the observation, and $\mathit{n}$ is the additive noise. When the PSF is given, this is known as the "non-blind" image deblurring problem; otherwise, it is called the "blind" deblurring.

Image deblurring is an ill-posed inverse problem, meaning that the solution does not exist and is not unique. Additionally, small perturbations in the data can result in significant errors during reconstruction. Therefore, regularization techniques are essential for transforming this problem into a well-posed one. Over the past three decades, numerous methods have been developed to address image deblurring. Due to the ill-posed nature of the image deblurring problem, it has been widely recognized that the prior information of images and the blurring process play a key role in enhancing the performance of image deblurring methods. A variety of image prior models have been developed, such as the statistical priors-based method^[1,2], the transform-based method^[3-7], the total variation (TV) based method ^[8-13], the dark channel-based method^[14-16], class-specific prior based method^[17], decoupled method^[18] etc. Image prior information is often enforced through regularization strategy in the deblurring process, which can be described in the following variational regularization model

$\underset{\mathit{f}}{\mathrm{m}\mathrm{i}\mathrm{n}}\left\{\frac{\mathrm{1}}{\mathrm{2}}{\Vert \mathit{g}-\mathit{H}\mathit{f}\Vert }_{\mathrm{2}}^{\mathrm{2}}+w\phi (\mathit{f})\right\}$(2)

where ${\Vert \cdot \Vert }_{\mathrm{2}}^{}$ denotes the Euclidean norm, $\phi (\mathit{f})$ is conventionally called a regularization functional, and $w$>0 is referred to as a regularization parameter, which controls the balance between fidelity term and regularization term in model (2).

Choosing a suitable regularization functional is an active area of research in imaging science. Total variation, a classical regularization term, has verified its great success in image deblurring. Indeed, the TV norm has played a crucial role in various image processing tasks such as denoising, deblurring, and inpainting. Wang and Yang et al ^[13] proposed a state-of-the-art algorithm for image deblurring called fast total variation deconvolution (FTVd), which utilizes splitting techniques. In comparison with FTVd, the alternating direction method of multipliers (ADMM) TV method can obtain a similar peak signal-to-noise ratio (PSNR), while the choice of the regularization parameter is not needed.

TV norm serves as a highly effective signal prior to piece-wise smooth images. However, TV is an isotropic function and is unsuitable for images with a dominant direction. The main disadvantage of TV-based regularization for natural images is that TV regularization terms do not preserve the details and textures very well on the regions of complex structures. In various applications such as seismic and medical imaging, image textures exhibit clear directionality, which is crucial in analyzing detailed information. Therefore, it is highly desirable to incorporate the directional characteristics of these textures to ensure high-quality image outcomes. In this paper, we propose a novel approach to enhance the directional information of images in the gradient domain by extending the standard TV norm to the STV norm.

Another regularization technique for deblurring involves incorporating sparsity-based priors to regularize images, leveraging the observation that images typically exhibit sparse representation in certain redundant transform domains, such as wavelet and framelet transforms. The framelet-based approach for super-resolution reconstruction was initially proposed in Refs. [19,20]. The framelet deblurring algorithm was initially examined in Ref.[21]. In Ref.[10], Jiao et al also proposed an efficient algorithm using the ADMM strategy to solve the frame-based deblurring problem.

The core contributions are summarized as follows:

1) The incorporation of enhanced directionality is achieved by proposing a novel shear gradient (SG) operator, which combines the gradient operator and shear operator. Through experiments, it is verified that the proposed SG operator can detect more abundant edge information than that detected by the gradient operator. We also explore the properties of the SG operator and demonstrate that its matrix form exhibits a Block-Circulant-with-Circulant-Blocks (BCCB) structure, rendering it diagonalizable through the 2D discrete Fourier transform.

2) Based on the SG operator, we extend the conventional total variation to the shear total variation to enhance the direction information of the image in the gradient domain.

3) In accordance with the STV norm, we propose a deblurring method that incorporates the STV norm as a regularization term. The proposed problem can be efficiently addressed by employing the ADMM, leveraging the BCCB structure of the shear matrices. We conduct various numerical experiments to demonstrate the potential of our method.

The remaining sections of this paper are structured as follows: Section 1 provides definitions for the shear gradient operator and shear TV norm, along with a brief introduction to the fundamentals of the shear operator, BCCB matrix, and ADMM algorithm. In Section 2, we propose the STV deblurring model and employ ADMM to solve it. The effectiveness of our proposed model is demonstrated through numerical results presented in Section 3. Finally, conclusions are drawn in Section 4.

1 Preliminaries

1.1 Shear Operator

Let ${\mathit{S}}_{{\theta }_d}$ be a shear operator with the shear angle ${\theta }_d$ (see Ref. [22]), and ${\theta }_d\in [\mathrm{0},\mathrm{\pi }]$, $d\in \{‘\mathrm{t}’,‘\mathrm{b}’,‘\mathrm{l}’,‘\mathrm{r}’\}$

${\mathit{S}}_{{\theta }_d}:{\mathbb{R}}^{\mathit{n}\times \mathit{n}}\mapsto {\mathbb{R}}^{\mathit{n}\times \mathit{n}}:X\mapsto X^{\mathrm{s}\mathrm{h}\mathrm{r}}$(3)

where the $(\widehat{i},\widehat{j})$-th entry of $X^{\mathrm{s}\mathrm{h}\mathrm{r}}$ is given, for $\widehat{i}=\mathrm{1},\cdots ,n$ and $\widehat{j}=\mathrm{1},\cdots ,n$, by $X_{\widehat{i},\widehat{j}}^{\mathrm{s}\mathrm{h}\mathrm{r}}:=X_{\overline{i},\overline{j}}$. The parameter $d$ indicates the axis of shear (top $\mathrm{‘}\mathrm{t}\mathrm{’}$, bottom $\mathrm{‘}\mathrm{b}\mathrm{’}$, left $\mathrm{‘}\mathrm{l}\mathrm{’}$, right $\mathrm{‘}\mathrm{r}\mathrm{’}$), and ${\theta }_d$ determines the degree of shear.

$\overline{i}:=\{\begin{array}{c}\left[\left(\widehat{i}-\mathrm{1}+[\frac{\mathrm{4}}{\mathrm{\pi }}\theta \widehat{j}]\right)\mathrm{m}\mathrm{o}\mathrm{d}n\right]+\mathrm{1},\mathrm{i}\mathrm{f}d=‘\mathrm{t}’\\ \begin{array}{l}\left[\left(\widehat{i}-\mathrm{1}+[\frac{\mathrm{4}}{\mathrm{\pi }}\theta \left(n-\widehat{j}\right)]\right)\mathrm{m}\mathrm{o}\mathrm{d}n\right]+\mathrm{1},\mathrm{i}\mathrm{f}d=‘\mathrm{b}’\\ \widehat{i},\mathrm{o}\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{w}\mathrm{i}\mathrm{s}\mathrm{e}\end{array}\end{array}$

$\overline{j}:=\{\begin{array}{c}\left[\left(\widehat{j}-\mathrm{1}+[\frac{\mathrm{4}}{\mathrm{\pi }}\theta \widehat{i}]\right)\mathrm{m}\mathrm{o}\mathrm{d}n\right]+\mathrm{1},\mathrm{i}\mathrm{f}d=‘\mathrm{r}’\\ \left[\left(\widehat{j}-\mathrm{1}+[\frac{\mathrm{4}}{\mathrm{\pi }}\theta \left(n-\widehat{i}\right)]\right)\mathrm{m}\mathrm{o}\mathrm{d}n\right]+\mathrm{1},\mathrm{i}\mathrm{f}d=‘\mathrm{l}’\\ \widehat{j},\mathrm{o}\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{w}\mathrm{i}\mathrm{s}\mathrm{e}\end{array}$

where $[\cdot ]$ represents the procedure rounding the real number to the nearest integer in the zero direction. For example, ${\mathit{S}}_{{\theta }_d}X$ means that the image $X$ is sheared with the top axis and ${\theta }_d$ degree. In particular, ${\theta }_d=\mathrm{0}$ means ${\mathit{S}}_{{\theta }_d}$ is the identity transformation, i.e., no shear. Figure 1 shows the effect of the shear operator in different angles and directions.

Fig. 1 Effect of the shear operator, ${\mathit{S}}_{{\mathit{\theta }}_{\mathit{d}}}$ with different angles and directions

1.2 Shear Gradient and Shear TV Norm

Traditional regularization methods include Tikhonov regularization ^[23] and TV regularization ^[8]. To avoid the Tikhonov regularization term being too smooth in the image restoration process and unable to preserve sharp edges, TV regularization was introduced by Rudin et al^[8] for the image denoising and deconvolution problems. For the grayscale image $\mathit{u}\in {\mathbf{R}}^{\mathit{n}\times \mathit{n}}$, its discrete total variation norm is:

$\mathrm{T}\mathrm{V}(\mathit{u})={{\displaystyle \sum _{i=\mathrm{1}}^{\mathrm{2}}}\Vert {\mathit{D}}^i\mathit{u}\Vert }_{\mathrm{1}}$(4)

where ${\mathit{D}}^{\mathrm{1}}={\mathit{D}}_x$ and ${\mathit{D}}^{\mathrm{2}}={\mathit{D}}_y$ are linear operators corresponding to the horizontal and vertical first-order differences. The TV norm only includes finite differences in the horizontal and vertical directions, which results in a loss of some detailed information in other directions.

By combining the first-order finite difference operator and the shear operator, we propose the shear gradient operator (SG) in the following way:

$\mathit{D}\triangleq ({\mathit{D}}^{\mathrm{3}},{\mathit{D}}^{\mathrm{4}})=({\mathit{S}}_{{\theta }_d}^{\mathrm{T}}{\mathit{D}}_x{\mathit{S}}_{{\theta }_d},{\mathit{S}}_{{\theta }_d}^{\mathrm{T}}{\mathit{D}}_y{\mathit{S}}_{{\theta }_d})$(5)

where the superscript "$\mathrm{T}$" means transpose. Then, for the grayscale image $\mathit{u}$, its discrete shear total variation norm is defined as:

$\mathrm{S}\mathrm{T}\mathrm{V}(\mathit{u})={{\displaystyle \sum _{i=\mathrm{1}}^{\mathrm{4}}}\Vert {\mathit{D}}^i\mathit{u}\Vert }_{\mathrm{1}}$(6)

The shear gradient has the following properties:

Remark 1   The shear operator ${\mathit{S}}_{{\theta }_d}$ and the horizontal gradient operator ${\mathit{D}}_x$ are commutative when $=‘\mathrm{l}’\mathrm{o}\mathrm{r}‘\mathrm{r}’$ ; Similarly, ${\mathit{S}}_{{\theta }_d}$ and the vertical operator ${\mathit{D}}_y$ are commutative when $d=‘\mathrm{t}’\mathrm{o}\mathrm{r}‘\mathrm{b}’$ . i.e.,

${\mathit{S}}_{{\theta }_d}^{\mathrm{T}}{\mathit{D}}_x{\mathit{S}}_{{\theta }_d}={\mathit{D}}_x,\mathrm{w}\mathrm{h}\mathrm{e}\mathrm{n}d=‘\mathrm{l}’\mathrm{o}\mathrm{r}‘\mathrm{r}’$(7)

${\mathit{S}}_{{\theta }_d}^{\mathrm{T}}{\mathit{D}}_y{\mathit{S}}_{{\theta }_d}={\mathit{D}}_y,\mathrm{w}\mathrm{h}\mathrm{e}\mathrm{n}d=‘\mathrm{t}’\mathrm{o}\mathrm{r}‘\mathrm{b}’$(8)

By Remark 1, the shear gradient operator is essentially a generalization of the gradient operator.

Remark 2   (Matrix form of ${\mathit{S}}_{{\theta }_d}$) The operator ${\mathit{S}}_{{\theta }_d}$ is linear so that it can be represented as a matrix denoted by ${\mathit{S}}_{{\theta }_d}$. Moreover, since ${\mathit{S}}_{{\theta }_d}$ is a permutation matrix, it is orthonormal.

The comparison between the gradient operator and the shear gradient operator applied to the Cameraman image is shown in Fig.2.

Fig. 2 Comparison between gradient operator and shear gradient operator applied on Cameraman image

It can be observed that the shear gradient operator captures more direction information and retains more detailed edge information.

1.3 Block-Circulant-with-Circulant-Blocks Matrices

Under the periodic boundary conditions, both the blurring matrix and gradient matrix are BCCB and thus are diagonalizable by the 2D discrete Fourier transform. This section will show that the shear gradient matrix is also a BCCB matrix.

Let ${\mathit{E}}_{ij}$ denote a matrix whose $(i,j)$ entry is 1 and zero elsewhere. Then $\{{\mathit{E}}_{ij}|i,j=\mathrm{1,2},\mathrm{3}\}$ forms a set of orthonormal basis in the linear space ${\mathbb{R}}^{\mathrm{3}\times \mathrm{3}}$. Let ${\mathit{S}}_{{\theta }_l}$ denote the shear matrix of the shear operator under the basis $\{{\mathit{E}}_{ij}|i,j=\mathrm{1,2},\cdots ,n\}$, i.e.

$\begin{array}{l}{\mathit{S}}_{{\theta }_l}\left({\mathit{E}}_{\mathrm{11}},{\mathit{E}}_{\mathrm{21}},{\mathit{E}}_{\mathrm{31}},{\mathit{E}}_{\mathrm{12}},{\mathit{E}}_{\mathrm{22}},{\mathit{E}}_{\mathrm{32}},{\mathit{E}}_{\mathrm{13}},{\mathit{E}}_{\mathrm{23}},{\mathit{E}}_{\mathrm{33}}\right)\\ =({\mathit{E}}_{\mathrm{11}},{\mathit{E}}_{\mathrm{21}},{\mathit{E}}_{\mathrm{31}},{\mathit{E}}_{\mathrm{12}},{\mathit{E}}_{\mathrm{22}},{\mathit{E}}_{\mathrm{32}},{\mathit{E}}_{\mathrm{13}},{\mathit{E}}_{\mathrm{23}},{\mathit{E}}_{\mathrm{33}}){\mathit{S}}_{{\theta }_l}\end{array}$(9)

where ${\theta }_l=\mathrm{45}°$. By Eq. (9), we can obtain that ${\mathit{S}}_{{\theta }_l}$ is a matrix with the following form

${\mathit{S}}_{{\theta }_l}=\left[\begin{array}{ccccccccc}\mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{1}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{1}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{1}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\end{array}\right]$

Since gradient matrices are BCCB under periodic boundary conditions, it can be verified that the ${\mathit{D}}^{\mathrm{3}}$ and ${\mathit{D}}^{\mathrm{4}}$ in formula (5) are also BCCB matrices. The shear gradient matrix ${\mathit{D}}^{\mathrm{4}}={\mathit{S}}_{{\theta }_l}^{\mathrm{T}}{\mathit{D}}_y{\mathit{S}}_{{\theta }_l}$ with ${\theta }_l=\mathrm{45}°$ has the following form:

${\mathit{D}}^{\mathrm{4}}=\left[\begin{array}{ccccccccc}-\mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{1}& \mathrm{0}\\ \mathrm{0}& -\mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{1}\\ \mathrm{0}& \mathrm{0}& -\mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{1}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{1}& \mathrm{0}& -\mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{1}& \mathrm{0}& -\mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& -\mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{1}& \mathrm{0}& -\mathrm{1}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{1}& \mathrm{0}& -\mathrm{1}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& -\mathrm{1}\end{array}\right]$

But unfortunately, when the shear angle ${\theta }_d\in (\mathrm{0,45}°)$,${\mathit{D}}^i(i=\mathrm{1,2},\mathrm{3,4})$ do not always have a BCCB structure. When the shear angle is ${\theta }_l=\mathrm{40}°$, the shear gradient matrix ${\mathit{D}}^{\mathrm{4}}={\mathit{S}}_{{\theta }_l}^{\mathrm{T}}{\mathit{D}}_y{\mathit{S}}_{{\theta }_l}$ has the following form:

${\mathit{D}}^{\mathrm{4}}=\left[\begin{array}{ccccccccc}-\mathrm{1}& \mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& -\mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{1}\\ \mathrm{0}& \mathrm{0}& -\mathrm{1}& \mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{1}& \mathrm{0}& -\mathrm{1}& \mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{1}& \mathrm{0}& -\mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& -\mathrm{1}& \mathrm{1}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{1}& \mathrm{0}& -\mathrm{1}& \mathrm{1}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{1}& \mathrm{0}& -\mathrm{1}& \mathrm{0}\\ \mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{1}& \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}& -\mathrm{1}\end{array}\right]$

It can be seen that the above matrix is not a BCCB matrix.

1.4 Alternating Direction Method of Multipliers

The ADMM was proposed in Ref.[24] and has been analyzed in various literatures^[25,26]. It is widely used to solve structured optimization problems due to its decomposability and superior flexibility. Recently, ADMM has gained popularity in modern signal and image processing, statistics, machine learning, and various other fields. ADMM is an algorithm utilizing proximal splitting techniques for solving the following convex optimization problem:

$\begin{array}{l}\underset{}{\mathrm{m}\mathrm{i}\mathrm{n}}{\theta }_{\mathrm{1}}(x_{\mathrm{1}})+{\theta }_{\mathrm{2}}(x_{\mathrm{2}})\\ \mathrm{s}.\mathrm{t}.{\mathit{A}}_{\mathrm{1}}{x}_{\mathrm{1}}+{\mathit{A}}_{\mathrm{2}}{x}_{\mathrm{2}}=\mathit{c},\\ x_i\in {\mathit{\chi }}_i,i=\mathrm{1,2}\end{array}$(10)

where ${\theta }_i:{\mathit{\chi }}_i\to \mathbb{R}$ are closed convex functions, ${\mathit{A}}_i$ are linear transforms, ${\mathit{\chi }}_i$ are nonempty closed convex sets, and $\mathit{c}$ is a given vec tor.

One of the advantages of ADMM is that the algorithm updates lend themselves to parallel implementations. The algorithm is given in Algorithm 1.

Algorithm 1: ADMM
Initialization: Starting point $(x_1^0,x_2^0,p^0)$, $\delta >0$. Iteration: 1. $x_1^{k+1}=\mathrm{a}\mathrm{r}\mathrm{g}\underset{x_1}{\mathrm{m}\mathrm{i}\mathrm{n}}{\theta }_1(x_1)+\frac{\delta }{2}{\Vert {\mathit{A}}_1{x}_1+{\mathit{A}}_2{x}_2^k-\mathit{c}+\frac{p^k}{\delta }\Vert }_2^2$ 2. $x_2^{k+1}=\mathrm{a}\mathrm{r}\mathrm{g}\underset{x_2}{\mathrm{m}\mathrm{i}\mathrm{n}}{\theta }_2(x_2)+\frac{\delta }{2}{\Vert {\mathit{A}}_1{x}_1^{k+1}+{\mathit{A}}_2{x}_2-\mathit{c}+\frac{p^k}{\delta }\Vert }_2^2$ 3. $p^{k+1}=p^k+\delta ({\mathit{A}}_1{x}_1^{k+1}+{\mathit{A}}_2{x}_2^{k+1}-\mathit{c})$ 4. $k=k+1$
Until a stopping criterion is satisfied.

2 STV Deblurring Model

The discussion in subsection 1.2 shows that the STV norm can capture sharp edges and more detailed information about the image. This leads us to consider the following problem:

$\{\begin{array}{c}\underset{\mathit{v}}{\mathrm{m}\mathrm{i}\mathrm{n}}\psi (\mathit{v})={\displaystyle \sum _{i=\mathrm{1}}^{\mathrm{4}}}{\Vert {\mathit{v}}_i\Vert }_{\mathrm{1}}\\ \mathrm{s}.\mathrm{t}.{\mathit{v}}_{\mathrm{1}}={\mathit{D}}_x\mathit{f}\\ {\mathit{v}}_{\mathrm{2}}={\mathit{D}}_y\mathit{f}\\ {\mathit{v}}_{\mathrm{3}}={\mathit{S}}_{{\theta }_t}^{\mathrm{T}}{\mathit{D}}_x{\mathit{S}}_{{\theta }_t}\mathit{f}\\ {\mathit{v}}_{\mathrm{4}}={\mathit{S}}_{{\theta }_l}^{\mathrm{T}}{\mathit{D}}_y{\mathit{S}}_{{\theta }_l}\mathit{f}\\ \mathit{H}\mathit{f}=\mathit{g}\end{array}$(11)

the corresponding augmented Lagrangian function of (11) can be written as:

$\begin{array}{l}{\mathrm{ℒ}}_{{\rho }_{\mathrm{1}},{\rho }_{\mathrm{2}}}(\mathit{f},\mathit{v};\mathit{\lambda },\mathit{\mu })=\psi (\mathit{v})+〈\mathit{\lambda },\mathit{H}\mathit{f}-\mathit{g}〉+{\displaystyle \sum _{i=\mathrm{1}}^{\mathrm{4}}}〈{\mathit{\mu }}_i,{\mathit{D}}^i\mathit{f}-{\mathit{v}}_i〉\\ +\frac{{\rho }_{\mathrm{1}}}{\mathrm{2}}{\Vert \mathit{H}\mathit{f}-\mathit{g}\Vert }^{\mathrm{2}}+\frac{{\rho }_{\mathrm{2}}}{\mathrm{2}}{{\displaystyle \sum _{i=\mathrm{1}}^{\mathrm{4}}}\Vert {\mathit{D}}^i\mathit{f}-{\mathit{v}}_i\Vert }^{\mathrm{2}}\end{array}$(12)

where ${\rho }_{\mathrm{1}}$ and ${\rho }_{\mathrm{2}}$ are two positive parameters. We decompose the optimization problem (12) into individual sub-problems and minimize them separately. Our ADMM algorithm then reads as follows:

$\{\begin{array}{l}{\mathit{f}}^{k+\mathrm{1}}=\mathrm{a}\mathrm{r}\mathrm{g}\underset{\mathit{f}}{\mathrm{m}\mathrm{i}\mathrm{n}}{\mathrm{ℒ}}_{{\rho }_{\mathrm{1}},{\rho }_{\mathrm{2}}}\left(\mathit{f},{\mathit{v}}^k;{\mathit{\lambda }}^k,{\mathit{\mu }}^k\right)\\ {\mathit{v}}_i^{k+\mathrm{1}}=\mathrm{a}\mathrm{r}\mathrm{g}\underset{{\mathit{v}}_i}{\mathrm{m}\mathrm{i}\mathrm{n}}{\mathrm{ℒ}}_{{\rho }_{\mathrm{1}},{\rho }_{\mathrm{2}}}\left({\mathit{f}}^{k+\mathrm{1}},\mathit{v};{\mathit{\lambda }}^k,{\mathit{\mu }}^k\right)\\ {\mathit{\lambda }}^{k+\mathrm{1}}={\mathit{\lambda }}^k+{\rho }_{\mathrm{1}}(\mathit{H}{\mathit{f}}^{k+\mathrm{1}}-\mathit{g})\\ {\mathit{\mu }}_i^{k+\mathrm{1}}={\mathit{\mu }}_i^k+{\rho }_{\mathrm{2}}({\mathit{D}}^i{\mathit{f}}^{k+\mathrm{1}}-{\mathit{v}}_i^{k+\mathrm{1}})\end{array}$(13)

Let us now examine each of these sub-problems one by one.

1) $\mathit{f}$ sub-problem

For fixed ${\mathit{\lambda }}^k$, ${\mathit{v}}^k$, and ${\mathit{\mu }}^k$, the minimization of ${\mathrm{ℒ}}_{{\rho }_{\mathrm{1}},{\rho }_{\mathrm{2}}}(\mathit{f},\mathit{v};\mathit{\lambda },\mathit{\mu })$ with respect to $\mathit{f}$ is a least squares problem, and the corresponding normal equations are:

$\begin{array}{l}({\rho }_{\mathrm{1}}{\mathit{H}}^{\mathrm{T}}\mathit{H}+{\rho }_{\mathrm{2}}{\displaystyle \sum _{i=\mathrm{1}}^{\mathrm{4}}}{{\mathit{D}}^i}^{\mathrm{T}}{\mathit{D}}^i)\mathit{f}\\ ={\rho }_{\mathrm{1}}{\mathit{H}}^{\mathrm{T}}\mathit{g}+{\displaystyle \sum _{i=\mathrm{1}}^{\mathrm{4}}}{{\mathit{D}}^i}^{\mathrm{T}}({\rho }_{\mathrm{2}}{\mathit{v}}_i^k-{\mathit{\mu }}_i^k)-{\mathit{H}}^{\mathrm{T}}{\mathit{\lambda }}^k\end{array}$(14)

Considering what has been discussed in subsection 1.3, both ${\mathit{H}}^{\mathrm{T}}\mathit{H}$ and ${{\mathit{D}}^i}^{\mathrm{T}}{\mathit{D}}^i(i=\mathrm{1,2},\mathrm{3,4})$ are BCCB matrices under the periodic boundary condition with a shear angle of $\mathrm{45}°$. Fast Fourier transforms (FFTs) can effectively calculate the BCCB matrix ^[28]. First, concerning formula (14), we compute the right-hand side vector and apply a forward FFT. Second, by component-wisely dividing the eigenvalues of ${\rho }_{\mathrm{1}}{\mathit{H}}^{\mathrm{T}}\mathit{H}+{\rho }_{\mathrm{2}}{\displaystyle \sum _{i=\mathrm{1}}^{\mathrm{4}}}{{\mathit{D}}^i}^{\mathrm{T}}{\mathit{D}}^i$ we get $\mathrm{ℱ}{\mathit{f}}^{k+\mathrm{1}}$, where $\mathrm{ℱ}$ represents the two-dimensional discrete Fourier transform. Third, we apply the inverse FFT to $\mathrm{ℱ}{\mathit{f}}^{k+\mathrm{1}}$ to get a new iteration of ${\mathit{f}}^{k+\mathrm{1}}$.

From the analysis of the shear gradient matrix in Section 1.3, we know that ${{\mathit{D}}^i}^{\mathrm{T}}{\mathit{D}}^i(i=\mathrm{1,2},\mathrm{3,4})$ does not have the BCCB matrix property when the shear angle is not $\mathrm{45}°$. Thus, minimizing ${\mathrm{ℒ}}_{{\rho }_{\mathrm{1}},{\rho }_{\mathrm{2}}}(\mathit{f},\mathit{v};\mathit{\lambda },\mathit{\mu })$ with $\mathit{f}$ does not have a simple closed form solution. Inspired by Ref. [27], one way to overcome this challenge is by adding an approximation term representing a semi-norm induced by a suitable positive semi-definite matrix. Then, the objective function of ${\mathit{f}}^{k+\mathrm{1}}$ in (13) can be written as follows:

${\mathit{f}}^{k+\mathrm{1}}=\mathrm{a}\mathrm{r}\mathrm{g}\underset{\mathit{f}}{\mathrm{m}\mathrm{i}\mathrm{n}}\{{\mathrm{ℒ}}_{{\rho }_{\mathrm{1}},{\rho }_{\mathrm{2}}}\left(\mathit{f},{\mathit{v}}^k;{\mathit{\lambda }}^k,{\mathit{\mu }}^k\right)+\frac{\mathrm{1}}{\mathrm{2}}{\Vert \mathit{f}-{\mathit{f}}^k\Vert }_{\mathit{S}}^{\mathrm{2}}\}$(15)

where $\mathit{S}$ is a positive semi-definite matrix and ${\Vert \cdot \Vert }_{\mathit{S}}$ is the semi-norm induced by the semi-inner product $\sqrt[]{{\mathit{f}}^{\mathrm{T}}\mathit{S}\mathit{f}}$. The solution to equation (15) can then be derived from the normal equation:

$\begin{array}{l}({\rho }_{\mathrm{1}}{\mathit{H}}^{\mathrm{T}}\mathit{H}+{\rho }_{\mathrm{2}}{\displaystyle \sum _{i=\mathrm{1}}^{\mathrm{4}}}{{\mathit{D}}^i}^{\mathrm{T}}{\mathit{D}}^i+\mathit{S})\mathit{f}\\ ={\rho }_{\mathrm{1}}{\mathit{H}}^{\mathrm{T}}\mathit{g}+{\displaystyle \sum _{i=\mathrm{1}}^{\mathrm{4}}}{{\mathit{D}}^i}^{\mathrm{T}}({\rho }_{\mathrm{2}}{\mathit{v}}_i^k-{\mathit{\mu }}_i^k)-{\mathit{H}}^{\mathrm{T}}{\mathit{\lambda }}^k+\mathit{S}{\mathit{f}}^k\end{array}$(16)

As a result,${\rho }_{\mathrm{1}}{\mathit{H}}^{\mathrm{T}}\mathit{H}+{\rho }_{\mathrm{2}}{\displaystyle \sum _{i=\mathrm{1}}^{\mathrm{4}}}{{\mathit{D}}^i}^{\mathrm{T}}{\mathit{D}}^i$ is positive semi-definite, we choose $\mathit{S}=\alpha I-({\rho }_{\mathrm{1}}{\mathit{H}}^{\mathrm{T}}\mathit{H}+{\rho }_{\mathrm{2}}{\displaystyle \sum _{i=\mathrm{1}}^{\mathrm{4}}}{{\mathit{D}}^i}^{\mathrm{T}}{\mathit{D}}^i)$ with $\alpha \ge {\sigma }_{\mathrm{m}\mathrm{a}\mathrm{x}}({\rho }_{\mathrm{1}}{\mathit{H}}^{\mathrm{T}}\mathit{H}+{\rho }_{\mathrm{2}}{\displaystyle \sum _{i=\mathrm{1}}^{\mathrm{4}}}{{\mathit{D}}^i}^{\mathrm{T}}{\mathit{D}}^i)$. Therefore, (16) can be rewritten as:

${\mathit{f}}^{k+\mathrm{1}}=\frac{{\rho }_{\mathrm{1}}{\mathit{H}}^{\mathrm{T}}\mathit{g}+{\displaystyle \sum _{i=\mathrm{1}}^{\mathrm{4}}}{{\mathit{D}}^i}^{\mathrm{T}}({\rho }_{\mathrm{2}}{\mathit{v}}_i^k-{\mathit{\mu }}_i^k)-{\mathit{H}}^{\mathrm{T}}{\mathit{\lambda }}^k+\mathit{S}{\mathit{f}}^k}{\alpha }$(17)

2) $\mathit{v}$ sub-problem

The $\mathit{v}$ sub-problem in (13) is decoupled so that it can be solved separately. We have:

${\mathit{v}}_i^{k+\mathrm{1}}=\mathrm{m}\mathrm{a}\mathrm{x}\left\{\Vert {\mathit{D}}^i{\mathit{f}}^{k+\mathrm{1}}+\frac{{\mathit{\mu }}_i^k}{{\rho }_{\mathrm{2}}}\Vert -\frac{\mathrm{1}}{{\rho }_{\mathrm{2}}},\mathrm{0}\right\}\circ \mathrm{s}\mathrm{g}\mathrm{n}({\mathit{D}}^i{\mathit{f}}^{k+\mathrm{1}}+\frac{{\mathit{\mu }}_i^k}{{\rho }_{\mathrm{2}}})$(18)

where "$\circ$" and "$\mathrm{s}\mathrm{g}\mathrm{n}$" represent the point-wise product and signum function, respectively.

3) $\mathit{\lambda }$ and $\mathit{\mu }$ update

We update $\mathit{\lambda }$ and $\mathit{\mu }$ by

${\mathit{\lambda }}^{k+\mathrm{1}}={\mathit{\lambda }}^k+{\rho }_{\mathrm{1}}(\mathit{H}{\mathit{f}}^{k+\mathrm{1}}-\mathit{g})$(19)

${\mathit{\mu }}_i^{k+\mathrm{1}}={\mathit{\mu }}_i^k+{\rho }_{\mathrm{2}}({\mathit{D}}^i{\mathit{f}}^{k+\mathrm{1}}-{\mathit{v}}_i^{k+\mathrm{1}})$(20)

where $i=\mathrm{1,2},\mathrm{3,4}$.

For simplicity, we terminate the algorithm with the relative change of $\mathit{f}$ in all experiments, i.e.,

$\frac{\Vert {\mathit{f}}^{k+\mathrm{1}}-{\mathit{f}}^k\Vert }{\mathrm{m}\mathrm{a}\mathrm{x}\left\{\Vert {\mathit{f}}^k\Vert ,\mathrm{1}\right\}}<ϵ$(21)

Based on the above discussion, we get the following algorithm to solve the problem (11): the ADMM STV algorithm.

Obviously, Algorithm 2 is an instance of ADMM. The convergence of the proposed algorithm is guaranteed if the sub-problem has closed-form solutions.

Algorithm 2: ADMM STV
Initialization: Starting point $({\mathit{f}}^0,{\mathit{v}}_i^0,{\mathit{\lambda }}^0,{\mathit{\mu }}_i^0)$, ${\rho }_1>0,{\rho }_2>0$. Iteration: 1. Compute ${\mathit{f}}^{k+1}$ according to (14) or (17). 2. Compute ${\mathit{v}}_i^{k+1}$ according to (18). 3. Update ${\mathit{\lambda }}^{k+1}$ according to (19). 4. Update ${\mathit{\mu }}_i^{k+1}$ according to (20). 5. $k=k+1$
Until the stopping criterion (21) is satisfied.

3 Experiments

In this section, we present various numerical results to illustrate the performance of the proposed algorithm. All experiments were conducted on a laptop equipped with an AMD Ryzen 7840H CPU and 16 GB of RAM, using Windows 11 (64-bit) and MATLAB 2022b.

3.1 Parameter Selection

The proposed algorithm requires careful adjustment of several parameters. The parameter ${\rho }_{\mathrm{1}}$ controls the fidelity of the data. We observe that the value of ${\rho }_{\mathrm{1}}$ plays an important role in removing the additive noise in the image. However, relying excessively on ${\rho }_{\mathrm{1}}$ to filter out noise can lead to the occurrence of blocking and speckle artifacts in the restored image. To this end, we apply our algorithm to the Monarch deblurring problem on the grid ${\rho }_{\mathrm{1}}$ = 50:50:1 300, ${\rho }_{\mathrm{2}}$ = 100:20:600. In Fig.3, we can observe that the values of ${\rho }_{\mathrm{1}}$ and ${\rho }_{\mathrm{2}}$ can be selected from a wide range while still maintaining stable reconstructions. For the shear angle ${\theta }_d$, through experiments, we find that better PSNR and SSIM values are often obtained when ${\theta }_d$ is set to 45° or its multiples. Therefore, we set ${\theta }_d=k\frac{\mathrm{\pi }}{\mathrm{4}}$ ($k=\mathrm{1,2},\mathrm{3,4}$) in the experiment and select the parameters that yielded the most optimal results.

Fig. 3 PSNR results of image Monarch by ADMM STV on the grid ${\mathit{\rho }}_{\mathrm{1}}$ = 50:50:1 300, ${\mathit{\rho }}_{\mathrm{2}}$ = 100:20:600

3.2 Comparisons with Other State-of-the-Art Methods

Figure 4 shows the test images and their corresponding blurred versions. The test images are blurred by various kernels, including motion, Gaussian, and average. Fig. 4(h)-(n) show images blurred by kernels generated using the MATLAB fspecial (̔motion',35,50), fspecial (̔average',10) and fspecial (̔Gaussian',[20 20], 10) , respectively.

Fig. 4 Some of the test images (the first row) and the corresponding blurred images (the second row)

In Fig.5 and Fig.6, we compare the deblurring performance of the proposed method with the FTVd, ADMM TV, and ADMM frame. One can observe that the images produced by the ADMM STV algorithm have fewer staircase artifacts and are visually much better than those obtained by FTVd, ADMM TV, and ADMM frame. From the local detail enlargement in Fig.5, one can observe that the image obtained by ADMM TV seems to be blockier compared with the proposed method. The restored image obtained by FTVd is too cartoon-like. The reconstructed image edges by the ADMM frame method are hampered by ring artifacts, also known as the Gibbs phenomenon. This indicates that the proposed model performs better for deblurring images while preserving edges and small details. The curves of PSNR value versus the number of iterations are illustrated in Fig. 7. It can be observed that as the number of iterations increases, the PSNR curve gradually increases and finally becomes flat and stable. In other words, our proposed algorithm shows good stability in terms of convergence.

Fig. 5 Visual comparison of Monarch by different methods with gaussian blur(the noise level is 0.256)

Fig. 6 Visual comparison of Pirate by different methods with average blur (the noise level is 0.512)

Fig. 7 PSNR results versus iteration number by ADMM STV algorithm

The results in terms of PSNR, SSIM, and CPU time are listed in Table 1. We report the PSNR and SSIM results for a collection of 9 test images for all competing methods. It is evident that the proposed ADMM STV outperforms other competing methods in terms of PSNR in most cases. ADMM TV only contains finite differences in two directions, while ADMM STV contains finite differences in four directions. Therefore, the cost of ADMM STV is slightly more expensive than that of ADMM TV, but it better preserves more details, such as directional texture, by the ADMM STV method.

4 Conclusion

The shear total variation norm is proposed in this paper to effectively preserve sharp image edges and enhance the direction information of the image in the gradient domain. It has been demonstrated as a highly efficient solution for non-blind image deblurring problems. Thanks to the BCCB structure of the shear gradient operator and blurring matrices, the proposed model can be efficiently solved using the ADMM algorithm. Comparative evaluations against three state-of-the-art non-blind image deblurring methods clearly demonstrate that our proposed method outperforms both subjectively and objectively.

References

All Tables

All Figures

	Fig. 1 Effect of the shear operator, ${\mathit{S}}_{{\mathit{\theta }}_{\mathit{d}}}$ with different angles and directions
In the text

	Fig. 2 Comparison between gradient operator and shear gradient operator applied on Cameraman image
In the text

	Fig. 3 PSNR results of image Monarch by ADMM STV on the grid ${\mathit{\rho }}_{\mathrm{1}}$ = 50:50:1 300, ${\mathit{\rho }}_{\mathrm{2}}$ = 100:20:600
In the text

	Fig. 4 Some of the test images (the first row) and the corresponding blurred images (the second row)
In the text

	Fig. 5 Visual comparison of Monarch by different methods with gaussian blur(the noise level is 0.256)
In the text

	Fig. 6 Visual comparison of Pirate by different methods with average blur (the noise level is 0.512)
In the text

	Fig. 7 PSNR results versus iteration number by ADMM STV algorithm
In the text