On a fast bilateral filtering formulation using functional rearrangements
aa r X i v : . [ c s . C V ] M a y On a fast bilateral filtering formulation using functionalrearrangements ∗ Gonzalo Galiano † Juli´an Velasco † Abstract
We introduce an exact reformulation of a broad class of neighborhood filters, amongwhich the bilateral filters, in terms of two functional rearrangements: the decreasing and therelative rearrangements.Independently of the image spatial dimension (one-dimensional signal, image, volume ofimages, etc.), we reformulate these filters as integral operators defined in a one-dimensionalspace corresponding to the level sets measures.We prove the equivalence between the usual pixel-based version and the rearranged ver-sion of the filter. When restricted to the discrete setting, our reformulation of bilateral filtersextends previous results for the so-called fast bilateral filtering. We, in addition, prove thatthe solution of the discrete setting, understood as constant-wise interpolators, converges tothe solution of the continuous setting.Finally, we numerically illustrate computational aspects concerning quality approxima-tion and execution time provided by the rearranged formulation.
Keywords:
Neighborhood filters, bilateral filter, decreasing rearrangement, relative rear-rangement, denoising.
To appear in Journal of Mathematical Imaging and Vision, DOI 10.1007/s10851-015-0583-y
Let Ω ⊂ R d ( d ≥
1) be an open and bounded set, u ∈ L ∞ (Ω) be an intensity image, and considerthe family of filters to which we shall refer broadly as to bilateral filters , defined byF u ( x ) = 1 C ( x ) Z Ω K h ( u ( x ) − u ( y )) w ρ ( | x − y | ) u ( y ) d y , (1)where h and ρ are positive constants, and C ( x ) = Z Ω K h ( u ( x ) − u ( y )) w ρ ( | x − y | ) d y is a normalization factor.Functions K h ( ξ ) = K ( ξ/h ) and w ρ are the range kernel and the spatial kernel of the filter,respectively, making reference to their type of interaction with the image domain. A usual choicefor K is the Gaussian K ( ξ ) = exp( − ξ ), while different choices of w ρ give rise to several wellknown neighborhood filters, e.g., ∗ Supported by Spanish MCI Project MTM2013-43671-P. † Dpt. of Mathematics, Universidad de Oviedo, c/ Calvo Sotelo, 33007-Oviedo, Spain ( [email protected],[email protected] )
1. Galiano and J. Velasco • The Neighborhood filter, see [10], for w ρ ≡ • The Yaroslavsky filter [51, 52], for w ρ ( | x − y | ) ≡ χ B ρ ( x ) ( y ), the characteristic function ofa ball centered at x of radios ρ . • The SUSAN [44] or Bilateral filters [46], for w ρ ( s ) = exp( − ( s/ρ ) ).These filters have been introduced in the last decades as alternatives to local methods suchas those expressed in terms of nonlinear diffusion partial differential equations (PDE), amongwhich the pioneering approaches of Perona and Malik [33], ´Alvarez, Lions and Morel [3] andRudin, Osher and Fatemi [42] are fundamental. We refer the reader to [12] for a review andcomparison of these methods. Neighborhood filters have been mathematically analyzed from different points of view. Forinstance, Barash [8], Elad [19], Barash et al. [9], and Buades et al. [11] investigate the asymptoticrelationship between the Yaroslavsky filter and the Perona-Malik equation. Gilboa et al. [24]study certain applications of nonlocal operators to image processing. In [34], Peyr´e establishes arelationship between nonlocal filtering schemes and thresholding in adapted orthogonal basis. Ina more recent paper, Singer et al. [43] interpret the Neighborhood filter as a stochastic diffusionprocess, explaining in this way the attenuation of high frequencies in the processed images.From the computational point of view, until the reformulation given by Porikli [37], theiractual implementation was of limited use due to the high computational demand of the directspace-range discretization. Only window-sliding optimization, like that introduced by Weiss [48]to avoid redundant kernel calculations, or filter approximations, like the introduced by Paris andDurand [32], were of computational use. In [32], the space and range domains are merged intoa single domain where the bilateral filter may be expressed as a linear convolution, followedby two simple nonlinearities. This allowed the authors to derive simple down-sampling criteriawhich were the key for filtering acceleration.However, in [37], the author introduced a new exact discrete formulation of the bilateral filterfor spatial box kernel (Yarsolavsky filter) using the local histograms of the image, h x = h | B ρ ( x ) ,where B ρ ( x ) is the box of radios ρ centered at pixel x , arriving to the formulaF u ( x ) = 1 C ( x ) n X i =1 q i h x ( q i ) K h ( u ( x ) − q i ) , (2)where the range of summation is over the quantized values of the image, q , . . . , q n , instead ofover the pixel spatial range. In addition, a zig-zag pixel scanning technique was used so thatthe local histogram is actualized only in the borders of the spatial kernel box.Formula (2) is an exact formulation of the box filter where all the terms but the localhistogram may be computed separately in constant time, and it is therefore referred to as a constant time O (1) method .Unfortunately, the use of local histograms is only valid for constant-wise spatial kernels, andsubsequent applications of the new formulation to general spatial kernels is, with the exception ofpolynomial and trigonometric polynomial kernels, only approximated. Thus, in [37] polynomialapproximation was used to deal with the usual spatial Gaussian kernel. This idea was improvedin [13] by using trigonometric expansions. 2n a fast bilateral filtering formulation using functional rearrangementsIn [49], Yang et al. introduced a new O (1) method capable of handling arbitrary spatial andrange kernels, as an extension of the ideas of Durand et al [17]. They use the so-called PrincipleBilateral Filtered Image Component J k , given by, for u ( x ) = q k , J q k ( x ) = P y ∈ N ( x ) K ( q k − u ( y )) w ρ ( | x − y | ) u ( y ) P y ∈ N ( x ) K ( q k − u ( y )) w ρ ( | x − y | ) , where N ( x ) is some neighborhood of x . Then, the bilateral filter may be expressed as F u ( x ) = J u ( x ) ( x ). In practice, only a subset of the range values is considered, and the final filtered imageis produced by linear interpolation. In this situation this filter is, thus, an approximation tothe bilateral filter. The same authors have recently extended and optimized [50] the method ofParis et al. [32] by solving cost volume aggregation problems.Other approaches include that of Kass et al. [26], who used a smoothed histogram toaccelerate median filtering and proposed mode-based filters. Adams et al. [1] proposed to useGaussian KD-trees for efficient high-dimensional Gaussian filtering, integrating this method withParis et al. method [32] for fast bilateral filtering. They later proposed to use permutohedrallattice [2] for bilateral filtering, which is faster than Gaussian KD-trees for relatively lowerdimensionality. However, both Gaussian KD-trees and permutohedral lattice are relatively lessefficient when applied to intensity images. Ram et al. [41] used a smooth patches reordering toimage denoising and inpainting which, when the patches are shrunk to pixels, contain similarideas to the decreasing rearrangement approach. A review of some of these methods may befound in [28]. In [22], we heuristically introduced a denoising algorithm based in the Neighborhood filter fordenoising images produced as time-frequency representations of one-dimensional signals, forwhich the large computational effort made useless other denoising methods based on PDE’s.The main observation was that the Neighborhood filter can be computed using only the levelsets of the image, instead of the pixel space, producing a fast denoising algorithm.In [23], this idea was rendered to a rigorous mathematical formulation through the notionof decreasing rearrangement of a function. In short, the decreasing rearrangement of u , denotedby u ∗ , is defined as the inverse of the distribution function m u ( q ) = |{ x ∈ Ω : u ( x ) > q }| , seeSection 2 for the precise definition and some of its properties.Realizing that the structure of level sets of u is invariant through the Neighborhood filteroperation as well as through the decreasing rearrangement of u allowed us to rewrite (1), for thehomogeneous kernel w ρ ≡
1, in terms of the one-dimensional integral expressionNF ∗ u ( x ) = R | Ω | K h ( u ( x ) − u ∗ ( s )) u ∗ ( s ) ds R | Ω | K h ( u ( x ) − u ∗ ( s )) ds , (3)which is computed jointly for all the pixels in each level set { x : u ( x ) = q } .Perhaps, the most important consequence of using the rearrangement was, apart from thelarge dimensional reduction, the reinterpretation of the Neighborhood filter as a local algorithm.Indeed, notice that since u ∗ is decreasing, we have | u ∗ ( t ) − u ∗ ( s ) | < h if | t − s | < c ( h ) ,
3. Galiano and J. Velascofor some continuous function c . Thus, the range kernel of the pixel-based filter (1), whose supportin the pixel space may be disconnected, is transformed into a local kernel in its rearranged version(3).Thanks to this we proved, among others, the following properties for the most usual nonlineariterative variant of the Neighborhood filter: • The asymptotic behavior of NF ∗ as a shock filter of the type introduced by ´Alvarez et al.[4], combined with a contrast loss effect. • The contrast change character of the Neighborhood filter, i.e. the existence of a continuousand increasing function g : R → R such thatNF u ( x ) ≡ NF ∗ u ( x ) = g ( u ( x )) . In this article we extend the use of functional rearranging techniques, as introduced in [23], tobilateral filters which, in the discrete case and for the spatial box window, coincides with thebilateral filter reformulation given by Porikli [37].We provide a rigorous mathematical ground which allows to interpret the discrete formulasof [37], and their extension to general spatial and range kernels, as the discrete counterpart ofcontinuous pixel-space formulations.The formulas we obtain for general bilateral filters allows to reinterpret these filters fromthe range space point of view, as a natural extension to Porikli’s discrete model. In addition,these formulas may be computationally competitive when a high degree of approximation to thepixel-based formulation is required.This work may be considered as a non trivial extension of the rearranging techniques in-troduced in [23] for the simpler Neighborhood filter. Thus, even if the spatial kernel w ρ isnon-homogeneous, we may still use our approach by introducing the relative rearrangement ofthe spatial kernel with respect to the image, see Section 2 for definitions.Using this tool, we may express the general bilateral filter (1) in terms of the one-dimensionalintegral expression F ∗ u ( x ) = R | Ω | K h ( u ( x ) − u ∗ ( s )) w ρ ( | x − ·| ) ∗ u ( s ) u ∗ ( s ) ds R | Ω | K h ( u ( x ) − u ∗ ( s )) w ρ ( | x − ·| ) ∗ u ( s ) ds , (4)where v ∗ u denotes the relative rearrangement of v with respect to u .To gain some insight into formula (4), let us consider a constant-wise interpolation of agiven image, u , quantized in n levels labeled by q i , with max( u ) = q > . . . > q n = 0. That is u ( x ) = P ni =1 q i χ E i ( x ) , where E i are the level sets of u , E i = { x ∈ Ω : u ( x ) = q i } , i = 1 , . . . , n. Similarly, let w ρ be a constant-wise interpolation of the spatial kernel quantized in m levels, r j ,with max( w ρ ) = r > . . . > r m = min( w ρ ) ≥
0. For each x ∈ Ω, consider the partition of E i given by F ij ( x ) = { y ∈ E i : w ρ ( | x − y | ) = r j } . Then, we show in Theorems 1 and 3 that for each x ∈ E k , k = 1 , . . . , n, F u ( x ) ≡ F ∗ u ( x ) = P ni =1 K h ( q k − q i ) W im ( x ) q i P ni =1 K h ( q k − q i ) W im ( x ) , (5)4n a fast bilateral filtering formulation using functional rearrangementswhere W im ( x ) = P mj =1 r j | F ij ( x ) | , with | · | denoting set measure (number of pixels). That is, | F ij ( x ) | is the number of pixels in the q i -level set of u that belongs to the r j -level set of w ρ .Observe that, like for the Neighborhood filter, the range kernel is transformed into a localkernel, independent of the pixel location, and it is therefore computed once and for all, explainingthe large gain in computational cost, as already observed in [37].In fact, notice that in the particular case in which w ρ is taken as a box window we have that W im ( x ) is just the local histogram, h x , used in Porikli’s formula, so (2) and (4) coincide in thiscase. However, the rearranged formula (4) is valid for any other shape of the spatial kernel.The results in this article may be summarized as follows:1. The rearranged formulation (4) is exact and has no restrictions on the spatial or rangekernels shape (Theorem 1). In fact, our results are not limited to the special functional form ofthe spatial kernel, w ρ ≡ w ρ ( | x − y | ), that we considered along this article. One can easily modifyour arguments to deal with more general kernels of the type w ρ ≡ w ρ ( x , y ), thus including otherwell known filters like the median filter or the cross-joint bilateral filter [18, 35]. We preferredto limit ourselves to bilateral type filters for the sake of clarity, and leave the more general caseto future developments of the work.2. When restricted to a spatial box kernel, the relative rearrangement term of formula (4) isjust the local histogram, like in Porikli’s discrete formulation. For more general spatial kernels,the relative rearrangement may be interpreted as an averaged local histogram evaluated on thekernel level sets, extending in this way the notion of local histogram and the discrete formulationof [37].3. We establish the relationship between the discrete and the continuous image setting of therearranged version of the bilateral filter, proving the convergence of constantwise interpolatorsto their continuous counterpart (Theorems 2, 3). The plan of the article is the following. In Section 2 we recall the notions of decreasing andrelative rearrangements and motivate the deduction of the rearranged formula (4) under somerestrictive regularity assumptions.In Section 3, we remove those assumptions and establish the equivalence between the usualpixel-based expression of the filter (1) and its rearranged formulation (4). In addition, we providea fully discrete algorithm to approximate by constant-wise functions the filter F ∗ u ( x ) given by(4), and thus the original equivalent filter F u ( x ) given by (1). We also prove the convergenceof this discretization to the solution of the continuous setting. In Section 4, we give the proofsof these results.In Section 5 we illustrate the performance of the rearranged filter formulation in comparisonto the brute force pixel-based implementation and to other state of the art filters: the O (1)method of Yang et al. [49], and the permutohedral lattice method of Adams et al. [1].Finally, in Section 6 we give the Summary. 5. Galiano and J. Velasco Let us denote by | E | the Lebesgue measure of any measurable set E . For a Lebesgue measurablefunction u : Ω → R , the function q ∈ R → m u ( q ) = |{ x ∈ Ω : u ( x ) > q }| is called the distributionfunction corresponding to u .Function m u is non-increasing and therefore admits a unique generalized inverse, called the decreasing rearrangement . This inverse takes the usual pointwise meaning when the function u has no flat regions, i.e. when |{ x ∈ Ω : u ( x ) = q }| = 0 for any q ∈ R . In general, the decreasingrearrangement u ∗ : [0 , | Ω | ] → R is given by: u ∗ ( s ) = ess sup { u ( x ) : x ∈ Ω } if s = 0 , inf { q ∈ R : m u ( q ) ≤ s } if s ∈ (0 , | Ω | ) , ess inf { u ( x ) : x ∈ Ω } if s = | Ω | . We shall also use the notation Ω ∗ = (0 , | Ω | ). Notice that since u ∗ is non-increasing in ¯Ω ∗ , it iscontinuous but at most a countable subset of ¯Ω ∗ . In particular, it is then right-continuous forall t ∈ (0 , | Ω | ].The notion of rearrangement of a function is classical and was introduced by Hardy, Little-wood and Polya [25]. Applications include the study of isoperimetric and variational inequalities[36, 7, 30], comparison of solutions of partial differential equations [45, 5, 47, 14, 6], and others.We refer the reader to the textbook [27] for the basic definitions.Two of the most remarkable properties of the decreasing rearrangement are the equi-measurabilityproperty Z Ω f ( u ( y )) d y = Z | Ω | f ( u ∗ ( s )) ds, for any Borel function f : R → R + , and the contractivity k u ∗ − v ∗ k L p (Ω ∗ ) ≤ k u − v k L p (Ω) , (6)for u, v ∈ L p (Ω), p ∈ [1 , ∞ ]. Apart from the pure mathematical interest, the reformulation of neighborhood filters in termsof functional rearrangements is useful for computational purposes, specially when the spatialkernel, w , is homogeneous, that is w ρ ≡
1. In this case, it may be proven [23] that the level setsof u are invariant through the filter, i.e. u ( x ) = u ( y ) implies F ( u )( x ) = F ( u )( y ), and thus, itis sufficient to compute the filter only for each (quantized) level set, instead of for each pixel,meaning a huge gain of computational effort.For non-homogeneous kernels the advantages of the filter rearranged version are kernel-dependent, and in any case, the gain is never comparable to that of homogeneous kernels. Themain reason is that the non-homogeneity of the spatial kernel breaks, in general, the invarianceof level sets through the filter application.In the following lines, we provide a heuristic derivation of the Bilateral filter rearrangedversion, as first noticed in [23]. We shall prove in Theorem 1 that the resulting formula is validin a very general setting. 6n a fast bilateral filtering formulation using functional rearrangementsUnder suitable regularity assumptions, the coarea formula states Z Ω g ( y ) |∇ u ( y ) | d y = Z ∞−∞ Z u = t g ( y ) d Γ( y ) dt. Taking g ( y ) = K h ( u ( x ) − u ( y )) w ρ ( | x − y | ) u ( y ) / |∇ u ( y ) | , and using u ( x ) ∈ [0 , Q ] for all x ∈ Ωwe get, for the numerator of the filter (1) , I ( x ) := Z Ω K h ( u ( x ) − u ( y )) w ρ ( | x − y | ) u ( y ) d y = Z Q K h ( u ( x ) − t ) t Z u = t w ρ ( | x − y | ) |∇ u ( y ) | d Γ( y ) dt. Assuming that u ∗ is strictly decreasing and introducing the change of variable t = u ∗ ( s ) we find I ( x ) = − Z | Ω | K h ( u ( x ) − u ∗ ( s )) u ∗ ( s ) du ∗ ( s ) ds × Z u = u ∗ ( s ) w ρ ( | x − y | ) |∇ u ( y ) | d Γ( y ) ds = Z | Ω | K h ( u ( x ) − u ∗ ( s )) u ∗ ( s ) w ρ ( | x − ·| ) ∗ u ( s ) ds. (7)Here, the notation v ∗ u stands for the relative rearrangement of v with respect to u which, underregularity conditions, may be expressed as v ∗ u ( s ) = Z u = u ∗ ( s ) v ( y ) |∇ u ( y ) | d Γ( y ) Z u = u ∗ ( s ) |∇ u ( y ) | d Γ( y ) , (8)see the next subsection for details. Transforming the denominator of the filter, C ( x ), in a similarway allows us to deduceF u ( x ) = R | Ω | K h ( u ( x ) − u ∗ ( s )) w ρ ( | x − ·| ) ∗ u ( s ) u ∗ ( s ) ds R | Ω | K h ( u ( x ) − u ∗ ( s )) w ρ ( | x − ·| ) ∗ u ( s ) ds . (9) The relative rearrangement was introduced by Mossino and Temam [29] as the directional deriva-tive of the decreasing rearrangement. Thus, if for u ∈ L (Ω) and v ∈ L p (Ω), with p ∈ [1 , ∞ ], weconsider the function ϕ ( s ) = Z u>u ∗ ( s ) v ( x ) d x + Z s −| u>u ∗ ( s ) | (cid:0) v | u = u ∗ ( s ) (cid:1) ∗ ( σ ) dσ, then the relative rearrangement of v with respect to u , v ∗ u , is defined as v ∗ u := dds ϕ ∈ L p (Ω ∗ ) .
7. Galiano and J. VelascoThis identity may also be understood as the weak L p (Ω ∗ ) directional derivative (weak* L ∞ (Ω ∗ ),if p = ∞ ), v ∗ u = lim t → ( u + tv ) ∗ − u ∗ t . (10)Under the additional assumptions u ∈ W , (Ω) and |{ y ∈ Ω : ∇ u ( y ) = 0 }| = 0, i.e. the non-existence of flat regions of u , the identity (8) is well defined and coincides with (10). In this case,the relative rearrangement represents an averaging procedure of the values of v on the level setsof u labeled by the superlevel sets measures, s .When formula (8) does not apply, that is, in flat regions of u , we may resort to identity (10)to interpret the relative rearrangement as the decreasing rearrangement of v restricted to suchflat regions of u .After the seminal work of Mossino and Temam [29], the relative rearrangement was furtherstudied by Mossino and Rakotoson [31] and applied to several types of problems, among whichthose related to variable exponent spaces and functional properties [20, 21, 39, 40], or nonlocalformulations of plasma physics problems [15, 16]. In the rest of the article, we shall make anextensive use of the results given in the monograph on the relative rearrangement by Rakotoson[38]. Consider the constant-wise functions u, v : [0 , → R given in Fig. 1 (a). Writing max( u ) =5 = q > . . . > q = 0 = min( u ), we may express u as u ( x ) = P i =1 q i χ E i ( x ), where E i arethe level sets of u , E = (10 , E = (8 , u isconstant-wise too, and given by u ∗ ( s ) = n X i =1 q i χ I i ( s ) , with I i = [ a i − , a i ) for i = 1 , . . . ,
6, and with a = 0, a = | E | = 1, a = | E | + | E | = 3, . . . , a = P i =1 | E i | = | Ω | = 13. The corresponding plot is shown in Fig. 1 (b).In Fig. 1 (c) we show the graphs { E i , v ( E i ) } i =1 transported as it was done in the step beforeto construct u ∗ . For instance, the highest level set of u , E = (10 , , E = (8 ,
10] to (1 , { E , v ( E ) } = { (10 , , { }} is transported to { [0 , , { }} ; { E , v ( E ) } = { (8 , , { }} is transported to { [1 , , { }} , etc.Finally, to obtain the decreasing rearrangement of v with respect to u , v ∗ u , we rearrangedecreasingly the restriction of v to the level sets of u , E i , as shown in Fig. 1 (d). In this section we provide analytical results showing that, given an image, u , a spatial kernel, w ρ , and a range kernel, K h , all bounded in L ∞ , we may approximate the filtered image v = F u by a sequence of constant-wise functions, v n , which are obtained by filtering a constant-wiseapproximation of u , u n , through a discrete filter, F m , involving a constant-wise approximationof w ρ , w ρ m .In addition, and fundamentally, we show that the filter and its approximations may beobtained using the rearranged version (9), and give explicit formulas for their computation.8n a fast bilateral filtering formulation using functional rearrangements (a) Functions u (blue) and v (red) (b) Decreasing rearrangement of u , u ∗ (c) u ∗ (blue) and transported v (red) (d) The relative rearrangement v ∗ u Figure 1:
Example of construction of the relative rearrangement. (c) shows u ∗ and v transported as bythe displacement of the level sets of u in the construction of u ∗ . (d) shows the decreasing rearrangementof v restricted to the level sets of u , that is v ∗ u . The main assumption we implicitly made for the heuristic deduction of formula (9) is thecondition |{ y ∈ Ω : ∇ u ( y ) = 0 }| = 0, i.e. the non-existence of flat regions of u , which gives senseon one hand to formula (8), and on the other hand, allows us to obtain the strictly decreasingbehavior of u ∗ , which justifies the change of variable in (7).Our first result is that formulas (1) and (9) for the pixel-based filter and its rearrangedversion are equivalent under weaker hypothesis. Due to the nature of our application, we keepthe assumption on the boundedness of u and w in L ∞ , although these can be also weakened.Theorems 1 and 2 are proven in Section 4. We use the notation R + = { t ∈ R : t ≥ } . Theorem 1
Let Ω ⊂ R d be an open and bounded set, d ≥ , K ∈ L ∞ ( R , R + ) and w ρ ∈ L ∞ ( R + , R + ) . Assume that u ∈ L ∞ (Ω) is, without loss of generality, non-negative. Consider F u ( x ) given by (1) and F ∗ u ( x ) = R | Ω | K h ( u ( x ) − u ∗ ( s )) w ρ ( | x − ·| ) ∗ u ( s ) u ∗ ( s ) ds R | Ω | K h ( u ( x ) − u ∗ ( s )) w ρ ( | x − ·| ) ∗ u ( s ) ds . (11)9. Galiano and J. Velasco Then F ∗ u ( x ) = F u ( x ) for a.e. x ∈ Ω . The next step towards the definition of a fully discrete algorithm to approximate the filterF ∗ given by (11), and thus the original equivalent filter F given by (1), is showing that we mayapproximate u and w ρ by constant-wise functions u n and w ρ m such that, as m, n → ∞ ,F m ∗ u n ( x ) → F ∗ u ( x )where F m ∗ is the discrete version of F ∗ , F m ∗ u ( x ) = R | Ω | K h ( u ( x ) − u ∗ ( s )) w ρ m ( | x − ·| ) ∗ u ( s ) u ∗ ( s ) ds R | Ω | K h ( u ( x ) − u ∗ ( s )) w ρ m ( | x − ·| ) ∗ u ( s ) ds . (12) Theorem 2
Assume the conditions of Theorem 1 and suppose that u has a finite number of flatregions. Then, there exist sequences of constant-wise functions u n , w ρ m , with a finite numberof flat regions, such that u n → u strongly in L ∞ (Ω) , w ρ m → w ρ strongly in L ∞ ( R + ) , and F m ∗ u n → F ∗ u a.e. in Ω and strongly in L ∞ (Ω) . (13) Remark 1
Theorem 2 may be extended to the case in which u has a countable number of flatregions. However, in this case the approximation sequence of constant-wise functions u n has alsoa countable number of levels. Since our aim is providing a finite discretization for numericalimplementation, such a sequence is not appropriate. Theorem 2 gives utility to Theorem 3, in which we produce the fully discrete formula actuallyused for the numerical implementation.The main difficulty in computing formula (12) is the determination of the relative rear-rangement. However, in the case of constant-wise functions with a finite number of levels thiscomputation is simplified thanks to identity (10), which may be easily applied to this situation,as shown in [38, Th.7.3.4].In few words, for the case of constant-wise functions u and v , the relative rearrangement v ∗ u may be computed as the decreasing rearrangement of v restricted to the level sets of u . Theorem 3
Let u ∈ L ∞ (Ω) be a constant-wise function quantized in n levels labeled by q i , with max( u ) = q > . . . > q n = 0 . That is u ( x ) = P ni =1 q i χ E i ( x ) , where E i are the level sets of u , E i = { x ∈ Ω : u ( x ) = q i } , i = 1 , . . . , n. Similarly, let w ρ ∈ L ∞ ( R + ) be constant-wise and quantized in m levels, r j , with max( w ρ ) = r > . . . > r m = min( w ρ ) ≥ . For each x ∈ Ω , consider the partition of E i given by F ij ( x ) = { y ∈ E i : w ρ ( | x − y | ) = r j } . Then, for each x ∈ E k , k = 1 , . . . , nF ∗ u ( x ) = P ni =1 K h ( q k − q i ) W im ( x ) q i P ni =1 K h ( q k − q i ) W im ( x ) , (14) with W im ( x ) = P mj =1 r j | F ij ( x ) | .
10n a fast bilateral filtering formulation using functional rearrangements
As it clear from formula (14), the main difficulty for its computation is the determination of themeasures of the sets F ij ( x ), which must be computed for each x ∈ Ω.The formula provides the algorithm complexity. If N is the number of pixels, then thecomplexity (without any further code optimization) is of the order O ( N nm ), where n is thenumber of levels of the image and m is the number of levels of the spatial kernel. Let usexamine some examples. The Neighborhood filter.
In this case, w ρ ≡
1, and therefore j = 1 and F ij ( x ) = E i is independentof x for all i = 1 , . . . , n . Thus, formula (14) is computed only on the level sets of u , that is, forall x ∈ E k F ∗ h u ( x ) = P ni =1 K h ( q k − q i ) | E i | q i P ni =1 K h ( q k − q i ) | E i | . In this case, the complexity is of order O ( n ). The Yaroslavsky filter.
In this case, w ρ ( | x − y | ) = χ B ρ ( x ) ( y ), and therefore there are only twolevels r = 1 , r = 0 of w ρ corresponding to the sets F i ( x ) = { y ∈ E i : | x − y ) | < ρ } ,F i ( x ) = { y ∈ E i : | x − y ) | ≥ ρ } , Thus, formula (14) reduces to: for each x ∈ E k , k = 1 , . . . , nF ∗ h u ( x ) = P ni =1 K h ( q k − q i ) | F i ( x ) | q i P ni =1 K h ( q k − q i ) | F i ( x ) | The complexity is then of order O ( N n ). The Bilateral filter.
In this case, w ρ ( | x − y | ) = exp( −| x − y | /ρ ) and therefore there is acontinuous range of levels of w ρ . However, for computational purposes the range of w ρ isquantized to some finite number of levels, determined by the size of ρ . Thus, the full formula(14) must be used in this case. The resulting complexity is of order O ( N nm ).The corresponding complexities for the pixel based formulation is, at least, of the order O ( N nµ ), where µ denote the discrete size of the spatial kernel. However, meanwhile for therearranged version the range kernel is computed only once, for the pixel-based version it mustbe computed for each pixel.Observe that more optimized code and a further reduction of algorithm complexity may bereached by standard zig-zag techniques of window scanning, see [37]. Proof of Theorem 1.
Let f ∈ L ∞ ( R ) and b ∈ L ∞ (Ω). We start showing Z Ω f ( u ( y )) b ( y ) d y = Z | Ω | f ( u ∗ ( s )) b ∗ u ( s ) ds. (15)Consider the sets of flat regions of u and u ∗ , P = [ i ∈ D P i , P i = { y ∈ Ω : u ( y ) = q i } , (16)11. Galiano and J. Velascoand P ∗ = ∪ i ∈ D P ∗ i , with , P ∗ i = { s ∈ Ω ∗ : u ∗ ( s ) = q i } , where the subindices set D is, at most,countable. According to [38, Lemma 2.5.2], we have Z | Ω | f ( u ∗ ( s )) b ∗ u ( s ) = Z Ω \ P f ( u ∗ ( m u ( u ( y )))) b ( y ) d y (17)+ X i ∈ D Z P i M b i ( h i )( y ) b ( y ) d y , (18)where b i = b | P i , h i ( s ) = f ( u ∗ ( s ′ i + s )) for s ∈ [ s ′ i , s ′′ i ) := P ∗ i , and M b i ( h i )( y ) = h i ( m b i ( b i ( y ))) if y ∈ P i \ Q i , | Q j | Z σ ′′ j σ ′ j h i ( s ) ds, if y ∈ Q ij , where Q i = ∪ j ∈ D ′ Q ij , with Q ij the flat regions of b i and [ σ ′ j , σ ′′ j ) := Q j ∗ i . In (17), since thefunctions u ∗ and m u are strictly decreasing and inverse of each other in the set Ω \ P , we obtain Z Ω \ P f ( u ∗ ( m u ( u ( y )))) b ( y ) d y = Z Ω \ P f ( u ( y )) b ( y ) d y . In the flat regions of u and u ∗ we have, on one hand, Z P f ( u ( y )) b ( y ) d y = X i ∈ D Z P i f ( u ( y )) b ( y ) d y = X i ∈ D f ( q i ) Z P i b ( y ) d y . And, on the other hand, since h i ( s ) = f ( u ∗ ( s ′ i + s )) = f ( q i ) for s ∈ P ∗ i , X i ∈ D Z P i M b i ( h i )( y ) b ( y ) d y = X i ∈ D Z P i \ Q i h i ( m b i ( b i ( y ))) b ( y ) d y + X j ∈ D ′ | Q j | Z σ ′′ j σ ′ j h i ( s ) ds Z Q ij b ( y ) d y = X i ∈ D f ( q i ) Z P \ Q i b ( y ) d y + f ( q i ) X j ∈ D ′ Z Q ij b ( y ) d y = X i ∈ D f ( q i ) Z P i b ( y ) d y . Therefore, both sides of (15) are equal.Finally, for fixed x ∈ Ω set b ( y ) = w ρ ( x , y ) and first, f ( t ) = K h ( u ( x ) − t ) t , for t ≥ f ( t ) = K h ( u ( x ) − t ) to obtain the equality between the denominators of those expressions. (cid:3)
12n a fast bilateral filtering formulation using functional rearrangements
Remark 2
As deduced in the proof, identity (15) follows from [38, Lemma 2.5.2]. In fact, alittle more than (15) may be obtained. Let f, b ∈ L ∞ (Ω) . Then, if f is constant in the flatregions of u , that is f | P i = f i = const , then Z | Ω | f ( s ) b ∗ u ( s ) ds = Z Ω f ( m u ( u ( y ))) b ( y ) d y + X i ∈ D f ( q i ) Z P i b ( y ) d y . (19) Proof of Theorem 2.
We split the proof in several steps.
Step 1.
We use the construction of the sequence of constant-wise functions u n given in [38, Th.7.2.1]. In our case, the construction is simpler because u has a finite number of flat regions,implying that u n has a finite number of levels.In any case, this construction is such that u n → u a.e. in Ω and strongly in L ∞ (Ω),and w ∗ u n → w ∗ u weakly* in L ∞ (Ω ∗ ). Besides, due to the strong continuity of the decreasingrearrangement (6), we also have ( u n ) ∗ → u ∗ strongly in L ∞ (Ω ∗ ). Therefore, we readily see firstthat F ∗ h u n ( x ) → F ∗ h u ( x ) for a.e. x ∈ Ω , and then, due to the dominated convergence theorem, F ∗ h u n → F ∗ h u strongly in L ∞ (Ω) . (20) Step 2.
We consider a sequence of constant-wise functions with a finite number of levels, w m ,such that w m → w strongly in L ∞ (Ω × Ω). Due to the contractivity property of the relativerearrangement, see [29], we also have, for a.e. x ∈ Ω, w m ( x , · ) ∗ v → w ρ ( x , · ) ∗ v strongly in L ∞ (Ω),for any v ∈ L ∞ (Ω). Thus, as m → ∞ , F ∗ h,m u n ( x ) = F ∗ h u n ( x ) for a.e. x ∈ Ω , and, again, the dominated convergence theorem implies F ∗ h,m u n → F ∗ h u n strongly in L ∞ (Ω) . (21) Step 3.
In view of (20) and (21), we have | F ∗ h,m u n − F ∗ h u | ≤| F ∗ h,m u n − F ∗ h u n | + | F ∗ h u n − F ∗ h u | → m → ∞ and n → ∞ , so (13) follows. (cid:3) Proof of Theorem 3.
Since u is constant-wise, the decreasing rearrangement of u is constant-wise too, and given by u ∗ ( s ) = n X i =1 q i χ I i ( s ) , with I i = [ a i − , a i ) for i = 1 , . . . , n , and a = 0, a = | E | , a = | E | + | E | , . . . , a n = P ni =1 | E i | = | Ω | . It is convenient to introduce here the cumulative sum of sets measurescum( E ◦ ,
0) = 0 , and cum( E ◦ , i ) = i X k =1 | E k | ,
13. Galiano and J. Velascofor i = 1 , ..., n , where the symbol ◦ denotes the summation variable. Thus, a i = cum( E ◦ , i ).We consider, for fixed x ∈ Ω and t >
0, the function H ( y ) := u ( y ) + tw ρ ( | x − y | ). Sinceboth u and w ρ are constant-wise with a finite number of levels, we have, for t small enough q i +1 < q i + tr j < q i − , implying that each level set of H is included in one and only one level set of u . Thus, we havethe orderings • For each j , if i > i and y ∈ F i j ( x ) , ¯ y ∈ F i j ( x ) then H ( y ) < H (¯ y ). • For each i , if j > j and y ∈ F ij ( x ) , ¯ y ∈ F ij ( x ) then H ( y ) < H (¯ y ).With these observations we may compute the decreasing rearrangement of H as follows. Forinstance, for s ∈ I = [0 , | E | ) we have (omitting x from the notation F ji ( x )) H ∗ ( s ) = q + tr if s ∈ [0 , | F | ) ,q + tr if s ∈ [ | F | , | F | + | F | ) ,. . . . . .q + tr m if s ∈ [cum( F ◦ , m − , cum( F ◦ , m )) . Notice that [0 , | F | ) = [cum( F ◦ , , cum( F ◦ , , [ | F | , | F | + | F | ) = [cum( F ◦ , , cum( F ◦ , ,. . . [cum( F ◦ , m − , cum( F ◦ , m )) = [cum( F ◦ , m − , | E | ) , so, in general, we may write for i = 1 , . . . , n and j = 1 , . . . , m , H ∗ ( s ) = q i + tr j if s ∈ J ij ( x ) , where J ij ( x ) := h b ij − ( x ) , b ij ( x ) (cid:17) , with b ij ( x ) := a i − + cum( F i ◦ ( x ) , j ). Observe that b im ( x ) = a i .Finally, since J ij ( x ) ⊂ E i we have for s ∈ J ij ( x ) H ∗ ( s ) − u ∗ ( s ) t = r j , implying w ρ ( | x − ·| ) ∗ u ( s ) = r j . We are now in disposition to compute formula (11). For x ∈ E k , Z | Ω | K h ( u ( x ) − u ∗ ( s )) u ∗ ( s ) w ρ ( | x − ·| ) ∗ u ( s ) ds = n X i =1 m X j =1 K h ( q k − q i ) q i r j | J ij ( x ) | In a similar way we obtain C ( x ) = Z | Ω | K h ( u ( x ) − u ∗ ( s )) w ρ ( | x − ·| ) ∗ u ( s ) ds = n X i =1 m X j =1 K h ( q k − q i ) r j | J ij ( x ) | and therefore, using the definition of the sets J ij ( x ) we deduce (14). (cid:3)
14n a fast bilateral filtering formulation using functional rearrangements
We conducted several experiments on standard natural images to check the performance of thediscrete formula (14) in comparison to well known state of the art denoising algorithms: Yang’set al. Real Time O (1) Bilateral Filter [49], and the Permutohedral Lattice Fast Filtering [2]. The aim of our experiments was to investigate the quality of the approximation of thesealgorithms to1. the ground truth, and2. the exact pixel-based bilateral filter.We also compared execution times as delivered straightly from the available codes. Noticethat execution time depends on code optimization and therefore a rigorous study of this aspectrequires some kind of code normalization which was out of the scope of our study.In both experiments we used four intensity images of different sizes to compare executiontimes and peak signal to noise ratio (PSNR): Clock (256 × Boat (512 × Airport (1024 × Still life (2144 × σ ( u ) /σ ( ν ), where σ is the empiricalstandard deviation, u is the original image, and ν is the noise. We repeated the experimentsdescribed in this section for higher levels of noise. Since the results were qualitatively similar,we omitted them for the sake of brevity.We considered different spatial window sizes determined by ρ , with ρ = 4 , , ,
32. Forthe Yaroslavsky filter, the spatial kernel is the characteristic function of a box of radius ρ , whilefor the bilateral filter a Gaussian function is taken in the experiments.The range size of the filter was taken, throughout all the experiments, as h = ρ which,according to [11], is the regime in which the corresponding iterative filter behaves asymptoticallyas a Perona-Malik type filter. The shape of the range filter is always a Gaussian.The discretization of the pixel-based and the rearranged version of both filters was imple-mented in non-optimized C ++ codes in which the main differences are those corresponding toformulas (23) and (24). Execution time was measured by means of function clock . All theexperiments were run on a standard laptop (Intel Core i7-2.80 GHz processor, 8GB RAM). The performance of formula (14) for filter computation strongly depends on the functional formof the spatial kernel w ρ . In the case in which w ρ is the characteristic function of a box of radius ρ , that is F u ( x ) = 1 C ( x ) Z B ρ ( x ) K h ( u ( x ) − u ( y )) u ( y ) d y , (22) C ++ code downloaded from ∼ qiyang C ++ code downloaded from http://graphics.stanford.edu/papers/permutohedral
15. Galiano and J. Velasco
Algorithm 1
Yaroslavsky filter (rearranged version) * Algorithm according to formula (24), with q k = k * for each pair ( i, k ) of image quantization levels do K h ( i, k ) = exp ( − (( k − i ) /h ) ) (* Main loop *) for each pixel, x , of the image do W i ← if x is the first pixel then for each y ∈ B ρ ( x ) do if u ( y ) = i then W i ( x )+ = 1 else for each y ∈ B ρ ( x ) do actualize W i ( x ) adding/deleting new/old bins w.r.t. the previous box for each quantization level, i do numerator + = K h ( i, u ( x )) ∗ W i ( x ) ∗ i denominator + = K h ( i, u ( x )) ∗ W i ( x ) uF iltered ( x ) = numerator/denominator the number of level sets of the spatial kernel is just two, and the rearranged version requiresonly the computation of the measure of the sets F i ( x ) for each pixel x , that is, counting amongall pixels in the i − level set of u how many of them belong to the spatial box.In the discretization of the pixel-based filter, for each pixel x we compute u in the squareneighborhood, N ( x ), centered at x of radius 2 h , containing (4 h + 1) pixels (close to the border,the image is extended by zero). Then, we approximate (22) asF u ( x ) ≈ P y ∈ N ( x ) K h ( u ( x ) − u ( y )) u ( y ) P y ∈ N ( x ) K h ( u ( x ) − u ( y )) . (23)For the rearranged version we compute, for each pixel x ∈ E k , k = 1 , . . . , n , the number of pixelsof the spatial box which belong to the different level sets of u , W i ( x ) = ( W ( x ) , . . . , W n ( x )).That is, W i ( x ) is the number of pixels of E i which are present in the box N ( x ). Then, we takeF ∗ u ( x ) ≈ P ni =1 K h ( q k − q i ) W i ( x ) q i P ni =1 K h ( q k − q i ) W i ( x ) , (24)where q = ( q , . . . , q n ) are the quantization levels of u . Observe that since the rearranged versiontransforms the range kernel K h ( u ( x ) − u ( y )) into, K h ( q k − q ), this term may be computed andstored outside the main loop running over all the pixels, see Algorithm 1.Thus, in (24), only M r ( x ) and q k must be actualized for each pixel, establishing an importantdifference with respect to the pixel-based version (23), in which both the spatial and range kernelsare actualized. In addition, following [37], we perform the actualization of the measures M r ( x )using previously computed information of neighbor boxes (zig-zag scanning technique).16n a fast bilateral filtering formulation using functional rearrangements Algorithm 2
Bilateral filter (rearranged version) * Algorithm according to formula (26), with q k = k * for each pair ( i, k ) of image quantization levels do K h ( i, k ) = exp ( − (( k − i ) /h ) ) for each j of spatial kernel levels do r ( j ) = exp ( − ( j/ρ ) ) for each pixel y in a generic Gaussian window do if exp ( − ( y /ρ ) ) ≈ r ( j ) then P ixelLevel ( y ) ← j (* Main loop *) for each pixel, x , of the image do countBins ← for each pixel, y , of Gauss. window centered at x do countBins ( u ( y ) , P ixelLevel ( x − y ))+ = 1 for each quantization level, i do W im ( x ) ← for each spatial kernel level, j do W im ( x )+ = r ( j ) ∗ countBins ( i, j ) numerator + = K h ( i, u ( x )) ∗ W im ( x ) ∗ i denominator + = K h ( i, u ( x )) ∗ W im ( x ) uF iltered ( x ) = numerator/denominator If the spatial kernel w ρ has a more complex profile than the characteristic function, e.g. theGaussian function, then the rearranged version requires the computation of, for each pixel x ,the measure of the sets F ij ( x ) for all the level setss, 1 ≤ j ≤ m , of some discretization of w ρ . Forinstance, for ρ = 4 , , ,
32, the corresponding number of the Gaussian function quantizationlevels arising from double precision arithmetic are m = m ( ρ ), with m (4) = 42, m (8) = 135, m (16) = 457, and m (32) = 1621.The approximations to the pixel based and the rearranged version are, respectively,F u ( x ) ≈ P y ∈ N ( x ) K h ( u ( x ) − u ( y )) w ρ ( | x − y | ) u ( y ) P y ∈ N ( x ) K h ( u ( x ) − u ( y )) w ρ ( | x − y | ) , (25)and F ∗ u ( x ) ≈ P ni =1 K h ( q k − q i ) W im ( x ) q i P ni =1 K h ( q k − q i ) W im ( x ) , (26)with W im ( x ) = P mj =1 r j | F ij ( x ) | , where, for x ∈ E k , | F ij ( x ) | is the number of pixels of E i whichbelong to the j − level set of w ρ .The bilateral filter may be accelerated by manipulating the quantization levels of the image,and of the spatial range. As shown in [49] for the Yaroslavsky filter, the reduction of the imagequantization levels leads to poor denoising results. However, we checked that a similar restriction17. Galiano and J. Velascoapplied to the spatial kernel reduces the execution time while conserving good denoising quality.In our experiments we tested this ansatz with a fixed number of levels, m = 20, for all theimages and ρ − values. The algorithm introduced by Yang et al. [49] is the following. Let { ˜ q , . . . , ˜ q ˜ n } ⊂ { q , . . . , q n } , ˜ n ≤ n, be a subset of the quantization levels of the image. Then, if u ( x ) ∈ [˜ q k , ˜ q k +1 ], Yang’s filter isgiven by the image interpolation F Y u ( x ) = (˜ q k +1 − u ( x )) J ˜ q k ( x ) + ( u ( x ) − ˜ q k ) J ˜ q k +1 ( x ) , (27)where J ˜ q k ( x ) = P y ∈ N ( x ) K h (˜ q k − y ) w ρ ( | x − y | ) u ( y ) P y ∈ N ( x ) K h (˜ q k − y ) w ρ ( | x − y | ) . (28)Thus, if ˜ n = n then F Y coincides with the exact discrete version of the bilateral filter, whereasif ˜ n < n an approximation of the bilateral filter is obtained by interpolation on the quantizationrange kernel levels. Observe that since x is no longer present in the range kernel of J ˜ q k ( x ),both factors in (28) may be computed by fast approximated convolution algorithms, which isthe main strength of the algorithm provided in [49]However, in the first experiment (Yaroslavsky filter), and in order to obtain exact resultswhen ˜ n = n , we implemented Yang’s algorithm by using our rearranged version (24) applied to J ˜ q k ( x ), explaining the increment of execution times when compared to the actual implementationof Yang et al.As already mentioned, we also use for comparison purposes the fast filtering using the per-mutohedral lattice introduced by Adams et al. This algorithm is based on resampling techniquesand specially useful for high dimensional filtering. We refer to [2] for a thoroughfull explanationof the method.We use the following notation for the algorithms to be compared: • YPB: Yaroslavsky pixel-based version, formula (23). • YRR: rearranged Yaroslavsky version, formula (24). • EYang˜ n : exact Yang’s algorithm, formula (27), computed with the rearranged strategy,for ˜ n = 256, 64 or 8 interpolators. • BPB: Bilateral pixel-based version, formula (25). • BRRM: rearranged bilateral version, formula (26), with the maximum number of spatialkernel discretization levels, m (4) = 42, m (8) = 135,... • BRR20: rearranged bilateral version, formula (26), with fixed number of spatial kerneldiscretization levels, m = 20. • Y˜ n : Yang’s code provided by the authors, see footnote 1 in page 15, for ˜ n = 256 or ˜ n = 8interpolators. 18n a fast bilateral filtering formulation using functional rearrangements • P: Permutohedral Lattice bilateral filter code provided by the authors, see footnote 2 inpage 15.
In the first experiment, we tested the approximation quality of the rearranged and the exactYang’s versions of the Yaroslavsky filter to the pixel-based version. In Table 1 we show thePSNR of YRR and EYang with respect to YPB. We recall that, according to [32], PSNR valuesabove 40dB often correspond to almost invisible differences. The high PSNR values obtainedfor YRR and EYang256 show that these algorithms give practically the same results than YPB,up to rounding errors. EYang64 still gives very good approximation results, whereas thesesignificantly diminishes for EYang8.In Table 1 we also show the execution times. YRR is up to 400 times faster than YPB.The execution times of EYang compared to YRR are always higher, but this has to be takencautiously since the code implementations were not optimized. In this experiment we were moreinterested in approximation quality.In the second experiment, we tested the approximation quality of the algorithms described inthe previous subsection to the ground truth image and to the bilateral pixel-based filter (BPB).In Table 2 we show the corresponding PSNR’s. We see that all the algorithms employedgive similar results when compared to the ground truth image. Thus, if this were the choicecriterium, the faster, that is Y
8, should be considered.However, when compared to the BPB, the PSNR’s are quite different. As expected, therearranged version with the maximum number of quantization levels for the spatial kernel,BRRM, gives the same result, up to rounding errors, than the BPB. However, the Yang’salgorithm with the maximum number of levels, Y256, although should give exact results, itdoes not, revealing other sources of error beyond rounding errors. Even the BRR20 and theYaroslavsky YRR give better approximations than Y256. In particular, BRR20 has alwaysvalues of PSNR around 40dB. The Permutohedral algorithm gives similar results than Y256 andY8.Let us mention two difficulties we found with the Yang code implementation. The firstis related to low h − values, for which Y8 gives poor results due to some artifact formation.The second is related to large images, for which huge memory allocation requirements result inrunning time execution errors.In Table 3 we collect the execution times obtained in this experiment. Only for the smallerimage sizes and h − values Y8 has a competitor in YRR. For large images and high h − valuesthe second fastest algorithm is, as expected, the permutohedral, although always taking aroundtwice the time of Y8. Both BRRM and BRR20 give execution times considerably higher thanthe other algorithms, for our non-optimized codes.Finally, in Figures 2, 3 and 4 we show some of the filters outcome. The columns give, fromleft to right, the denoised image, a detail of the image, the contour plot of the detail, and thehistogram of the denoised image. The rows give, from above to below, the noisy image, theresults of BPB, BRR20, Y8, YRR and P, and in the last row, the ground truth clean image.It may be observed that, in general, the permutohedral algorithm performs a hardest denois-ing than the other algorithms. This is also reflected in the formation of picks in the histogram.The Y8 gives the softest denoising and performs some smoothing in the histogram of the denoisedimage, unlike the other algorithms. Finally, BPB and BRR20 are almost indistinguishable, while19. Galiano and J. VelascoYRR is very close to them, although with a harder denoising effect. The staircasing effect dueto the reduction of level sets is present in all the algorithms, although not specially relevant innone of them. In this paper we introduced the use of functional rearrangements to express bilateral type filtersin terms of integral operators in the one-dimensional space [0 , | Ω | ].We have proved the equivalence between the pixel-based formulation of the original filterand its rearranged version. In addition, we proved the convergence of discrete finite step-wiseapproximations of image and filter to their corresponding continuous limits. Although this isa property easy to obtain in the pixel-based formulation, it is far from trivial in its rearrangedversion.In the case in which the spatial kernel, w ρ , is homogeneous (e.g. the Neighborhood filter), itwas proven in [23] that the level set structure of the image is left invariant through the filteringprocess, allowing to compute the filter jointly for all the pixels in each level set, instead ofpixel-wise.However, if the spatial kernel is not homogeneous the invariance of the u − level sets throughthe filter is, in general, broken. Despite this fact, there still remains an important property ofthe rearranged version: the range kernel K h ( u ( x ) − u ( y )) is transformed into a pixel-independentkernel K ( q k − q i ), implying a large gain in computational effort, as already observed for particularcases in [37]. We have illustrated numerically this property.The present state of the rearranging technique has an inherent restriction: the need of arelation of order in the range space. While this may be a minor issue for vector-valued colorimages, for which a map to some one-dimensional color space may be used, it is unclear how toextend the method to patch-based algorithms such as the Nonlocal Means. Whether or not thepatch reordering techniques introduced, among others, by Ram et al. [41] may be adapted tothe rearranging approach will be the focus of future research. The authors are partially supported by the Spanish DGI Project MTM2013-43671.The authors thank to the anonymous reviewers for their interesting comments and sugges-tions, that highly contributed to the improvement of our work.
References [1] Adams A, Gelfand N, Dolson J, Levoy M (2009) Gaussian kd-trees for fast high-dimensionalfiltering. ACM Trans Graph 28:21:1–21:12.[2] Adams A, Baek J, Davis MA (2010) Fast high-dimensional filtering using the permutohedrallattice. Computer Graphics Forum,29(2):753–762[3] ´Alvarez L, Lions PL, Morel JM (1992) Image selective smoothing and edge detection bynonlinear diffusion, II. Siam J Numer Anal, 29:845–86620n a fast bilateral filtering formulation using functional rearrangements[4] Alvarez L, Mazorra L (1994) Signal and image restoration using shock filters and anisotropicdiffusion. Siam J Numer Anal 31(2):590–605[5] Alvino A, Trombetti G (1978) Sulle migliori costanti di maggiorazione per una classe diequationi ellittiche degeneri. Ricerche Mat 27:413–428[6] Alvino A, D´ıaz JI, Lions PL, Trombetti G (1996) Elliptic Equations and Steiner Sym-metrization. Comm Pure Appl Math XLIX:217–236[7] Bandle C (1980) Isoperimetric inequalities and applications. Pitman.[8] Barash D (2002) Fundamental relationship between bilateral filtering, adaptive smoothing,and the nonlinear diffusion equation. IEEE T Pattern Anal 24(6):844–847[9] Barash D, Comaniciu D (2004) A common framework for nonlinear diffusion, adaptivesmoothing, bilateral filtering and mean shift. Image Vision Comput 22(1):73–81[10] Buades A, Coll B, Morel JM (2005) A review of image denoising algorithms, with a newone. Multiscale Model Sim 4(2):490–530[11] Buades A, Coll B, Morel JM (2006) Neighborhood filters and pde’s. Numer Math 105(1):1–34[12] Buades A, Coll B, Morel JM (2010) Image denoising methods. A new nonlocal principle.Siam Rev 52(1):113–147[13] Chaudhury KN, Sage D, Unser M (2011) Fast O(1) bilateral filtering using trigonometricrange kernels. TIP 2011.[14] D´ıaz JI, Nagai T (1995) Symmetrization in a parabolic-elliptic system related to chemotaxis.Adv Math Sci Appl, 5:659–680[15] D´ıaz JI, Rakotoson JM (1996) On a nonlocal stationary free boundary problem arising inthe confinement of a plasma in a Stellarator geometry. Arch Rat Mech Anal, 134(1):53–95[16] D´ıaz JI, Padial JF, Rakotoson JM (1998) Mathematical treatement of the magnetic con-finement in a current carrying Stellerator. Nonlinear Anal. TMA 34:857–887[17] Durand F, Dorsey J (2002) Fast bilateral filtering for the display of high-dynamic-rangeimages. ACM Siggraph 21:257–266[18] Eisemann E, Durand F (2004) Flash photography enhancement via intrinsic relighting.Siggraph 23(3):673–678[19] Elad M (2002) On the origin of the bilateral filter and ways to improve it. IEEE T ImageProcess 11(10):1141–1151[20] Fiorenza A, Rakotoson JM (2007) Relative rearrangement and Lebesgue spaces L p ( · ) withvariable exponent. J Math Pures Appl 88(6):506–521[21] Fiorenza A, Rakotoson JM (2009) Relative rearrangement method for estimating dualnorms. Indiana Univ Math J 58(3):1127–115021. Galiano and J. Velasco[22] Galiano G, Velasco J (2014) On a non-local spectrogram for denoising one-dimensionalsignals. Appl Math Comput, DOI:10.1016/j.amc.2014.07.003.[23] Galiano G, Velasco J (2014) Neighborhood filters and the decreasing rearrangement. J MathImaging Vision, DOI: 10.1007/s10851-014-0522-3.[24] Gilboa G, Osher S (2008) Nonlocal operators with applications to image processing. Mul-tiscale Model Sim 7(3):1005–1028[25] Hardy GH, Littlewood JE, Polya G (1964) Inequalities. Cambridge University Press[26] Kass M, Solomon J (2010) Smoothed local histogram filters. ACM TOG 29(4):100:1–100:10[27] Lieb EH, Loss M (2001) Analysis, vol 4. American Mathematical Soc.[28] Milanfar P (2012) A tour of modern image filtering: new insights and methods, bothpractical and theoretical. IEEE Signal Process Magazine (30):106–128[29] Mossino J, Temam R (1981) Directional derivative of the increasing rearrangement mappingand application to a queer differential equation in plasma physics. Duke Math J, 48(3):475–495[30] Mossino J (1984) In´egalit´es Isop´erm´etriques et applications en physique. Hermann.[31] Mossino J, Rakotoson JM (1986) Isoperimetric inequalities in parabolic equations. Ann ScNorm Super Pisa Sci(4) 13(1):51–73[32] Paris S, Durand F (2006) A fast approximation of the bilateral filter using a signal processingapproach. ECCV 2006[33] Perona P, Malik J (1990) Scale-space and edge detection using anisotropic diffusion. Ieee TPattern Anal 12(7):629–639[34] Peyr´e G (2008) Image processing with nonlocal spectral bases. Multiscale Model Sim7(2):703–730[35] Petschnigg G, Agrawala M, Hoppe H, Szeliski R, Cohen M, Toyama K (2004) Digitalphotography with flash and no-flash image pairs. ACM Siggraph, 23(3):664–672[36] P´olya G Szeg¨o G (1951) Isoperimetric inequalities in mathematical physics. PrincentonUniversity Press[37] Porikli F (2008) Constant time O(1) bilateral filtering. CVPR 2008[38] Rakotoson JM (2008) R´earrangement Relatif: Un instrument destimations dans lesproblˇcmes aux limites, vol 64. Springer[39] Rakotoson JM (2010) Lipschitz properties in variable exponent problems via relative rear-rangement. Chin Ann Math, 31B(6):991–1006.[40] Rakotoson JM (2008) New Hardy inequalities and behaviour of linear elliptic equations. JFunct Anal, 263(9):2893–2920. 22n a fast bilateral filtering formulation using functional rearrangements[41] Ram I, Elad M, Cohen I (2013) Image processing using smooth ordering of its patches.IEEE T Image Process 22(7):2764–2774[42] Rudin LI, Osher S, Fatemi E (1992) Nonlinear total variation based noise removal algo-rithms. Physica D 60(1):259–268[43] Singer A, Shkolnisky Y, Nadler B (2009) Diffusion interpretation of nonlocal neighborhoodfilters for signal denoising. SIAM J Imaging Sci 2(1):118–139.[44] Smith SM, Brady JM (1997) Susan. a new approach to low level image processing. Int JComput Vision 23(1):45–78[45] Talenti G (1976) Best constant in Sobolev inequality. Ann Mat Pura Appli, (4)110:353–372[46] Tomasi C, Manduchi R (1998) Bilateral filtering for gray and color images. In: Sixth Inter-national Conference on Computer Vision, 1998, IEEE, pp 839–846[47] V´azquez JL (1982) Sym´etrization pour u t = ∆ φ ( u ) et applications. C R Acad Paris, 295:71–74[48] Weiss B (2006) Fast median and bilateral filtering. ACM Siggraph 25:519–526[49] Yang Q, Tan KH, Ahuja N (2009) Real-time O(1) bilateral filtering. CVPR 2009.[50] Yang Q, Ahuja N, Tan KH (2014) Constant time median and bilateral filtering. Int JComput Vis DOI 10.1007/s11263-014-0764-y[51] Yaroslavsky LP (1985) Digital picture processing. An introduction. Springer Verlag, Berlin[52] Yaroslavsky LP, Eden M (2003) Fundamentals of Digital Optics. Birkh¨auser, Boston23. Galiano and J. Velasco h YPB YRR EYang256 EYang64 EYang8 Clock (256 × Boat (512 × Airport (1024 × Still life (2144 × Table 1: Comparison between the pixel-based Yaroslavsky (YPB) filter algorithm, the rear-ranged version introduced in this article (YRR), and several instances of Yang’s exact algorithm(˜ n = 256 , , Compared to ground truth Compared to BPB
Clock (256 × Boat (512 × Airport (1024 × Still life (2144 × Table 2: Left box: PSNR between the ground truth image and the pixel-based Bilateralfilter (BPB), its rearranged version for the maximum number of meaningful levels (BRRM), therearranged version for 20 levels (BRR20), two instances of Yang’s algorithm, (Y256 and Y8), thePermutohedral algorithm (P), and the Yaroslavsky filter in its rearranged version (YRR). Rightbox: PSNR between the pixel-based Bilateral filter (BPB) and the other algorithms. “ND”means “no data available“ due to too huge Y256 memory requirements for these images. Imagesizes enclosed in parentheses. 25. Galiano and J. Velasco
Execution time Execution time relative to BPB
Clock (256 × Boat (512 × Airport (1024 × Still life (2144 × Table 3: Execution times for the images used in the experiments (sizes enclosed in parentheses).Left box: Execution times ot the pixel-based Bilateral filter (BPB), its rearranged version for themaximum number of meaningful levels (BRRM), the rearranged version for 20 levels (BRR20),two instances of Yang’s algorithm, (Y256 and Y8), the Permutohedral algorithm (P), and theYaroslavsky filter in its rearranged version (YRR). Right box: Execution times relative to BPB,that is, number of times the other algorithms are faster than the BPB algorithm. “ND” means“no data available“ due to too huge Y256 memory requirements for these images.26n a fast bilateral filtering formulation using functional rearrangements
NoisyBPB0 50 100 150 200 25005001000150020002500
NoisyBRR200 50 100 150 200 25005001000150020002500
NoisyY80 50 100 150 200 25005001000150020002500
NoisyYRR0 50 100 150 200 25005001000150020002500
NoisyP
Figure 2: Image
Clock , h = 8. The columns are, from left to right, the denoised image, a detailof the image, the contour plot of the detail, and the histogram of the denoised image. The rowsgive, from above to below, the noisy image, the results of BPB, BRR20, Y8, YRR and P, andin the last row, the ground truth clean image. 27. Galiano and J. Velasco NoisyBPB0 50 100 150 200 25001000200030004000500060007000
NoisyBRR200 50 100 150 200 25001000200030004000500060007000
NoisyY80 50 100 150 200 25001000200030004000500060007000
NoisyYRR0 50 100 150 200 25001000200030004000500060007000
NoisyP
Figure 3: Image
Boat , h = 16. The columns are, from left to right, the denoised image, a detailof the image, the contour plot of the detail, and the histogram of the denoised image. The rowsgive, from above to below, the noisy image, the results of BPB, BRR20, Y8, YRR and P, andin the last row, the ground truth clean image. 28n a fast bilateral filtering formulation using functional rearrangements NoisyBPB0 50 100 150 200 25000.511.52x 10 NoisyBRR200 50 100 150 200 25000.511.52x 10 NoisyY80 50 100 150 200 25000.511.52x 10 NoisyYRR0 50 100 150 200 25000.511.52x 10 NoisyP
Figure 4: Image
Still life , hh