A corner-weighted bounded Hessian model for image denoising

Abstract

Image denoising is an essential step in the image processing task. The first-order variational model can remove noise while preserving edges, but it also generates the staircase effect. Although the bounded Hessian regulariser can reduce this side effect, it tends to blur object edges. In this paper, we propose a corner-weighted bounded Hessian model (CWBH) for image denoising, which has capability of removing noise without causing blurring object edges and artifacts. The bounded Hessian regularization at each pixel is controlled by a weight function which has an exponential form and depends on the corner response of the pixel. The split Bregman algorithm is adapted to decompose the proposed minimization problem into several subproblems which are solved directly using fast Fourier transform and the shrinkage operators. The proposed model is evaluated on synthetic and real noisy images for both spatially invariant and variant additive white Gaussian noise (AWGN). Extensive experiments demonstrate that our proposed model outperforms some state-of-the-art variational models for various types of noise and images. For uniform AWGN, CWBH surpasses other models on average by 0.014 for SSIM and by 0.77dB for PSNR; for spatially variant AWGN, these figures are 0.033 and 0.89dB, respectively.

Appendix: Discretisation scheme

In this section, we describe the discretisation schemes for differential operators and norms. For more details, we refer the reader to [20, 23, 24]. In this paper, we use the periodic boundary condition for operators so that the FFT can be applied to solve the u-subproblem (17) efficiently. To define the discrete forms of differential operators, we first provide the forward and the backward difference operators.

$$ \begin{array}{l} D_{x}^{+} u_{i,j} = \begin{cases} u_{i,j+1}-u_{i,j} &\text{if} 1 \leq i \leq N, 1\leq j<M, \\ u_{i,1}-u_{i,j} &\text{if} 1\leq i\leq N, j=M, \end{cases}, \\ D_{y}^{+} u_{i,j} = \begin{cases} u_{i+1,j}-u_{i,j} &\text{if} 1\leq i<N, 1\leq j\leq M, \\ u_{1,j}-u_{i,j} &\text{if} i=N, 1\leq j\leq M. \end{cases}, \end{array} $$

(35)

$$ \begin{array}{l} D_{x}^{-} u_{i,j} = \begin{cases} u_{i,j}-u_{i,j-1} &\text{if} 1 \leq i \leq N, 1 < j\leq M, \\ u_{i,j}-u_{i,M} &\text{if} 1\leq i\leq N, j=1, \end{cases}, \\ D_{y}^{-} u_{i,j} = \begin{cases} u_{i,j}-u_{i-1,j} &\text{if} 1 < i\leq N, 1\leq j\leq M, \\ u_{i,j}-u_{N,j} &\text{if} i=1, 1\leq j\leq M. \end{cases}, \end{array} $$

(36)

where $u \in \mathbb {R}^{N\times M}$.

The first, second, and fourth order differential operators are expressed as follows. For simplicity, we do not list the boundary conditions.

• the gradient operator:

  $$ \nabla u = \left( D_{x}^{+} u_{i,j}, D_{y}^{+} u_{i,j} \right), $$

  (37)

• the divergence operator is defined by analogy with the continuous case by −div(p) ⋅ u = p ⋅∇u, where $p \in \left (\mathbb {R}^{N\times M}\right )^{2}$ and “⋅” denotes the Euclidean inner product:

  $$ \text{div}(p)_{i,j}= D_{x}^{-} \left( p_{1} \right)_{i,j}+ D_{y}^{-} \left( p_{2} \right)_{i,j}. $$

  (38)

• the Laplace operator:

  $$ {\varDelta} u_{i,j} = \text{div}(\nabla u)_{i,j} = D_{x}^{+} D_{x}^{-} u_{i,j} + D_{y}^{-} D_{y}^{+} u_{i,j}. $$

  (39)

• the second-order operator $\text {div}^{2}: \left (\mathbb {R}^{N\times M} \right )^{2} \rightarrow \mathbb {R}^{N\times M}$ is related to the Hessian operator via the adjointness property as div²(q) ⋅ u = q ⋅∇²u, where $q \in \left (\mathbb {R}^{N\times M}\right )^{4}$:

  $$ \text{div}^{2} (q)_{i,j} = D_{x}^{-} D_{x}^{+} q_{1_{i,j}} + D_{y}^{-} D_{x}^{-} q_{2_{i,j}} + D_{x}^{-} D_{y}^{-} q_{3_{i,j}} + D_{y}^{+} D_{y}^{-} q_{4_{i,j}}. $$

  (40)

• the Hessian operator:

  $$ \nabla^{2} u_{i,j} = \begin{pmatrix} D_{x}^{-} D_{x}^{+} u_{i,j} & D_{y}^{+} D_{x}^{+} u_{i,j} \\ D_{x}^{+} D_{y}^{+} u_{i,j} & D_{y}^{-} D_{y}^{+} u_{i,j} \end{pmatrix}. $$

  (41)

• the fourth-order differential operator:

  $$ \begin{array}{ll} \text{div}^{2} \left( \nabla^{2} u\right)_{i,j} = & D_{x}^{+} D_{x}^{-} D_{x}^{-} D_{x}^{+} u_{i,j} + D_{y}^{-} D_{x}^{-} D_{y}^{+} D_{x}^{+} u_{i,j}\\ & + D_{x}^{-} D_{y}^{-} D_{x}^{+} D_{y}^{+} u_{i,j} + D_{y}^{+} D_{y}^{-} D_{y}^{-} D_{y}^{+} u_{i,j}. \end{array} $$

  (42)

By applying the forward and backward difference operators (35), (36) to differential operators (37)–(42), one can obtain final formulas.

The L¹ and L² norms are defined as

$$ \begin{array}{@{}rcl@{}} \Vert \omega \Vert_{1} = \sum\limits_{i=1}^{N} \sum\limits_{j=1}^{M} \left( \sum\limits_{k=1}^{K} \left( \omega_{k} \right)_{i,j}^{2} \right)^{1/2}, \end{array} $$

(43)

$$ \begin{array}{@{}rcl@{}} \Vert \omega \Vert_{2} = \left( \sum\limits_{i=1}^{N} \sum\limits_{j=1}^{M} \sum\limits_{k=1}^{K} \left( \omega_{k} \right)_{i,j}^{2} \right)^{1/2}, \end{array} $$

(44)

where $\omega \in (\mathbb {R}^{N\times M})^{K}$.