Ares Lagae

Sylvain Lefebvre

Philip Dutré

TRANSACTIONS ON VISUALIZATION & COMPUTER GRAPHICS 1 Improving Gabor Noise

Keywords: procedural noise, sparse convolution noise, Gabor noise, isotropic Gabor kernel, circular Gabor filter, Hankel transform, circularly symmetric functions, Gabor noise error analysis, spatially varying Gabor noise

We have recently proposed a new procedural noise function, Gabor noise, which offers a combination of properties not found in existing noise functions. In this paper, we present three significant improvements to Gabor noise: (1) an isotropic kernel for Gabor noise, which speeds up isotropic Gabor noise with a factor of roughly two, (2) an error analysis of Gabor noise, which relates the kernel truncation radius to the relative error of the noise, and (3) spatially varying Gabor noise, which enables spatial variation of all noise parameters. These improvements make Gabor noise an even more attractive alternative for existing noise functions.

I. INTRODUCTION

Since its introduction by Perlin in 1985, procedural noise has become an essential component in computer graphics [START_REF] Ebert | Texturing and Modeling: A Procedural Approach[END_REF]. Several noise functions have been proposed, for example Perlin noise [START_REF] Perlin | Improving noise[END_REF], sparse convolution noise [START_REF] Lewis | Algorithms for solid noise synthesis[END_REF], wavelet noise [START_REF] Cook | Wavelet noise[END_REF] and anisotropic noise [START_REF] Goldberg | Anisotropic noise[END_REF]. For a recent survey of procedural noise functions we refer the reader to Lagae et al.. We have recently proposed a new procedural noise function, sparse Gabor convolution noise, or, in short, Gabor noise [Lagae et al., 2009a] (section II). Gabor noise has several interesting properties: it is procedural, it offers significant spectral control, it supports anisotropy, it can be mapped onto surfaces without using a parametrization, it can be filtered, and it is interactive. This combination of properties is not found in existing noise functions.

In this paper, we present three significant improvements to Gabor noise. As a first improvement, we present an isotropic kernel for Gabor noise (section III). This improvement speeds up isotropic Gabor noise with a factor of roughly two. This can result in significant savings of 3D rendering time. Indeed, in the 1990's it was informally observed that "90% of 3D rendering time is spent in shading, and 90% of that time is spent computing Perlin noise" 1 . To our knowledge, the n-dimensional real and even isotropic or circularly symmetric Gabor kernel we derive is not known in literature. [START_REF] Zhang | Invariant texture segmentation via circular Gabor filters[END_REF] have presented a circular Gabor filter in the context of rotation-invariant texture segmentation, but their filter is complex (that is, it has an imaginary part), and is therefore not usable in the context of Gabor noise. As a second improvement, we present an error analysis of Gabor noise (section IV). This improvement relates the kernel truncation radius, an important quality parameter of Gabor noise, to the relative error of the noise, and replaces the Ares Lagae is with the Katholieke Universiteit Leuven and REVES/INRIA Sophia-Antipolis.

Sylvain Lefebvre is with REVES/INRIA Sophia-Antipolis and AL-ICE/INRIA Nancy.

Philip Dutré is with the Katholieke Universiteit Leuven. 1 industry lore, as relayed by J.P. Lewis ad hoc method to choose this parameter with a more principled one. As a third improvement, we present spatially varying Gabor noise (section V). This improvement enables spatial variation of all noise parameters, a property not found in existing procedural noise functions. This allows an artist to create spatially varying procedural textures. Indeed, spatial variation of noise parameters was shown to be useful in the context of non-procedural methods related to Gabor noise [START_REF] Jarke | Spot noise texture synthesis for data visualization[END_REF][START_REF] Ware | Using visual texture for information display[END_REF][START_REF] Holten | A perceptually based spectral model for isotropic textures[END_REF]. These three orthogonal improvements build upon and augment the strong theoretical foundation of Gabor noise, and make Gabor noise an even more attractive alternative for existing noise functions. Although this paper focuses on procedural noise for computer graphics, the family of methods that Gabor noise belongs to is also relevant to visualization [START_REF] Jarke | Spot noise texture synthesis for data visualization[END_REF].

II. GABOR NOISE

In this section, we briefly review Gabor noise [Lagae et al., 2009a]. We focus on its procedural nature and spectral control, the two properties most relevant to this paper. Since this paper addresses improvements to Gabor noise, we assume the reader is generally familiar with Gabor noise.

Anisotropic Gabor noise is a sum of randomly weighted and positioned Gabor kernels, nK,F 0 ,a,ω 0 (x, y) = X i wi gK,F 0 ,a,ω 0 (xxi, yyi),

where the magnitude K, the frequency F 0 , the bandwidth a and the orientation ω 0 are the noise parameters, g is the Gabor kernel, {w i } are the random weights, and {(x i , y i )} are the random positions. The Gabor kernel is the product of a radially symmetric Gaussian and a 2D cosine,

gK,F 0 ,a,ω 0 (x, y) = K exp[-πa 2 (x 2 + y 2 )] cos[2πF0(x cos ω0 + y sin ω0)], (2)
where K and a control the magnitude and width of the Gaussian, and (F 0 , ω 0 ) is the frequency of the cosine. Anisotropic Gabor noise is the convolution of sparse white noise and the Gabor kernel,

nK,F 0 ,a,ω 0 (x, y) = " X i wiδ (x i ,y i ) * gK,F 0 ,a,ω 0 # (x, y), (3)
where the random weights {w i } are distributed according to a random variable W with a uniform distribution on the interval [-1, +1], and the random positions {(x i , y i )} are distributed according to a Poisson distribution with impulse density λ.

Because sparse white noise has a constant power spectrum, the power spectrum of anisotropic Gabor noise is a scaled version of the power spectrum of the Gabor kernel. The Gabor kernel in the frequency domain is a pair of Gaussians,

GK,F 0 ,a,ω 0 (fx, fy) = K 2a 2 exp n - π a 2 ˆ(fx ± F0 cos ω0) 2 + (fy ± F0 sin ω0) 2 ˜o , (4)
where K and a control the magnitude and width of the Gaussians, and (F 0 , ω 0 ) and its symmetrical counterpart are the locations of the Gaussians. Because the power spectrum of anisotropic Gabor noise is a scaled version of the power spectrum of the Gabor kernel, the parameters K, F 0 , a and ω 0 directly control the power spectrum of the noise. Isotropic Gabor noise is a sum of randomly weighted, positioned and oriented Gabor kernels,

nK,F 0 ,a(x, y) = X i wi gK,F 0 ,a(ω0i, x -xi, y -yi), (5)
where the magnitude K, the frequency F 0 and the bandwidth a are the noise parameters, and the random orientations {ω 0i } are distributed according to a random variable Ω with a uniform distribution on the interval [0, 2π). Similar to anisotropic Gabor noise, the parameters K, F 0 and a directly control the power spectrum of the noise.

The procedural evaluation of Gabor noise is similar to that of [START_REF] Lewis | Algorithms for solid noise synthesis[END_REF] sparse convolution noise and [START_REF] Worley | A cellular texture basis function[END_REF] cellular texture basis function. Gabor noise is evaluated procedurally by truncating the Gabor kernel and introducing a grid with a cell size equal to the radius of the truncated kernel. This restricts the evaluation of the noise to the grid cell containing the point of evaluation and the eight neighboring grid cells. The Gabor kernels in each cell are generated on-the-fly using a pseudorandom number generator.

Next to its procedural nature and spectral control, Gabor noise has several other interesting properties for computer graphics: it supports anisotropy, it can be mapped onto surfaces without using a parametrization, it can be filtered, and it is interactive. This is the major difference between Gabor noise and related methods in computer graphics, such as sparse convolution noise, and related methods in visualization, such as spot noise [van Wijk, 1991, Ware and[START_REF] Ware | Using visual texture for information display[END_REF].

III. AN ISOTROPIC KERNEL FOR GABOR NOISE

Isotropic Gabor noise is defined using an anisotropic Gabor kernel (see equation 5). In this section, we show that isotropic Gabor noise can also be defined using an isotropic Gabor kernel, and that the isotropic kernel has several advantages over the anisotropic kernel. Most importantly, we show that isotropic noise using the isotropic kernel is roughly two times faster than isotropic noise using the anisotropic kernel.

We assume that the reader is familiar with circularly symmetric functions (see appendix I), more specifically, with hyperspherical coordinates (see appendix I-A), the integration (see appendix I-B) and convolution (see appendix I-C) of circularly symmetric functions, and the Hankel transform (see appendix I-D).

A. The Isotropic Gabor Kernel

We define the n-dimensional isotropic Gabor kernel similar in spirit as other kinds of Gabor kernels: using a Gaussian and a harmonic, which are related by multiplication in the spatial domain, and by convolution in the frequency domain. We use the Hankel transform (see appendix I-D), which is the method of choice for working with Fourier transforms of isotropic or circularly symmetric functions. We denote the n-dimensional isotropic Gabor kernel, Gaussian, and harmonic as n

n I g (r) = n I gG (r) n I gH (r) n H ⇐⇒ [ n I GG * n I GH ] (fr) = n I G (fr) , (6)
where n H

⇐⇒ denotes an order-n Hankel transform pair. The Gaussian is a Gaussian in both the spatial domain and the frequency domain. We define the Gaussian in the spatial domain, n I g G (r), as the circularly symmetric Gaussian,

n I gG (r) = Ke -πa 2 r 2 , (7)
where K and a are the magnitude and width of the Gaussian.

The Gaussian in the spatial domain is illustrated in figure 1(a) and in figure 2(a). We obtain the Gaussian in the frequency domain, n I G G (fr), as the order-n Hankel transform of n

I g G (r), n I GG (fr) = n H [ n I gG (r)] = 2π f 1 2 n-1 r Z ∞ 0 Ke -πa 2 r 2 J 1 2 n-1 (2πfrr) r 1 2 n dr = K a n e -π a 2 f 2 r , ( 8
)
where Jn is the order-n Bessel function of the first kind. The Gaussian in the frequency domain is illustrated in figure 1(b) and in figure 2(b).

The harmonic is typically defined as an impulse in the frequency domain, located at the principal frequency, F 0 , of the Gabor kernel. We therefore define the harmonic in the frequency domain, n I G H (fr), as a circularly symmetric impulse,

n I GH (fr) = δ (fr -F0) , (9)
where F 0 is the frequency of the harmonic. The harmonic in the frequency domain is illustrated in figure 1(d) and in figure 2(d).

We obtain the harmonic in the spatial domain, n I g H (r), as the order-n Hankel transform of n I G H (fr),

n I gH (r) = n H [ n I GH (fr)] = 2π r 1 2 n-1 Z ∞ 0 δ (fr -F0) J 1 2 n-1 (2πrfr) f 1 2 n r dfr = 2π r 1 2 n-1 J 1 2 n-1 (2πF0r) F 1 2 n 0 . (10)
The harmonic in the spatial domain is illustrated in figure 1(c) and in figure 2(c).

We obtain the n-dimensional isotropic Gabor kernel in the spatial domain, n I g (r), as the multiplication of n I g G (r) and

n I g H (r), n I g (r) = n I gG (r) n I gH (r) = Ke -πa 2 r 2 2πF 1 2 n 0 r 1 2 n-1 J 1 2 n-1 (2πF0r) . (11)
The kernel in the spatial domain is illustrated in figure 1(e) and in figure 2(e). We obtain the n-dimensional isotropic Gabor kernel in the frequency domain, n I G (fr), as the convolution of n I G G (fr) and n I G H (fr). We simplify the convolution of isotropic or circularly symmetric functions by exploiting their symmetry (see appendix I-C). First, we convolve n I G G (fr) and n I G H (fr) using this simplification, where Γ is the Gamma function. Then, we integrate over f ′ r . Next, we simplify the convolution using the integral

n I G (fr) = [ n I GG * n I GH ] (fr) = 2π n-1 2 Γ `n-1 2 ´Z ∞ f ′ r =0 Z π f φ =0 δ `f ′ r -F0 ´K a n e -π a 2 (f 2 r +f ′ r 2 -2fr f ′ r cos f φ) f ′ r n-1 sin n-2 f φ df ′ r df φ , (12)
Z π 0 sin n θ e x cos θ dθ = 2 n 2 √ πx -n 2 I n 2 (x) Γ " n + 1 2 « , ( 13
)
where In is the order-n modified Bessel function of the first kind. Finally, we obtain

n I G (fr) = 2πKF 1 2 n 0 a 2 f 1 2 n-1 r e -π a 2 (f 2 r +F 2 0 ) I 1 2 n-1 " 2πF0 a 2 fr « . ( 14
)
The kernel in the frequency domain is illustrated in figure 1(f) and in figure 2(f).

B. Isotropic Gabor Noise using the Isotropic Gabor Kernel

We define n-dimensional isotropic Gabor noise using the isotropic Gabor kernel similar to anisotropic (see equation 1) and isotropic (see equation 5) Gabor noise using the anisotropic Gabor kernel,

n I n (x1, . . . , xn) = X i wi n I g " q (x1 -xi,1) 2 + . . . + (xn -xi,n) 2 « . (15)
Note that, in contrast to the anisotropic kernel, the isotropic kernel does not need to be randomly oriented. The variance of the noise n

I σ 2 n is n I σ 2 n = λE ˆW 2 ˜Sn Z ∞ r=0 n I g 2 (r) r n-1 dr. ( 16
)
The power spectrum of the noise n I Snn is

n I Snn (fr) = λE ˆW 2 ˜|n I G (fr)| 2 . ( 17
)
We provide equations for working with one-, two-, three-and four-dimensional isotropic Gabor noise using the isotropic Gabor kernel (see appendix II), including the isotropic Gabor kernel in the spatial domain (see appendix II-A) and in the frequency domain (see appendix II-B), the integral of the isotropic Gabor kernel squared (see appendix II-C), the envelope of the isotropic Gabor kernel (see appendix II-D), and the radius of the truncated isotropic Gabor kernel (see appendix II-E).

C. Implementation, Results and Comparison

We have implemented isotropic Gabor noise using the isotropic Gabor kernel, and we have verified most equations experimentally. We evaluate the Bessel functions using code based on Press et al. [2002, 6.5, 6.6] (polynomial approximations are also available in Abramowitz and Stegun [1972, 9.4,9.8]), and the Lambert-W function (see appendix II-E) using code based on [START_REF] Keith | Lambert W function[END_REF]. We verify the noise by comparing the estimated and expected power spectrum and the actual and expected intensity distribution [START_REF] Ares Lagae | Procedural noise using sparse Gabor convolution -auxiliary material[END_REF].

We illustrate one-, two-and three-dimensional isotropic Gabor noise using the isotropic Gabor kernel in figure 3, figure 4 and figure 5. Note the similarity between isotropic noise using the isotropic kernel and isotropic noise using the anisotropic kernel [Lagae et al., 2009a, figure 4]. Also note how closely the estimated and expected power spectrum and the actual and expected intensity distribution match.

We have found that isotropic noise using the isotropic kernel is significantly faster than isotropic noise using the anisotropic kernel. This is because of two reasons. First, the isotropic kernel is more compact than the anisotropic kernel (see figure 1(e) and figure 6). This results in a lower number of impulses per kernel for the same impulse density, which in turn results in a shorter time to evaluate the noise, since this time is directly proportional to the number of impulses. Second, in contrast to the anisotropic kernel, the isotropic kernel does not need to be randomly oriented (see equation 5 and equation 15). This avoids the generation of random orientations, which also results in a shorter time to evaluate the noise. For example, for two-dimensional isotropic noise with parameters K = 0.709645, a = 0.0443528 and F 0 = 0.0625, where the kernel was truncated at 5% of its maximum value, the radius of the kernel is 22.0169 for the anisotropic kernel, but only 15.6741 for the isotropic kernel, a reduction of 28.8086%, and for an impulse density of λ = 0.0414605, the number of impulses per kernel is 63.1386 for the anisotropic kernel, but only 32 for the isotropic kernel, a reduction of 49.3179%. The time to evaluate the noise using our CPU implementation is 4.11625 s for the anisotropic kernel, but only 1.64481 s for the isotropic kernel, a speedup of 2.50257 (512 × 512, Intel(R) Xeon(R) CPU 5160 at 3.00GHz). With similar parameter settings, the performance of the noise using our GPU implementation is 77 FPS (frames per second) for the anisotropic kernel, and 144 FPS for the isotropic kernel, a speedup of 1.87 (512 × 512, NVIDIA Quadro FX 5800 GPU). We have found that the speedup is roughly independent of the noise parameters. The speedup for three-and four-dimensional isotropic noise should be even larger, since the isotropic kernel gets more compact with increasing dimension (see figure 1(e) and figure 6), and the number of random orientations that has to be sampled gets larger with increasing dimension (see appendix I-A).

We have compared isotropic Gabor noise using the isotropic kernel with Perlin noise [START_REF] Perlin | Improving noise[END_REF] and wavelet noise [START_REF] Cook | Wavelet noise[END_REF], two other isotropic noise functions. The visual quality of these noise functions is illustrated in figure 7. Perlin noise, and to a lesser degree also wavelet noise, both exhibit an undesired axis-aligned anisotropy, which is not the case for Gabor noise (also see [Lagae et al.]). We have measured the performance of these noise functions using our GPU implementations. The performance of Perlin noise is 3388 FPS, and that of wavelet noise is 730 FPS (512 × 512, NVIDIA Quadro FX 5800 GPU). The performance of Gabor noise is dependent on the parameter settings and the approximations used in the implementation. With similar parameter settings as above, the performance of Gabor noise is 144 FPS, and 255 FPS for half the number of impulses per kernel. When approximating the Poisson distribution by its mean, the performance is 198 FPS and 357 FPS. The filtering quality of wavelet noise and Gabor noise is illustrated in figure 8 and in video 1. We filter wavelet noise by deriving a filtering weight by evaluating the filtering Gaussian in the frequency domain at the principal frequency of the noise. We filter Gabor noise by deriving new kernel parameters K, F 0 and a by multiplying the filtering Gaussian in the frequency domain with a Gaussian approximation of the isotropic kernel in the frequency domain. In both cases, the anisotropic filtering Gaussian in the frequency domain is approximated with an isotropic one. Gabor noise enables a better filtering quality than wavelet noise since Gabor noise is better than wavelet noise at preserving detail at high spatial frequencies in the far distance. This is clearly visible in video 1. This is because wavelet noise only allows to apply a single filtering weight w to an entire noise octave, while Gabor noise allows to adapt the frequency f and bandwidth a of individual kernels. We have not included Perlin noise in this comparison, since Cook and DeRose [2005, figure 1] already showed that it is worse than wavelet noise in terms of filtering quality, and because filtering Perlin noise in a principled way is difficult since Perlin noise is not band-pass. We would like to note that there are many other criteria that could also be taken into account in a comparison [Lagae et al.], which might or might not be relevant depending on the application.

Our results and comparison show that isotropic Gabor noise using the isotropic Gabor kernel should be used whenever an isotropic noise function is required with a higher quality than Perlin noise and wavelet noise and a lower cost than isotropic Gabor noise using the anisotropic Gabor kernel.

IV. AN ERROR ANALYSIS OF GABOR NOISE

The procedural evaluation of Gabor noise requires that the Gabor kernel is truncated. This is typically done at a radius where the envelope of the kernel reaches a sufficiently small value, for example 5% of its maximum value [Lagae et al., 2009a, 4] [Lagae et al., 2009b, 1]. However, this is an ad hoc approach, since truncating the kernel introduces an error in the noise, and this error is not quantified. In this section, we quantify the effect of truncating the Gabor kernel, and we present a more principled approach to truncate the Gabor kernel.

A. Relation of Kernel Truncation Radius to Noise Error

We relate the noise error resulting from truncating the Gabor kernel g to the truncation radius r t . In this analysis, we use the n-dimensional isotropic Gabor kernel (see section III), but the analysis is also valid for other kernels, such as the anisotropic Gabor kernel. We define a truncated Gabor kernel g t , where

gt(r) =  g(r) 0 ≤ rt < r 0 rt ≥ r , (18)
and an error kernel ∆g, where

∆g(r) =  0 0 ≤ rt < r g(r) rt ≥ r . ( 19
)
Note that g(r) = g t (r) + ∆g(r). We obtain the noise error ∆n by subtracting the noise using the truncated kernel n t from the " .

(20)

Our key insight is that the noise error ∆n is, similar to Gabor noise, a random pulse process [Lagae et al., 2009a, 2.2] [van Etten, 2005, 8] [ Papoulis and Pillai, 2002, 10.2]. The variance σ 2 ∆n of the noise error ∆n is therefore

σ 2 ∆n = λ E ˆW 2 ˜Sn Z ∞ 0 ∆g 2 (r)dr. ( 21
)
We define the root mean square error and the relative error using the variances σ 2 n and σ 2 ∆n of the noise and the noise error. We define the root mean square error e RM S as the square root of the variance of the noise error,

eRMS = q σ 2 ∆n . ( 22
)
We define the relative error e as the root mean square error e RM S over the root mean square amplitude

p σ 2 n , e = eRMS √ σ 2 n = s R ∞ 0 ∆g 2 (r) dr R ∞ 0 g 2 (r) dr = s 1 - R r t 0 g 2 (r) dr R ∞ 0 g 2 (r) dr . ( 23
)
Note that the relative error e only depends on the kernel g and the truncation radius r t , and not on the parameters of the sparse white noise λ and W . Our analysis allows to determine the relative noise error for a given kernel truncation radius, but also allows to determine the kernel truncation radius for a given relative error, by solving equation 23 for r t . The usage of relative error in this context is motivated by Weber's Law [START_REF] Blackwell | Luminance difference thresholds[END_REF]. This is a much more principled approach than truncating the kernel at an arbitrarily chosen value.

B. Implementation, Results and Discussion

We have implemented the noise error analysis and we have verified the equations experimentally. For most kernels, a closed-form expression for R ∞ 0 g 2 (r) dr is available (see appendix II-C), but a closed-form expression for R rt 0 g 2 (r) dr is not. Therefore, we generally solve equation 23 for r t numerically, using bracketing and bisection [Press et al., 2002, 9.1], the closed-form expression to evaluate R ∞ 0 g 2 (r) dr, and Simpson's rule [Press et al., 2002, 4.2] to evaluate R rt 0 g 2 (r) dr. Note that for isotropic noise using the isotropic kernel, all integrals are one-dimensional.

We illustrate the error analysis of one-dimensional isotropic noise using the isotropic kernel in figure 9. Note how closely the actual and expected root mean square error match. We illustrate the error analysis of two-dimensional isotropic noise using the isotropic kernel in figure 10. We plot the relative error versus the kernel radius for one-, two-, three-and four-dimensional isotropic noise using the isotropic kernel in figure 11. Note that the relative error quickly decreases with increasing kernel truncation radius. We now revisit the example of subsection III-C. When truncating the kernel using a relative error of 2%, the radius of the kernel is 25.25 for the anisotropic kernel, but only 20.8984 for the isotropic kernel, a reduction of 17.2341%, the number of impulses per kernel is 46.7139 for the anisotropic kernel, but only 32 for the isotropic kernel, a reduction of 31.498%, and the time to evaluate the noise using our CPU implementation is 2.50428 s for the anisotropic kernel, but only 1.02537 s for the isotropic kernel, a speedup of 2.44233. We generalize this example and Each subfigure shows a graph and three images. The graph shows the noise using the untruncated kernel (n) (reference implementation), the noise using the truncated kernel (nt) (procedural implementation), the noise error (∆n), the actual root mean square noise error, and the estimated root mean square noise error (e RM S ). The three images show the noise using the untruncated kernel (n), the noise using the truncated kernel (nt), and the noise error (∆n). plot the relative error versus the kernel radius for two-dimensional isotropic noise using the anisotropic and the isotropic kernel in figure 12. This figure shows that for the same truncated kernel radius, the relative error is always smaller for the isotropic kernel than for the anisotropic kernel, and for the same relative error, the isotropic kernel is always smaller and therefore faster than the anisotropic kernel.

We believe that in general equation 23 should be used rather than the ad hoc approach to determine the radius of the truncated kernel, except for reasons of computational expense, for example in the case of spatially varying noise (see section V), or simplicity.

V. SPATIALLY VARYING GABOR NOISE

The procedural evaluation of Gabor noise introduces a grid with a cell size equal to the radius of the truncated kernel [Lagae et al., 2009a, 4]. This restricts the evaluation of the noise to the grid cell containing the point of evaluation and the eight neighboring grid cells. However, this also prohibits spatial variation of the noise parameters, since the kernel radius, and therefore the cell size, is in general dependent on these parameters. In this section, we present a procedural evaluation of Gabor noise that enables noise with spatially varying parameters.

A. Procedural Evaluation of Gabor Noise with Spatially Varying Parameters

Evaluating Gabor noise using a single grid is optimal for a single kernel radius corresponding to a fixed set of noise parameters. In order to handle arbitrary kernel radii, corresponding to arbitrary or noise parameters, we use a hierarchy of grids, where the cell size of consecutive grids differs by a factor of two. More specifically, the grid at level l of the hierarchy has a cell size of 2 l . This implies that the grid at level l is optimal for a kernel radius r l = 2 l . We have observed that, for different noise parameters, an equal noise quality is obtained by maintaining a constant number of impulses N per kernel area. Therefore, we associate the grid at level l with an impulse density of λ l = N/πr 2 l . When evaluating the noise using an arbitrary kernel radius r, we use the two grids with kernel radii r 0 and r 1 that bracket r. These are the grids with level l 0 = ⌊log 2 r⌋ and l 1 = l 0 + 1, for which r 0 ≤ r < r 1 . We parametrize the radius r in terms of r 0 and r 1 by introducing a parameter α in [0, 1) such that r = (1α)r 0 + αr 1 . When evaluating the noise, we combine the impulses of both grids by weighing the contribution due to each grid,

n(x, y) =w0(α) X i 0 wi 0 g(x -xi 0 , y -yi 0 ) +w1(α) X i 1 wi 1 g(x -xi 1 , y -yi 1 ), (24)
where the first and second term correspond to the grid with level l 0 and l 1 , and w 0 and w 1 are weighting functions. This is illustrated in figure 13. It is important to note that this is a linear combination of two noises with the same parameters, and that the goal of the interpolation is to transition between grids with different impulse densities, and not to obtain a noise with interpolated parameters. Interpolating between two noises with frequencies F 0 and F 1 (for example two octaves of Perlin noise) cannot produce a noise with an intermediate frequency F , due to the linearity of the Fourier transform, while interpolating between two noises that already have frequency F obviously can.

For obvious reasons, we require that w 0 (0) = 1, that w 0 (1) = 0, and that w 0 is monotonically decreasing, and similarly that w 1 (0) = 0, that w 1 (1) = 1, and that w 1 is monotonically increasing. However, several choices for w 0 and w 1 remain. A simple option is to choose w 0 (α) = 1w(α) and w 1 (α) = w(α), where w(α) is the linear weighting function w(α) = α or the cubic weighting function w(α) = 3α 2 -2α 3 . However, the resulting weighting functions do not result in visually pleasing transitions. This is illustrated in video 2. The weighting functions can be determined by imposing additional constraints. For example, imposing the constraint that the power spectrum of both sides of equation 24 is the same leads to w 0 (α)(1 + α) + w 1 (α)(1+α)/2 = 1. We have noticed that the visually unpleasing transitions are caused by the fact that the two grids have a different impulse density. The number of impulses per kernel area N 0 due to the grid with level l 0 equals N 0 = (1 + α) 2 N , while the number of impulses per kernel area N 1 due to the grid with level l 1 equals N 1 = N 0 /4. Therefore, we impose the constraint that the weighted number of impulses per kernel area remains constant, more specifically such that w 0 (α)N 0 + w 1 (α)N 1 = N or equivalently such that w 0 (α

)(1 + α) 2 + w 1 (α)(1 + α) 2 /4 = 1.
We now choose w 0 (α)(1+α) 2 = 1-w(α) and w 1 (α)(1+α) 2 /4 = w(α), where w(α) is a weighting function. The simplest option for w(α) is the linear weighting function w(α) = α, which also satisfies the power spectrum constraint, and results in the weighting functions

w 0 (α) = 1 -α (1 + α) 2 (25) w 1 (α) = 4α (1 + α) 2 .
(26) This is illustrated in figure 14. These weighting functions result in visually pleasing transitions. This is illustrated in video 2.

B. Implementation, Results, Comparison and Discussion

We have implemented spatially varying Gabor noise and we have verified the equations experimentally. Note that in contrast with regular Gabor noise, the evaluation of spatially varying Gabor noise includes the computation of the kernel truncation radius. Because of computational expense, we use the ad hoc approach (see section IV). Also note that the noise evaluation is not restricted to exactly nine grid cells, and that the seeding strategy has to take into account the grid level.

We illustrate spatially varying 2D Gabor noise in figure 15 and video 2. We show procedural textures generated with spatially varying surface Gabor noise in figure 16. Note that the simple alternative of warping the domain variables of a regular Gabor noise to obtain spatially varying Gabor noise might work in the case of a 2D noise, as in figure 15, but would probably not work in the case of spatially varying surface noise, as in figure 16, and might be problematic for filtering.

We have compared the performance of spatially varying Gabor noise and regular Gabor noise. The hierarchical solution of spatially varying Gabor noise is only roughly a factor two slower than regular Gabor noise. We have also compared the performance of Perlin noise, regular Gabor noise and spatially varying Gabor noise in the specific case of figure 16. The performance of Perlin noise was 763 FPS, that of regular Gabor noise 22.8 FPS, and that of spatially varying Gabor noise 20.6 FPS (2 octaves, 1024×1024, NVIDIA Quadro FX 3800 GPU). In this specific case, spatially varying Gabor noise is almost as fast as regular Gabor noise. Perlin noise is significantly faster than Gabor noise, but also lacks several attractive features of Gabor noise such as spatial variation and filtering.

A somewhat similar mechanism as the one used for spatially varying Gabor noise was recently used by [START_REF] Benard | A dynamic noise primitive for coherent stylization[END_REF] to construct a dynamic noise primitive for coherent stylization in expressive rendering.

VI. CONCLUSION

In this paper, we have presented three significant improvements to Gabor noise: an isotropic kernel for Gabor noise, which speeds up isotropic Gabor noise with a factor of roughly two, an error analysis of Gabor noise, which relates the kernel truncation radius to the relative error of the noise, and spatially varying Gabor noise, which enables spatial variation of all parameters. These improvements build upon and augment the strong theoretical foundation of Gabor noise, and make Gabor noise an even more attractive alternative for existing noise functions.

to a one-dimensional integral,

Z +∞ x 1 =-∞ . . . Z +∞ xn=-∞ f (x1, . . . , xn) dx1 . . . dxn = Sn Z ∞ r=0 f (r) r n-1 dr, (29)
where Sn is the hyper-surface area of an n-sphere of unit radius,

Sn = 2π n 2 Γ `n 2 ´, ( 30
)
where Γ is the Gamma function.

C. Convolution of Circularly Symmetric Functions

We simplify the convolution of circularly symmetric functions by exploiting their symmetry. First, we formulate the convolution of two n-dimensional circularly symmetric functions f and g using hyperspherical coordinates (see subsection I-A) as

[f * g] (r) = Z ∞ r ′ =0 Z π φ ′ 1 =0 . . . Z π φ ′ n-2 =0 Z 2π φ ′ n-1 =0 f `r′ ´g (R) r ′n-1 sin n-2 φ ′ 1 sin n-3 φ ′ 2 . . . sin φ ′ n-2 dr ′ dφ ′ 1 . . . dφ ′ n-1 , (31)
where R 2 = r 2 +r ′2 -2rr ′ cos φ ′ 1 . We set the angular coordinates

φ ′ 1 . . . φ ′
n-1 to zero in the expression for R 2 , since all functions are circularly symmetric, leaving only the angular coordinate φ ′ 1 . Next, we simplify the convolution using the integral

Z π 0 sin n θ dθ = √ π Γ `n+1 2 Γ `n+2 2 ´, (32)
where Γ is the Gamma function. Finally, we obtain

[f * g] (r) = 2π n-1 2 Γ `n-1 2 ´Z ∞ r ′ =0 Z π φ=0 f `r′ ´g (R) r ′n-1 sin n-2 φ dr ′ dφ, (33)
with R 2 = r 2 + r ′2 -2rr ′ cos φ. For n = 2, equation 33 corresponds to the equation given in Bracewell [2000, table 13.3].

D. The Hankel Transform

The Hankel transform [START_REF] Bracewell | The Fourier Transform and its Applications[END_REF], 13] [Bracewell, 2004, 9] is the method of choice for working with Fourier transforms of circularly symmetric functions.

When a function f is circularly symmetric in n dimensions, that is, when f (x 1 , . . . , xn) = f (r), where r 2 = x 2 1 + . . . + x 2 n , then F , the n-dimensional Fourier transform of f , is also circularly symmetric, that is, F (fx 1 , . . . , fx n ) = F (fr), where

f 2 r = f 2 x1 + . . . + f 2 xn .
The relation between the one-dimensional functions f (r) and F (fr) is given by the Hankel transform of order n. More specifically, the n-dimensional Fourier transform of f (x 1 , . . . , xn) is given by the order-n Hankel transform of f (r), that is, n F [f (x 1 , . . . , xn)] = n H [f (r)]. Note that, somewhat counterintuitively, the n-dimensional Fourier transform of f (x 1 , . . . , xn) is not equal to the one-dimensional Fourier transform of f (r), that is, n F [f (x 1 , . . . , xn)] = 1 F [f (r)].

The Hankel transform, also called the Fourier-Bessel transform, is a one-dimensional integral transform with a Bessel function kernel. The Hankel transform of order n is n H [f (r)] = F (fr) = 2π

f 1 2 n-1 r Z ∞ 0 f (r) J 1 2 n-1 (2πfrr) r 1 2 n dr, ( 34
)
where Jn is the order-n Bessel function of the first kind. The Hankel transform is strictly reciprocal. An order-n Hankel transform pair is denoted as f (r) n H ⇐⇒ F (fr). For n = 1, the Hankel transform corresponds to the Fourier transform, since circularly symmetric functions are real and even, and J 1/2 (x) = (2/πx) 1/2 sin x and J -1/2 (x) = (2/πx) 1/2 cos x.

APPENDIX II EQUATIONS FOR ISOTROPIC GABOR NOISE USING THE ISOTROPIC GABOR KERNEL

In this section, we provide equations for working with one-, two-, three-and four-dimensional isotropic Gabor noise using the isotropic Gabor kernel. We provide equations for the isotropic Gabor kernel in the spatial domain (subsection II-A) and in the frequency domain (subsection II-B), the integral of the isotropic Gabor kernel squared (subsection II-C), the envelope of the isotropic Gabor kernel (subsection II-D), and the radius of the truncated isotropic Gabor kernel (subsection II-E).

A. The Isotropic Gabor Kernel in the Spatial Domain

We obtain the one-, two-, three-and four-dimensional isotropic Gabor kernel in the spatial domain from equation 11, keeping into account that that J 1/2 (x) = (2/πx) 1/2 sin x and J -1/2 (x) = (2/πx) 1/2 cos x. The one-, two-, three-and four-dimensional isotropic Gabor kernel in the spatial domain is

1 I g (r) = Ke -πa 2 r 2 2 cos (2πF0r) , (35)
We illustrate the isotropic Gabor kernel in the spatial domain in figure 1(e) and in figure 2(e).

B. The Isotropic Gabor Kernel in the Frequency Domain

We obtain the one-, two-, three-and four-dimensional isotropic Gabor kernel in the frequency domain from equation 14, keeping into account that I 1/2 (x) = (2πx) -1/2 `ex -e -x ´and I -1/2 (x) = (2πx) -1/2 `ex + e -x ´. The one-, two-, three-and four-dimensional isotropic Gabor kernel in the frequency domain is We illustrate the isotropic Gabor kernel in the frequency domain in figure 1(f) and in figure 2(f).

C. The Integral of the Isotropic Gabor Kernel Squared

We obtain the integral of the isotropic Gabor kernel squared by integrating the Gabor kernel squared in the spatial domain or in the frequency domain. The integral of the one-, two-, threeand four-dimensional isotropic Gabor kernel squared is

I

  g (r), n I g G (r), and n I g H (r) in the spatial domain, and as n I G (fr), n I G G (fr), and n I G H (fr) in the frequency domain. We summarize their relations

Fig. 1 .

 1 Fig. 1. The isotropic Gabor kernel. (a) Gaussian, spatial domain. (b) Gaussian, frequency domain. (c) Harmonic, spatial domain. (d) Harmonic, frequency domain. (e) Kernel, spatial domain. (f) Kernel, frequency domain.

  (a) Gaussian, spat. dom. (b) Gaussian, freq. dom. (c) Harmonic, spat. dom. (d) Harmonic, freq. dom. (e) Kernel, spat. dom. (f) Kernel, freq. dom.

Fig. 2 .

 2 Fig. 2. The 2D isotropic Gabor kernel. (a) Gaussian, spatial domain. (b) Gaussian, frequency domain. (c) Harmonic, spatial domain. (d) Harmonic, frequency domain. (e) Kernel, spatial domain. (f) Kernel, frequency domain.

  Intensity distribution (histogram and expected).

  Spectrum estimated power spectrum expected power spectrum (f) Radial power spectrum (estimate and expected).

Fig. 4 .

 4 Fig. 4. 2D isotropic Gabor noise using the isotropic kernel. (a) Noise. (b) Actual and expected intensity distribution. (c) Fourier transform (magnitude). (d) Estimated power spectrum. (e) Expected power spectrum. (f) Estimated and expected radial power spectrum.

  Intensity distribution (histogram and expected).

  Spectrum estimated power spectrum actual power spectrum (f) Radial power spectrum (estimate and expected).

Fig. 5 .

 5 Fig. 5. 3D isotropic Gabor noise using the isotropic kernel. (a) Noise. (b) Actual and expected intensity distribution. (c) Fourier transform (magnitude). (d) Estimated power spectrum. (e) Expected power spectrum. (f) Estimated and expected radial power spectrum.

  Fourier transform (magnitude).

  (d) Power spectrum (estimate and expected).

Fig. 3 .

 3 Fig. 3. 1D isotropic Gabor noise using the isotropic kernel. (a) Noise. (b) Actual and expected intensity distribution. (c) Fourier transform (magnitude). (d) Estimated and expected power spectrum.

Fig. 6 .

 6 Fig. 6. The envelope of the isotropic Gabor kernel. (a) Envelope of the harmonic. (b) Envelope of the kernel.

Fig. 7 .

 7 Fig. 7. A comparison of visual quality of different noise functions. (a) Perlin noise. (b) Wavelet noise. (c) Gabor noise using the isotropic kernel.Note that Perlin noise, and to a lesser degree also wavelet noise, both exhibit an undesired axis-aligned anisotropy.

Fig. 8 .

 8 Fig. 8. A comparison of filtering quality of different noise functions. (a) Wavelet noise. (b) Gabor noise using the isotropic kernel. Each subfigure shows a tilted plane with unfiltered noise on the left half and filtered noise on the right half. Note that Gabor noise is better than wavelet noise at preserving detail at high spatial frequencies in the far distance. Please also refer to video 1.

  noise using the untruncated kernel n, ∆n (x1, . . . , xn) = n (x1, . . . , xn)nt (x1, . . . , xn) x1 i ) 2 + . . . + (xnxn i ) 2

Fig. 9 .

 9 Fig. 9. An error analysis of 1D isotropic Gabor noise. (a) Noise with a high relative error (e = 25%). (b) Noise with a low relative error (e = 1%).Each subfigure shows a graph and three images. The graph shows the noise using the untruncated kernel (n) (reference implementation), the noise using the truncated kernel (nt) (procedural implementation), the noise error (∆n), the actual root mean square noise error, and the estimated root mean square noise error (e RM S ). The three images show the noise using the untruncated kernel (n), the noise using the truncated kernel (nt), and the noise error (∆n).

  (a) High relative error. (b) Low relative error.

Fig. 10 .

 10 Fig. 10. An error analysis of 2D isotropic Gabor noise. (a) Noise with a high relative error (e = 25%). (b) Noise with a low relative error (e = 1%).Each subfigure shows the noise using the untruncated kernel (n) (reference implementation), the noise using the truncated kernel (nt) (procedural implementation), and the noise error (∆n).

  Kernel radius vs. relative error.

Fig. 11 .

 11 Fig. 11. Relative error versus kernel radius for isotropic Gabor noise using the isotropic kernel. (a) Relative error versus kernel radius. (b) Kernel radius versus relative error.

Fig. 12 .

 12 Fig. 12. Relative error versus kernel radius for 2D isotropic Gabor noise using the anisotropic and the isotropic kernel (a) Relative error versus kernel radius. (b) Kernel radius versus relative error.

Fig. 13 .Fig. 14 .

 1314 Fig.13. The hierarchy of grids used in spatially varying Gabor noise. (a) The grids at level l 0 (red) and l 1 (blue). (b-f) The weighted combined impulses of both grids for different kernel radii. Note that the red points disappear faster than the blue points appear.

Fig. 15 .

 15 Fig. 15. Spatially varying 2D Gabor noise. The noise parameters vary from K = 0.7096, F 0 = 0.03125 and a = 0.02218 at the left to K = 0.7096, F 0 = 0.25 and a = 0.1774 at the right. The value of α is visualized at the bottom of the image. Note that this image spans four grid levels.

ACKNOWLEDGMENT

We would like to thank the anonymous reviewers. We are grateful to the NVIDIA Professor Partnership Program for equipment donations. Ares Lagae is a Postdoctoral Fellow of the Research Foundation -Flanders (FWO), and acknowledges K.U.Leuven CREA funding (CREA/08/017). Sylvain Lefebvre receives support from the Agence Nationale de la Recherche (SIMILAR-CITIES ANR-2008-COORD-021).

APPENDIX I CIRCULARLY SYMMETRIC FUNCTIONS

In this section we cover circularly symmetric functions. We review hyperspherical coordinates (subsection I-A) and the integration of circularly symmetric functions (subsection I-B), we simplify the convolution of circularly symmetric functions (subsection I-C), and we review the Hankel transform (subsection I-D).

A. Hyperspherical Coordinates

The hyperspherical coordinate system, the generalization of two-dimensional polar coordinates and three-dimensional spherical coordinates, is a natural coordinate system for working with circularly symmetric functions. The hyperspherical coordinates of a point in n-dimensional space with Cartesian coordinates (x 1 , . . . , xn) consist of a radial coordinate, r, and n -1 angular

where r ∈ [0, ∞), φ 1 . . . φ n-2 ∈ [0, π) and φ n-1 ∈ [0, 2π). The corresponding volume element is r n-1 sin n-2 φ1 sin n-3 φ2 . . . sin φn-2 dr dφ1 . . . dφn-1. (28)

B. Integration of Circularly Symmetric Functions

The integration of circularly symmetric functions can be simplified by exploiting their symmetry. When a function f is circularly symmetric in n dimensions, then the integral of f over R n reduces

D. The Envelope of the Isotropic Gabor Kernel

We define the envelope of the isotropic Gabor kernel as the product of the Gaussian with the envelope of the harmonic. The envelope of cos (x) and sin (x) is 1, and the envelope of Jn (x) is

, where Yn is the order-n Bessel function of the second kind. The envelope of the one-, two-, threeand four-dimensional isotropic Gabor kernel is

The envelope of the isotropic Gabor kernel is illustrated in figure 6.

E. The Radius of the Truncated Isotropic Gabor Kernel

We obtain the radius of the truncated isotropic Gabor kernel by solving the envelope of the isotropic Gabor kernel for the radius, and evaluating the result for a specific value of the envelope. The envelope of the one-dimensional isotropic Gabor kernel solved for the radius is

We solve the envelope of the two-dimensional isotropic Gabor kernel for the radius by approximating Mn (x) by Mn (x) ≈ (2/πx) 1/2 , using the asymptotic expansions for large arguments [Abramowitz and Stegun, 1972, 9.2.1,9.2.2]. The envelope of the two-dimensional isotropic Gabor kernel solved for the radius is

where W is the Lambert W-function, that is, the inverse of f (W ) = W e W . The envelope of the three-dimensional isotropic Gabor kernel solved for the radius is

The envelope of the four-dimensional isotropic Gabor kernel cannot be solved for the radius.


Lina Fahed
email: lina.fahed@loria.fr

Armelle Brun
email: armelle.brun@loria.fr

Anne Boyer
email: anne.boyer@loria.fr

Episode Rules Mining Algorithm for Distant Event Prediction

Keywords: Data Mining, Episode Rules Mining, Minimal Rules, Distant Event Prediction

This paper focuses on event prediction in an event sequence, where we aim at predicting distant events. We propose an algorithm that mines episode rules, which are minimal and have a consequent temporally distant from the antecedent. As traditional algorithms are not able to mine directly rules with such characteristics, we propose an original way to mine these rules. Our algorithm, which has a complexity similar to that of state of the art algorithms, determines the consequent of an episode rule at an early stage in the mining process, it applies a span constraint on the antecedent and a gap constraint between the antecedent and the consequent. A new confidence measure, the temporal confidence, is proposed, which evaluates the confidence of a rule in relation to the predefined gap. The algorithm is validated on an event sequence of social networks messages. We show that minimal rules with a distant consequent are actually formed and that they can be used to accurately predict distant events.

INTRODUCTION

The flow of messages posted in blogs and social networks is an important and valuable source of information that can be analyzed, modeled (through the extraction of hidden relationships) and from which information can be predicted, which is the focus of our work. For example, companies may be interested in the prediction of what will be said about them in social networks. Similarly, this prediction can be a way to recommend items. We consider that the sooner an event is predicted, the more useful this prediction is for the company or the person concerned, since this allows to have enough time to act before the occurrence of the event. Predicting distant events is thus the focus of our work.

Temporal data mining is related to the mining of sequential patterns ordered by a given criterion such as time or position [START_REF] Laxman | A survey of temporal data mining[END_REF]. Episode mining is the appropriate pattern discovery task related to the case the data is made up of a single long sequence. An episode is a temporal pattern made up of "relatively close" partially ordered items (or events), which often appears throughout the sequence or in a part of it [START_REF] Mannila | Discovery of frequent episodes in event sequences[END_REF]. When the order of items is total, the episode is said to be serial. Similarly to the extraction of association rules from itemsets, episode rules can be extracted from episodes to predict events [START_REF] Cho | On-line rule matching for event prediction[END_REF]. The rule mining task is usually decomposed into two subproblems. The first one is the discovery of frequent itemsets or episodes that have a support higher than a predefined threshold. The second one is the generation of rules from those frequent itemsets or episodes, with the constraint of a minimal confidence threshold [START_REF] Agrawal | Mining association rules between sets of items in large databases[END_REF]. In general, rules are generated by considering some items in the itemset (or the last items in the episode) as the consequent of the rule, and the rest of the items as the antecedent. Since the second sub-problem is quite straightforward, most of the researches focus on the first one: the extraction of itemsets or episodes.

Episode and episode rules mining are used in many areas, such as telecommunication alarm management [START_REF] Mannila | Discovery of frequent episodes in event sequences[END_REF], intrusion detection [START_REF] Luo | Mining fuzzy association rules and fuzzy frequency episodes for intrusion detection[END_REF], discovery of relation between financial events [START_REF] Ng | Mining frequent episodes for relating financial events and stock trends[END_REF], etc.

Our goal is to reliably predict events that will occur after a predetermined temporal distance, in order to have enough time to act before the occurrence of events. Therefore, serial episode rules with a consequent distant from the antecedent will be mined. Traditional episode rules mining algorithms form episode rules with a consequent close to the antecedent. To mine rules with a distant consequent, these algorithms have to perform a post-processing step: the extracted rules are filtered to keep only rules with a consequent that may occur far from the antecedent. This postprocessing is not only time consuming, but also the time required for mining rules which will be filtered out is useless. Thus, we propose a new algorithm for mining serial episode rules, with a consequent temporally distant from the antecedent, to predict distant events, and with a small antecedent (in number of events and in time), to be able to predict events as soon as possible. This algorithm has a complexity equal to that of traditional algorithms.

An example of a required rule is presented below (from a sequence of annotated messages of blogs about finance issues, where each event includes a sentiment polarity): R: (interest rate, neutral) (credit, negative) (waiting, neutral) → (concurrence, negative); the antecedent occurs within 5 days, the gap between the antecedent and the consequent is 15 days.

The rest of this paper is organized as follows: section 2 presents related works about episode rules mining. Our algorithm is introduced in section 3, followed by experimental results in section 4. We conclude and provide some perspectives in section 5.

RELATED WORKS

We first start by introducing few concepts. Let I = {i 1 , i 2 , ..., i m } be a finite set of items. I t is the set of items that occur at a timestamp t referred to as an event. An event sequence S is an ordered list of events, S =< (t 1 , I t 1 ), (t 2 , I t 2 ), ..., (t n , I t n ) > with t 1 < t 2 < ... < t n , (see Figure 1). The serial episode P =< p 1 , p 2 , ..., p k > on I k is an ordered list of events. Its support, denoted to as supp(P), represents the number of occurrences of P, according to a frequency measure. P is said to be a frequent episode if supp(P) ≥ minsupp where minsupp is the predefined minimal threshold. An occurrence window of the episode P is a segment < I t s , ..., I t e > of the sequence, denoted as OW (S,t s ,t e ) which starts at timestamp t s and ends at timestamp t e , where P ⊆< I t s , ..., I t e >, p 1 ⊆ I t s and p k ⊆ I t e . It represents the interval that bounds the episode. Let P, Q be two episodes. An episode rule R : P → Q means that Q appears after P. The confidence of this episode rule is the probability to find Q after P: con f (P → Q) = supp(P.Q)/supp(P). The rule is said to be confident if its confidence exceeds a predefined threshold mincon f .

Winepi and Minepi are seminal episode mining algorithms [START_REF] Mannila | Discovery of frequent episodes in event sequences[END_REF], and are the basis of many other algorithms. To extract episodes, both algorithms start by extracting 1-tuple episodes (made up of one item), then iteratively seek larger ones by merging items on their right side. This approach is the one used by many other algorithms [START_REF] Laxman | A fast algorithm for finding frequent episodes in event streams[END_REF] [START_REF] Huang | Efficient mining of frequent episodes from complex sequences[END_REF].

In the frequent itemsets mining task (which includes episode mining), it is suggested that antimonotonicity is a common property that has to be respected by any frequency measure [START_REF] Agrawal | Mining association rules between sets of items in large databases[END_REF]. Several frequency measures for episode mining have been proposed. In [START_REF] Mannila | Discovery of frequent episodes in event sequences[END_REF], window-based and minimal occurrence-based frequency measures have been introduced through both Winepi and Minepi. Winepi evaluates the frequency of an episode as the number of windows of length w that contain the episode. Minepi evaluates the frequency of an episode as the number of minimal windows that contain the episode. A minimal window is a window such that no sub-window contains the episode. The non-overlapped occurrencebased frequency measure is proposed [START_REF] Laxman | A fast algorithm for finding frequent episodes in event streams[END_REF], where two occurrences of an episode are nonoverlapped if no item of one occurrence appears in between items of the other. It is shown that the nonoverlapped occurrence-based algorithms are much more efficient in terms of space and time needed. When mining serial episodes, additional constraints on the episodes may be imposed. Such as the span constraint [START_REF] Achar | Pattern-growth based frequent serial episode discovery[END_REF], that imposes an upper bound (of distance or time) between the first and last event in the occurrence of an episode. This constraint has been mainly introduced for complexity reasons. Another constraint is the gap constraint [START_REF] Méger | Constraint-based mining of episode rules and optimal window sizes[END_REF], which imposes an upper bound between successive events in the occurrence of an episode. If the extracted serial episodes have to represent causative chains then such constraints are important.

Traditional episode rules mining algorithms construct episode rules with a large antecedent (made up of many events) [START_REF] Pasquier | Discovering frequent closed itemsets for association rules[END_REF]. Discovering rules with a small antecedent was introduced for association rules, called "minimal association rules discovery" [START_REF] Rahal | Mining confident minimal rules with fixedconsequents[END_REF]. Minimal rules are considered as sufficient (as no new knowledge is given by larger ones). The constraint is that the consequents are fixed in advance by the user. Minimal rules have also been studied with the aim to reduce time and space complexity of the mining task, as well as to avoid redundancy in the resulting set of rules [START_REF] Neeraj | Overview of nonredundant association rule mining[END_REF]. These works focus on association rules; recall we want to form episode rules.

Mining episodes in an event sequence is a task which has received much attention [START_REF] Gan | Fast mining of nonderivable episode rules in complex sequences[END_REF]. In an event sequence each data element may contain several items (an event). In (Huang and 1 2 3 4 5 6 7 8 9 10 11 12 where the event sequence is encoded with frequent itemsets then serial episodes are mined. In [START_REF] Gan | Fast mining of nonderivable episode rules in complex sequences[END_REF], episodes are first extracted, then nonderivable episodes rules are formed (where no rule can be derived from another).

H A A B C I D F E A B C F E G F E t

THE PROPOSED ALGORITHM

Principle

Our goal is to form episode rules that can also be used to efficiently predict distant events.

To achieve this goal, the episode rules formed have to hold the characteristic that the consequent is temporally distant from the antecedent. Traditional algorithms are not designed to form such rules. Recall that they first form episodes from left to right by iteratively appending events temporally close to the episode being formed (the minimal occurrence), in the limit of the predefined span. Then, the episode rule is built by considering the last element(s) of the episode as consequent of the rule. In addition, when forming these episodes, it is impossible to know if the event being appended will be part of the consequent or not, so it is impossible to constrain its distance to other events while forming the episode. The only way to mine rules with a consequent distant from the antecedent is to mine all rules and then filter the occurrences that do respect this distance. Due to the limited span, this distance cannot be large.Thus, we propose to mine episode rules without any episode mining phase. To be able to constrain the distance between the antecedent and the consequent, we propose to determine the consequent early in the mining process. We think that, by determining the consequent at an early stage, the occurrence windows of an episode rule will be filtered early, thus the search space will be pruned, and no post-processing is required.

We aim also at predicting events early. We assume that the more the antecedent of a rule is small in number of events and in time, the earliest it ends, and the earliest the consequent can be predicted. Therefore, we propose to extract episode rules that have an additional characteristic: an antecedent as small as possible (in number of events and in time), which we call "minimal episode rules". For that, we apply the traditional minimal occurrence-based frequency measure.

The required characteristics force us to propose an approach totally different from traditional ones. We propose an episode rules mining algorithm, where the prefix (the first event) of a rule is first fixed. Then, the consequent is determined to constrain its distance from the prefix. Finally, the antecedent is completed.

Before explaining our algorithm, we present new concepts on which our algorithm relies.

Sub-windows of Win(S,t s , w): Let Win(S,t s , w) be a window in the sequence S of length w that starts at t s , with its first element containing the prefix of an episode rule (the first event) to be built. In order to mine episode rules with a distant consequent and a minimal antecedent, we propose to split this window into three sub-windows as follows (see Figure 2):

Win begin is a segment of Win(S,t s , w) of length w begin < w, starting at t s . Win begin can be viewed as an expiry time for the antecedent of an episode rule. It represents the span of the antecedent of an episode rule to guarantee that the antecedent occurs within a determined time.

Win end is a segment of Win(S,t s , w) of length w end < w, that ends at t s + w. Win end represents the time window of occurrence of the consequent.

Win between is the remaining sub-window of length w between , in which neither the antecedent nor the consequent can appear. Win between guarantees the temporal distance between the antecedent and the consequent of an episode rule. It represents a minimal gap between the antecedent and the consequent to guarantee that the consequent is far from the antecedent.

Steps of the algorithm

We present now the different steps of our algorithm.

Initialization

Our algorithm starts by an initialization phase, which reads the event sequence to extract all frequent events and their associated occurrence timestamps. An event represents a 1-tuple episode and will be denoted as P. Table 1 presents the list of 1-tuple episodes of the sequence S (Figure 1) and their associated occurrence windows for minsupp = 2 (see Figure 3 line 2). 1).

Consequent identification

A candidate consequent of an antecedent ant (here, ant corresponds to a unique element, the prefix) is chosen in the windows Win(S,t s , w) where: ant ⊆ I t s .

Recall we want to form episode rules with a consequent as far as possible from the antecedent. Thus, the candidate consequents are not searched in the entire window, they are searched only in Win end : where the farthest candidates are. We construct P end (ant), the ordered list of 1-tuple episodes that occur frequently in Win end from Win(S,t s , w) that contains ant.

Let p j ∈ P end (ant) be a candidate consequent of ant. The episode rule R : ant → p j is formed and its support is computed (see Figure 3 line 7). At this stage, the occurrence windows of the episode rule R are filtered to get minimal occurrences as well as to preserve the anti-monotonicity property. This filtering is done by counting only once the occurrence windows containing the same occurrence of the consequent p j . For example, let w begin = 2, w end = 2 and w = 6. The episode rule R : A → E has three occurrence windows: 7,12]). In the two first occurrence windows, the consequent E is a common 1-tuple episode which occurs at timestamp t 6 . Therefore, we consider only the interval [2, 6]. However, all occurrence windows are kept in memory, to be used to complete the antecedent. This enables not losing any of the interesting episode rules which could be missed if we kept only the minimal occurrences in memory.

([1, 6], [2, 6], [
Next, we compute the support of the corresponding episode rules (see Figure 3 line 7).

We define the support of a rule P → Q, referred to as supp end (P.Q), as the number of minimal occurrence windows computed above. It is different from the traditional supp(P.Q) as it considers only occurrence windows where P occurs in Win begin and Q occurs in Win end .

If R : ant → p j is not frequent, we consider that p j cannot be a consequent of ant. This iteration is stopped and the rule is discarded. There is no need to complete the antecedent of the rule R, as whatever the events that complete the antecedent are, the resulting rule will not be frequent. The algorithm will iterate on another consequent. If R : ant → p j is frequent, its confidence is computed. its confidence is computed.

We define the confidence of a rule P → Q (see Equation ( 1)) as the probability that the consequent occurs in Win end , given that P appears in Win begin .

con f (P → Q) = supp end (P.Q) supp(P) (1)
If the rule R : ant → p j is confident, this rule is added to the set of rules formed by the algorithm. It is minimal and has a consequent far from the antecedent; it fulfills our goal. If the rule is frequent but not confident, the antecedent of the rule R : ant → p j is completed (as in the next subsection).

For example, let w = 6, w begin = 2 and w end = 2.

For the episode rule R with prefix

A. P end (A) = [E, F, EF, A].
We first try to construct the episode R with the consequent E. Thus, for R : A → E, supp(R) = 2 and con f (R) = 2/3 = 0.67. For minsupp = 2 and mincon f = 0.7, R is frequent but not confident, so its antecedent has to be completed.

Antecedent completion

In this step, the antecedent, referred to as ant, is iteratively completed with 1-tuple episodes, placed on its right side in the limit of the predefined sub-window Win begin . At the first iteration: ant is a unique element (the prefix) (see Figure 3 line 15). Recall that we aim at forming rules having the last event of the antecedent as far as possible from the consequent, so as close as possible of the prefix. Thus, we construct P begin (ant), the ordered list of 1-tuple episodes that occur frequently after ant in the windows that starts with it: the 1-tuple episodes that occur in Win begin .

Similarly to the consequent identification step, the occurrence windows of the episode rule R are filtered to get minimal occurrences and to preserve the anti-monotonicity property. In addition, we apply the same support, confidence verifications as for the consequent identification.

To speed up the episode rules mining process we use a heuristic. We propose to order the list of candi-dates P begin (ant) in descending order of the number of Win begin in which the candidates appear. We assume that this number is highly correlated with the support of the corresponding episode rules. So, in the traversal of this list, when we observe that candidates tend to form infrequent episode rules (several consecutive 1-tuple episodes lead to infrequent episode rules), we stop the traversal. We consider that the remaining candidates in this list will lead to infrequent episode rules. This heuristic is used to reduce the number of iterations. Although this heuristic may discard interesting rules, it allows to reduce the iterations thus to increase the size of the span of the rule (Win).

For example, let minsupp = 2 and mincon f = 0.7, for R :

A → E, P begin (A) = [B,C, A].
The antecedent of R is completed with B and forms the episode rule R : A, B → E. Thus, supp(R) = 2 and con f (R) = 2/2 = 1. The episode rule R is now confident, the phase of completing its antecedent is stopped.

Temporal confidence

Win between has been introduced so as to guarantee that the consequent of an episode rule occurs in Win end , after a w begin + w between temporal distance from the prefix of the episode rule. However, given an occurrence of the antecedent, the consequent may also occur closer to the antecedent, in the window Win between , which should affect the confidence of the rule. This information may be important in some applications.

In the example of social networks, the prediction of a negative event allows the company to act so as to prevent its occurrence. So, it is important to mine rules with a consequent that never occurs in Win between . Indeed, predicting a consequent at a given distance, which may appear closer is useless, even dangerous. Consequently, we have to take into consideration the occurrence of the consequent in Win between . We introduce a new measure, the temporal confidence, which represents the probability that the consequent occurs in Win end and only in it. For an episode rule R : P → Q, this measure takes into account the support of P.Q when Q occurs only in Win end (and not in Win between ), denoted as supp end (P.Q). The temporal confidence is defined as follows:

con f t (P → Q) = supp end (P.Q) supp(P.Q) (2)
con f t (R) = 1 if for each occurrence of the consequent of R in Win end , no occurrence of the consequent is found in Win between . The temporal confidence of rules from the previous section is computed and the rules with a temporal confidence above mincon f t are kept. For example, let w = 6, w begin = 2 and w end = 2, the temporal confidence of the frequent confident episode rule R : A → E depends on the occurrences of E in Win between which is equal to 1 (E appears in the timestamp t 10 in Win between ). Thus, con f t (R : A → E) = 1/2 = 0.5. For mincon f t = 0.5, R is temporally confident and is a rule formed by our algorithm.

EXPERIMENTAL RESULTS

In this section, we evaluate our algorithm through the study of the characteristics of the episode rules formed, as well as its performance in a prediction task.

Dataset

The dataset we use is made up of 27, 612 messages extracted from blogs about finance. Messages are an- In this dataset, messages are annotated with 4.8 items on average, ranging from 1 to 50 and the median is 4. There are about 4, 000 distinct annotations (items), with an average frequency of 88.5. 1, 981 of these items have a frequency equal to 1. These items will be automatically filtered out in the initialization phase of our algorithm.

Characteristics of the resulting rules

Initialization phase

In this phase, frequent 1-tuple episodes are extracted (1-tuple episode is made up of one or more items). We fix minsupp for these 1-tuple episodes to 30. This phase results in 652 frequent 1-tuple episodes. Table 2 shows that the 1-tuple episodes are made up of one to three items only. 1-tuple episodes are mainly made up of one item (76% of them), which also have a high support (on average 149.7).

In order to study the episode rules formed, we make vary minsupp, mincon f and w between one at a time, while fixing others.

Making vary minsupp

In Figure 4, we make vary minsupp from 10 to 35 to study the number of rules. Other parameters remain fixed: mincon f = 0.4, w = 40 (with w begin = 20, w between = 10 and w end = 10). As expected, the smaller minsupp, the higher the number of rules. When minsupp is fixed to 10, the number of rules is high (about 10 4 ). Recall that the number of 1tuple episodes is only 652. However, the number of rules is dramatically decreased when minsupp is increased: only 1, 200 when minsupp = 15 and 250 when minsupp = 20. These low values are due to the small average frequency of 1-tuple episodes (about 88). In addition, these rules represent a temporal dependence between the antecedent and the consequent which, in these experiments, is at least w between = 10. This may explain the low number of rules. A thorough study shows that their average confidence increases with minsupp.

Making vary mincon f

We make vary mincon f from 0.1 to 0.5. Figure 5 presents the number of episode rules according to the value of mincon f (minsupp = 20, w = 40 (w begin = 20, w between = 10 and w end = 10)). The number of rules is particularly high with mincon f = 0.2. This is explained by the way our rules are formed. When the confidence of a rule does not exceed mincon f , its antecedent is extended. Given an antecedent ant, events in P begin (ant) (up to 652) are appended to it, resulting in a large number of candidate rules. Some of them are confident, which explains the increase of the number of rules. When mincon f exceeds 0.2, the number of rules decreases.

Table 3 presents the length of the antecedent of the rules (in number of events), according to mincon f . The maximum length of an antecedent is three. Thus, our algorithm, forms rules with a small antecedent, which was one of its goals. The average length of the antecedent increases with mincon f : when mincon f = 0.1, most of the rules have an antecedent of length 1, whereas when mincon f = 0.3, most of the rules have an antecedent of length 2. This was expected, as minimal antecedents are searched. Indeed, when a frequent rule has a confidence below mincon f , its antecedent is extended, till it is confident or not frequent. So, the higher mincon f , the larger the antecedents. A thorough study shows that the average length of occurrence windows of the antecedents is 8 timestamps (for antecedents of length 2 and 3), which is smaller than the span of the antecedent (w begin ). We conclude that our algorithm succeeds in forming rules with a distant consequent, a small antecedent (in length and time) and a relatively high confidence.

Making vary w between

We now focus on the number of rules formed, according to w between (w and w begin remain fixed), presented in Figure 6 where minsupp = 20 and mincon f = 0.4. Two values of w are studied: w = 40 (with w begin = 20) and w = 100 (w begin = 20). Notice that the cases w between = 0 represent similar cases than state of the art. We note that the larger w between , the smaller the number of rules. Two reasons may explain this decrease. First, when w between increases, w end (the window in which the consequent is searched) decreases, as well as the number of consequents studied. Second, the larger w between , the more distant the consequent, thus the lower the probability of having a dependence between the antecedent and the consequent. However, even with a large value of w between , some rules are formed: 210 rules when w between = 70. We conclude that there is actually a temporal dependence between messages in blogs. When the minimal distance between the antecedent and the consequence is 50, more than 140k confident rules are formed: there is a strong dependence between messages with such a distance. For example, when w = 100 an episode rule: (price, positive), (information, positive) → (buy, positive), means that when someone talks about the price of an article then asks for information, he/she will buy this article after some time. Thus, we have time to recommend him similar articles or to propose to him a credit to buy it.

Influence of mincon f t : In this section we study the temporal confidence of the resulting rules. Table 4 presents the evolution of the temporal confidence according to w between . We remark that the smaller w between , the higher the temporal confidence: the consequent does rarely occur between the antecedent and the consequent when the gap between them is small, which was expected. When w between = 30, the average temporal confidence is 0.6, which is quite high. A thorough study shows that among 2.7 • 10 5 rules (see Figure 6), about 40 ones have a temporal confidence equal to 1 (the consequent never occurs in Win between ) and 1,400 rules have a temporal confidence higher than 0.9 (the consequent appears in Win between in less than 10% of the cases). This shows that in this dataset there is a strong temporal dependence between events, and that some events are interdependent at a distance of 30.

So, when exploiting the temporal confidence as a filter, a great number of rules remain.

Table 4: w between vs. Temporal confidence (con f t ).

w between min-con f t max-con f t mean-con f t median-con f t 70 0 0.5 0.2 0.2 50 0.2 0.9 0.4 0.4 30 0.3 1 0.6 0.6 10 0.2 1 0.9 0.9

Performance

In this section we focus on the accuracy of the rules formed when they are used to predict events, and we perform a comparison of the rules with those of a traditional algorithm.

We evaluate the accuracy of our algorithm with the traditional recall and precision measures. The episode rules are trained on the first 75% messages and are tested on the 25% of messages left. Figure 7 presents the resulting precision and recall at 20. We fix minsupp = 20, mincon f = 0.4, w = 100 and w begin = 20. We make vary w between from 10 to 70. First of all, mention that two precision and recall values with two different values of w between are not directly comparable as they are not computed on the same data (the windows Win end , on which they are computed vary in size). Both precision and recall curves decrease as w between increases. This was expected as the number of rules decreases. When w between = 70 (and w end = 10), both precision and recall values are quite low. This was expected as the rules aim at predicting events distant to at least 70, in an occurrence window of length w end =10. The prior probability of predicting events accurately is low. Let us now consider w between = 30, as in the previous section. We can see that both precision and recall values are quite high. When an event is predicted, in 37% of the cases, it actually occurs and events that occur in the sequence are predicted by our rules in 70% of the cases. Comparison with Minepi: Contrary to traditional algorithms, our algorithm forces a minimum gap of length w between between the antecedent and the consequent of an episode rule. Our algorithm can be compared to traditional algorithms when w between = 0. Since we apply the minimal occurrence-based frequency, we choose to compare it to the well-known Minepi [START_REF] Mannila | Discovery of frequent episodes in event sequences[END_REF]. For minsupp = 20, mincon f = 0.4, w = 40, Minepi forms more than 136, 000 episode rules, whereas our algorithm (when w between = 0) extracts about 40, 000 episode rules (70% less). This decrease is due to two reasons. First, the constraint about the position of consequent of the episode rules from our algorithm (in this case the distance between the antecedent and the consequent is at least 20, even if w between = 0), makes the number of rules resulting from our algorithm lower (also their support is lower). Second, our algorithm aims at forming minimal rules, thus few rules with a large antecedent are formed. We remark that 25% of the rules extracted by Minepi have an antecedent larger or equal to 3, whereas this rate is only 1.8% for our algorithm.

Here an example of an episode rule extracted by both our algorithm (w between = 0) and by Minepi: (credit, positive), (consultant, positive) → (loan subscription, positive).

Here is a rule that has not been extracted by our algorithm, as it does not satisfy the desired characteristics of episode rule (minimal antecedent): (consultant, neutral), (interest rate, positive) → (request interest rate 0, positive), where the antecedent occurs in 5 timetstamps and consequent occurs in the 7 th timestamp. This rule is useful in traditional cases of event prediction (prediction of close events). However, it does not fit our objective of early prediction of distant events, as the antecedent is so long both in time and in number and the consequent is too close to the antecedent.

Concerning the running time, our algorithm runs 5 times faster than Minepi when w = 100, and 4 times faster when w = 40. This decrease is due to two factors. The first one is related to the consequent, which is fixed at an early stage of the algorithm and which allows to filter infrequent rules early in the process. The second one is due to the fact that our algorithm mines rules with a minimal antecedent, which avoids some iterations once a confident rule is found.

CONCLUSION

In this paper, we propose an algorithm that mines episode rules, in order to predict distant events. To achieve our goal, the algorithm mines serial episode rules with distant consequent. We determine several characteristics of the episode rules formed: minimal antecedent, and a consequent temporally distant from the antecedent. A new confidence measure, the temporal confidence, is proposed to evaluate the confidence on distant consequents. Our algorithm is evaluated on an event sequence of annotated social networks messages. We show that our algorithm is efficient in extracting episode rules with the desired characteristics and in predicting distant events.

Since we use data from social networks, we aim to use multi-thread sequences. This means that we construct a sequence for each thread of messages: user messages thread, topic messages thread and discussion thread, etc. and the algorithm is run on each one. Using multi-thread sequences allows to build more di-verse episode rules which are all together more significant. The presence of a rule in several threads will increase its confidence.

Acknowledgements.

This research is supported by Crédit Agricole S.A.

Figure 1 :

 1 Figure 1: Example of an event sequence S

Figure 3 :

 3 Figure 3: Episode rules mining algorithm

FigureFigure 5 :

 5 Figure 4: number of rules (10 4 ) vs. minsupp

Figure 6 :

 6 Figure 6: number of rules (10 5 ) vs. w between

Figure 7

 7 Figure 7: precision, recall vs. w between

Table 1 :

 1 1-tuple episodes of S.Episode rules are built iteratively by first fixing the prefix (the first event of the antecedent). The antecedent is denoted as ant. In the encoded sequence, each 1-tuple episode p is viewed as a prefix of the antecedent of an episode rule to be built. Once the prefix of an episode rule R is fixed, its occurrence windows OW (S,t s ,t

	1-tuple episode p List of occurrence windows
	A	[1,1], [2,2], [7,7]
	B	[3,3], [8,8]
	C	[4,4], [9,9]
	E	[6,6], [10,10], [12,12]
	F	[6,6], [10,10], [12,12]
	EF	[6,6], [10,10], [12,12]
	3.2.2 Prefix identification

e ) are known. For example, let minsupp = 2, A can be considered as a prefix of an episode rule. The list of occurrence windows of

A = ([1, 1], [2, 2], [7, 7]) (see Table

Table 2 :

 2 1-tuple episodes : length and support.

	1-tuple episode length (#items)	number(%)	support min max mean median
	1	498(76.4) 30 797 149.7 89.5
	2	147(22, 5) 30 376 71.5	55
	3	7(1)	32 54 44.4	46
	notated using the Temis 1 software. Each message is
	represented by its corresponding set of annotations
	(items). For example, the message: "Not only my
	bank propose the best, but also online banks, so how
	to optimize my savings? accounts in other banks? life
	insurance? i need more info." is annotated with the
	following items : {(online banks, negative), (savings,
	neutral), (life insurance, neutral ), (needs, neutral)} ,
	where each item is associated with an opinion degree.

Table 3 :

 3 Antecedent length when making vary mincon f .

	mincon f #Rules	%Rules ant 1 ant 2 ant 3
	0.1	16,850 57.2 42.84 0
	0.2	31,972 4.2 95.6 0.2
	0.3	5,072 0.9 97.5 1.6
	0.4	251	0.8 90.9 8.3
	0.5	9	0	100	0

http://www.temis.com