CC BY
22
Fast L1 Gauss Transforms for Edge-Aware Image Filtering
1,2Dina Bashkirova <dina.bashkirova@riken.jp>
1Shin Yoshizawa <shin@riken.jp> 2Roustam Latypov <roustam.latypov@kpfu.ru> iHideo Yokota <hyokota@riken.jp> 1 Image Processing Research Team, RIKEN Center for Advanced Photonics, RIKEN, 2-1, Hirosawa, Wako, Saitama, 351-0198, Japan 2Institute of Computational Mathematics and Information Technologies, Kazan Federal University, 35 Kremlyovskaya, Kazan, Russia, 420008
Abstract. 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 LI 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 LI Gauss transform give higher quality of the result and are faster than the original filters that use box kernel for Gaussian convolution approximation.
Keywords: Gaussian smoothing, Laplace distribution, fast approximation algorithms.
DOI: 10.15514/ISPRAS-2017-29(4)-4
For citation: Bashkirova D.R., Yoshizawa S., Latypov R.H., Yokota H. Fast LI Gauss Transforms for Edge-Aware Filtering. Trudy ISP RAN/Proc. ISP PAS, vol. 29, issue 4, 2017, pp. 55-72. DOI: 10.15514/ISPRAS-2017-29(4)-4
1. 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 Ll Gauss transform approximation [11] based on using 1} distance instead of conventional L2 Euclidean metric. This Ll metric preserves most of the properties of the L2 Gaussian, and is separable, hence it allows to perform computations along each dimension separately, which is very beneficial in terms of computational complexity. Also, L1 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 (ID) 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 non-uniformly 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 two-dimensional (2D) image transform. Our method is based on extending the fast Ll 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 Ll 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.
2. Fast L1 Gauss Transform
In this section, we briefly describe the ID domain splitting algorithm [11] employed for fast 1} Gauss transforms.
Consider the ordered point set X = Xi € R, ж» Vi = 2, JV. Each
point Xi has a corresponding value Ii € M, e.g. pixel intensity in case of images. The 1} Gauss transform for each point in set X is given by
N | |
J(xj) = ^G(xj-Xi)Ii, G(x) = exp(——), (1) i = 1
where G(x), x 6 R, is a Ll Gaussian function (also called Laplace distribution in statistics) with its standard deviation a. It is convenient to decompose L1 norm by splitting its domain by using the point xi such that
(2)
. . J \xj — x-i | — \xi — x\ | if a?! < Xi < Xj,
3 1 \ \'Xi — III — I Xj — Iij if < Xj < Xi.
Thus, Gauss transform (1) using the equation (2) becomes
J-1 i
Jixj) = Ii + G(xj-x1)Y,G{xi*_Xi) + 1
xi> i=j+i (3)
Such representation (3) allows to reduce the amount of computational operations, since values G(xj — x{), G( r x_ r | j and the sums hG(xi — x{) and
2i=i G(xj-xi) can precomputed in linear time. However, using the equation (3) may imply some numerical issues, such as overflow, if the distance between and Xl, I £ {hj} is relatively large. To avoid such issues, this algorithm introduced certain representative points (poles) {»fc <= R} instead of using the single point Xi, where the distance between at and xi 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 G(atk — Xj) 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.
3. Two-Dimensional Algorithm
For a given 2D point set X = {Xj}^.,, x,; = {<Ci,yi) € M2, L1 distance between two points in is given by |X; -x,:| = \xj-Xi\ + \yj-y,\, thus the Gauss transform (1) is represented by the formula:
N
J(xi) = >>xp(-
- Xi\ + I У] - Уг
i=1
)h
Domain splitting (2) for 2D points is given by
'\Xj -Xl\ - \xt - Zi| + [ijj -yi I - \yi - Jill
\xi-X\\ - \xj X\| + \yj -3/11 - \yi -yi\ |Xj - III - \xi - III + \yi - î/î I - \yj - yi\ Jxj - III - I Xj - III + \vi -yi| - \yj - yi| see Fig. 1 a for geometric illustration of the domains.
I Xj -Xi\ + I ]jj -Vi\ = <
if X| G Di if Xi € D2 if x, € D:i if х^ e D\,
■xi ......
Di d2
d4 d3
d9 d d8 d9
d3 d d2 d3
d6 d5 I Xi j D4 d6
d9 Ok+i d7 p i °8 <*k+l,p+l d9
(a) Single pole Xj case (b) Multipole (at! case
Fig. 1. Illustration of 2D domain splliting. Using the above decomposition, Gauss transform is represented similar to (3):
./(X.) = i(xj} + Fixi) sp F(Vi) lh: ) . F(Vi) V4 F(xi) u.. л
xieDj(j) v
where F(xj) = G(xj — Xi) and F(yf) = G(yj — yi).
Precomputation and storage of values p^.] and require 0(AN) operations
and 0(AN) space, and all the subsequent sums F(xj)F(yj), can be iteratively computed in 0(N) operations. Gauss transform for all points using the formula (4) requires 0(1 CW) as opposed to employing the separable implementation of equation 58
(3) for 0(6N) operations. Since computing the Gauss transform using the equation (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 {ak}jf= i> = (<ik,bk) € K2- The distance between points using poles at is given by
\Xi - dk\ - \Xj - dk + lift — Ьлг I — I yj -bk \ f X, eDj
- ak - U-i - ak + ~h\~ I Vj -bk\ f x, €D2
l^i - a*| + -ak + \Vi ~h\~ |У, -4 f Xj G D3
\xi ~ ak\ _ \Xj - ak + 1 Уз - bfc| - \ш ~bk\ f X, G йл
\xj - ak - - a,k + \Уз - Ь/bl - Ы ~bk\ f x, € Д,
- ak\ + \Xj - ak + I Vj - bfc| - |yi -bk 1 f Xj € Da
Is* - a,k\ - Xj - ik + 1 Vi -h\ + \Vj -bk\ f X e Du
1 Xj - <ik - \Xi - dk + \ы -bk I f Xj € £>s
~ o-k\ 1 Xj - o-k + I Vi - h\ + 1 yj -bk 1 f X, G D9
where
Dy = {xilxi € D^ Vi e D\}.D2 = {xi|xi e Dl,yi € D\}.
£»3 = {x^ e € = {x.l^ 6 Di.yi €
Ds = {Xilii e Dlvi e D\),D6 = {x^ 6 DJ.j/i € D'i),
D7 = {xilxi e Df,y< e D^.Dg = {Xi|i< e € £>?},
A> = {x,|x, € £>f?i/i e Z).?},
= {£i|afc ^ ii or xj < xi < afc}>
= {xi|afc < Xj < Xi or Xi < Xj < at},
£>3 = {xi\xi <ak< Xj or xj <ak< Xi},
Di = {Vi\bk <Vi< Vj or yj <yi< 6fc},
Dl = {yi\bk < Vj < Vi or yi < yj < bk}.
D3 = {Vi\Vi <hk< Vj or yj <bk< //, },
see Fig. lb for geometric illustration of the domains with their poles. The point x j is
assigned for one representative pole defined by
ak(xj) = max{at|at < Xj, bk < yj}, k
which is the closest pole to Xj that has absolute values of coordinate smaller than Xj. For each point Xj, the multipole L1 Gauss transform is given by the equation (5),
+ E 4 + E E ci+ E E ts)
afcej?£. ofc^Dr ttfc^a
4 = E GteWvJii, Bi = e(Xj)c(K) y = E G(*0<?(vOA,
x;-A(fe) Xi-A(fc) ^ ^ Xi-A(fe)
A(fc+1)-1 A(fc+1)-1
* = __ E ( fgf" = Igf xE( «(*<№(*)/,
where g^-) = G(xf-ak), G(yj) = G(Vj - &*),
and A(-) is an index function defined by
A(fc) = min (x, |afc < Xj < and bk < yj < bk+i). l<j<N
For the sake of simplicity, we assume that the numbers of poles in 2D are same M. Following [11], M and the poles {a*..}- are given by
r i ri \ {0: f; 2, —, (A/ — !)}•»;
W} = {&*} = --jj-—, (6)
w
w = max(l*i - xN\, \n - yN\), M = l^log(MAX)]
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 user-specified parameter (0.5 is employed in our numerical experiments). The above pole selection scheme leads to m.ax(G(a£+i — a^), G(bk+i — bk)) < MAX which theoretically guarantees numerical stability in our method.
When the distance between poles is determined by the equation (6) and G(«fc — Xj) becomes numerically zero if — ' we can efficiently truncate Gauss
transform by approximating the values:
£ £ 4, £ Bi~ £ H,
otkGD,, afe€M (£>9) «(.eDj ak€ii(D7)
£ eg" £ Cl Y, £ Dl
akGDg otk€p(Dg) ak£D3 ak€fi{D3)
£ £ Ei
ak€Dc, ak£iJ.(De)
where /-t(A.) = {*i € D* I ~ afc(*i)| < ^jj}-
In other words, instead of computing terms Ak, Bk, Ck, Dk, Ek 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 ID algorithm [11], the terms can be iteratively computed in linear time. Assume that an image consists of Vn x Vn pixels and the number of poles along each dimension is M, total complexity of our method is 0(16^+2^+4^ ) which is
a little bit slower than the separable implementation employed in [11] that requires 0{\2N + 2s/N + M) operations.
(a) Input image 1 (b) Input image 2
Fig. 2. Input images.
4. 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 ID 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 1282to 51202and 10 various
values of a = 5,10.....50. The radius for the Box method was chosen equal to a. The
timing results (see Fig. 5) show that our method is slightly slower than the ID 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 Emux and PSNR measures. Consider Ic is the exact result of J_} Gauss transform, I" is the approximation achieved by a given algorithm, and mux is calculated using formula
£niax =,51.aJCr <k-1 <t<N
We also use peak signal-to-noise ratio (PSNR) [2] to measure the performance of our algorithm according to the equation
PSNR = -101og(f;( * )').
fri max(/|, If )
We performed linear image smoothing by the following normalized convolutions for each color channel:
fG(x-y)I(y)dy J(Xj)
fG(x - y)dy EfG(Xj-x,)
where the denominator is also obtained by our method convolving L' 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, andFDCT 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.
Table 1. Precision and speed evaluation results (speed measured in Mpix/sec).
Our YY14 FDCT Box
A'max 1.8x10 11 3.8* 1CT10 0.44 3.73
PSNR 291.05 281.81 58.98 41.45
Speed 7.19 9.76 3.37 8.58
(a) Exact
Я
(Ъ) Our
(c) Box (d) FDCT
Fig. 3. Results of smoothing (a = 20), where the input image is given by Fig. 2a.
i? >
(a) Exact (b) Our
v$ > M. iirVH
(c) Box (e) Exact (d) FDCT (ft Our
(?) Box (h) FDCT
Fig. 4: Visualisation of\ VL\for comparison of artifacts (a = 20).
10 15
Image sixe, Mpix
Fig. 5: Timing with respect to image size (averaged by a).
Гцувлнц
(a) Exact
(b) Our
(c) FDCT
(a) Exact
(b) Our
(c) FDCT
Fig. 6: Visualisation of\ VL\for comparison of artifacts of FDCT (a = 20), where the input
image is given by Fig. 2b.
5. 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].
Guided filter is categorized into a joint image filtering technique consisting of two input images where one of them is called guidance image, and reflects guidance colors into the other input. One of the most popular joint image filters is the joint bilateral filter [14] which averages the neighbouring colors using the weights that depend on the guidance image. Guided filter is an approach for joint image filtering that allows to overcome a problem with the undesirable gradient reversal artifacts that joint bilateral filter suffers from. Besides edge-aware filtering, it has various image processing applications such as matting, flash-noflash synthesis, HDR-compression, and haze removal.
Consider a point set X = {x,}^!, x,; = (xt,yi) e M2, a guidance image g = g(x) e K, an input image /{xj) el, a desired output image //(x;) e M, and an image region fi(x) centered at x. The guided filter is defined as the following linear transformation:
H(y)=aS(y) + fr,yen(x),
where a, b e M are the coefficients constant in S7(x) that depend on the input image I. Such representation is very useful for image processing tasks, since it preserves the gradient extrema V// — «V;/, and hence the edges of the guidance image. The coefficients a and b are obtained using the linear ridge regression model [15]:
yeii(x)
where W(x-y) is the weight that determines the importance of the point y in il(x) and f is the regularization parameter. One can obtain values a and b by minimizing K(a, b): S-K{a, b) = 0 and b) = 0. This leads to the following representation:
Here /(») is an averaging function. Since a point y is included in many overlapping regions fi(x) and values a and b for y are different for each region, the final coefficients are found by averaging over all possible values of y: H(x) = f(a)g(x) + f(b). (8)
Guided filtering of color images involves inversion of 3 x 3 coefficient matrix in order to solve the equation (7) (see [12] for further details). If we set I = g, then the guided filter preserves salient edges while smoothing the flat regions (edge-aware filtering). In the simplest case of I = g and I being is a grayscale image, computing guided filter involves performing 4 smoothing operations (e.g. f(I), f(I2), f(a), f(b)), and it takes 33 smoothing operations for a color image if I ^ g. Which is why the choice of the smoothing operator /(*) is crucial, since it determines the overall speed and quality of filtering. Authors of the guided filter [12] suggested employing classic L~ Gauss transform or box kernel method but prefer the latter due to its simplicity and speed despite the fact that box kernel produces undesired artifacts discussed above.
We introduce the new approach for computing guided filtering where our Ll 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 /(a) and f(b) in the equation (8):
H{x.) = ag{x) + b (9)
Thus, using our algorithm involves 2 operations of /(*) compared to 4 operations in the original method if I = g (grayscale case), and 21 operations compared to 33 operations if I ^ g 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:
Dfx) = I(x) + t(/(x) - H(x)), where r is the enhancement parameter. The experiments show that applying our approach for detail enhancement filtering gives high quality results (see Fig. 8).
Fig. 10: Edge-aware filtering results (a=8, s=0.04). a: input image, b-d: visualization of gradients \ VH \ of edge-aware filtering via our approach, eq. (9) and box kernel using eqs.
(9) and (8) respectively.
6. Conclusion
In this paper1 we presented a novel and fast approximation method for L1 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 ID 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 L1 Gauss transform and showed its advantages in terms of quality and speed over [12].
1 It is an extension of our previous work [16]. The main difference from [16] is the novel approach to joint image filtering and its numerical experiments.
Our method is applicable only to uniformly distributed structures, such as images. Hence our future work includes extending the proposed method to higher-dimensional nonuniform cases which can be done for example by using treelike structures. We also would like to investigate possible applications of the proposed method to various machine learning and image processing tasks, such as regression, segmentation, and registration.
(a) Input (b) L1 GT (#f: 10)
Fig. 7: Edge-aware filtering results (a
(c) Box (#f: 21) =8. e=0.0016).
(a) Input
(b) Edge-aware filtering (c) Detail enhancement
Fig. 8: Our results of edge-aware filtering and detail enhancement (a=8, s=0.04, z=3).
H HI
(a) L*GT (#f: 10) (b)Box(#f:9) (с) L} GT (#f: 22) (d) Box (#f: 21)
\ W //
\
(e) | m\L'GTof(a)
(ft I VH\ Box of(b)
VH\ L'GT of(c) (h) I VH\ Box of(cl)
Fig. 9: Edge-aware filtering results (a=8, e=0.0016). a: L1 Gauss transform with eq. (9), b: using box kernel with eq. (9), c: L1 Gauss transform with eq. (8), d: box kernel with eq. (8). e-li: visualization of \ VH \ of the corresponding images.
(a) Input (b) | VH| Our (#f: 10) (c) \ VH\ Box (#f: 9) (d) \ VH\ Box (#f: 21)
Acknowledgements
We would like to thank the anonymous reviewers of this paper for their valuable comments. This work was supported in part by Grants-in-Aid for Scientific Research of Japan (15H05954 and 16K15019).
References
[1]. A. Elgammal, R. Duraiswami, and L. Davis, "Efficient kernel density estimation using the fast Gauss transform with applications to color modeling and tracking," IEEE Transactions on Pattern Analysis and Machine Intelligence (TPAM1), vol. 25, no. 11, pp. 1499-1504,2003.
[2]. S. Paris and F. Durand, "A fast approximation of the bilateral filter using a signal processing approach," in Proc. of European Conference on Computer Vision (ECCV). Springer, 2006, pp. 568-580.
[3]. L. Greengard and J. Strain, "The fast Gauss transform," SI4AI Journal on Scientific and Statistical Computing, vol. 12, no. 1, pp. 79-94, 1991.
[4]. D. Lee, A. Gray, and A. Moore, "Dual-tree fast Gauss transforms," Advances in Neural Information Processing Systems (NIPS), vol. 18, pp. 747-754, 2006.
[5]. C. Yang, R. Duraiswami, N. A. Gumerov, and L. Davis, "Improved fast Gauss transform and efficient kernel density estimation." in Proc. of International Conference on Computer Vision (ICCV), vol. 1, 2003, pp. 464^171.
[6]. A. Adams, N. Gelfand, J. Dolson, and M. Levoy, "Gaussian kd-trees for fast high-dimensional filtering," in ACM Transactions on Graphics (TOG), vol. 28, no. 3, ACM, 2009.
[7]. E. Dougherty, Digital Image Processing Methods. CRC Press, 1994.
[8]. R. Deriche, "Fast algorithms for low-level vision," IEEE Transactions on Pattern Analysis and Machine Intelligence (TPAMI), vol. 12, no. 1, pp. 78-87, 1990.
[9]. D. Lang, M. Klaas, and N. de Freitas, "Empirical testing of fast kernel density estimation algorithms," University of British Columbia, Technical Report UBC TR-2005-03, 2005.
[10]. P. Getreuer, "A survey of Gaussian convolution algorithms," Image Process. On Line, vol. 3, pp. 276-300,2013.
[11]. S. Yoshizawa and H. Yokota, "Fast L1 Gaussian convolution via domain splitting," in Proc. of IEEE International Conference on Image Processing (ICIP). IEEE, 2014, pp. 2908-2912.
[12]. He, Kaiming, Jian Sun, and Xiaoou Tang. "Guided image filtering", IEEE Transactions on Pattern Analysis and Machine Intelligence (TPAMI), vol. 35, no. 6, pp. 1397-1409, 2013.
[13]. T. Ooura, General Purpose FFT (Fast Fourier/Cosine/Sine Transform) Package. www.kurims.kyoto-u.ac.jp/~ooura/fft.html, 2006.
[14]. Kopf, Johannes, et al. "Joint bilateral upsampling." ACM Transactions on Graphics (TOG), vol. 26. no. 3. ACM, 2007.
[15]. Tikhonov, Andrey. "Solution of incorrectly formulated problems and the regularization method." SovietMeth. Dokl. 1965, Vol. 163, No. 3.
[16]. D. Bashkirova, S. Yoshizawa, R. Latypov and H. Yokota. "Fast LI Gauss 2D Image Transforms", in Proc. of Spring/Summer Young Researchers' Colloquium on Software Engineering (SYRCoSE), Institute for System Programming, RAS, pp. 145-149, 2017. Available at http://syrcose.ispras.ru/2017/SYRCoSE2017_Proceedings.pdf, accessed June 10, 2017.
Быстрое L1-преобразование Гаусса для сглаживания изображений с сохранением границ
1,2Дина Башкирова <dina.bashkirova@riken.jp>
1Шин Йошидзава <shin@riken.jp> 2Рустам Латыпов <roustamJatypov@kpfu.ru> 1Хидео Йокота <hyokota@riken.jp> 1 Image Processing Research Team, RIKEN Center for Advanced Photonics, RIKEN 2-1, Hirosawa, Wako, Saitama, 351-0198, Japan 2Институт вычислительной математики и информационных технологий, Казанский (Приволжский) Федеральный Университет, 420008 Россия, г. Казань, Кремлевская 35
Аннотация. Преобразование Гаусса, также как и его дискретный аналог, является важнейшим инструментом во множестве математических дисциплин и находит свое применение во многих научных и инженерных областях, таких как математическая статистика и теория вероятностей, физика, математическое моделирование, машинное обучение и обработка изображений и прочие. Ввиду высокой вычислительной сложности преобразования Гаусса (квадратичная сложность относительно количества точек и экспоненциальная — относительно размерности точек), необходимы эффективные и быстрые методы его аппроксимации, обладающие большей точностью по сравнению с существующими ныне методами, такими как Быстрое Преобразование Фурье или оконное преобразование. В данной статье предложен новый метод аппроксимации преобразования Гаусса для равномерно распределенный множеств точек (например, двумерных изображений), основанный на использовании L2 метрики и метода разделения доменов. Такой подход позволяет значительно сократить количество вычислительных операций путем выполнения предварительных вычислений, и снизить вычислительную сложность метода до линейной. Результаты ряда численных экспериментов показали, что разработанный алгоритм позволяет получить более высокую точность аппроксимации без потери скорости вычисления в сравнении со стандартными методами. Также в качестве примера применения предлагаемого алгоритма была разработана новая схема смежной фильтрации изображения. Было показано, что новый фильтр на основе быстрого L1 преобразования Гаусса позволяет получить результат более высокого качества при сопоставимой скорости вычисления и при этом избежать появления нежелательных артефактов в результате обработки, таких как эффект ореола.
Ключевые слова: фильтр Гаусса, распреледение Лапласа, быстрый метод аппроксимации
DOI: 10.15514/ISPRAS-2017-29(4)-4
Для цитирования: Башкирова Д.Р., Йошидзава Ш., Латыпов Р.Х., Йокота X. Быстрое Ll-преобразование Гаусса для сглаживания изображений с сохранением границ. Труды ИСП РАН, том 29, вып. 4, 2017 г., стр. 55-72 (на английском). DOI: 10.15514/ISPRAS-2017-29(4)-4
Список литературы
[1]. A. Elgammal, R. Duraiswami, and L. Davis, "Efficient kernel density estimation using the fast Gauss transform with applications to color modeling and tracking," IEEE Transactions on Pattern Analysis and Machine Intelligence (TPAhO), vol. 25, no. 11, pp. 1499-1504,2003.
[2]. S. Paris and F. Durand, "A fast approximation of the bilateral filter using a signal processing approach," in Proc. of European Conference on Computer Vision (ECCV). Springer, 2006, pp. 568-580.
[3]. L. Greengard and J. Strain, "The fast Gauss transform," SI4A1 Journal on Scientific and Statistical Computing, vol. 12, no. 1, pp. 79-94, 1991.
[4]. D. Lee, A. Gray, and A. Moore, "Dual-tree fast Gauss transforms," Advances in Neural Information Processing Systems (NIPS), vol. 18, pp. 747-754, 2006.
[5]. C. Yang, R. Duraiswami, N. A. Gumerov, and L. Davis, "Improved fast Gauss transform and efficient kernel density estimation." in Proc. of International Conference on Computer Vision (ICCV), vol. 1, 2003, pp. 464^171.
[6]. A. Adams, N. Gelfand, J. Dolson, and M. Levoy, "Gaussian kd-trees for fast high-dimensional filtering," in ACM Transactions on Graphics (TOG), vol. 28, no. 3, ACM, 2009.
[7]. E. Dougherty, Digital Image Processing Methods. CRC Press, 1994.
[8]. R. Deriche, "Fast algorithms for low-level vision," IEEE Transactions on Pattern Analysis and Machine Intelligence (TP,4A1I), vol. 12, no. 1, pp. 78-87, 1990.
[9]. D. Lang, M. Klaas, and N. de Freitas, "Empirical testing of fast kernel density estimation algorithms," University of British Columbia, Technical Report UBC TR-2005-03, 2005.
[10]. P. Getreuer, "A survey of Gaussian convolution algorithms,"/;«<7ge Process. On Line, vol. 3, pp. 276-300,2013.
[11]. S. Yoshizawa and H. Yokota, "Fast L1 Gaussian convolution via domain splitting," in Proc. of IEEE International Conference on Image Processing (ICIP). IEEE, 2014, pp. 2908-2912.
[12]. He, Kaiming, Jian Sun, and Xiaoou Tang. "Guided image filtering." IEEE Transactions on Pattern Analysis and Machine Intelligence (TPAMI), vol. 35, no. 6, pp. 1397-1409, 2013.
[13]. T. Ooura, General Purpose FFT (Fast Fourier/Cosine/Sine Transform) Package. www.kurims.kyoto-u.ac.jp/~ooura/fft.html, 2006.
[14]. Kopf, Johannes, et al. "Joint bilateral upsampling." ACM Transactions on Graphics (TOG), vol. 26. no. 3. ACM, 2007.
[15]. Тихонов, A. H. "О некорректных задачах линейной алгебры и устойчивом методе их решения. "ДАН СССР, 1965, 163.3.
[16]. D. Bashkirova, S. Yoshizawa, R. Latypov and H. Yokota. "Fast LI Gauss 2D Image Transforms", in Proc. of Spring/Summer Young Researchers' Colloquium on Software Engineering (SYRCoSE), Institute for System Programming, RAS, 2017, pp. 145-149. Доступно по ссылке http://syrcose.ispras.ru/2017/SYRCoSE2017_Proceedings.pdf, дата обращения 10.06.2017.
