Fast L 1 Gauss Transforms for Edge-Aware Image Filtering

. Gaussian convolution and its discrete analogue, Gauss transform, have many science and engineering applications, such as mathematical statistics, thermodynamics and machine learning, and are widely applied to computer vision and image processing tasks. Due to its computational expense (quadratic and exponential complexities with respect to the number of points and dimensionality, respectively) and rapid spreading of high quality data (bit depth/dynamic range), accurate approximation has become important in practice compared with conventional fast methods, such as recursive or box kernel methods. In this paper, we propose a novel approximation method for fast Gaussian convolution of two-dimensional uniform point sets, such as 2D images. Our method employs L1 distance metric for Gaussian function and domain splitting approach to achieve fast computation (linear computational complexity) while preserving high accuracy. Our numerical experiments show the advantages over conventional methods in terms of speed and precision. We also introduce a novel and effective joint image filtering approach based on the proposed method, and demonstrate its capability on edge-aware smoothing and detail enhancement. The experiments show that filters based on the proposed L1 Gauss transform give higher quality of the result and are faster than the original filters that use box kernel for Gaussian convolution approximation.


Introduction
Gaussian convolution is a core tool in mathematics and many related research areas, such as probability theory, physics, and signal processing.Gauss transform is a discrete analogue to the Gaussian convolution, and has been widely used for many applications including kernel density estimation [1] and image filtering [2].Despite its reliable performance and solid theoretical foundations, Gauss transform in its exact form along with other kernel-based methods has a drawback -it is very computationally expensive (has quadratic computational complexity w.r.t. the number of points) and hard to scale to higher dimensions.Which is why there have been many attempts to overcome these problems by creating approximation algorithms, such as fast Gauss transform [3], dualtree fast Gauss transforms [4], fast KDE [5], and Gaussian kd-trees [6].Also, box kernel averaging [7] and recursive filtering [8] have been popular in computer graphics and image processing because of their simplicity, see the surveys [9], [10] for numerical comparisons of these approximation methods.Since high bit depth (also dynamic range) images have become popular in both digital entertainment and scientific/engineering applications, it is very important to acquire high approximation precision and to reduce artefacts caused by drastic truncation employed in many conventional methods focused on computational speed.One of the highly accurate methods is called fast L 1 Gauss transform approximation [11] based on using L 1 distance instead of conventional L 2 Euclidean metric.This L 1 metric preserves most of the properties of the L 2 Gaussian, and is separable, hence it allows to perform computations along each dimension separately, which is very beneficial in terms of computational complexity.Also, L 1 Gaussian has only one peak in Fourier domain at the coordinate origin, and therefore its convolution does not have some undesirable artefacts that box kernels and truncation methods usually have.However, this algorithm works only on one-dimensional (1D) point sets, although it can be extended to uniformly distributed points in higher dimensions by performing it separately in each dimension.In order to be able to acquire Gauss transform for nonuniformly distributed two-dimensional points and to further generalize it to higher dimensional cases, we need to extend existing method [11] to the 2D uniform case.In this paper we propose a novel approximation method for fast Gauss twodimensional (2D) image transform.Our method is based on extending the fast L 1 Gauss transform approximation on uniformly distributed 2D points that allows to perform Gaussian convolution quickly while preserving high accuracy.We demonstrate that efficiency of the proposed method in terms of computational complexity, numerical timing, and approximation precision.We also successfully applied our method in the novel filtering approach based on combining the approximated L 1 Gauss transformations into the so-called guided filter [12] (joint image filtering via ridge regression).Our approach reduces computational costs while providing higher quality results compared to the conventional one.We show the application to edge-aware smoothing and image detail enhancement.57

Fast L 1 Gauss Transform
In this section, we briefly describe the 1D domain splitting algorithm [11] employed for fast L 1 Gauss transforms.Consider the ordered point set Each point has a corresponding value e.g.pixel intensity in case of images.The L 1 Gauss transform for each point in set is given by where G(x), x ∈ R, is a L 1 Gaussian function (also called Laplace distribution in statistics) with its standard deviation σ.It is convenient to decompose L 1 norm by splitting its domain by using the point x1 such that Thus, Gauss transform (1) using the equation ( 2) becomes (3) Such representation (3) allows to reduce the amount of computational operations, since values and the sums and can be precomputed in linear time.However, using the equation (3) may imply some numerical issues, such as overflow, if the distance between and is relatively large.To avoid such issues, this algorithm introduced certain representative points (poles) instead of using the single point where the distance between is smaller than the length that causes the numerical instability.Hence the equation (3) becomes more complex form, a highly accurate truncation can be applied where is numerically equal to zero, see [11] for further technical details.Although this algorithm can be used in case of multidimensional images by applying it separately in each dimension, this separable implementation approach is not applicable to nonuniformly distributed high-dimensional point sets.Therefore, we present a novel and natural extension of the domain splitting concept on 2D cases (images) in the following sections.

Two-Dimensional Algorithm
For a given 2D point set L 1 distance between two points in is given by thus the Gauss transform ( 1) is represented by the formula: Domain splitting (2) for 2D points is given by see Fig. 1a for geometric illustration of the domains.4) is numerically troublesome, it is reasonable to divide the space into smaller groups and perform computations separately, as it was proposed in [11].Let us introduce a novel 2D multipole approach for solving this problem.

Consider a set of poles
The distance between points using poles αk is given by where see Fig. 1b for geometric illustration of the domains with their poles.The point is assigned for one representative pole defined by which is the closest pole to that has absolute values of coordinate smaller than For each point the multipole L 1 Gauss transform is given by the equation ( 5), For the sake of simplicity, we assume that the numbers of poles in 2D are same M.
Following [11], M and the poles are given by where [•] is the ceiling function, MAX is the maximum value of precision (e.g., double floating point: DBL_MAX in C programming language), and ϕ is a userspecified parameter (0.5 is employed in our numerical experiments).The above pole selection scheme leads to which theoretically guarantees numerical stability in our method.When the distance between poles is determined by the equation ( 6) and becomes numerically zero if we can efficiently truncate Gauss transform by approximating the values: where In other words, instead of computing terms across all the corresponding point sets, we consider only the neighbouring points, which allows to avoid nested loop structures in our implementation and speed up the computational process.As in the 1D algorithm [11], the terms can be iteratively computed in linear time.Assume that an image consists of pixels and the number of poles along each dimension is M, total complexity of our method is which is a little bit slower than the separable implementation employed in [11] that requires operations.

Numerical Experiments
We held all the experiments on Intel Core i7-6600U 2.60 GHz dual core computer with 16GB RAM and a 64-bit operating system.We compared the multipole version of our algorithm with box kernel (Box) using moving average method [7], the 1D domain splitting (YY14) with separable implementations [11], and Fast Discrete Cosine Transform (FDCT) via the FFT package [13] well-known for its efficiency.
To evaluate the performance of the methods mentioned above we used randomly generated 2D point sets with 10 different sizes from 128 2 to 5120 2 and 10 various values of σ = 5,10,...,50.The radius for the Box method was chosen equal to σ.The timing results (see Fig. 5) show that our method is slightly slower than the 1D domain splitting (YY14) despite its theoretical complexity is much larger.It is worth noticing that the implementation of our method can be further improved by using GPU-based or parallel computing techniques.However, the accuracy evaluation results (see Table 1) show that our method achieves best approximation quality among the discussed methods.We evaluate the precision using and PSNR measures.Consider is the exact result of Gauss transform, is the approximation achieved by a given algorithm, and is calculated using formula We also use peak signal-to-noise ratio (PSNR) [2] to measure the performance of our algorithm according to the equation We performed linear image smoothing by the following normalized convolutions for each color channel:

62
where the denominator is also obtained by our method convolving Gaussian with the image whose intensity is equal to one everywhere.Fig. 3 illustrates the smoothing results using naive implementation (Exact), our method, Box kernel, and FDCT algorithms.The gradient magnitude of smoothed images on Figs. 4 and 6 show that, in contrast to FDCT and box kernel, our method does not produce some undesirable artifacts and is extremely close to the exact implementation.

Edge-Aware Filtering
The proposed algorithm for Gauss transform approach can be applied in various computer vision tasks.We present one of the possible applications of our method by introducing the novel approach for improving the so-called guided filter [12].
We introduce the new approach for computing guided filtering where our Gauss transform algorithm is employed for instead of the box kernel method.As it was shown before, our algorithm gives a much higher quality of smoothing, and this allows us to eliminate smoothing of in the equation ( 8): Thus, using our algorithm involves 2 operations of compared to 4 operations in the original method if (grayscale case), and 21 operations compared to 33 operations if and both of them are color images.We examined edge-aware filtering on color images, where the number of is equal to 21 for the box kernel method and 10 for our approach (9 operations for smoothing of the coefficients and one operation for normalization).As seen on the Figs.7 and 9, our approach with the reduced amount of smoothing operations gives quality of edge-aware filtering higher than [12] with the box kernel method, and is faster (0.24 and 0.28 sec for Figs.9a and 9d respectively).We examine the differences of equations ( 8) and ( 9) in terms of filtering quality on Figs. 9 and 10, which show us that the box kernel method causes artifacts similar to linear filtering case.We also applied our approach for the detail enhancement filter defined by: where is the enhancement parameter.The experiments show that applying our approach for detail enhancement filtering gives high quality results (see Fig. 8).9) and box kernel using eqs.( 9) and ( 8) respectively.

Conclusion
In this paper1 we presented a novel and fast approximation method for L 1 Gauss 2D image transforms.Series of numerical experiments have shown that our method is generally more accurate than the conventional methods and faster than the widely used FFT.We also demonstrated capability of the proposed method in image smoothing application where the conventional box kernel averaging and FFT both suffer from undesirable artifacts.Despite our method is slightly slower than the separable implementations of 1D algorithm [11], this approach can be efficiently used for non-uniformly distributed points.
We have also proposed a novel approach for improving the guided filtering [12] via our Gauss transform and showed its advantages in terms of quality and speed over [12].
Our method is applicable only to uniformly distributed structures, such as images.

Table 1 .
Precision and speed evaluation results (speed measured in Mpix/sec).