Considering the general form of bilateral filtering Tomasi and Manduchi, 1998, for the input image f, the output result is obtained by scale adaptive bilateral filtering, written as:

where u (p) is the output value at the pixel p, and w (l) and φ (t) represent the spatial kernel and the range kernel, respectively. We use a box function for the spatial kernel in this paper, that is, the window Ω_p of the spatial kernel centered at the pixel p. We assume that W_p represents the scale of the spatial kernel at the pixel p, then Ω_p can be expressed as |Ωp|=(2Wp+1)2 and q is the pixel that belongs to Ω_p. The Gaussian range kernel φ (t) in Tomasi and Manduchi (1998) is defined as:

where t is the intensity difference between the pixels p and q. The parameter σ_r is the standard deviation of the Gaussian kernel, which determines the width of the range kernel, that is, the smoothing parameter, and σ_r is fixed at each pixel. A small σ_r gives rise to superior structure-preservation and inferior texture smoothing, and on the contrary, a large σ_r gives rise to better texture smoothing but the undesired blurring of the structure edges. Hence, it is significant to find an appropriate parameter σ_r for achieving better structure-preserving texture filtering results.

In our proposal, we apply the large spatial kernel sizes in the homogeneous regions for texture elimination and the small sizes near structures for edge-preserving. And the kernel size at each pixel is adaptively optimized by structure measurement. So we calculate the structure information as follows.

First, we blur the input image f using a Gaussian filter to get image f_σ. The gradient of the image f_σ is calculated by:

where G_p is the gradient value at the pixel p, and ∂_xf_σ and ∂_yf_σ are the partial derivatives of f_σ in x and y directions.

The gradient map G is calculated by Equation (3), as shown in Figure 1B. It is clear that the textures with strong gradients are also preserved when only gradient information is considered, that is, the structures and textures with similar gradients cannot be completely distinguished. Therefore, in our proposal, we further conduct four-directional structure detection relying on gradient information. For each pixel, the detection neighborhood is a (2m + 1) × (2m + 1) neighborhood centered at it. To determine a more accurate structural inspection value for each pixel, we detect the (m + 1) × (m + 1) sub-neighborhoods located in four directions. Figure 2 shows the sub-neighborhoods in four directions for structure detection. To be more specific, taking into account the distance from the pixel to the exampled pixel, we computer the weighted average of the gradient values in each sub-neighborhood to obtain Ap(j), j = {NW, NE, SW, SE}, which is used to evaluate the appearance of structure edges.

In the four detecting neighborhoods of each pixel, a strong structure edge corresponds to a large G_p value while a weak structure edge corresponds to a small G_p value. For this reason, we adopt the Gaussian function as the weight to calculate the Ap(j) in the four directions, and the maximum value of Ap(j) is selected as the result of structure measurement for each pixel.

where Ψp(j) represents jth sub-neighborhood, b is the pixel that belongs to Ψp(j), g_m (p, b) is the Gaussian function of the distance between the pixel p and b, and max {•} represents the maximum value of the elements in the bracket. Ap(j) is the comprehensive manifestation of the structure edges in the jth sub-neighborhood for the pixel p, whose maximum value S_p denotes more likelihood of the edges occurring. Therefore, the larger value of S_p implies less smoothing and the smaller value of S_p implies more smoothing around the pixel p.

From the analysis of structure measurement, a large value of S_p suggests that the pixel is in the vicinity of the structure edges, where the scale of the spatial kernel should be adjusted as small as possible. Conversely, the spatial kernel scale should be adjusted as large as possible in textural regions. To estimate the scale W_p in terms of S_p, we establish an inverse mapping function from S_p to W_p, so that the function satisfies the above conditions. The mapping can be expressed as:

where Sp2 is the square of S_p. λ is the denominator of the base of an exponential function, whose value must be greater than 1 to ensure (1/λ)Sp2 ranges in [0, 1]. η is the upper limit of the scale of filtering windows, so η(1/λ)Sp2 denotes the estimated value of the spatial kernel scale. The introduction of δ is to keep the size of windows from approaching 0 so that it prevents the filtering result from over-sharpening or aliasing (δ = 1 by default). Therefore, W_p ranges in [ δ, η ].

The brute force computation of Equation (1) requires O(Wp2) operations for each pixel, which is time-consuming in practical applications, especially in textural regions, where the scale W_p is usually large. For the computational limitation of traditional bilateral filtering, various acceleration algorithms have been proposed to approximate the bilateral filter (Chaudhury, 2011, 2015; Chaudhury et al., 2011), whose computational complexity is decreased to O(1), that is, the complexity no longer depends on the scale W_p. However, some of these algorithms cannot guarantee that the error of the approximate value of the discrete points is within the tolerance range, and the poor approximated estimation may result in color distortion in the filtering image.

In this paper, we adopt the Fourier expansion of the range kernel in Ghosh and Chaudhury (2016) to approximate the scale adaptive bilateral filter. Specifically, Equation (2) can be approximated in another manner:

where τ² = −1, v = π/T, φ^(t) is an approximate estimate of φ (t), N denotes the order of Fourier expansion, c_n is the corresponding coefficient, t is the pixel intensity differences in Ω_p, and the range of t is {−T, ⋯ , 0, ⋯T}, where T can be calculated by:

For all t ∈ [−T, T], the following constraint must be satisfied:

where ε is the tolerance of the approximation for the Gaussian range kernel (ε = 0.01 by default).

For the given range kernel φ (t) and tolerance ε, the specific solution of the approximation order N, and the corresponding coefficients c_n is provided in Ghosh and Chaudhury (2016). By using Equation (6) to approximate Equation (2), we can reformulate Equation (1) as:

where û (p) is an approximation of u (p), E (p) and H (p) represent the approximate value of numerator and denominator of Equation (1), respectively, which can be expressed as:

We can further express Equations (10) and (11) as:

where e_n (p) and h_n (p) are expressed as follows:

Since a box function is employed for the spatial kernel, in conclusion, the adaptive bilateral filtering can be decomposed into a series of box filtering. Therefore, Equations (14) and (15) can be simplified as follows:

It can be added point-by-point in the neighborhood of the pixel p to obtain e_n (p) and h_n (p), whose computation is expensive. Hence, in our proposal, we compute Equations (16) and (17) by the recursive algorithm in Crow (1984). We assume that the integrated element of pixel p in e_n (p) is r (q):

First, we compute the integral image R (p) at the pixel p:

where (x, y) is the coordinate of pixel p and (k₁, k₂) is the coordinate of the pixel in the integral region.

By using recursive theory, the integral image R (x + 1, y + 1) at the pixel (x + 1, y + 1) can be expressed as:

For any scale W_p, e_n (p) can be computed as follows:

Similarly, h_n (p) can be obtained.

In conclusion, we can calculate the Equations (10) and (11) according to e_n (p) and h_n (p), and instead of directly computing scale adaptive bilateral filtering, we replace each convolution with pointwise operation through Fourier expansion of the range kernel, as shown in Equations (12) and (13). Furthermore, we can compute e_n (p) and h_n (p) at O(1) complexity with a recursive algorithm, that is, e_n (p) and h_n (p) require a fixed number of operations for any scale W_p.

To be specific, since Equation (21) requires three additions, this means that it takes three additions to compute both Equations (16) and (17). In summary, we can compute Equations (16) and (17) using addition operations, then we can compute Equations (12) and (13) in terms of Equations (16) and (17) by pointwise operations. It is quite clear that the computation of Equation (9) is based on Equations (12) and (13), which proves that the scale adaptive bilateral filter can be computed at O(1) complexity.

We compute the approximation of the output value in this paper, particularly, we consider the error to be:

which provides the largest difference between the exact and approximate scale-adaptive bilateral filtering pixelwise.

According to the Equation (6), the error comes from the approximation of the range kernel, meanwhile, for all t ∈ [−T, T], |φ(t)-φ^(t)|≤ε. From the conclusion of Ghosh and Chaudhury (2016), we can ensure that Equation (22) is within some tolerance:

Since there are complex and irregular textures in many natural images, generally, single iterative filtering cannot completely smooth out the textures. Considering this limitation, we adopt the multiple iteration operation of adaptive bilateral filtering in this paper. Algorithm 1 summarizes the overall process of our method.

Our method is implemented using MATLAB. In our algorithm, the relevant parameters are σ_r, m, λ, η, and I. σ_r determines the scale of the Gaussian range kernel, as we adopt the suggested setting by Ghosh et al. (2019): σ_r ranges in [20, 40]. While evaluating structure measurement, generally, we manually set the radius of detection neighborhood m = 4 to practice a majority of cases. λ is used to normalize the values of structure measurement to the interval [0, 1], so we fix λ = 10 throughout.

The relatively important parameters are the upper limit of the spatial kernel scale η and iteration number I. The value of η depends on the roughness of the textures and the sharpness of the structures and by the parameter recommendation of Ghosh et al. (2019), we set η ranges in [8, 16]. In most situations, setting I = 2 can achieve the desired filtering visual effect. Figure 3 shows the filtering results in different combinations of η and I.

For subjective evaluation, we compute our algorithm with the state-of-the-art texture smoothing techniques, including relative total variation (RTV) (Xu et al., 2012), structure gradient and texture decorrelation regularization (SGTD) (Liu et al., 2013), rolling guidance filter (RGF) (Zhang et al., 2014), bilateral texture filtering (BTF) (Cho et al., 2014), scale-aware structure-preserving texture filtering (SATF) (Jeon et al., 2016), and relativity-of-Gaussian (ROG) (Cai et al., 2017). Generally, we use the suggested parameters to obtain optimal filtering results for previous methods. In Figures 4–6, we display the visual effect comparison for three images containing various textures and structures. The reason why we choose these three images is that they contain different types of textures and different shaped structures, which can illustrate the superiority of our algorithm from many aspects.

Figure 4 shows the filtered results of different methods on the mosaic art “Pompeii fish mosaic,” where the image contains coarse textures and highlighted small-scale structure edges. It is observed that all methods can eliminate fine details in homogenous regions; however, the approaches of SGTD and ROG perform better in removing high-contrast textures effectively. Moreover, for the preservation of small structures highlighted in the image, the methods of SGTD, RGF, BTF, and ROG can hardly preserve the fine structures of fish's eyes which are overly smoothed since the size of the windows is oversize. In the enlarged box, we can clear that the methods of RTV, SGTD, and ROG may result in excessive sharpness near structure edges, which appears as an unwanted jaggy artifact.

Compared with these existing advanced methods, our algorithm works better in eliminating coarse textures while preserving main structures as much as possible in Figure 4H. Particularly, our method can completely preserve the structure of fish's eyes.

Figure 5 shows the smoothing effect on a face image. Especially, we focus on the region highlighted with the red box, whose meaningful structures and textures on the left and right sides of the nose bridge are very similar in appearance. Since the previous methods apply texture filtering with the fix-scale kernel to remove textures, the visual effect is not always well. The results of RTV, BTF, and ROG exist unwished artifacts in the bridge of the nose. On the side, the methods of RGF and SATF perform poorly when removing high-contrast textures.

Our algorithm handles the pixel around structures with a small scale and the pixel in the textural region with a large scale. In Figure 5H, we obtain a better filtering result than the state-of-the-art methods and our method can remove coarse textures without creating artifacts.

Figure 6 shows small structures comparison of different filtering results on the mosaic art “fish.” All the existing methods blur the fine structures and cause artifacts near edges. Relatively serious are the results of RTV, SGTD, RGF, and ROG, and the whiskers and teeth of fish even became sticky. Meanwhile, in methods of BTF and SATF, the teeth of fish are barely preserved.

In contrast, our method achieves the superior property of preserving multi-scale structures, as shown in Figure 6H. The edges and details can maintain the original structure as much as possible.

Figure 7 shows the comparison of denoising effects on a gray image. Intuitively, it can be seen that the effects of RGF, BTF, and SATF methods are not ideal when removing gray image noise, and cannot completely smooth the noise in the background. The methods of RTV and ROG cause edge sharpening in smoothing results.

In comparison, our proposed algorithm can remove the noise of gray images and retain the edge features of people in the image, as shown in Figure 7H.

The widely used image objective quantitative evaluation methods include Peak Signal-to-Noise Ratio (PSNR) and Structural Similarity Index (SSIM). PSNR is an image quality evaluation based on error sensitivity. SSIM comprehensively measures image similarity from the aspects of brightness, contrast, and structure. In our evaluation, we also take the Feature Similarity index (FSIM) (Zhang et al., 2011) and Blind Image Spatial Quality Evaluator (BRISQUE) (Chen et al., 2018) as the evaluation indexes. We selected four ground truth images in Dong et al. (2015), Abiko and Ikehara (2019), and Shen et al. (2015), and then added salt and pepper noise along with periodic noise to these four images as the texture images, as shown in Figure 8, using ground truth images as the reference to calculate PSNR, SSIM, and FSIM. In contrast, the BRISQUE is obtained only by the filtered result.

Table 1 shows the statistics of the mean values of the objective evaluation indexes of four images in Figure 8. First, on the metric of PSNR, our method performs best among these seven methods, which suggests that our results have less image distortion. The methods of RTV and SATF get higher PSNR results that are only inferior to ours. Concerning SSIM results, our method also achieves the highest result. In contrast, we only obtain the third-highest FSIM value, which is inferior to BTF and SATF. In general, the similarity between the results filtered by our approach and ground truth images is relatively good. Finally, we take a look at BRISQUE results, whose smaller score implies better perceptual quality. It just so happens that our method has the smallest BRISQUE value.

To verify that the complexity of our algorithm does not depend on the size of the spatial kernel, we set different parameters for the upper limit of the spatial kernel scale, that is, we change the width of scale to smooth the 400 × 324 image and recode the timings required for a single iteration.

Table 2 shows the statistics of timings for a single iteration of the image in Figure 9A, Figures 9B–D show the different filtering results for different η. It can be seen that the timings required for the filtering process are not much different when the values of η are different, which illustrates that the complexity of the adaptive bilateral filtering does not depend on the size of scales. This result verifies our algorithm that the complexity is decreased to O(1) by our approximation of bilateral filtering.

Our approach can be applied to image detail enhancement (Fei et al., 2017). It aims to highlight image details and improve the visual effects of the images. Figure 10 displays the application of our method in detail enhancement. We first subtract the filtered image from the input image to generate the textures, which are magnified three times and superimposed on the input image, so that we can achieve the purpose of detail enhancement.

The existence of high-contrast textures will keep some irrelevant information and produce false edges in edge detection. Due to the severe influence of textures, we execute our method for texture removal before edge detection. As shown in Figure 11, compared to the edge detection of the original image, the edge map of the filtered image extracted by the canny detection (Canny, 1986) operator is clearer.

The texture smoothing method proposed in this paper can also be applied to image abstraction and pencil sketching. Following (Winnemöller et al., 2006), our method is employed in replacing the bilateral filter to generate abstraction results. Furthermore, we obtain pencil sketching results based on image abstraction. The results are shown in Figure 12.