Beyond Staircasing Effect: Robust Image Smoothing via ℓ₀ Gradient Minimization and Novel Gradient Constraints ^†

Abstract

In this paper, we propose robust image-smoothing methods based on ${\mathit{\ell }}_0$ gradient minimization with novel gradient constraints to effectively suppress pseudo-edges. Simultaneously minimizing the ${\mathit{\ell }}_0$ gradient, i.e., the number of nonzero gradients in an image, and the ${\mathit{\ell }}_2$ data fidelity results in a smooth image. However, this optimization often leads to undesirable artifacts, such as pseudo-edges, known as the “staircasing effect”, and halos, which become more visible in image enhancement tasks, like detail enhancement and tone mapping. To address these issues, we introduce two types of gradient constraints: box and ball. These constraints are applied using a reference image (e.g., the input image is used as a reference for image smoothing) to suppress pseudo-edges in homogeneous regions and the blurring effect around strong edges. We also present an ${\mathit{\ell }}_0$ gradient minimization problem based on the box-/ball-type gradient constraints using an alternating direction method of multipliers (ADMM). Experimental results on important applications of ${\mathit{\ell }}_0$ gradient minimization demonstrate the advantages of our proposed methods compared to existing ${\mathit{\ell }}_0$ gradient-based approaches.

1. Introduction

Image smoothing is a key technique in image processing and is used in applications such as deblurring [1,2,3,4], detail enhancement [5,6], tone mapping [7,8,9], and so on. Furthermore, it is used in various fields, such as in medical imaging [10,11,12], computer graphics [13,14,15], remote sensing [16,17,18], and character recognition [4,19,20].

Filtering-based smoothing methods, e.g., a bilateral filter [6,21], a guided filter [22,23,24,25], and a nonlocal mean filter [26,27], have been actively studied for a long time and are often used in practical situations, as they can easily obtain smooth images that roughly maintain the structural gradients of input images with a relatively low computational cost.

Optimization-based smoothing methods were recently proposed in [1,28,29,30,31,32,33,34,35,36,37,38,39,40,41], and they can flexibly incorporate a priori information into minimization problems. The most standard a priori knowledge of natural images is local smoothness. The popular total variation (TV) [1] is designed as the total magnitude of the vertical and horizontal discrete gradients of an image and promotes this smoothness property on optimization. This method can be employed in more advanced image restoration problems, e.g., reflection removal [42,43], rain streak removal [44,45,46], intrinsic image decomposition [36,47,48], image blending [49], and multispectral pansharpening [50,51,52,53], which are complex problems to solve using filtering-based techniques; this is because optimization-based methods can explicitly and quantitatively design observation models by modeling using some norms and specific functions approximately. For example, we can flexibly define not only the commonly used ${\mathit{\ell }}_2$-norm but also the ${\mathit{\ell }}_1$-norm [54,55] and Huber loss function [49,56,57], which are robust to outliers, and so on, as a data-fidelity term in minimization problems.

Among these optimization-based methods, the ${\mathit{\ell }}_0$ gradient minimization proposed in [31] is known to provide the best approximation of the properties of the local piece-wise flatness of natural images, which is among the most essential prior knowledge in image processing. In this method, the ${\mathit{\ell }}_0$ gradient, which is the number of nonzero gradients of an image, is minimized together with the ${\mathit{\ell }}_2$ data fidelity to an input image. Recently, the applications of this method were actively studied, such as algorithm acceleration [58], deblurring [3,4], multi-dimensional data smoothing [58], and layer separation [48,59]. Ono investigated the ${\mathit{\ell }}_0$ gradient projection in [37], which allows us to control the degree of smoothing of the output image intuitively. Several extensions of the ${\mathit{\ell }}_0$ gradient minimization were recently proposed to flexibly separate the structure and texture components in smoothing texture-rich images [39,40,41].

However, the ${\mathit{\ell }}_0$ gradient-based smoothing methods have inherent issues, resulting in undesirable artifacts, such as pseudo-edges and halos in the output images. When an input image contains regions with gradations, the output image often exhibits a piecewise constant appearance, commonly referred to as the staircasing effect. These gradation regions produce strong edges in the output, which are not present in the input image. These issues arise because characteristics related to the intensity and sign of gradients within the observed images are not explicitly considered constraints in optimization problems. Furthermore, these artifacts become more noticeable, particularly in detail enhancement and tone-mapping applications. For example, although these applications are frequently used for diagnosing X-ray images, the staircasing effect introduces a potential risk of misdiagnosis. Similarly, tone mapping is often applied in low-light scenes, such as autonomous driving applications. However, this effect may lead to a reduction in object recognition accuracy. Our primary objective is to mitigate those unexpected artifacts caused by the staircasing effect by incorporating appropriate constraints into the gradient domain within optimization problems.

In this paper, we propose effective smoothing methods based on minimizing the ${\mathit{\ell }}_0$ gradient with novel gradient constraints that suppress those pseudo-edge artifacts. In particular, we introduce two types of gradient constraints: a box-type gradient constraint and a ball-type constraint. We incorporate them in the ${\mathit{\ell }}_0$ gradient minimization problem. These constraints are imposed using an appropriate reference image (e.g., an input image is used as a reference for image smoothing) to suppress pseudo-edges in homogeneous regions and the blurring effect of strong edges. We find a solution to the proposed minimization problems via the alternating direction method of multipliers (ADMM) [37,60,61] for nonconvex optimization.

The contributions of this paper are as follows.

• Image smoothing while maintaining gradient characteristics of reference image: Existing smoothing methods based on the ${\mathit{\ell }}_0$ gradient and TV do not explicitly consider any constraints in the gradient domain. In the proposed method, the local smoothness properties of a reference image can be explicitly considered constraints in the gradient domain. Therefore, we can suppress artifacts, including the staircasing effect, through image smoothing.
• Strict or flexible gradient constraints on the sign of gradients: Since the box-type gradient constraint is strict with respect to the sign of gradients, gradient reversals are well suppressed. In contrast, the ball-type constraint is flexible with respect to the sign of gradients, allowing robust image smoothing, even when a reference image is degraded by noise and has different shading characteristics, including gradient reversals.

The remainder of this paper is organized as follows. In Section 2, we present several mathematical preliminaries and the ADMM algorithm. In Section 3, we introduce the novel box- and ball-type gradient constraints and the minimization problems for image smoothing. In Section 4, several examples are shown to confirm the validity of the proposed methods compared with the existing smoothing methods based on the ${\mathit{\ell }}_0$ gradient. Finally, Section 5 concludes the paper.

In the previous paper [62], we introduced only the box-type gradient constraint to extend the ${\mathit{\ell }}_0$ gradient minimization. Additionally, robustness to noise and performance when using some reference images with different luminance and colors with/without reversed edges were not sufficiently evaluated. In this paper, we show the limitation of the box-type gradient constraint and propose a novel ball-type gradient constraint that can be applied flexibly, even when using references with noise or edge reversals. We also extend the proposed gradient constraints to the ${\mathit{\ell }}_0$ gradient projection [37] and the details of the algorithm derivation. Furthermore, we examine an experiment on JPEG artifact removal.

2. Preliminaries

Throughout this paper, bold-faced lowercase and uppercase letters indicate vectors and matrices, respectively, and ${\mathbb{R}}^N$ denotes N-dimensional vector space. We define the set of $M\times N$ real-valued matrices as ${\mathbb{R}}^{M\times N}$, and the block diagonal matrix of ${\mathbf{A}}_1,$$\dots ,{\mathbf{A}}_M$ by $\mathrm{diag}({\mathbf{A}}_1,\dots ,{\mathbf{A}}_M)$. The transpose operation of a vector and matrix is denoted as ${(·)}^{\top }$.

2.1. ${\mathit{\ell }}_0$ Gradient

Let $\mathbf{x}\in {\mathbb{R}}^N$ be an input vector. The ${\mathit{\ell }}_0$ pseudo-norm (note that this norm does not satisfy the properties of the norm; for simplicity in this paper, we refer to the ${\mathit{\ell }}_0$ pseudo-norm as the ${\mathit{\ell }}_0$-norm) counts the number of nonzero elements, and it is denoted by

$${\Vert \mathbf{x}\Vert }_0:=\sum _{n=1}^{N}C\left(x_n\right),C\left(x_n\right):=\left\{\begin{array}{cc}0,\hfill & \mathrm{if}{x}_n=0,\hfill \\ 1,\hfill & \mathrm{otherwise}.\hfill \end{array}\right.$$

(1)

Let $\mathbf{x}:={\left[{\mathbf{x}}_R^{\top }{\mathbf{x}}_G^{\top }{\mathbf{x}}_B^{\top }\right]}^{\top }(\in {\mathbb{R}}^{3N})$ be a vectorized color image, and let N be the number of pixels. Further, let ${\mathbf{D}}_v$ and ${\mathbf{D}}_h\in {\mathbb{R}}^{N\times N}$ be the vertical and horizontal first-order differential operators with a Neumann boundary, respectively; the differential operator is then defined by ${\mathbf{D}}_1:={\left[{\mathbf{D}}_v^{\top }{\mathbf{D}}_h^{\top }\right]}^{\top }(\in {\mathbb{R}}^{2N\times N})$ for a gray image and $\mathbf{D}\hspace{0.17em}:=\mathrm{diag}\left({\mathbf{D}}_1,{\mathbf{D}}_1,{\mathbf{D}}_1\right)(\in {\mathbb{R}}^{6N\times 3N})$ for a color image. In the ${\mathit{\ell }}_0$ gradient minimization, the group sparsity of the $RGB$ gradients is considered by concatenating the ${\mathit{\ell }}_1$-norm [31,37,48], and it is defined as

$${\Vert \mathbf{D}\mathbf{x}\Vert }_{0,1}:=\sum _{n=1}^{N}C\left(\left|{\left({\mathbf{D}}_v{\mathbf{x}}_R\right)}_n\right|+\left|{\left({\mathbf{D}}_h{\mathbf{x}}_R\right)}_n\right|+\left|{\left({\mathbf{D}}_v{\mathbf{x}}_G\right)}_n\right|+\left|{\left({\mathbf{D}}_h{\mathbf{x}}_G\right)}_n\right|+\left|{\left({\mathbf{D}}_v{\mathbf{x}}_B\right)}_n\right|+\left|{\left({\mathbf{D}}_h{\mathbf{x}}_B\right)}_n\right|\right).$$

(2)

A piece-wise smooth image is obtained by minimizing an optimization problem based on (2). In [31,48], a half quadratic optimization method is used to solve this optimization problem.

2.2. Alternating Direction Method of Multipliers

Alternating the direction method of multipliers (ADMM) [60] is a method of the proximal splitting algorithm that can treat convex optimization problems of the form

$$\underset{\mathbf{x}\in {\mathbb{R}}^{N_1},\mathbf{z}\in {\mathbb{R}}^{N_2}}{min}F\left(\mathbf{x}\right)+G\left(\mathbf{z}\right)\mathrm{s}.\mathrm{t}.\mathbf{z}=\mathbf{L}\mathbf{x},$$

(3)

where F and G are usually assumed to be quadratic and proximable functions, respectively, and $\mathbf{L}\in {\mathbb{R}}^{N_2\times N_1}$ is a matrix with full-column rank. For any ${\mathbf{x}}^{\left(0\right)}\in {\mathbb{R}}^{N_1},{\mathbf{z}}^{\left(0\right)}\in {\mathbb{R}}^{N_2},{\mathbf{b}}^{\left(0\right)}\in {\mathbb{R}}^{N_2}$ and $\gamma >0$, the ADMM algorithm is given by

$$⌊\begin{array}{cc}\hfill {\mathbf{x}}^{(\tau +1)}& =arg\underset{\mathbf{x}}{min}\left\{F\left(\mathbf{x}\right)+\frac{1}{2\gamma }{\Vert {\mathbf{z}}^{\left(\tau \right)}-\mathbf{L}\mathbf{x}-{\mathbf{b}}^{\left(\tau \right)}\Vert }_2^2\right\},\hfill \\ \hfill {\mathbf{z}}^{(\tau +1)}& =arg\underset{\mathbf{z}}{min}\left\{G\left(\mathbf{z}\right)+\frac{1}{2\gamma }{\Vert \mathbf{z}-\mathbf{L}{\mathbf{x}}^{(\tau +1)}-{\mathbf{b}}^{\left(\tau \right)}\Vert }_2^2\right\},\hfill \\ \hfill {\mathbf{b}}^{(\tau +1)}& ={\mathbf{b}}^{\left(\tau \right)}+\mathbf{L}{\mathbf{x}}^{(\tau +1)}-{\mathbf{z}}^{(\tau +1)}.\hfill \end{array}$$

(4)

When the optimization problem is convex, the sequence generated by (4) converges to an optimal solution of (3) under the existence of a saddle point.

2.3. Proximal Tools

The proximity operator [63] is a key tool of proximal splitting techniques. Let $\mathbf{x}\in {\mathbb{R}}^N$ be an input vector. For any $\gamma >0$, the primary operator of f over ${\mathbb{R}}^N$ is defined by

$${\mathrm{prox}}_{\gamma f}\left(\mathbf{x}\right):=arg\underset{\mathbf{y}\in {\mathbb{R}}^N}{min}f\left(\mathbf{y}\right)+\frac{1}{2\gamma }{\Vert \mathbf{x}-\mathbf{y}\Vert }^2.$$

(5)

For a given nonempty closed convex set $\mathcal{C}$, the indicator function of $\mathcal{C}$ is defined by

$${\iota }_{\mathcal{C}}\left(\mathbf{x}\right):=\left\{\begin{array}{cc}0,\hfill & \mathrm{if}\mathbf{x}\in \mathcal{C},\hfill \\ +\infty ,\hfill & \mathrm{otherwise}.\hfill \end{array}\right.$$

(6)

The proximity operator of ${\iota }_{\mathcal{C}}$ is expressed as

$${\mathrm{prox}}_{\gamma {\iota }_{\mathcal{C}}}\left(\mathbf{x}\right):=arg\underset{\mathbf{y}\in {\mathbb{R}}^N}{min}{\iota }_{\mathcal{C}}\left(\mathbf{x}\right)+\frac{1}{2\gamma }{\Vert \mathbf{x}-\mathbf{y}\Vert }_2^2.$$

(7)

The solution of ${\mathrm{prox}}_{\gamma {\iota }_{\mathcal{C}}}$ should be in the set $\mathcal{C}$ and minimize ${\Vert \mathbf{x}-\mathbf{y}\Vert }_2^2$. Thus, for any index $\gamma >0$, this proximity operator is equivalent to the metric projection onto $\mathcal{C}$, i.e., ${\mathcal{P}}_{\mathcal{C}}\left(\mathbf{x}\right)={\mathrm{prox}}_{\gamma {\iota }_{\mathcal{C}}}\left(\mathbf{x}\right)$.

Let $\mathbf{l}$ and $\mathbf{u}$$\in {\mathbb{R}}^N$ be the lower and upper bounds, respectively. The box constraint forces each element of $\mathbf{x}$ into the dynamic range $[l_i,u_i]$ for $i=1,\dots ,N$, and its closed convex set is defined as

$${\mathcal{C}}_{[\mathbf{l},\mathbf{u}]}:=\left\{\mathbf{x}\in {\mathbb{R}}^N\mid l_i\le x_i\le u_i(i=1,\dots ,N)\right\}.$$

(8)

The computation of the metric projection onto ${\mathcal{C}}_{[\mathbf{l},\mathbf{u}]}$ for $i=1,\dots ,N$ is given by

$${\left[{\mathcal{P}}_{{\mathcal{C}}_{[\mathbf{l},\mathbf{u}]}}\left(\mathbf{x}\right)\right]}_i=\left\{\begin{array}{cc}{l}_i,\hfill & \mathrm{if}{x}_i<l_i,\hfill \\ u_i,\hfill & \mathrm{if}{x}_i>u_i,\hfill \\ x_i,\hfill & \mathrm{if}{l}_i\le x_i\le u_i.\hfill \end{array}\right.$$

(9)

The ${\mathit{\ell }}_2$ ball constraint forces the Euclidean distance between the vector $\mathbf{x}$ and the centered vector $\mathbf{v}$ to be less than the radius $ϵ$, and its closed convex set is defined as

$${\mathcal{B}}_{\mathbf{v},ϵ}:=\left\{\mathbf{x}\in {\mathbb{R}}^N\mid {\Vert \mathbf{x}-\mathbf{v}\Vert }_2\le ϵ\right\}.$$

(10)

The computation of the metric projection onto ${\mathcal{B}}_{\mathbf{v},ϵ}$ is given by

$${\mathcal{P}}_{{\mathcal{B}}_{\mathbf{v},ϵ}}\left(\mathbf{x}\right)=\left\{\begin{array}{cc}\mathbf{x},\hfill & {\mathrm{if}\Vert \mathbf{x}-\mathbf{v}\Vert }_2\le ϵ,\hfill \\ \mathbf{v}+ϵ\frac{\mathbf{x}-\mathbf{v}}{{\Vert \mathbf{x}-\mathbf{v}\Vert }_2},\hfill & \mathrm{otherwise}.\hfill \end{array}\right.$$

(11)

3. Proposed Methods

To suppress the undesirable artifacts mentioned in Section 1, we introduce the following two types of gradient constraints.

3.1. Gradient Constraints

3.1.1. Box-Type Gradient Constraint

Let $\mathbf{x}:={\left[{\mathbf{x}}_R^{\top }{\mathbf{x}}_G^{\top }{\mathbf{x}}_B^{\top }\right]}^{\top }$$(\in {\mathbb{R}}^{3N})$ be a vectorized color image. We introduce the box-type gradient constraints in the positive and negative ranges:

$$0\le {\left(\mathbf{D}\mathbf{x}\right)}_i\le {\left(\mathbf{D}{\mathbf{x}}_{\mathrm{ref}}\right)}_i,i \in {\mathsf{\Omega }}_{+},$$

(12)

$${\left(\mathbf{D}{\mathbf{x}}_{\mathrm{ref}}\right)}_j\le {\left(\mathbf{D}\mathbf{x}\right)}_j<0,j \in {\mathsf{\Omega }}_{-},$$

(13)

where ${\mathbf{x}}_{\mathrm{ref}}\in {\mathbb{R}}^{3N}$ is a vectorized reference color image, $\mathbf{D}$ is the differential operator for a color image introduced in (2), and $\mathbf{D}\mathbf{x}$, $\mathbf{D}{\mathbf{x}}_{\mathrm{ref}}(\in {\mathbb{R}}^{6N})$ are gradient images. The notations ${\mathsf{\Omega }}_{+}$ and ${\mathsf{\Omega }}_{-}$ are the index sets of positive and negative elements in $\mathbf{D}{\mathbf{x}}_{\mathrm{ref}}$, respectively, where these sets satisfy ${\mathsf{\Omega }}_{+}\cup {\mathsf{\Omega }}_{-}=\{1,\dots ,6N\}$ and ${\mathsf{\Omega }}_{+}\cap {\mathsf{\Omega }}_{-}=\varnothing$. The gradient intensity of an output image obtained by solving a minimization problem with (12) and (13) is limited to the range of the gradient values of the reference image ${\mathbf{x}}_{\mathrm{ref}}$. Thus, these constraints suppress pseudo-edges that do not exist in the reference image.

The major drawback of minimizing the total variation regularization based on the ${\mathit{\ell }}_0$- or ${\mathit{\ell }}_1$-norms is that it yields pseudo-edge artifacts, such as an undesirable staircasing effect on the homogeneous regions in many cases. These artifacts are more severe in the case of the ${\mathit{\ell }}_0$-norm, as it has significant flattening effects. Figure 1 plots an example of the vertical first-order gradient signals of an input image and the smoothing results obtained by the two conventional methods—one is the ${\mathit{\ell }}_0$ gradient minimization [31], and the other is the ${\mathit{\ell }}_0$ gradient projection [37]—and the proposed methods with the box-type gradient constraint. The gradient signals of the ${\mathit{\ell }}_0$ gradient minimization and the ${\mathit{\ell }}_0$ gradient projection contain the gradients denoted by arrows that do not originally exist in the input image. Owing to the proposed gradient constraint, the gradient signal of the proposed methods is lower than that of the input image. When a reference image has local smoothness properties, artifacts such as halos and pseudo-edges can be suppressed by the proposed constraint. Similar results are obtained with the other channels and the horizontal direction.

Figure 2a illustrates an example of the box-type gradient constraint in the case where the vertical gradient is positive, and the horizontal gradient is negative in the reference image. We see from the figure that the ideal gradient is not included in the box-type gradient constraint. This often occurs when there is noise in the input image. Because the sign of gradients is sensitive to noise, our method generates better results for noise-free images; however, the method may generate unexpected artifacts for noisy images. To overcome this limitation, we introduce a ball-type gradient constraint.

3.1.2. Ball-Type Gradient Constraint

Let the sets of indexes ${\mathcal{G}}_1,\dots ,{\mathcal{G}}_N\in \mathcal{G}$ be ${\mathcal{G}}_n:=\{n,N+n,\dots ,5N+n\}$ for $n=1,\dots ,N$. We introduce the ball-type gradient constraint with respect to the n-th pixel as

$${\Vert \left(\mathbf{D}\mathbf{x}\right)}^{{\mathcal{G}}_n}{\Vert }_2 \le {\Vert {\left(\mathbf{D}{\mathbf{x}}_{\mathrm{ref}}\right)}^{{\mathcal{G}}_n}\Vert }_2,$$

(14)

where ${\left(\mathbf{D}\mathbf{x}\right)}^{{\mathcal{G}}_n}\in {\mathbb{R}}^6$ is the n-th sub-vector of $\mathbf{D}\mathbf{x}$ with the entries specified by ${\mathcal{G}}_n$, such as

$$\begin{array}{cc}\hfill {\left(\mathbf{D}\mathbf{x}\right)}^{{\mathcal{G}}_n}& :={\left[{\left(\mathbf{D}\mathbf{x}\right)}_n{\left(\mathbf{D}\mathbf{x}\right)}_{N+n}{\left(\mathbf{D}\mathbf{x}\right)}_{2N+n}{\left(\mathbf{D}\mathbf{x}\right)}_{3N+n}{\left(\mathbf{D}\mathbf{x}\right)}_{4N+n}{\left(\mathbf{D}\mathbf{x}\right)}_{5N+n}\right]}^{\top }\hfill \\ \hfill & ={\left[{\left({\mathbf{D}}_v{\mathbf{x}}_R\right)}_n{\left({\mathbf{D}}_h{\mathbf{x}}_R\right)}_n{\left({\mathbf{D}}_v{\mathbf{x}}_G\right)}_n{\left({\mathbf{D}}_h{\mathbf{x}}_G\right)}_n{\left({\mathbf{D}}_v{\mathbf{x}}_B\right)}_n{\left({\mathbf{D}}_h{\mathbf{x}}_B\right)}_n\right]}^{\top }.\hfill \end{array}$$

(15)

The ball-type gradient constraint limits the magnitude of the gradients of an output image to less than that of a reference image. Note that the signs of the gradients are not considered explicitly in these constraints.

Figure 2b illustrates an example of the ball-type gradient constraint with the same reference image as in Figure 2a. We see from this figure that the ball-type constraint is defined by a circle whose radius is the gradient magnitude of the reference image (see the dotted blue circle). Furthermore, the range of the ball-type constraint is more relaxed compared with the box-type one in Figure 2a, even with the same reference image. Therefore, even when we use a noisy image as a reference, in which the sign of the gradients varies locally, the ball-type constraint enables us to generate a smooth image that is locally smoother and more natural than the box-type one.

As described above, the ball-type constraint can obtain better results than the box-type one for a noisy image. However, when it is applied to a noise-free image, gradient reversals may occur, especially in the case of small-magnitude texture.

3.2. ${\mathit{\ell }}_0$-Smoothing Based on Box-Type Gradient Constraint

3.2.1. Minimization Problem

Let ${\mathbf{s}}_{\mathrm{in}}$ and ${\mathbf{s}}_{\mathrm{ref}}$ be a vectorized input color image and a vectorized reference image for the proposed box-type gradient constraint, respectively. The image-smoothing problem with the proposed gradient constraints (12) and (13) is defined by

$$\begin{array}{cc}\hfill & \underset{\mathbf{s}}{min}{\Vert \mathbf{D}\mathbf{s}\Vert }_{0,1}+\frac{\lambda }{2}{\Vert \mathbf{s}-{\mathbf{s}}_{\mathrm{in}}\Vert }_2^2\hfill \\ \hfill & \mathrm{s}.\mathrm{t}.0\le {\left(\mathbf{D}\mathbf{s}\right)}_i\le {\left(\mathbf{D}{\mathbf{s}}_{\mathrm{ref}}\right)}_i,i \in {\mathsf{\Omega }}_{+},{\left(\mathbf{D}{\mathbf{s}}_{\mathrm{ref}}\right)}_j\le {\left(\mathbf{D}\mathbf{s}\right)}_j<0,j \in {\mathsf{\Omega }}_{-},\hfill \end{array}$$

(16)

where $\lambda$ is a balancing weight of the two costs. The first term in (16) is defined by (2). The second term is a data-fidelity term calculating the squared error between the input and output images. Note that the objective function, excluding the gradient constraints from (16), is equivalent to that of the ${\mathit{\ell }}_0$ gradient minimization introduced in [31].

To find a solution of (16), we use ADMM described in Section 2.2. Since the first term of (16) is a nonconvex function, there is no guarantee that the solution obtained by ADMM is a global minimum solution of (16). However, as stated in [37,61], we experimentally confirmed that a reasonable solution satisfying the constraints can be obtained by carefully decreasing the step size of the algorithm in each iteration.

3.2.2. Optimization

Let ${\mathbf{g}}^{+}:=\mathbf{P}\mathbf{D}{\mathbf{s}}_{\mathrm{ref}}$, ${\mathbf{g}}^{-}:=\mathbf{M}\mathbf{D}{\mathbf{s}}_{\mathrm{ref}}$, where $\mathbf{P}$ and $\mathbf{M}$ are the matrices that extract the elements from $\mathbf{D}{\mathbf{s}}_{\mathrm{ref}}$ corresponding to the indexes included in ${\mathsf{\Omega }}_{+}$ and ${\mathsf{\Omega }}_{-}$, respectively. The convex sets ${\mathcal{C}}_{[\mathbf{0},{\mathbf{g}}^{+}]}$ and ${\mathcal{C}}_{[{\mathbf{g}}^{-},\mathbf{0}]}$ are defined as

$${\mathcal{C}}_{[\mathbf{0},{\mathbf{g}}^{+}]}:=\{\mathbf{x}\in {\mathbb{R}}^{m_{+}}\mid 0\le x_i\le g_i^{+},i=1,\dots ,m_{+}\},$$

$${\mathcal{C}}_{[{\mathbf{g}}^{-},\mathbf{0}]}:=\{\mathbf{x}\in {\mathbb{R}}^{m_{-}}\mid g_j^{-}\le x_j<0,j=1,\dots ,m_{-}\}.$$

where $m_{+}$ and $m_{-}$ denote the number of the indices in their sets, which satisfy $6N=m_{+}+m_{-}$, and $\mathbf{0}$ is a vector filled with zeros. To apply ADMM, we first reformulate (16) into the following unconstrained problem:

$$\underset{\mathbf{s}}{min}{\Vert \mathbf{D}\mathbf{s}\Vert }_{0,1}+\frac{\lambda }{2}{\Vert \mathbf{s}-{\mathbf{s}}_{\mathrm{in}}\Vert }_2^2+{\iota }_{{\mathcal{C}}_{[\mathbf{0},{\mathbf{g}}^{+}]}}\left(\mathbf{P}\mathbf{D}\mathbf{s}\right)+{\iota }_{{\mathcal{C}}_{[{\mathbf{g}}^{-},\mathbf{0}]}}\left(\mathbf{M}\mathbf{D}\mathbf{s}\right)$$

(17)

where ${\iota }_{{\mathcal{C}}_{[\mathbf{0},{\mathbf{g}}^{+}]}}(·)$ and ${\iota }_{{\mathcal{C}}_{[{\mathbf{g}}^{-},\mathbf{0}]}}(·)$ are the indicator functions (let $\mathbf{x}\in {\mathbb{R}}^N$ be an input vector; for a given non-empty closed convex set $\mathcal{C}$, the indicator function of $\mathcal{C}$ is defined by ${\iota }_C\left(\mathbf{x}\right)$, which returns 0 if $\mathbf{x}\in \mathcal{C}$, and $+\infty$ otherwise) of ${\mathcal{C}}_{[\mathbf{0},{\mathbf{g}}^{+}]}$ and ${\mathcal{C}}_{[{\mathbf{g}}^{-},\mathbf{0}]}$. These functions guarantee that the positive and negative gradient intensities of the optimal image ${\mathbf{s}}^{*}$ fall in the ranges $[\mathbf{0},{\mathbf{g}}^{+}]$ and $[{\mathbf{g}}^{-},\mathbf{0}]$. The role of the third and fourth terms of (17) correspond to the role of the first and second constraints of (16), respectively.

By introducing auxiliary variables ${\mathbf{z}}_1,{\mathbf{z}}_2$, and ${\mathbf{z}}_3$, we rewrite the minimization problem (17) into the following equivalent expression:

$$\begin{array}{cc}\hfill & \underset{\mathbf{s},{\mathbf{z}}_1,{\mathbf{z}}_2,{\mathbf{z}}_3}{min}F\left(\mathbf{s}\right)+G\left(\mathbf{z}\right)\mathrm{s}.\mathrm{t}.\mathbf{z}=\mathbf{L}\mathbf{s},\hfill \\ \hfill & F\left(\mathbf{s}\right):=\frac{\lambda }{2}\Vert \mathbf{s}-{\mathbf{s}}_{\mathrm{in}}{\Vert }_2^2,G\left(\mathbf{z}\right):={\Vert {\mathbf{z}}_1\Vert }_{0,1}+{\iota }_{{\mathcal{C}}_{[\mathbf{0},{\mathbf{g}}^{+}]}}\left({\mathbf{z}}_2\right)+{\iota }_{{\mathcal{C}}_{[{\mathbf{g}}^{-},\mathbf{0}]}}\left({\mathbf{z}}_3\right),\hfill \\ \hfill & \mathbf{z}={\left[{\mathbf{z}}_1^{\top }{\mathbf{z}}_2^{\top }{\mathbf{z}}_3^{\top }\right]}^{\top }({\mathbf{z}}_1\in {\mathbb{R}}^{6N},{\mathbf{z}}_2\in {\mathbb{R}}^{m_{+}},{\mathbf{z}}_3\in {\mathbb{R}}^{m_{-}}),\hfill \\ \hfill & \mathbf{L}=\left[\begin{array}{c}\mathbf{D}\\ \mathbf{P}\mathbf{D}\\ \mathbf{M}\mathbf{D}\end{array}\right]\left(\in {\mathbb{R}}^{12N\times 3N}\right).\hfill \end{array}$$

(18)

The algorithm for solving Equation (18) with ${\gamma }_i(i=1,2,3)$ is summarized in Algorithm 1.

Algorithm 1 Proposed algorithm for (18).

1: Input    : ${\mathbf{z}}_i^{\left(0\right)}$, ${\mathbf{b}}_i^{\left(0\right)},{\gamma }_i$$(i=1,2,3)$, $\lambda$, $\mu$, $\eta$$(0<\eta <1)$ 2: Output: ${\mathbf{s}}^{\left(\tau \right)}$ 3: while A stopping criterion is not satisfied do 4:    ${\mathbf{s}}^{\left(\tau +1\right)}\leftarrow arg{min}_{\mathbf{s}}\frac{\lambda }{2}\Vert \mathbf{s}-{\mathbf{s}}_{\mathrm{in}}{\Vert }_2^2+\frac{1}{2{\gamma }_1}\Vert {\mathbf{z}}_1^{\left(\tau \right)}-\mathbf{D}\mathbf{s}-{\mathbf{b}}_1^{\left(\tau \right)}{\Vert }_2^2+\frac{1}{2{\gamma }_2}\Vert {\mathbf{z}}_2^{\left(\tau \right)}-\mathbf{P}\mathbf{D}\mathbf{s}-{\mathbf{b}}_2^{\left(\tau \right)}{\Vert }_2^2+\frac{1}{2{\gamma }_3}{\Vert {\mathbf{z}}_3^{\left(\tau \right)}-\mathbf{M}\mathbf{D}\mathbf{s}-{\mathbf{b}}_3^{\left(\tau \right)}\Vert }_2^2$; 5:    for $n=1$ to N do 6:        ${\mathbf{z}}_1^{{\mathcal{G}}_n(\tau +1)}\leftarrow {\mathrm{prox}}_{{\gamma }_1{\Vert ·\Vert }_{0,1}}\left({(\mathbf{D}{\mathbf{s}}^{(\tau +1)}+{\mathbf{b}}_1^{\left(\tau \right)})}^{{\mathcal{G}}_n}\right)$; 7:    end for 8:    ${\mathbf{z}}_2^{(\tau +1)}\leftarrow {\mathcal{P}}_{[\mathbf{0},{\mathbf{g}}^{+}]}\left(\mathbf{P}\mathbf{D}{\mathbf{s}}^{(\tau +1)}+{\mathbf{b}}_2^{\left(\tau \right)}\right)$; 9:    ${\mathbf{z}}_3^{(\tau +1)}\leftarrow {\mathcal{P}}_{[{\mathbf{g}}^{-},\mathbf{0}]}\left(\mathbf{M}\mathbf{D}{\mathbf{s}}^{(\tau +1)}+{\mathbf{b}}_3^{\left(\tau \right)}\right)$; 10:    ${\mathbf{b}}_1^{(\tau +1)}\leftarrow {\mathbf{b}}_1^{\left(\tau \right)}+\mathbf{D}{\mathbf{s}}^{(\tau +1)}-{\mathbf{z}}_1^{(\tau +1)}$; 11:    ${\mathbf{b}}_2^{(\tau +1)}\leftarrow {\mathbf{b}}_2^{\left(\tau \right)}+\mathbf{P}\mathbf{D}{\mathbf{s}}^{(\tau +1)}-{\mathbf{z}}_2^{(\tau +1)}$; 12:    ${\mathbf{b}}_3^{(\tau +1)}\leftarrow {\mathbf{b}}_3^{\left(\tau \right)}+\mathbf{M}\mathbf{D}{\mathbf{s}}^{(\tau +1)}-{\mathbf{z}}_3^{(\tau +1)}$; 13:    ${\gamma }_1\leftarrow \eta {\gamma }_1$; 14:    ${\gamma }_2\leftarrow \eta {\gamma }_2$; 15:    ${\gamma }_3\leftarrow \eta {\gamma }_3$; 16:    $\tau \leftarrow \tau +1$; 17: end while

The update of $\mathbf{s}$ in Algorithm 1 is achieved by solving the following quadratic minimization problem (hereafter, the superscript ($\tau$) is omitted for simplicity):

$$\underset{\mathbf{s}}{min}\frac{\lambda }{2}\Vert \mathbf{s}-{\mathbf{s}}_{\mathrm{in}}{\Vert }_2^2+\frac{1}{2{\gamma }_1}\Vert {\mathbf{z}}_1-\mathbf{D}\mathbf{s}-{\mathbf{b}}_1{\Vert }_2^2+\frac{1}{2{\gamma }_2}\Vert {\mathbf{z}}_2-\mathbf{P}\mathbf{D}\mathbf{s}-{\mathbf{b}}_2{\Vert }_2^2+\frac{1}{2{\gamma }_3}{\Vert {\mathbf{z}}_3-\mathbf{M}\mathbf{D}\mathbf{s}-{\mathbf{b}}_3\Vert }_2^2.$$

(19)

Thus, by setting the first-order derivative to zero, the optimal solution is characterized by the system of linear equations:

$$\begin{array}{cc}\hfill & \mathbf{A}\mathbf{s}=\mathbf{q},\hfill \\ \hfill & \left\{\begin{array}{c}\mathbf{A}=\lambda \mathbf{I}+{\gamma }_1^{-1}{\mathbf{D}}^{\top }\mathbf{D}+{\gamma }_2^{-1}{\mathbf{D}}^{\top }{\mathbf{P}}^{\top }\mathbf{P}\mathbf{D}+{\gamma }_3^{-1}{\mathbf{D}}^{\top }{\mathbf{M}}^{\top }\mathbf{M}\mathbf{D},\hfill \\ \mathbf{q}=\lambda {\mathbf{s}}_{\mathrm{in}}+{\gamma }_1^{-1}{\mathbf{D}}^{\top }\left({\mathbf{z}}_1-{\mathbf{b}}_1\right)+{\gamma }_2^{-1}{\mathbf{D}}^{\top }{\mathbf{P}}^{\top }\left({\mathbf{z}}_2-{\mathbf{b}}_2\right)+{\gamma }_3^{-1}{\mathbf{D}}^{\top }{\mathbf{M}}^{\top }\left({\mathbf{z}}_3-{\mathbf{b}}_3\right),\hfill \end{array}\right.\hfill \end{array}$$

(20)

where $\mathbf{A}$ is a matrix of size $\left\{3N\times 3N\right\}$, $\mathbf{q}$ is a $3N$-dimensional vector, and $\mathbf{I}\in {\mathbb{R}}^{3N\times 3N}$ is an identity matrix. The optimal solution is obtained by the inverse problem ${\mathbf{s}}^{*}={\mathbf{A}}^{-1}\mathbf{q}$. If we set ${\gamma }_3={\gamma }_2$, the third term of the right side in the formula $\mathbf{A}$ can be rewritten as ${\gamma }_2^{-1}{\mathbf{D}}^{\top }\mathbf{D}={\gamma }_2^{-1}{\mathbf{D}}^{\top }{\mathbf{P}}^{\top }\mathbf{P}\mathbf{D}+{\gamma }_3^{-1}{\mathbf{D}}^{\top }{\mathbf{M}}^{\top }\mathbf{M}\mathbf{D}$. Since the boundary condition of $\mathbf{D}$ is periodic and the matrix $\mathbf{A}$ is a block circulant matrix with circulant blocks (BCCB), it is diagonalized by the 2D fast discrete Fourier transform (2DFFT). Thus, the inverse problem of (20) can be calculated as

$${\mathbf{s}}^{*}={\mathbf{F}}^{*}{\left\{\lambda \mathbf{I}+({\gamma }_1^{-1}+{\gamma }_2^{-1})\mathsf{\Sigma }\right\}}^{-1}\mathbf{F}\mathbf{q},$$

(21)

where $\mathbf{F}$ and ${\mathbf{F}}^{*}$ are the 2D-FFT and its inverse matrices, respectively, and $\mathsf{\Sigma }$ is the diagonal matrix with its entries being the Fourier-transformed Laplacian filter kernel. Therefore, its inversion is reduced to entry-wise division.

For the update of ${\mathbf{z}}_1$, we need to compute the pseudo-proximity operator of the ${\mathit{\ell }}_{0,1}$-norm. By dividing ${\mathbf{z}}_1$ to sub-vectors as ${\mathbf{z}}_1:={[{\mathbf{z}}_1^{{\mathcal{G}}_1\top }\dots {\mathbf{z}}_1^{{\mathcal{G}}_N\top }]}^{\top }$, in which ${\mathcal{G}}_n(n=1,\dots ,N)$ is introduced in Section 3.1.2 and ${\mathbf{z}}_1^{{\mathcal{G}}_n}$ can be regarded as the auxiliary variable for ${\left(\mathbf{D}\mathbf{s}\right)}^{{\mathcal{G}}_n}$ (see (15)), the pseudo-proximity operator equals to a group hard-thresholding operation for each sub-vector [48,61]. Thus, the optimal solution of the n-th sub-vector ${\mathbf{z}}_1^{{\mathcal{G}}_n}$ is obtained by

$${\mathbf{z}}_1^{{\mathcal{G}}_n*}\leftarrow {\mathrm{prox}}_{{\gamma }_1{\Vert ·\Vert }_{0,1}}\left({\mathbf{d}}_n\right)=\left\{\begin{array}{cc}{\mathbf{d}}_n,\hfill & \mathrm{if}\sqrt{{\sum }_{m=1}^6{d}_{n,m}^2}\ge \sqrt{2{\gamma }_1},\hfill \\ \mathbf{0},\hfill & \mathrm{otherwise},\hfill \end{array}\right.$$

(22)

where $d_{n,m}$ is the m-th element of ${\mathbf{d}}_n:={(\mathbf{D}\mathbf{s}+{\mathbf{b}}_1)}^{{\mathcal{G}}_n}$$(\in {\mathbb{R}}^6)$.

To update ${\mathbf{z}}_2$ and ${\mathbf{z}}_3$, we need to compute the proximity operator of the indicator functions ${\iota }_{{\mathcal{C}}_{[\mathbf{0},{\mathbf{g}}^{+}]}}$ and ${\iota }_{{\mathcal{C}}_{[{\mathbf{g}}^{-},\mathbf{0}]}}$. By (9), these are simply projections onto each set, respectively, i.e., for $i=1,\dots ,m_{+}$ by

$${\mathbf{z}}_2\leftarrow {\left[{\mathcal{P}}_{{\mathcal{C}}_{[\mathbf{0},{\mathbf{g}}^{+}]}}(\mathbf{P}\mathbf{D}\mathbf{s}+{\mathbf{b}}_2)\right]}_i=\left\{\begin{array}{cc}0,\hfill & \mathrm{if}{(\mathbf{P}\mathbf{D}\mathbf{s}+{\mathbf{b}}_2)}_i<0,\hfill \\ g_i^{+},\hfill & \mathrm{if}{(\mathbf{P}\mathbf{D}\mathbf{s}+{\mathbf{b}}_2)}_i>g_i^{+},\hfill \\ {(\mathbf{P}\mathbf{D}\mathbf{s}+{\mathbf{b}}_2)}_i,\hfill & \mathrm{if}0\le {(\mathbf{P}\mathbf{D}\mathbf{s}+{\mathbf{b}}_2)}_i\le g_i^{+},\hfill \end{array}\right.$$

(23)

and for $j=1,\dots ,m_{-}$ by

$${\mathbf{z}}_3\leftarrow {\left[{\mathcal{P}}_{{\mathcal{C}}_{[{\mathbf{g}}^{-},\mathbf{0}]}}(\mathbf{M}\mathbf{D}\mathbf{s}+{\mathbf{b}}_3)\right]}_j=\left\{\begin{array}{cc}{g}_j^{-},\hfill & \mathrm{if}{(\mathbf{M}\mathbf{D}\mathbf{s}+{\mathbf{b}}_3)}_j<g_j^{-},\hfill \\ 0,\hfill & \mathrm{if}{(\mathbf{M}\mathbf{D}\mathbf{s}+{\mathbf{b}}_3)}_j>0,\hfill \\ {(\mathbf{M}\mathbf{D}\mathbf{s}+{\mathbf{b}}_3)}_j,\hfill & \mathrm{if}{g}_j^{-}\le {(\mathbf{M}\mathbf{D}\mathbf{s}+{\mathbf{b}}_3)}_j\le 0.\hfill \end{array}\right.$$

(24)

To stabilize ADMM for nonconvex optimization, a scalar $0<\eta <1$ is introduced to gradually decrease the value of ${\gamma }_i(i=1,2,3)$ as in the 11th–13th lines of Algorithm 1. The algorithm continues until the following stopping criterion is satisfied:

$$\frac{\Vert {\mathbf{s}}^{(\tau +1)}-{\mathbf{s}}^{\left(\tau \right)}{\Vert }_2}{\Vert {\mathbf{s}}^{(\tau +1)}{\Vert }_2}\le \mu ,$$

(25)

where we set a small value to $\mu$ (e.g., $\mu \le {10}^{-4}$). Although there is no theoretical guarantee of convergence for our algorithm, we show that the sequence generated by this algorithm experimentally converges to a reasonable solution with sufficiently small values for ${\gamma }_i(i=1,2,3)$; similar strategies have been employed in the existing minimization algorithms for the ${\mathit{\ell }}_0$-norm [37,48,61].

3.3. ${\mathit{\ell }}_0$-Smoothing Based on Ball-Type Gradient Constraint

The image-smoothing problem with the proposed ball-type gradient constraint (14), which is the modified version of (16), is defined by

$$\underset{\mathbf{s}}{min}{\Vert \mathbf{D}\mathbf{s}\Vert }_{0,1}+\frac{\lambda }{2}\Vert \mathbf{s}-{\mathbf{s}}_{\mathrm{in}}{\Vert }_2^2{\mathrm{s}.\mathrm{t}.\Vert \left(\mathbf{D}\mathbf{s}\right)}^{{\mathcal{G}}^n}{\Vert }_2\le {\Vert {\left(\mathbf{D}{\mathbf{s}}_{\mathrm{ref}}\right)}^{{\mathcal{G}}_n}\Vert }_2,n=1,\dots ,N.$$

(26)

Similar to the preceding section, we solve the minimization problem (26) by applying ADMM to it. Now, the convex set ${\mathcal{B}}_{\mathbf{0},{ϵ}_n}$ is defined for $n=1,\dots ,N$ by

$${\mathcal{B}}_{\mathbf{0},{ϵ}_n}:=\{\mathbf{x}\in {\mathbb{R}}^6\mid {\Vert \mathbf{x}\Vert }_2\le {ϵ}_n\},$$

(27)

where ${ϵ}_n:={\Vert {\left(\mathbf{D}{\mathbf{s}}_{\mathrm{ref}}\right)}^{{\mathcal{G}}_n}\Vert }_2$. We reformulate (26) into the following unconstrained problem by introducing the indicator functions of ${\mathcal{B}}_{\mathbf{0},{ϵ}_n}(n=1,\dots ,N)$, and then, we also introduce auxiliary variables ${\mathbf{z}}_1$ and ${\mathbf{z}}_2$ to non-differential functions

$$\underset{\mathbf{s},{\mathbf{z}}_1,{\mathbf{z}}_2}{min}\frac{\lambda }{2}\Vert \mathbf{s}-{\mathbf{s}}_{\mathrm{in}}{\Vert }_2^2+{\Vert {\mathbf{z}}_1\Vert }_{0,1}+\sum _{n=1}^N{\iota }_{{\mathcal{B}}_{\mathbf{0},{ϵ}_n}}({\mathbf{z}}_2^{{\mathcal{G}}_n})\mathrm{s}.\mathrm{t}.{\mathbf{z}}_1=\mathbf{D}\mathbf{s},{\mathbf{z}}_2=\mathbf{D}\mathbf{s},$$

(28)

where ${\mathbf{z}}_2^{{\mathcal{G}}_n}$ is the n-th sub-vector of ${\mathbf{z}}_2$ specified by ${\mathcal{G}}_n$.

The algorithm for solving (28) is obtained by slightly modifying Algorithm 1. In particular, the update of $\mathbf{s}$ in line 4 of Algorithm 1 is obtained by

$$\begin{array}{cc}\hfill & {\mathbf{s}}^{*}={\mathbf{F}}^{*}{\left\{\lambda \mathbf{I}+({\gamma }_1^{-1}+{\gamma }_2^{-1})\mathsf{\Sigma }\right\}}^{-1}\mathbf{F}{\mathbf{q}}^{\prime },\hfill \\ \hfill & \mathrm{where}{\mathbf{q}}^{\prime }=\lambda {\mathbf{s}}_{\mathrm{in}}+{\gamma }_1^{-1}{\mathbf{D}}^{\top }({\mathbf{z}}_1-{\mathbf{b}}_1)+{\gamma }_2^{-1}{\mathbf{D}}^{\top }({\mathbf{z}}_2-{\mathbf{b}}_2).\hfill \end{array}$$

(29)

In the update of ${\mathbf{z}}_2$, we calculate the following metric projection onto the ${\mathit{\ell }}_2$ ball according to (11) instead of the 8th and 9th lines of Algorithm 1, i.e., for $n=1,\dots ,N$, by

$${\mathbf{z}}_2^{{\mathcal{G}}_n}\leftarrow {\mathcal{P}}_{{\mathcal{B}}_{\mathbf{0},{ϵ}_n}}\left({(\mathbf{D}\mathbf{s}-{\mathbf{b}}_2)}^{{\mathcal{G}}_n}\right)=\left\{\begin{array}{cc}{\left(\mathbf{D}\mathbf{s}-{\mathbf{b}}_2\right)}^{{\mathcal{G}}_n},\hfill & \mathrm{if}\Vert {(\mathbf{D}\mathbf{s}-{\mathbf{b}}_2)}^{{\mathcal{G}}_n}{\Vert }_2\le {ϵ}_n,\hfill \\ {ϵ}_n\frac{{(\mathbf{D}\mathbf{s}-{\mathbf{b}}_2)}^{{\mathcal{G}}_n}}{\Vert {(\mathbf{D}\mathbf{s}-{\mathbf{b}}_2)}^{{\mathcal{G}}_n}{\Vert }_2},\hfill & \mathrm{otherwise}.\hfill \end{array}\right.$$

(30)

3.4. ${\mathit{\ell }}_0$ Gradient Projection with Gradient Constraint

Our proposed box-/ball-type gradient constraints can be incorporated with not only the minimization problem of the ${\mathit{\ell }}_0$ gradient minimization [31] but also with that of the ${\mathit{\ell }}_0$ gradient projection [37]. An image-smoothing problem based on the ${\mathit{\ell }}_0$ gradient projection (for simplicity, we only discuss using the ball-type gradient constraint; note that the box-type gradient constraints can be also used as well) with the proposed gradient constraint is defined by

$$\underset{\mathbf{s}}{min}\frac{1}{2}\Vert \mathbf{s}-{\mathbf{s}}_{\mathrm{in}}{\Vert }_2^2{\mathrm{s}.\mathrm{t}.\Vert \mathbf{D}\mathbf{s}\Vert }_{0,1}{ \le \alpha ,\Vert \left(\mathbf{D}\mathbf{s}\right)}^{{\mathcal{G}}_n}{\Vert }_2 \le {\Vert {\left(\mathbf{D}{\mathbf{s}}_{\mathrm{ref}}\right)}^{{\mathcal{G}}_n}\Vert }_2,n=1,\dots ,N.$$

(31)

where $\alpha$ is the parameter specified by a user that is the least upper bound of ${\Vert \mathbf{D}\mathbf{s}\Vert }_{0,1}$. It allows us to intuitively control the degree of smoothing of the output image ${\mathbf{s}}^{*}$.

Similar to the preceding section, a solution of (31) can be estimated based on Algorithm 1 after transforming it into an equation applicable to ADMM (the indicator function of the inequality constraint on ${\Vert \mathbf{D}\mathbf{s}\Vert }_{0,1}$ is defined in (9) of [37], and its pseudo-metric projection is derived in (16) of [37]).

4. Experiments

To clarify the differences in the characteristics of the proposed methods with the box- and ball-type gradient constraints, we first conducted noise removal experiments using some reference images. We then applied the proposed methods to various applications, such as image smoothing, detail enhancement, tone mapping, and JPEG artifact removal, and the results were compared with those of two existing methods based on the ${\mathit{\ell }}_0$ gradient: ${\mathit{\ell }}_0$ gradient minimization [31], and ${\mathit{\ell }}_0$ gradient projection [37] (for the implementation of the existing methods, we used the source code distributed by the authors of each paper).

All experiments were performed using MATLAB R2021a on a computer with an AMD EPYC 7402P 2.80 GHz processor and 128 GB RAM. To accelerate the computation time of the proposed methods, we used an NVIDIA GeForce RTX 3090 GPU. For the parameters of our methods, we set ${\gamma }_i=5(i=1,2,3)$, $\eta =0.97$, and $\mu ={10}^{-4}$ in all experiments. For the reference image ${\mathbf{s}}_{\mathrm{ref}}$, we set ${\mathbf{s}}_{\mathrm{in}}$ in the three experiments described in Section 4.2, Section 4.3 and Section 4.4, respectively. Note that the dynamic range of input images was normalized in $[0,1]$.

4.1. Box-Type vs. Ball-Type Gradient Constraint

The performance of our methods depend on the variation of the luminance and color in a reference image. Accordingly, we studied the sensitivity of our methods with respect to a reference image by conducting an experiment on Gaussian noise removal.

We generated a noisy image, which is shown in Figure 3b, by adding zero-mean Gaussian noise with standard deviation $\sigma =5.0\times {10}^{-3}$ to the ground truth image shown in Figure 3a. Then, four types of images, shown in Figure 4(a-1) ground truth image, a-2 horizontally flipped image of a-1, a-3 hue shift image of a-1, and a-4 horizontally flipped and hue shift image of a-1, were prepared, each of which was used as a reference image of our proposed gradient constraints. For the quality metric, we used the peak signal-to-noise ratio (PSNR). For the smoothing parameter setting, the balancing weight $\lambda$ was set to maximize PSNR. Note that the purpose of this experiment is not to evaluate the pure performance of noise removal but to clarify the sensitivity of the method.

Figure 4 shows the noise removal results and their PSNR values. One can see from Figure 4b,c that our images with the ball-type constraint have higher PSNR and visually better results than the box-type image for all cases. The results of Figure 4(b-2,b-4) are over-smoothed because the reference images have luminance variations that are considerably different from those of the input image (there are flipped relationships between the luminance variations of the input and reference images in the horizontal direction). This is because the sign of the gradients of the reference images is often reversed from the ideal one around the edges. The box-type constraint does not prevent the smoothed image from having the gradients whose sign is reversed from that of the reference image, thus yielding blurred images. In Figure 4(b-1,b-3), this constraint generates a better smooth result with sharper edges similar to those of the ground truth image; however, pseudo-color artifacts often occur around edges as shown in Figure 4(b-3). Thus, this constraint is very sensitive to the sign of image gradients unless an input image is used as a reference. In contrast, the ball-type constraint robustly generates a better smooth result with sharper edges that is the same as that in the input image, even if the reference images have luminance and color variations considerably different from those in the input image.

4.2. Detail Enhancement

Based on the methods [5,6], we performed the detail enhancement as follows: (i) obtain a base layer by applying each method to an input image; (ii) calculate a detail layer by subtracting the base layer from the input image; (iii) amplify the detail layer by scaling factor $s>1$; and (iv) reconstruct a detail-enhanced image by summing the base layer and the enhanced detail layer.

Figure 5 and Figure 6 show the smoothing and detail enhancement results obtained by the existing methods and the proposed methods (we took input images from an extensive database of royalty-free images https://pixabay.com/en/ accessed on 8 January 2018). In Section 4.2 and Section 4.3, since there are visually no differences between the results of the proposed methods with the box- and ball-type gradient constraints, we only show the results of the proposed method with the box-type gradient constraint. For a fair comparison, we adjusted the smoothing parameters of each method so that the ${\mathit{\ell }}_2$-norms are the same in the input and output images. The scaling factor s was set to $2.5$. Figure 5 and Figure 6 indicate that the contrast in the entire image is maintained for all methods. However, the results of [31] have halo-like artifacts in the region with gradation indicated by the red square window in (b), and the color of the clouds is changed. Pseudo-edge artifacts, such as the staircasing effect, occur in the sky region of the results of [37] (indicated by the yellow and red square windows in (c)). In contrast, our methods do not have these undesirable artifacts because the gradient constraints suppress pseudo-edges in the homogeneous regions.

4.3. Tone Mapping

Next, we applied these methods for layer separation in the tone-mapping framework [7].

Figure 7 shows the tone-mapping results, where we also show the result of Reinhard et al.’s global tone-mapping operator [64] for reference. Figure 7b,c show that the halo and pseudo-edge artifacts occurred, especially at the boundary between the bright and dark regions. In contrast, the proposed methods sufficiently suppressed these artifacts.

To show the robustness of the degree of smoothing, we applied the ${\mathit{\ell }}_0$ gradient minimization [31] and the proposed methods to an input HDR image with different smoothing parameters. The smoothing parameters were adjusted so that the mean squared error (MSE) of the smoothing results was $1.2\times {10}^{-5},1.9\times {10}^{-5},$ and $2.5\times {10}^{-5}$, and these results are shown in Figure 8. One sees from Figure 8(b-1,c-1,d-1) that the staircasing effects and pseudo-color edges occur. In addition, these artifacts are more noticeable as the smoothing degree increases. In contrast, the results of the proposed methods have few artifacts. Even in Figure 8(d-2), which is the smoothest result, noticeable artifacts do not appear. In Section 4.2 and Section 4.3, we applied the proposed methods to dozens of images and obtained similar results.

4.4. JPEG Artifact Removal in Clip-Art Images

Input images were generated from clip-art images (we took ten test images from the dataset provided in https://google.github.io/cartoonset/ accessed on 20 July 2019), which were compressed by standard JPEG with quality values in the range of $[10,90]$. The resultant images were evaluated by using two quality metrics: PSNR and structural similarity (SSIM) [65]. In each method, we adjusted the degree of noise removal to obtain the visually best restoration results, i.e., maximizing the smoothness while preserving the edges of the images as much as possible (we only showed the results of the proposed methods obtained with the ball-type gradient constraint because it is more robust than the box-type constraint for noise removal).

Figure 9a,b plot the PSNR [dB] and SSIM values averaged over ten test images, respectively. It can be seen that the PSNR and SSIM of the proposed methods are higher than those of the existing methods in most cases.

Figure 10 shows the results and their PSNR values [dB] of two images in the cases of the quality value with 20. Block-wise artifacts in the JPEG images shown in Figure 10(b-1,b-2) are sufficiently removed by the proposed methods, while the existing methods retain them (see close-up images indicated by the red and blue boxes). Our methods have the highest PSNR in both images.

5. Conclusions

We proposed effective smoothing methods based on minimizing the ${\mathit{\ell }}_0$ gradient using novel gradient constraints. To suppress the pseudo-edge artifacts such as the staircasing effect, box- and ball-type gradient constraints were introduced. These constraints limited the gradient intensity of an output image to the range of a reference image. By solving the ${\mathit{\ell }}_0$ gradient minimization problem based on the proposed gradient constraints via ADMM for nonconvex optimization, we effectively preserved the strong edges and removed weak edges without yielding pseudo-edges and halo artifacts. It was confirmed that the smoothed image obtained by the proposed method can generate better contrast-enhanced images and tone map images without pseudo-edges or halos. These results suggest the potential for delivering more dependable image diagnosis and recognition outcomes in medical imaging and autonomous driving applications, where image enhancement applications are essential. Throughout all our experiments, we consistently observed the effectiveness of the proposed methods in comparison to existing ${\mathit{\ell }}_0$ gradient-based methods.

In future work, we will apply our proposed gradient constraints for other optimization problems, e.g., intrinsic image decomposition, reflection removal, and flash/no-flash image blending.