Embedding Bilateral Filter in Least Squares for Efficient Edge-preserving Image Smoothing

TL;DR

This paper introduces a new global edge-preserving smoothing method that combines bilateral filtering with least squares, achieving high performance with significantly improved efficiency over existing global methods.

Contribution

The paper presents a novel approach embedding bilateral filter into the least squares model, balancing the accuracy of global methods with the speed of local methods.

Findings

Comparable performance to state-of-the-art global methods

Significantly faster processing times

Flexible framework adaptable with bilateral filter variants

Abstract

Edge-preserving smoothing is a fundamental procedure for many computer vision and graphic applications. This can be achieved with either local methods or global methods. In most cases, global methods can yield superior performance over local ones. However, local methods usually run much faster than global ones. In this paper, we propose a new global method that embeds the bilateral filter in the least squares model for efficient edge-preserving smoothing. The proposed method can show comparable performance with the state-of-the-art global method. Meanwhile, since the proposed method can take advantages of the efficiency of the bilateral filter and least squares model, it runs much faster. In addition, we show the flexibility of our method which can be easily extended by replacing the bilateral filter with its variants. They can be further modified to handle more applications. We validate the effectiveness and efficiency of the proposed method through comprehensive experiments in a range of applications.

Tables

Table 1. TABLE I: Property comparison summary of different methods. The running time is measured on processing a 1024 × 1024 1024 1024 1024\times 1024 RGB color image.

	GF [20]	$L_0$ norm [25]	WLS [2]	BLF [1]	AMF [17]	NC filter [16]	BLF-LS	AMF-LS	NC-LS
runtime efficiency	0.32s	4.76s	6.83s	0.16s	0.2s	0.13s	0.48s	0.56s	0.42s
halos	$✓$	$\times$	$\times$	$✓$	$✓$	$✓$	$\times$	$\times$	$\times$
gradient reversals	$\times$	$✓$	$\times$	$✓$	$✓$	$✓$	$\times$	$\times$	$\times$
local/global	local	global	global	local	local	local	global	global	global

Equations

$$\begin{array}{rl}& u_s=\frac{1}{Z_s}\sum _{t\in N(s)}{G}_{{\sigma }_s}(\Vert s-t\Vert )G_{{\sigma }_r}(\Vert g_s-g_t\Vert )g_t,\\ & Z_s=\sum _{t\in N(s)}{G}_{{\sigma }_s}(\Vert s-t\Vert )G_{{\sigma }_r}(\Vert g_s-g_t\Vert ),\end{array}$$

$$\underset{u}{\min }\sum _s\left((u_s-g_s{)}^2+\lambda \left({\omega }_{x,s}(\nabla u_{x,s}{)}^2+{\omega }_{y,s}(\nabla u_{y,s}{)}^2\right)\right),$$

$${\omega }_{x,s}={\left(|\nabla {\ell }_{x,s}{|}^{\alpha }+\epsilon \right)}^{-1},{\omega }_{y,s}={\left(|\nabla {\ell }_{y,s}{|}^{\alpha }+\epsilon \right)}^{-1},$$

$$(I+\lambda L_g)u=g\to u=(I+\lambda L_g{)}^{-1}g,$$

$$\underset{u}{\min }\sum _s\left((u_s-g_s{)}^2+\lambda \left((\nabla u_{x,s}{)}^2+(\nabla u_{y,s}{)}^2\right)\right).$$

$$u={\mathcal{F}}^{-1}\left(\frac{\mathcal{F}(g)}{\mathcal{F}(1)+\lambda (\overline{\mathcal{F}({\partial }_x)}\cdot \mathcal{F}({\partial }_x)+\overline{\mathcal{F}({\partial }_y)}\cdot \mathcal{F}({\partial }_y))}\right),$$

$$\underset{u}{\min }\sum _s\left((u_s-g_s{)}^2+\lambda \sum _{\ast \in \{x,y\}}{\left(\nabla u_{\ast ,s}-{\left(f_{BLF}(\nabla g_{\ast })\right)}_s\right)}^2\right),$$

$$u={\mathcal{F}}^{-1}\left(\frac{\mathcal{F}(g)+\lambda \sum _{\ast \in \{x,y\}}\overline{\mathcal{F}({\partial }_{\ast })}\cdot \mathcal{F}(f_{BLF}(\nabla g_{\ast }))}{\mathcal{F}(1)+\lambda \sum _{\ast \in \{x,y\}}\overline{\mathcal{F}({\partial }_{\ast })}\cdot \mathcal{F}({\partial }_{\ast })}\right).$$

$$u_s=c_k,s\in N(k),$$

$$u_s=a_k\cdot s+b_k,s\in N(k),$$

$$\frac{\partial u_s}{\partial s}=a_k,s\in N(k).$$

$$\nabla u=f_{BLF}(\nabla g),$$

$$\displaystyle\begin{array}[]{c}\underset{u}{\min}\|u-g\|^{2}\\ s.t.~{}~{}\nabla u=f_{BLF}(\nabla g).\end{array}$$

$$\displaystyle\begin{array}[]{c}\underset{u}{\min}\|u-g\|^{2}\\ s.t.~{}~{}\|\nabla u-f_{BLF}(\nabla g)\|^{2}\leq\epsilon.\end{array}$$

Full text

I Introduction

Edge-preserving smoothing (EPS) has attracted increasing research interests in the fields of both computer vision and computational photography for decades. The main aim of EPS is to smooth out small details or perturbations in images and preserve the major edges and structures. To this end, many approaches have been proposed in the literature which can be roughly classified into two categories: local methods and global methods. Local methods usually calculate output pixel values as a weighted average of input pixel values inside a local window. Bilateral filter (BLF) [1] is one of the well-known local methods and has been widely used in a number of applications such as HDR tone mapping [3], image detail enhancement [4], mesh smoothing [5], noise smoothing [1], artifact removal [6], and etc. As its variant, joint bilateral filter [7] has also been used in flash/no flash filtering [7, 8] and depth map upsampling [9]. There are also methods that built on the BLF for advanced applications such as the rolling guidance filter [10] and the bilateral texture filtering [11]. Since the brute-force implementation of the BLF is computationally expensive, a number of approaches have been proposed to either accelerate the BLF [12, 3, 13, 14, 15] or introduce fast alternatives [16, 17]. There are also other local methods such as weighted median filter [18, 19], guided filter (GF) [20] and its variants [21, 22]. Generally, local methods are computationally efficient. Some of them can even achieve real-time or near real-time image smoothing [16, 15]. However, these approaches can cause artifacts such as gradient reversals and halos in their results [2, 20]. Besides, they are less appropriate in preserving image details at arbitrary spatial scales [2].

Global methods usually formulate the image smoothing as an optimization framework, consisting of a data term and a regularization term that takes global properties of output images into account. The smoothed output image is the solution to the object function. The total variation smoothing [23] is the very early work which regularizes gradients with the $L_{1}$ norm penalty. Recently, Farbman et al. [2] proposed to smooth images with a weighted least squares (WLS) framework which imposed a weighted $L_{2}$ norm penalty on gradients. Their method shows superior performance over BLF [1] and methods based on BLF [24, 4] in producing results free of gradient reversals and halos. Gradient $L_{0}$ norm smoothing was adopted by Xu et al. [25] to sharpen salient edges while smooth out weak edges. There are also other global methods that are designed for specific tasks. For example, Xu et al. [26] proposed the relative total variation smoothing which was specifically efficient for image texture removal. Ham et al. [27] proposed the static dynamic filter which was especially effective for guided depth map restoration [28]. Generally, global methods can have state-of-the-art performance when compared with local methods. However, the superior performance is achieved with much higher computational cost. Even with the recent acceleration techniques [29, 30], the global methods are still an order of magnitude slower than the local methods.

In this paper, we propose a new global method for efficient edge-preserving smoothing. It can show comparable performance with the state-of-the-art global method but run much faster. Its computational cost is comparable with that of the state-of-the-art local methods. The main contributions of this paper are as follows:

• •

We propose a new global method that embeds the BLF [1] in the least squares model for efficient edge-preserving smoothing. We show that the proposed method can produce results free of gradient reversals and halos which are comparable with the WLS [2]. Meanwhile, it can take advantages of the efficiency of the BLF and least squares model, which makes it more than $10\times$ faster than the WLS.

• •

The proposed method is flexible which can be extended by replacing the BLF with its variants [16, 17] to yield different smoothing methods. In addition, the extended methods can be modified to handle more complex applications in an efficient manner.

The rest of this paper are organized as follows. In Sec. II, we briefly present three smoothing methods that are highly related to the proposed method. In Sec. III, we first show our motivation and then present the proposed method together with its possible extensions to handle more applications. Experimental results of our method and comparison with other state-of-the-art methods are shown in Sec. IV. We conclude our work in Sec. V.

II Related Work

In this section, we briefly present three smoothing methods that are highly related to our method: the BLF [1], the WLS framework [2] and the least squares (LS) framework which is an un-weighted version of the WLS framework.

II-A The Bilateral Filter

Given an input image $g$, each pixel in the output image $u$ of BLF [1] is a weighted average of its neighborhoods in $g$:

[TABLE]

where $s$ and $t$ denote pixel positions. The spatial kernel $G_{\sigma_{s}}$ and the range kernel $G_{\sigma_{r}}$ are typically Gaussians, where ${\sigma_{s}}$ determines the spatial support and ${\sigma_{r}}$ controls the sensitivity to edges.

The advantage of BLF is that it has fast approximations [12, 3, 13, 14, 15] and alternatives [16, 17] which can achieve real-time or near real-time image processing. The disadvantage is that it can produce results with gradient reversals and halos which are also shared by its variants [16, 17]. The reason of gradient reversals is that edges are sharpened in the smoothed image and these edges are boosted in a reverse direction in the enhanced image. This usually happens when small $\sigma_{r}$ is adopted. Halos may occur when large $\sigma_{s}$ and $\sigma_{r}$ are adopted. Fig. 1(b) shows examples of these two kinds of artifacts with image detail enhancement.

II-B The WLS Model

For an input image $g$, the output image $u$ of the WLS [2] smoothing is the optimum of the following objective function:

[TABLE]

where $s$ denotes the pixel position, $\nabla u_{x,s}$ and $\nabla u_{y,s}$ represent the gradients of $u$ at $s$ along $x$-axis and $y$-axis, respectively. The weights $\omega_{x,s}$ and $\omega_{y,s}$ are defined as:

[TABLE]

where $\mathcal{\ell}$ is the log-luminance channel of the input image $g$. The exponent $\alpha$ determines the sensitivity to the gradients of $\mathcal{\ell}$, $\varepsilon$ is a small constant that prevents division by zero in areas where $g$ is constant.

The objective function in Eq. (2) is quadratic and has closed-form solution which can be obtained by solving the following large sparse linear system:

[TABLE]

where $L_{g}=D_{x}^{T}W_{x}D_{x}+D_{y}^{T}W_{y}D_{y}$. $W_{x}$ and $W_{y}$ are diagonal matrixes containing $\omega_{x,s}$ and $\omega_{y,s}$. $D_{x}$ and $D_{y}$ are forward discrete difference operators along $x$-axis and $y$-axis direction. $D_{x}^{T}$ and $D_{y}^{T}$ are backward difference operators. Hence, $L_{g}$ is a five-point non-homogenous Laplacian matrix.

The advantage of WLS smoothing is that it can produce results free of gradient reversals and halos. Fig. 1(c) shows examples of image detail enhancement produced by the WLS. The disadvantage is that solving Eq. (4) is time-consuming. Even with recent acceleration techniques [29, 30], it is still an order of magnitude slower than the state-of-the-art local methods [16, 17, 20]. Note that the recently proposed bilateral solver [31] is reported to enable a fast solution to the WLS problem like Eq. (2). However, it is designed to only handle Gaussian guidance weights in Eq. (2) which are in exponential form. It thus cannot solve the problem in Eq. (2) of which the guidance weights are in fractional form as defined in Eq. (3).

II-C The LS model

In fact, the high computational cost of Eq. (4) results from the non-homogeneous Laplacian matrix $L_{g}$. This makes solving Eq. (4) have to be performed in the image domain which requires an inverse of a very large matrix. However, if $L_{g}$ is a homogeneous Laplacian matrix, Eq. (4) can be solved in the Fourier domain, which is much more efficient. By setting $\omega_{x,s}=\omega_{y,s}=1$ for all the pixel position $s$, we can get a homogeneous Laplacian matrix $L_{g}$. This leads Eq. (2) into an un-weighted one:

[TABLE]

Since Eq. (5) is un-weighted, we denote it as least squares (LS) model in this paper. Solving the LS model in Eq. (5) can be achieved as follows:

[TABLE]

where $\mathcal{F}(\cdot)$ and $\mathcal{F}^{-1}(\cdot)$ are the fast Fourier transform (FFT) and inverse fast Fourier transform (IFFT) operators. $\overline{\mathcal{F}(\cdot)}$ denotes the complex conjugate of $\mathcal{F}(\cdot)$. $\mathcal{F}(1)$ is the FFT of the delta function. The plus, multiplication and division are all point-wise operations.

The advantage of Eq. (5) is that it can be solved efficiently. This is because Eq. (6) transforms the matrix inverse in image domain into a point-wise division in the Fourier domain. And since both FFT and IFFT can be implemented very efficiently, solving Eq. (5) with Eq. (6) is much more efficient than solving Eq. (4) in the image domain. For example, for a $1024\times 1024$ color image, solving Eq. (5) with Eq. (6) in the Fourier domain only needs $\sim 0.12$ seconds while solving Eq. (4) in the image domain needs $\sim 6.8$ seconds. The disadvantage is that Eq. (5) is not an edge-preserving smoothing operator which can cause halos in the results. Since the LS model penalizes all the gradients, the gradients in the smoothed image cannot be larger than the corresponding ones in the input image. It thus does not cause gradient reversals. We show the results produced by the LS model in Fig. 1(d) where clear halos exist.

III The Proposed Method

III-A The BLF-LS Model

All the three approaches mentioned in Sec. II directly perform the smoothing on the intensities of each image. Each of them has its advantage and disadvantage where a tradeoff between the smoothing quality and processing speed exists. In contrast, we propose to first smooth the image gradients with the BLF and then embed the smoothed gradients into the LS framework for efficient edge-preserving smoothing. This can enable our method to maintain both the smoothing quality and processing efficiency.

Our method is motivated by the following two facts. First, we illustrate the gradients (along $y$-axis direction) of the image smoothed by the WLS in Fig. 2(b), which correspond to the smoothed image in the first row of Fig. 1(c). Compared with the gradients of the input image in Fig. 2(a), the small gradients in Fig. 2(b) have been removed while the large gradients remain almost un-changed and sharp. This motivates us that we can directly smooth the gradients of the input image with an edge-preserving smoothing operator, for example, the BLF [1], to yield similar gradients. Fig. 2(d) shows the result of smoothing the gradients in Fig. 2(a) with the BLF. Clearly, large gradients are properly preserved with small gradients being removed, which is similar to the one in Fig. 2(b). Second, the LS model in Eq. (5) is in fact a framework that penalizes all the gradients close to zero. It thus cannot preserve the large gradients, as shown in Fig. 2(c). To solve this problem, we propose to modify the LS model in Eq. (5) to make the gradients of the smoothed image close to that in Fig. 2(d) rather than zero. This procedure can be achieved with the following optimization framework:

[TABLE]

where $f_{BLF}(\nabla g_{x})$ and $f_{BLF}(\nabla g_{y})$ denote smoothing the $x$-axis and $y$-axis direction gradients of the input image $g$ with the BLF, respectively. With a sufficient large value of $\lambda$, the gradients of the output $u$, i.e., $\nabla u_{\ast}$, will be close to $f_{BLF}(\nabla g_{\ast})(\ast\in\{x,y\})$. In this way, the smoothed image $u$ in Eq. (7) will be similar to the one produced by the WLS. Fig. 1(e) shows the results obtained with our method. Specially, the image in the first row of Fig. 1(e) is obtained with the smoothed gradients in Fig. 2(d). Our method can produce results free of gradient reversals and halos.

Similar to Eq. (6), solving Eq. (7) can also be implemented with a few FFTs and IFFTs. The output is obtained through:

[TABLE]

Compared with Eq. (6), only two additional FFTs are needed in Eq. (8). Thus, Eq. (8) can also maintain the efficiency of Eq. (6). For example, it only takes $\sim 0.16$ seconds to process a $1024\times 1024$ color image.

Since the BLF smoothing can also be achieved with its fast approximations [12, 3, 13, 14, 15], our method can be implemented in an efficient manner. In our experiments, we adopt the fast approximation proposed by Paris et al. [13] for BLF smoothing. In this way, our method is able to achieve high smoothing quality and high processing efficiency as well.

In summary, our method can be implemented with the following two steps: (1) smoothing the $x$-axis and $y$-axis gradients of the input image with the BLF and (2) solving Eq. (7) with Eq. (8) to obtain the smoothed image. Since our method embeds the BLF into the LS model, we denote our method in Eq. (7) as BLF-LS.

III-B Mathematical Explanation

To have a better understanding of our model in Eq. (7), we present its mathematical explanation in this subsection. It is well-known that in the BLF [1] smoothing, output is assumed to be piecewise constant inside a local support, which can be formulated as:

[TABLE]

where $c_{k}$ is the expected constant value of output pixel values inside the local support $N(k)$.

We assume that the output image is spatially piecewise linear, which is formulated as:

[TABLE]

where $a_{k}$ and $b_{k}$ are the linear coefficients that are assumed to remain constant inside $N(k)$. $N(k)$ is defined the same as that in Eq. (9).

By taking the derivative of $u_{s}$ with respect to $s$, we have:

[TABLE]

Note that $\frac{\partial u_{s}}{\partial s}$ is in fact the gradient of $u$ at $s$, i.e. $\frac{\partial u_{s}}{\partial s}=\nabla u_{s}$. Based on the above facts, we can have the following conclusion: when the image is spatially piecewise linear as formulated in Eq. (10), its gradients are piecewise constant as formulated in Eq. (11).

The comparison between Eq. (11) and Eq. (9) can motivate us with the following facts: When the image is piecewise constant as formulated in Eq. (9), we can adopt the BLF to smooth the intensities of the input image to get the intensities of the output image. When the image is spatially piecewise linear as formulated in Eq. (10), its gradients are piecewise constant as formulated in Eq. (11). Similarly, we can also adopt the BLF to smooth the gradients of the input image to get the gradients of the output image. By denoting the input image as $g$, the output image as $u$, their gradients as $\nabla g$ and $\nabla u$, respectively, we can formulate the above statement as:

[TABLE]

where both $\nabla u$ and $\nabla g$ represent gradients along either $x$-axis or $y$-axis.

Having stated the key conclusion, we also have a problem at hand: we cannot get the final smoothed image $u$ only from its gradients $\nabla u$ which are equal to $f_{BLF}(\nabla g)$. Note that the input image $g$ is also available at the same time. Hence, the smoothed image $u$ can be obtained through the following optimization problem:

[TABLE]

However, the optimization problem in Eq. (15) is essentially intractable because of the constraint. To make the problem tractable, we relax the constraint $\nabla u=f_{BLF}(\nabla g)$ with $\|\nabla u-f_{BLF}(\nabla g)\|^{2}\leq\epsilon$. Then for a sufficient small $\epsilon$, we have $\nabla u\approx f_{BLF}(\nabla g)$. The smoothed image $u$ now can be obtained via the following optimization problem:

[TABLE]

Then for a proper value of $\lambda\propto\frac{1}{\epsilon}$, Eq. (18) can be exactly transformed into our model in Eq. (7). In this way, we show our model in Eq. (7) is in fact an approximation of the smoothing operator that assumes the intensities of the output image are spatially piecewise linear as formulated in Eq. (10). Note that the trilateral filter in [32] also adopted a piecewise linear assumption, however, our method differs from theirs in both the motivation and the mathematical formulations of the approach.

III-C Parameter Discussion

As our method is highly related to the BLF [1], we first show a simple comparison of the smoothing ability between our method and the BLF. Fig. 3 shows two smoothed images and the corresponding detail enhanced images obtained by the BLF and our method with the same $\sigma_{s}$ and $\sigma_{r}$. As observed from Fig. 3(b) and Fig. 3(c), our method can have stronger smoothing on the input image than the BLF. We can explain this phenomenon with Fig. 3(d) and show that this results from the change in the range kernel of the BLF. More specifically, for two distinct pixels $a$ and $b$ in the intensity domain, the corresponding pixels in the gradient domain are denoted as $a^{\prime}$ and $b^{\prime}$. It is clear that the difference between $a$ and $b$ is much larger than the difference between $a^{\prime}$ and $b^{\prime}$. Note that the squared difference is used to calculate the weight between two pixels (see Eq. (1) for more details). Thus, for the same $\sigma_{r}$, the weight between $a^{\prime}$ and $b^{\prime}$ will be larger than the one between $a$ and $b$. This results in a stronger smoothing between $a^{\prime}$ and $b^{\prime}$ than that between $a$ and $b$. The above analysis also shows that we need to set a smaller value of $\sigma_{r}$ in our method than its usually required value in the BLF.

We then elaborate how to control the smoothing strength of our method. There are three parameters in our method: $\lambda$, $\sigma_{s}$ and $\sigma_{r}$. For a fixed $\lambda$, the smoothing strength of our method can be controlled by $\sigma_{s}$ and $\sigma_{r}$ in a similar way as that in BLF. Fig. 4 shows smoothed examples of our method with different values of $\sigma_{s}$ and $\sigma_{r}$. The visual comparison indicates that larger $\sigma_{s}$ and $\sigma_{r}$ result in stronger smoothing on the input image. The reason is that larger $\sigma_{s}$ and $\sigma_{r}$ can have stronger smoothing on the gradients which accordingly results in stronger smoothing on the input image. Fig. 4 also shows the $y$-axis gradients corresponding to the smoothed images. Clearly, small gradients are gradually smoothed out as $\sigma_{s}$ and $\sigma_{r}$ get larger while large gradients remain sharp. This allows our method to avoid halos when the smoothing strength is enlarged.

For the value of $\lambda$, we generally have $\lambda\propto\frac{1}{\epsilon}$ as stated in Sec. III-B. Fig. 5(b) shows the mean square difference between $\nabla u$ (gradients of the smoothed image) and $f_{BLF}(\nabla g)$ (the result of smoothing the gradients of the input image with the BLF). As the value of $\lambda$ increases, the difference between $\nabla u$ and $f_{BLF}(\nabla g)$ becomes smaller, which is consistent with the fact $\lambda\propto\frac{1}{\epsilon}$. According to our experimental results, for small $\sigma_{s}$ and $\sigma_{r}$ (i.e., weak smoothing on the image), our model with a smaller $\lambda=32$ can produce similar results as it with a larger $\lambda=1024$. We show two examples in Fig. 5(c) and Fig. 5(d). However, when $\sigma_{s}$ and $\sigma_{r}$ become large (i.e., strong smoothing on the image), our model with a smaller $\lambda=32$ can produce results with halos along salient edges as illustrated in Fig. 5(e). When a larger $\lambda=1024$ is adopted, the halo artifacts are properly eliminated as shown in Fig. 5(f). Thus, if not specified, we fix $\lambda=1024$ in all the experiments in this paper. This is different from the one in the WLS model where a quite small value of $\lambda$ is usually adopted.

III-D Model Extensions

Our method described in the above subsections can be extended in several ways. First, as there are many variants of the BLF such as the adaptive manifold filter (AMF) [17] and the domain transform filter [16], we can replace the BLF in our method with these alternatives. Especially, we find the normalized convolution (NC) of the domain transform filter [16] can well suit our method and produce promising results. We use AMF-LS and NC-LS to denote our method embedding the AMF and the NC filter, respectively. Fig. 6 shows examples of detail enhancement produced by our AMF-LS and NC-LS. When compared with the results produced by the AMF and the NC filter, our AMF-LS and NC-LS can properly eliminate the gradient reversals. We will further show our method in handling the halos in the experimental section.

In addition, as can be observed in Fig. 6, our AMF-LS can produce results quite similar to our BLF-LS, while the results of our NC-LS are different from those of AMF-LS and BLF-LS. This is because the weighting scheme of the AMF is quite similar to that of the BLF but the weighting scheme of the NC filter is different from that of the other two. This means embedding different methods in the proposed framework can result in methods with different smoothing properties.

Besides the above extension, we can also use our method to perform joint image filtering other than the single image filtering. One straightforward extension can be achieved in a similar way as the joint bilateral filter [7]. In this case, the guidance weights of the BLF, AMF and NC filter in our method are based on the gradients of another guidance image. They can be applied to tasks such as flash/no flash image filtering [7], which will be shown in our experimental part.

In addition, inspired by the rolling guidance filter (RGF) [10], our NC-LS can also be used to handle more complex applications such as image texture removal and clip-art compression artifacts removal. We can use our NC-LS to iteratively jointly smooth the input image with the guidance of the updated image at each iteration. The initial guidance image is set as the output of smoothing the input image with a Gaussian filter of standard deviation $\sigma$. We show our NC-LS can achieve better performance than the RGF. Fig. 7(c) illustrates one texture removal example of our method. When compared with the result of the RGF in Fig. 7(b), our method can better smooth out textures while properly preserve image structures. In a similar way, we denote the above modification as rolling guidance NC-LS. We use $n$ to denote the iteration number of our rolling guidance NC-LS. In general, we find our rolling guidance NC-LS can achieve promising results within three iterations, i.e., $n\leq 3$.

IV Applications and Experimental Results

To validate the effectiveness and efficiency of our method, we apply our method to several tasks including image detail enhancement, HDR tone mapping, flash/no flash image filtering, image texture removal and clip-art compression artifacts removal. For the BLF [1], we adopt its fast implementation in [13]. In all the applications, input images and the gradients, are firstly normalized into range $[0,1]$ before smoothing and then normalized back to their original range after smoothing.

Image detail enhancement is a fundamental application of image smoothing. By subtracting the smoothed image from the original input image, one can obtain the detail layer which contains the high-frequency details of the original input image. Then the detail layer is magnified and added back to the original input image to get the detail enhanced image. If edges in the smoothed image are sharper than those in the original input image, then these edges will be boosted in a reverse direction which causes the gradient reversals [2, 20]. When the edges are heavily blurred, there will be halos along the blurred edges in the detail enhanced image. This usually occurs along salient edges. The GF [20] does not cause gradient reversals where a theoretical guarantee exists. However, it can cause halos in the enhanced images as shown in the region highlighted with the red box in Fig. 8(b). The gradient $L_{0}$ norm [25] can cause clear gradient reversals in the result as shown in the region highlighted with the blue box in Fig. 8(c). This is due to its smoothing properties which eliminate low-amplitude structures while sharpen salient edges. The WLS [2] can produce results free of gradient reversals and halos as illustrated in Fig. 8(d).

As we have stated in Sec. II-A, the BLF can cause both gradient reversals and halos in the results. Small $\sigma_{s}$ and $\sigma_{r}$ could cause gradient reversals while large ones may result in halos. Fig. 8(e) shows the result of the BLF obtained with $\sigma_{s}=12,\sigma_{r}=0.08$. Clear gradient reversals exist in the result as shown in the region highlighted with the blue box. When larger $\sigma_{s}=12,\sigma_{r}=0.3$ are adopted, halos occur along salient edges as shown in the region highlighted with the red box in Fig. 8(g). We show the results of our BLF-LS in Fig. 8(f) and (h) which have similar smoothing strength on the input image as those in Fig. 8(e) and (g). However, our method can properly eliminate the gradient reversals and halos. Similar comparison between the results of the AMF [17] and our AMF-LS is shown in Fig. 8(i)$\sim$(l).

HDR tone mapping is another popular application that needs edge-preserving smoothing. It can be achieved by decomposing the HDR image into a piecewise smooth base layer and a detail layer [3, 2, 35]. The base layer is the smoothed output of input HDR image. The base layer is then nonlinearly mapped to a low dynamic range and re-combined with the detail layer to get the final tone mapped image. We adopt the framework in [3] where layer decomposition is applied to the logarithmic HDR images. Similar to image detail enhancement, gradient reversals and halos are also the challenging issues in HDR tone mapping. Fig. 9 shows results of different methods. Similar to the results in image detail enhancement, except for the WLS, all the other compared methods have either gradient reversals (highlighted with the blue box) or halos (highlighted with the red box) in their results. The AMF, BLF and NC filter have both artifacts in their results. On the contrary, our AMF-LS, BLF-LS and NC-LS can properly eliminate these artifacts in the results which are comparable with that of the WLS.

Flash/no-flash filtering was proposed by Petschnigg et al. [7] to smooth a noisy no-flash image with the guidance of the flash image. In [20], He et al. showed that gradient reversal artifacts also existed in the smoothed results of joint bilateral filter [7]. Fig. 10(d) shows one example. In our experiments, we find that the result of AMF [17] also suffers from this problem as shown in Fig. 10(c). Our method can also be used to perform this kind of joint image filtering in a similar way. We show the results of our AMF-LS and BLF-LS in Fig. 10(e) and (f) which are comparable to that of GF [20] in Fig. 10(b).

In the above paragraphs, we compare our method with other state-of-the-art local and global methods through several applications. Besides the effectiveness in handling gradient reversals and halos, our method can also run in an efficient manner. We summarize the detailed comparison between our method and the other methods in Table I. Especially, our method has the same properties as the WLS [2] but our method is more than $10\times$ faster. As a global method, the computational cost of our method is even close to that of the GF [20] which is a local method. The running time of $L_{0}$ norm smoothing and WLS is based on the implementations provided by the authors. For GF, AMF and NC, we adopt their official released version in OpenCV. Our method is implemented in Python, and we import the “cv2” package to make use of the implementation of AMF and NC in OpenCV. We re-implement the BLF with Python based on its fast approximation proposed by Paris et al. [13]111Their MATLAB implementation can be download here: http://people.csail.mit.edu/jiawen/software/bilateralFilter.m..

Besides the above applications, our rolling guidance NC-LS can also be applied to more complex tasks: image texture removal and clip-art compression artifacts removal. We adopt the implementations, which are provided by the authors, of all the compared methods for comparison.

Image texture removal aims at extracting semantically meaningful structures under the complex texture patterns, which could be regular, near-regular or irregular [26]. Fig. 11 shows texture removal results of different methods. The relative total variation (RTV) [26] can efficiently smooth the textures while properly preserve edges. However, it can also smooth out some small structures as labeled with the red arrow in Fig. 11(c). Besides, it is prone to smooth out the shading on the object surface [11]. The RGF [10] and bilateral texture filtering (BTF) [11] perform better in restoring the shading on the object surface but they cannot completely remove the image textures as shown in Fig. 11(b) and (d). In contrast, our rolling guidance NC-LS can properly smooth out the image textures while preserve small structures. It also performs better than the RTV in restoring the shading on the object surface as shown in Fig. 11(e). In addition, our method are more efficient than the RTV and BTF. They take 3.2 seconds and 5.3 seconds, respectively, to process the image. In contrast, our method only needs 0.27 seconds. The RGF needs 0.21 seconds which is a little faster than ours, however, our method shows better performance.

Clip-art compression artifacts removal aims to remove compression artifacts in compressed clip-art/carton images. This kind of images are piecewise constant with sharp edges. Clear compression artifacts exist around edges when they are compressed. General denoising approaches do not suit this task as the compression artifacts are strongly correlated with edges. The early work would be the image analogies approach proposed by Wang et al. [33], but their method cannot handle images with a high compression rate as shown in Fig. 12(b). The gradient $L_{0}$ norm smoothing [25] cannot completely remove the strong noise near edges as shown in the region highlighted with the blue box in Fig. 12(d). The $L_{0}$ norm region fusion [34] also shares the same problem. In addition, it can also cause clear artifacts along salient edges as shown in the highlighted region in the red box in Fig. 12(e). The BTF [11] was also applied to remove the compression artifacts. As shown in Fig. 12(c), it can properly eliminate the strong noise near edges but it cannot preserve sharp edges. Fig. 12(f) shows the result of our rolling guidance NC-LS. As can be observed from the figure, strong noise near edges are well smoothed while sharp edges are preserved. In term of processing time, the gradient $L_{0}$ norm smoothing takes 0.51 seconds, the $L_{0}$ norm region fusion takes 0.12 seconds, the BTF takes 1.8 seconds and our rolling guidance NC-LS only takes 0.1 seconds.

V Conclusion

In this paper, we propose a new global method that embeds the bilateral filter [1] and its variants [17] [16] in the least squares framework. The proposed method can have comparable performance with the state-of-the-art global method in producing results free of gradient reversals and halos. In addition, since the proposed method can take advantages of the efficiency of the bilateral filter, its variants and the least squares model, it can run much faster. The computational cost of our method is even comparable with the state-of-the-art local methods. We validate the effectiveness and efficiency of the proposed method in a number of applications with comprehensive experimental results. We also show the flexibility of our method with several extensions in handling complex applications in an efficient manner.