Herwig Wendt

Stéphane G Roux

Stéphane Jaffard
email: jaffard@univ-paris12.fr

Patrice Abry

† P Abry

Wavelet Leaders and Bootstrap for Multifractal Analysis of Images

Keywords: Multifractal, Wavelet Leaders, Image regularity, Bootstrap, Confidence intervals

Multifractal analysis is considered a promising tool for image processing, notably for texture characterization. However, practical operational estimation procedures based on a theoretically well established multifractal analysis are still lacking for image (as opposed to signal) processing. In the present contribution, a wavelet Leader based multifractal analysis, known to be theoretically strongly grounded, is described and assessed for 2D functions (images). By means of Monte Carlo simulations conducted over both self-similar and multiplicative cascade synthetic images, it is shown here to benefit from much better practical estimation performance than those obtained from a 2D discrete wavelet transform coefficient analysis. Furthermore, this is complemented by the original analysis and design of procedures aiming at practically assessing and handling the theoretical function space embedding requirements faced by multifractal analysis and certain image analysis methods. In addition, a bootstrap based statistical approach developed in the wavelet domain is proposed and shown to enable the practical computation of accurate confidence intervals for multifractal attributes from a given image. It is based on an original joint time and scale block non parametric bootstrap scheme. Performance are assessed by Monte Carlo simulations. Finally, the use and relevance of the proposed wavelet Leader and bootstrap based tools are illustrated at work on real-world images.

Introduction: Multifractal analysis for image processing

Nowadays, in a large number of applications of very different natures, the data collected by sensors for analysis consist of images, i.e., they are naturally bi-dimensional signals. This is mostly due to the recent and significant progresses achieved in digital sensor, fast rate and high resolution camera and video camera design. For a number of these applications the corresponding statistical analysis of the images amounts to performing texture characterization. This is the case notably for clouds or rainfalls analyses in geophysics [START_REF] Gupta | A statistical analysis of mesoscale rainfall as a random cascade[END_REF][START_REF] Schertzer | Multifractal cascade dynamics and turbulent intermittency[END_REF][START_REF] Roux | A wavelet-based method for multifractal image analysis. III. Applications to high-resolution satellite images of cloud structure[END_REF], bio-medical diagnosis for human body rhythms or structure (bones, tissues, mammography,. . . ) [START_REF] Jennane | Fractal analysis of bone x-ray computed microscopy projections[END_REF][START_REF] Lowen | Fractal-Based Point Processes[END_REF][START_REF] Richard | A statistical methodology for testing the anisotropy of brownian textures with an application to full-field digital mammography[END_REF], universe or galaxy structures in astronomy [START_REF] Starck | Image Processing and Data Analysis: The Multiscale Approach[END_REF], growth phenomena in physics [START_REF] Decoster | A wavelet-based method for multifractal image analysis. II. Applications to synthetic multifractal rough surfaces[END_REF][START_REF] Ponson | Anisotropic self-affine properties of experimental fracture surface[END_REF] or texture classification in computer vision [10,[START_REF] Xu | A projective invariant for texture[END_REF], to name but a few examples.

Texture characterization is now often envisaged by measuring the fluctuations (with respect to space) of the regularity of the amplitude of the image. There is an increasing number of research articles suggesting that regularity characterization should be conducted within the mathematical framework of multifractal analysis. In practice, multifractal analysis is quasisystematically performed using the coefficients of continuous [START_REF] Bacry | Singularity spectrum of fractal signals from wavelet analysis: Exact results[END_REF] or discrete [START_REF] Abry | Multiscale nature of network traffic[END_REF][START_REF] Jaffard | Wavelet leaders in multifractal analysis[END_REF] wavelet transforms. The regularity fluctuations are then inferred from the behavior in the limit of fine scales of the (sample mean estimators of the) moments of order q. However, surveying the literature related to applications, one realizes that image multifractal analysis suffers from two major limitations: It remains either incomplete, as it is simplified to the computation of Fourier spectra (via classical standard spectral estimation) or of wavelet coefficients variograms (regression of the variance against the scales), i.e., for both cases, of a second order statistical quantity only (e.g. [START_REF] Sharifi-Salamantian | A robust index for spatial heterogeneity in breast cancer[END_REF]); or it is reduced to the analyses of collections of 1D slices of data instead of real 2D images (e.g. [START_REF] Benhamou | Fractal analysis of radiographic trabecular bone texture and bone mineral density: two complementary parameters related to osteoporotic fractures[END_REF]).

The reasons for this are mostly associated to the theoretical requirement that a complete multifractal analysis involves, namely a variety of negative statistical orders q as well as positive ones [START_REF] Jaffard | Wavelet techniques in multifractal analysis[END_REF][START_REF] Riedi | Multifractal processes[END_REF]. However, wavelet coefficients, by nature, consist of quantities that mostly concentrate around 0, rendering the numerical computation of negative q moments extremely unstable or even theoretically infinite. To overcome this limitation, using the Modulus Maxima of a Continuous Wavelet Transform (MMWT) has been proposed for the 1D case [START_REF] Bacry | Singularity spectrum of fractal signals from wavelet analysis: Exact results[END_REF]. This method relies on the determination of the coefficients of a Continuous Wavelet Transform (CWT), from which a skeleton, consisting of maxima along scales lines, is extracted [START_REF] Mallat | A Wavelet Tour of Signal Processing[END_REF]. The wavelet coefficients living on this skeleton are then involved in the computation of the fluctuations of the image regularity. For details, the reader is referred to [START_REF] Bacry | Singularity spectrum of fractal signals from wavelet analysis: Exact results[END_REF]. This technique has been further extended to image analysis, at the price though of significant computational (2D-CWT + 2D-Skeleton) and conceptual complexities (maxima lines become maxima manifolds) [START_REF] Arneodo | A wavelet-based method for multifractal image analysis: from theoretical concepts to experimental applications[END_REF]. Therefore, it remains hardly ever used for images (see, a contrario, [START_REF] Arneodo | A wavelet-based method for multifractal image analysis: from theoretical concepts to experimental applications[END_REF]). In addition to the practical difficulties related to its implementation, the MMWT approach, be it 1D or 2D, despite its showing satisfactory experimental performance, is still lacking a theoretical mathematical support. Recently, an alternative approach has been proposed [START_REF] Jaffard | Wavelet leaders in multifractal analysis[END_REF][START_REF] Jaffard | Wavelet techniques in multifractal analysis[END_REF]: The Wavelet Leader (WL) Multifractal Analysis (MA). This method is theoretically backed up by a strong mathematical framework [START_REF] Jaffard | Wavelet leaders in multifractal analysis[END_REF][START_REF] Jaffard | Wavelet techniques in multifractal analysis[END_REF][START_REF] Lashermes | Wavelet leader based multifractal analysis[END_REF][START_REF] Jaffard | The contribution of wavelets in multifractal analysis[END_REF]. Also, its being defined from an (orthonormal) Discrete Wavelet Transform (DWT) enables reasonably easily its theoretical and practical extensions to higher dimension [START_REF] Jaffard | Wavelet leaders in multifractal analysis[END_REF][START_REF] Lashermes | Wavelet leader based multifractal analysis[END_REF]. It has been shown to be one of the best available tools for (1D) signals, enabling a relevant and general multifractal analysis of empirical data, with both excellent theoretical properties and statistical performance [START_REF] Wendt | Bootstrap for empirical multifractal analysis[END_REF].

Therefore, elaborating on theoretical results proven for functions in R n and recently published in [START_REF] Jaffard | Wavelet leaders in multifractal analysis[END_REF][START_REF] Jaffard | Wavelet techniques in multifractal analysis[END_REF][START_REF] Lashermes | Wavelet leader based multifractal analysis[END_REF][START_REF] Jaffard | The contribution of wavelets in multifractal analysis[END_REF], the first goal of the present contribution is to provide practitioners with an explicit formulation of a 2D wavelet Leader based multifractal analysis that can actually be applied to real-world images. To this end, Section 2 reviews the key notions underlying multifractal analysis, introduces the 2D wavelet Leader multifractal analysis, and details the corresponding practical estimation procedures. In addition, the validity and practical relevance of these 2D multifractal analysis procedures are assessed in the present work by numerical simulation studies of their estimation performance in Sections 6.1 and 6.2. The performance of the proposed procedures are shown to compare very favorably to those obtained from wavelet coefficients. The synthetic selfsimilar and multiplicative cascade images used in Monte Carlo simulations are described in Section 5.

The second goal of the present contribution is to further complement the theoretical understanding of a comprehensive use of the wavelet Leader multifractal analysis. For this, we address a number of theoretical issues, related to which function space data are embedded in, and constituting mandatory prerequisites in the derivation of the wavelet Leader based multifractal analysis, as well as for any other procedure aiming at performing a multifractal analysis of images or signals. Notably, the wavelet Leader multifractal analysis is well-defined for bounded functions only. However, a digital image consisting of an intensity local average, can naturally be seen as the approximation -at a given resolution level -of a positive measure, and there is hence no a priori guarantee for it to belong to the space of bounded functions. Section 3 addresses such issues and proposes tools for an a priori analysis aiming at enabling the characterization of which function space the data belong to, and solutions to handle the bounded function requirement. This part has pratical implications for the validation of many image processing models, where an a priori assumption is made on the function spaces that contain the image. The proposed solutions are supported by mathematical proofs, whose founding arguments are reported in Section 3 (yet in Section 3.3 detailed versions of some proofs are beyond the scope of the present article and will be published independently in a theoretically oriented journal, see [START_REF] Jaffard | Wavelet analysis of multifractal measure[END_REF]). The performance of the proposed procedures are assessed numerically. Also, their practical relevance and usefulness are supported through the examination of the large reference texture image database used in [START_REF] Xu | A projective invariant for texture[END_REF].

In the use of multifractal analysis in actual applications, practitioners are potentially as much interested in the confidence that should be granted to the estimates as in the value of the estimates themselves. However, the (asymptotic) statistical performance of the estimation procedures for multifractal attributes -and hence the design of confidence intervals -is beyond analytic derivation, mostly because of the involved statistical properties of multifractal images. To overcome such limitations, we recently proposed, for 1D signal, the use of non parametric bootstrap techniques [START_REF] Lahiri | Resampling Methods for Dependent Data[END_REF][START_REF] Zoubir | On confidence intervals for the coherence function[END_REF], applied in the wavelet domain [START_REF] Wendt | Bootstrap for empirical multifractal analysis[END_REF][START_REF] Wendt | Multifractality tests using bootstrapped wavelet leaders[END_REF]. Therefore, the third goal of the present contribution is to define an original bootstrap scheme that matches the 2D wavelet Leader multifractal analysis and to validate its performance in producing confidence intervals for multifractal attributes from a single image. This 2D scheme differs from the 1D one proposed earlier insofar as it consists of a joint time and scale block bootstrap procedure, designed to better account for the naturally joint time and scale correlation structure of wavelet coefficients and wavelet Leaders. This bootstrap approach is detailed in Section 4. The statistical performance of this bootstrap scheme, its relevance and usefulness are assessed numerically and discussed in Section 6.3.

The relevance and practical use of the multifractal Leader and bootstrap based analysis framework proposed here is illustrated at work on a real-world image (cf. Section 6.4). Also, comments and discussions related to the general applicability to images conditions of the proposed Leader and bootstrap based multifractal analysis procedures are proposed in Section 6.5.

Multifractal formalism and wavelet Leaders

This section aims at providing to practitioners in image processing actual procedures to conduct relevant multifractal analysis of real-world images. It is based on general theoretical results published in [START_REF] Jaffard | Wavelet leaders in multifractal analysis[END_REF][START_REF] Jaffard | Wavelet techniques in multifractal analysis[END_REF][START_REF] Lashermes | Wavelet leader based multifractal analysis[END_REF][START_REF] Jaffard | The contribution of wavelets in multifractal analysis[END_REF] and makes explicit their formulation and use for images.

Multifractal spectrum

Let X(t) : t ∈ R d → X ∈ R denote the function to be analyzed, assumed to belong to the space of bounded functions. For an image, d = 2 and the variable t = (x 1 , x 2 ) is commonly referred to as the space variable. Multifractal analysis consists in describing the local regularity of X(t), around t 0 by comparing X(t 0 ) against a local power law behavior: Hausdorff dimension, the multifractal spectrum is confined to 0 ≤ D(h) ≤ d. By convention, D(h) = -∞ for the Hölder exponents that are not present in X. For theoretical introductions to multifractal analysis, the reader is referred to e.g., [START_REF] Jaffard | Wavelet techniques in multifractal analysis[END_REF][START_REF] Riedi | Multifractal processes[END_REF]. If X is a bounded function, its Hölder exponent is necessarily nonnegative, so that D(h) = -∞, ∀h < 0. For ease of exposition, we hereafter only consider images that are uniformly singular, which implies that the Hölder exponent is bounded from above, see [START_REF] Jaffard | The contribution of wavelets in multifractal analysis[END_REF]. From a practical point of view, multifractal analysis consists in estimating D(h) from the finite size image under analysis. This is commonly achieved via the socalled multifractal formalisms. The wavelet Leader multifractal analysis is defined in the next section.

|X(t) -P t0 (t)| ≤ C|t -t 0 | α ,

Wavelet Leader Multifractal Analysis

Wavelet coefficients.

Let H 0 (k), G 0 (k) denote the (lowpass and high-pass, respectively) quadrature mirror filters defining a 1D orthonormal DWT. The underlying mother wavelet possesses N ψ ≥ 1 vanishing moment(s). A 2D orthonormal DWT can be defined via the use of 4 bi-dimensionel filters G (m) (k 1 , k 2 ), m = 0, 1, 2, 3 obtained as products of H 0 and G 0 . By convention, G (0) (k 1 , k 2 ) = H 0 (k 1 )H 0 (k 2 ) corresponds to the 2D low pass filter providing a lower approximation, while G (m) , m = 1, 2, 3, correspond to the high pass filters yielding the wavelet coefficients. Let {X(k 1 , k 2 ), k 1 = 1, ..., N 1 , k 2 = 1, ..., N 2 } denote the 2D grey level digitized image to be analyzed, regarded as a sampled version of X(x 1 , x 2 ). The 2D wavelet coefficients D (m) X (j, k 1 , k 2 ), m = 1, 2, 3, are obtained by, first (j = 1), convolution of G (m) (k 1 , k 2 ), m = 0, 1, 2, 3 with {X(k 1 , k 2 )} and down-sampling, and then (j ≥ 2) iterative convolution of G (m) (k 1 , k 2 ), m = 0, 1, 2, 3 with D (0) X (j -1, k 1 , k 2 ) and down-sampling. For detailed introduction to continuous time and discrete time 1D or 2D WT, the reader is referred to, e.g., the tutorial books [START_REF] Mallat | A Wavelet Tour of Signal Processing[END_REF][START_REF] Antoine | Two-Dimensional Wavelets and their Relatives[END_REF].

For scaling analysis, it is more suitable (cf. [START_REF] Bacry | Singularity spectrum of fractal signals from wavelet analysis: Exact results[END_REF][START_REF] Abry | Multiscale nature of network traffic[END_REF]) to renormalize the standard L 2 -norm wavelet coefficients according to a L 1 -norm:

d (m) X (j, k 1 , k 2 ) = 2 -jd/2 D (m) X (j, k 1 , k 2 ).
Wavelet Leaders. Let us now further assume that the H 0 (k), G 0 (k) have finite impulse responses, and introduce a dyadic indexing of squares as:

λ j,k1,k2 = [k 1 2 j , (k 1 + 1)2 j ), [k 2 2 j , (k 2 + 1)2 j ) .
The union of nine such neighbor intervals is denoted as:

3 2 λ j,k1,k2 = n1,n2={-1,0,1} λ j,k1+n1,k2+n2 .
Wavelet Leaders are theoretically defined by

L X (j, k 1 , k 2 ) = sup m=1,2,3, λ ′ ⊂3 2 λ j,k 1 ,k 2 |d (m) X (λ ′ )|, (1)
(see [START_REF] Jaffard | Wavelet techniques in multifractal analysis[END_REF]). This definition is illustrated in Fig. 1 and practically means that the Leader L X (j, k 1 , k 2 ) is obtained as the

X (j ′ , k ′ 1 , k ′ 2 )|, m = 1, • • • , 3 (
red, green and blue dots) within a spatial neighborhood of (k 1 , k 2 ) and within all finer scale 2 j ′ ≤ 2 j (red volume). The wavelet coefficients over which the supremum is taken are marked by fat dots. largest wavelet coefficient amongst those, |d

(m) X (j ′ , k ′ 1 , k ′ 2 )|, m = 1, 2, 3, existing in a (narrow) spatial neighborhood of (k 1 , k 2 ), at any finer scale 2 j ′ ≤ 2 j .
The central result underlying the use of wavelet Leaders for multifractal analysis lies in them accurately measuring local Hölder exponents (cf. [START_REF] Jaffard | Wavelet techniques in multifractal analysis[END_REF] for the theoretical proof): If X has Hölder exponent h(t 0 ) ≥ 0 in t 0 , then, on condition that N ψ > h(t 0 ) and, if 2 -j k is the dyadic point of scale j which is closest to t 0 , then

L X (j, k) ∼ 2 j →0 2 jh(t0) , (2)
where X a ∼ a→0 Y a means that both quantities share the same lim inf on a log-log scale.

Structure functions and scaling. Let us now form the structure functions, i.e., spatial averages of (the q-th order of) the Leaders at a given scale 2 j :

S(2 j , q) = 1 n j k1,k2 L X (j, k 1 , k 2 ) q . (3)
The scaling function ζ(q) is then defined as

ζ(q) = lim inf 2 j →0 ln S(2 j , q) ln 2 j . (4)
In practice, this definition means that the S(2 j , q) exhibit power law behaviors with respect to the analysis scale 2 j , in the limit of small scales 2 j → 0 [14]: S(2 j , q) ∼ λ q 2 jζ ( q) when j → -∞.

The ζ(q) are hence often termed the scaling exponents. This power law behavior establishes a clear and deep connection between the concepts of scale invariance and multifractal analysis. Eq. ( 4) also implies that the scaling function necessarily is a concave function of q [START_REF] Jaffard | Wavelet leaders in multifractal analysis[END_REF].

Multifractal formalism. Let L(h) be the Legendre spectrum defined from the scaling function through a Legendre transform:

L(h) = inf q∈R (d + qh -ζ(q)). (6)
Under a mild uniform regularity condition for X (uniform Hölder, cf. Section 3.2 below), it can be shown [START_REF] Jaffard | Wavelet leaders in multifractal analysis[END_REF][START_REF] Jaffard | Wavelet techniques in multifractal analysis[END_REF] that the Legendre spectrum and the multifractal spectrum are simply related : D(h) ≤ L(h). There exist (large) classes of functions and processes X(t) for which this inequality turns to an equality, referred to as the wavelet Leader based multifractal formalism:

D(h) = L(h). (7)
Notably, for the self-similar and multiplicative cascade processes used in the present contribution to evaluate the statistical properties of the proposed estimation procedures (cf. Section 5), the wavelet Leader multifractal formalism holds. An obvious necessary condition for it to hold is that the spectrum D(h) of X(t) is a concave function. In practice, one cannot have access to the spectrum of singularities of a real-life signal; therefore the purpose of practical multifractal analysis focuses on the estimation of the Legendre spectrum L(h).

Cumulant expansion.

Let us now assume that the scaling exponents ζ(q) are a smooth function of q around 0; then, we can consider a Taylor expansion of this scaling exponents at q = 0:

ζ(q) ≃ p≥1 c p q p /p ! (8)
It has been shown in [START_REF] Castaing | Log-similarity for turbulent flows[END_REF][START_REF] Delour | Intermittency of 1d velocity spatial profiles in turbulence: A magnitude cumulant analysis, The Euro[END_REF] that the log-cumulants c p are related to the cumulant of order p, C(2 j , p), of ln L X (j, k) as:

C(2 j , p) = c 0,p + c p ln 2 j . (9)
The concavity of ζ(q) implies that c 2 ≤ 0. Recently, we showed (detailed calculations reported in [START_REF] Wendt | Contributions of wavelet Leaders and bootstrap to multifractal analysis: Images, estimation performance, dependence structure and vanishing moments. Confidence intervals and hypothesis tests[END_REF]), that Eq. ( 8) can be cast into a polynomial expansion of the spectrum as:

L(h) = d + c 2 2! h -c 1 c 2 2 + -c 3 3! h -c 1 c 2 3 + + -c 4 + 3c 2 3 /c 2 4! h -c 1 c 2 4 + • • • (10)
Also, it indicates that c 1 corresponds to the location of the maximum of L(h), c 2 to its width, while c 3 is an asymmetry parameter.

In applications, practitioners cannot compute the entire functions ζ(q) or L(h) and instead often concentrate on a small number of well chosen representative multifractal attributes. The c p are candidates of particular interest and one usually computes c 1 and c 2 and, when data permit, c 3 . This is further discussed in Section 6 for the specific context of image analysis.

Estimation procedures

For images, estimation of the scaling or multifractal attributes can be conducted as in the 1D case, mutatis mutandis. Estimation procedures have been described at length in [START_REF] Wendt | Bootstrap for empirical multifractal analysis[END_REF][START_REF] Wendt | Multifractality tests using bootstrapped wavelet leaders[END_REF] and are only briefly outlined here. As suggested by Eqs. ( 5) and ( 8) above, the estimation of the ζ(q) and of the c p can be performed by means of linear regressions:

ζ(q) = j2 j=j1 w j log 2 S(2 j , q), (11)
ĉp = (log 2 e) • j2 j=j1 w j Ĉ(2 j , p), (12)
where the estimates Ĉ(2 j , p) for the cumulants of ln L X (j, •) are obtained from standard sample cumulant estimators. The weights w j have to satisfy the constraints j2 j1 jw j ≡ 1 and j2 j1 w j ≡ 0 and can be expressed as

w j = b j V0j-V1 V0V2-V 2 1 with V i = j2 j1 j i b j , i = 0, 1, 2.
The freely selectable positive numbers b j reflect the confidence granted to each Ĉ(2 j , p) or log 2 S(2 j , q). To estimate L(h), it has been proposed [START_REF] Chhabra | Direct determination of the singularity spectrum and its application to fully developed turbulence[END_REF] to handle practically the Legendre transform via its parametric form, h(q), L(q). Such a formulation efficiently matches the bootstrap requirements for confidence interval derivation (cf. Section 4) and is hence used here. Tedious yet simple calculations lead to the following estimators: [START_REF] Abry | Multiscale nature of network traffic[END_REF] where

ĥ(q) = j2 j=j1 w j V (2 j , q), L(q) = j2 j=j1 w j U (2 j , q),
V (j, q) = nj k=1 L q X (j, k) log 2 L X (j, k) n j S(j, q) , and
U (j, q) = nj k=1 L q X (j, k)(log 2 L q X (j, k) -log 2 n j S(j, q))
n j S(j, q) + log 2 n j .

Wavelet coefficients versus Leaders

Empirical scaling or multifractal analysis for images has so far been mostly conducted using a multifractal analysis based directly on the wavelet coefficients of a 2D DWT. The corresponding structure functions are computed as:

S d (2 j , q) = 1 3n j 3 m=1 k1,k2 |d (m) X (j, k 1 , k 2 )| q . ( 14
)
The remainder of the formalism (Eqs. ( 5) to (10)) and of the estimation procedures (Eqs. [START_REF] Xu | A projective invariant for texture[END_REF] to ( 13)) can be obtained by replacing L X (j, k) with |d X (j, k)| (and similarly for the bootstrap procedures defined in Section 4 below).

As pointed out earlier in Section 1, the wavelet coefficient multifractal analysis suffers from a major theoretical drawback: A complete multifractal analysis cannot be conducted, as negative q structure functions are numerically unstable and diverge because wavelet coefficients are mostly concentrated around the 0 value, see [START_REF] Jaffard | Wavelet leaders in multifractal analysis[END_REF] where a precise mathematical analysis of this phenomenon is carried out for some classes of stochastic processes. The wavelet Leader multifractal analysis is designed to overcome this difficulty and benefits from two other key theoretical properties: It is proven to be valid for selfsimilar and multiplicative cascade processes (while the wavelet coefficient multifractal analysis holds only partially, and for a restricted class), and the scaling exponents ζ(q) are shown not to depend on the precise choice of the wavelet1 H 0 , G 0 [START_REF] Jaffard | Wavelet leaders in multifractal analysis[END_REF][START_REF] Jaffard | Wavelet techniques in multifractal analysis[END_REF]. In Section 6, we further compare (very favorably) the Leader based multifractal analysis to the wavelet coefficient based one both in terms of practical estimation performance and bootstrap based confidence interval derivation.

However, an important limitation of the wavelet Leader based multifractal analysis results from the fact that it can be directly applied only to bounded functions, while the wavelet coefficient multifractal analysis is not restricted to this class. This question, together with a number of related theoretical issues of importance for a practical use of multifractal analysis on real images, is further addressed in the next section.

Advanced considerations on multifractal analysis of images

As mentioned in Section 1, there is a priori no guarantee that actual digital images are the discretization of a function that belongs to an L q (R 2 ) space, q ≥ 1, (or, a fortiori, belongs to the space of bounded variation (BV) functions, i.e. functions whose gradient is a bounded measure). Along the same line, one can hardly decide a priori whether the images under analysis fall within the class of bounded functions, permitting a straightforward application of the Leader multifractal analysis. Indeed, scanning large image databases (such as the one reported in [START_REF] Xu | A projective invariant for texture[END_REF]) confirms that there are roughly as many images which are bounded as images which are unbounded (see in Section 3.2 how this can be can be practically determined on discretized signals). Hence, the Leader based multifractal analysis cannot be applied to images without prior checking. Note that the violation of the bounded function requirement can be related to the existence of negative Hölder exponents in the data (see Section 3.2 below), and it is indeed commonly reported in the literature dedicated to empirical multifractal analysis that images exhibit negative Hölder exponents [START_REF] Ossiander | Statistical estimation for multiplicative cascades[END_REF]. There is no general consensus on a precise mathematical definition of negative Hölder exponents, and a discussion of this notion requires theoretical developments beyond the scope of this contribution (this is addressed in [START_REF] Wendt | Analyse multifractale d'images : l'apport des coefficients dominants[END_REF][START_REF] Jaffard | Wavelet analysis of multifractal measure[END_REF]); we can however indicate that the underlying heuristic remains yet the same: in some averaged sense, |X(t) -

X(t 0 )| ∼ |t-t0|→0 |t -t 0 | h .
These involved issues have been partially addressed in [START_REF] Wendt | Analyse multifractale d'images : l'apport des coefficients dominants[END_REF] and are further elaborated in this section. Furthermore, the analyses reported here indicate that, though wavelet coefficients offer a restricted analysis of the multifractal properties of an image, a number of useful information can still be extracted from the wavelet coefficient based structure functions.

Images and function space models

Let us define, for q > 0,

ζ d (q) = lim inf 2 j →0 ln S d (2 j , q) ln 2 j . ( 15
)
Then, ζ d (q) > 0 indicates that the image X belongs to L q (R 2 ). Indeed, the wavelet characterization of the Sobolev spaces L q,s implies that if q ≥ 1, then ζ d (q) = sup{s : f ∈ L q,s/q }, see [START_REF] Jaffard | The contribution of wavelets in multifractal analysis[END_REF]; therefore ζ d (q) > 0 implies that f belongs to a Sobolev space L q,s for an s > 0; a fortiori, it belongs to L q,0 , which coincides with the space L q . For q = 2, this provides the practitioners with tools helping to decide whether or not X fits the assumptions of the Osher-Rudin-Fatemi model [START_REF] Rudin | Nonlinear total variation based noise removal algorithms[END_REF]. This remark is not only pertinent for this model: All other models which have been later proposed as alternatives (see for instance the book of Y. Meyer [START_REF] Meyer | Oscillating patterns in image processing and nonlinear evolution equations[END_REF] and references therein, or [START_REF] Abry | Wavelet decomposition of measures: Application to multifractal analysis of images[END_REF]) are based on the hypothesis that the image considered belongs to a given function space (or a sum of function spaces). Since all the function spaces which have been proposed in this context either are Sobolev spaces, or satisfy sharp embeddings with Sobolev spaces, it follows that the determination of the wavelet scaling function of a given image allows to determine if it satisfies a particular function space assumption, and, therefore, if the corresponding algorithm can be reasonably applied in this case. When q = 1, the same argument as above yields that, if ζ d (1) > 0, then X belongs to L 1 (R 2 ) and additionally, ζ d (1) > 1 implies that X belongs to the class of bounded variations images. Furthermore, X is a measure necessarily yields

ζ d (1) ≥ 0.
Scanning image databases reveals that a non negligible proportion of images (e.g. ≈ 10% for the database in [START_REF] Xu | A projective invariant for texture[END_REF]) are characterized with ζ d (2) < 0, with confidence intervals (computed with the technique proposed in Section 4) clearly validating the negativity for the estimate. Also, ζ d (1) is positive for most images (as expected for positive measures), yet that 0 < ζ d (1) < 1 for a large proportion of images (e.g. ≈ 90% for the database in [START_REF] Xu | A projective invariant for texture[END_REF]), which are hence not within the bounded variation class. For illustration purpose, the image used in Section 6.4 (cf. Fig. 7) is characterized by ζ d (1) = -0.08 ± 0.03 and ζ d (2) = -0.17 ± 0.08, confidence intervals being obtained with the bootstrap approach described in Section 4.

Uniform regularity

A measure or a function X(t) is said to belong to C ε , ε ∈ R, when its wavelet coefficients satisfy:

∃C > 0 : ∀j, k 1 , k 2 , m |d (m) X (j, k 1 , k 2 )| ≤ C2 jε . (16)
A uniform regularity exponent can hence be defined by

h min = sup{ε : X ∈ C ε }. (17)
A function or a measure is said to be uniformly Hölder if h min > 0. In turns, h min > 0 implies that

∀t 0 , h(t 0 ) ≥ h min ,
and also that X is a continuous function, hence bounded, hence possesses finite wavelet Leaders in the limit of fine scales (all implications are strict).

As stated in Section 2.2 above, the computation of wavelet Leaders, hence the computation of the Legendre spectrum is possible only for functions which are uniformly Hölder. Applied to images that do not respect this restriction, the wavelet Leader based multifractal analysis produces results that are meaningless. However, one can meet the following pitfall: In practice, any image will yield empirical wavelet Leaders with finite value (because the supremum in the definition of wavelet leaders (1) is in practice always taken on a finite set), independently of the fact that h min is positive or not. Hence, practitioners have little or no means to decide a posteriori whether Leaders are meaningful or not and this needs to be checked a priori. A sufficient condition is that the image is uniformly Hölder and hence that h min > 0. This can be evaluated via the formulation:

h min = lim inf 2 j →0 ln sup m,k1,k2 |d (m) X (j, k 1 , k 2 )| ln 2 j , (18)
which is a direct consequence of Eq. ( 16) and practically amounts to performing linear regressions of the log of the magnitudes of the largest wavelet coefficients at scales 2 j versus log of scales. The performance of this estimator, assessed over 500 realizations of a synthetic multifractal stochastic images of size 1024 × 1024 (cf. Section 5.3), are summarized in Tab. 1. They clearly demonstrate that estimation of the uniform regularity exponent is feasible from a standard size image.

Examination of image databases shows that images are found as often with negative h min as with positive one (e.g. ≈ 50% for the database in [START_REF] Xu | A projective invariant for texture[END_REF]). For a simple illustration purpose, the image used in Section 6.4 (cf. Fig. 7 and Fig. 8) is characterized by h min = -0.24. This calls for a modification of the Leader based multifractal analysis proposed in Section 2.2 for bounded functions. This is the subject of the next section.

h min CMC-LP theo mean std mse -0.391 -0.395 0.117 0.117 Table 1: Estimation of h min . Mean, standard deviation and root mean squared error (cf. Eq. ( 27)) of wavelet coefficient based estimation of h min (Eq. ( 18)) for CMC-LP, obtained by numerical simulation as described in Section 5.3 below with linear regressions over the finest scales 2 1 ≤ 2 j ≤ 2 3 . The column on the left states the theoretical value of h min .

Fractional integration

To overcome the negative h min issue, it has been proposed to fractionally integrate with an order larger than -h min , which implies that the uniform regularity exponent of the new image is positive, and thus insures that all Hölder exponents are positive. This has been for instance abundantly used in [START_REF] Decoster | A wavelet-based method for multifractal image analysis. II. Applications to synthetic multifractal rough surfaces[END_REF][START_REF] Arneodo | A wavelet-based method for multifractal image analysis: from theoretical concepts to experimental applications[END_REF][START_REF] Arneodo | A wavelet-based method for multifractal image analysis. I. Methodology and test applications on isotropic and anisotropic random rough surfaces[END_REF].

The fractional integration (of order η) of a function or measure X is defined in the Fourier domain as:

( I η X)(ξ) = (1 + |ξ| 2 ) η/2 X(ξ). ( 19
)
If X is such that h min ≤ 0, then I η X is a uniformly Hölder function as soon as η > -h min ;

this result is a direct consequence of the interpretation of h min in terms of the Lipschitz spaces C α (R n ):

h min = sup{α : X ∈ C α (R n )},
and of the fact that if X ∈ C α , then I η X ∈ C α+η . Indeed, it follows that if η > -h min , then I η X belongs to C ε for an ε > 0, hence is a bounded function, see [START_REF] Jaffard | Wavelet analysis of multifractal measure[END_REF] for details.

Instead of actually computing the fractionally integrated version of X and then applying the wavelet Leader multifractal formalism to it, both operations can be combined into a single one, as follows.

i) First compute the 2D WT coefficients d (m) X (j, k 1 , k 2 ) and replace them with:

d (m),η X (j, k 1 , k 2 ) = 2 ηj d (m) X (j, k 1 , k 2 ).
This amounts to computing the wavelet coefficients of Ĩη X, a pseudo-fractionally integrated version of X, whose local and global regularity properties are identical to that of I η X, as soon as η > -h min , see [START_REF] Jaffard | Wavelet analysis of multifractal measure[END_REF].

ii) Wavelet Leaders are computed from these new wavelet coefficients :

L η X (j, k 1 , k 2 ) = sup m, λ ′ ⊂3λ j,k 1 ,k 2 |d (m),η X (λ ′ )|. (21)
Such modified Leader coefficients are equivalent to those computed from I η X, in the sense that if I η X has Hölder exponent h at t 0 then, when 2 -j k is the closest dyadic point of t 0 ,

L η X (j, k 1 , k 2 ) ∼ 2 jh .
iii) Compute new structure functions:

S η (j, q) = 1 n j k1,k2 L η X (j, k 1 , k 2 ) q . ( 22
)
They behave as power laws with respect to the analyzing scale 2 j , in the limit of fine scales 2 j → 0 :

S η (j, q) ≈ G η q 2 jζη(q) . ( 23
)
yielding the multifractal spectrum of I η X:

L η (h) = min q =0 (1 + qh -ζ η (q)). (24)
One can prove that the scaling function (and therefore the Legendre spectrum) thus obtained is the same one as if a "true" fractional integration had been performed, see [START_REF] Jaffard | Wavelet analysis of multifractal measure[END_REF]. Note that Condition [START_REF] Arneodo | A wavelet-based method for multifractal image analysis: from theoretical concepts to experimental applications[END_REF] does not prescribe a unique value for η. Therefore one meets the following alternative:

-A thorough possibility consists in computing the Legendre spectra L η (h) for a wide range of values of η, and use this collection of spectrums for classification, or in order to understand the nature of the pointwise singularities of X. Preliminary results in this direction are given in [START_REF] Jaffard | Wavelet analysis of multifractal measure[END_REF].

-A simpler option consists in noticing that, for large classes of mathematical functions and stochastic processes, one can prove that all the functions L η (h+η) are actually the same, and therefore define a unique function L(h):

∀η > -h min , L(h) = L η (h + η). (25)
A precise analysis of the validity of Eq. ( 25) is beyond the scope of this paper and can be found in [START_REF] Jaffard | Wavelet analysis of multifractal measure[END_REF][START_REF] Abry | Wavelet decomposition of measures: Application to multifractal analysis of images[END_REF]; let us just mention that it is related to the absence of strongly oscillating features in the image, such as the ones provided by chirps of the form |x -x 0 | α sin(|x -x 0 | -β ). Note however that some classical models do not satisfy Eq. [START_REF] Wendt | Multifractality tests using bootstrapped wavelet leaders[END_REF]; typical examples are supplied by lacunary wavelet series, see [START_REF] Jaffard | Lacunary wavelet series[END_REF]. If Eq. ( 25) holds, then the particular value of η which is picked is irrelevant (as soon as it is large enough): In practice, one picks one value of η and defines the Legendre spectrum by Eq. ( 25); the multifractal analysis is performed using this function L(h) for spectrum. Though this second option is successfully used for classification in a number of applications, it is important to be aware that this heuristic faces limitations: Besides the problem of the validity of Eq. ( 25) which has already been mentioned, the function L(h) can not generally be related to the multifractal properties of X. Nonetheless, it is deeply related to the intrinsic properties of the original image X and can be used as such (for classification, for instance).

Finally, note that fractional integration also increases the ζ d (q) (the counterpart of the heuristic translation in Eq. ( 25) reads ζ d η (q) = ζ d (q) + ηq), hence ensuring, for large enough η, that ζ d η (1) > 1 and/or ζ d η (2) > 0.

Validity of the multifractal formalism.

It is also of importance to be aware that theoretically Eq. ( 7) does not hold in general. More precisely, let us assume that X belongs to the class of mild functions, defined as:

∃C 1 , C 2 , A, B > 0 such that ∀j, k, C 1 2 jA ≤ L X (j, k) ≤ C 2 2 jB (
this implies that X is nowhere C ∞ , but satisfies a minimal uniform regularity assumption). When X is a mild function, its ζ(q) is an admissible function, i.e., it satisfies:

i) ζ(0) = 0, ii) ζ(q) is concave and increasing on R, iii) 0 ≤ ζ(q) -qζ(q)/dq ≤ d.
Conversely, we can show [START_REF] Jaffard | The contribution of wavelets in multifractal analysis[END_REF] that any admissible function corresponds to the ζ(q) of a mild function X. However, this does not ensure that X satisfies the wavelet Leader multifractal formalism. Let, for instance, λ(q) be an admissible function. We can design the following trifractal function whose ζ(q) = λ(q): It is constructed such that its multifractal spectrum reads D(h min ) = D(h max ) = 0, D(h med ) = d and D(h) = -∞ elsewhere, with h med = ζ ′ (0), h min = ζ ′ (+∞), and h max = ζ ′ (-∞). As soon as these three values h min , h med , h max differ, the multifractal formalism (as in Eq. ( 7)) obviously ceases to hold, since the Legendre transform of ζ(q) is a concave function. Another counterexample is detailed in [START_REF] Jaffard | The contribution of wavelets in multifractal analysis[END_REF] where the generality of the validity of the multifractal formalism is further investigated.

Wavelet coefficients versus Leaders

This section leads to the conclusion that wavelet coefficients are providing preliminary information regarding the regularity properties of X and should hence be used prior to applying the wavelet Leader multifractal analysis, and in a complementary manner, rather than with the usual competition perspective.

Space-scale block bootstrap

We now extend the line of work proposed in [START_REF] Wendt | Bootstrap for empirical multifractal analysis[END_REF][START_REF] Wendt | Multifractality tests using bootstrapped wavelet leaders[END_REF] for the design of bootstrap confidence limits for 1D multifractal processes (signals) to 2D ones (images). While most scale invariant processes are characterized by difficult statistical properties (non stationary, intricate dependence), their wavelet coefficients often form stationary sequences with less involved dependence structure [START_REF] Wendt | Bootstrap for empirical multifractal analysis[END_REF]. Therefore, it is more natural to apply bootstrap schemes in the wavelet domain rather than in the space one. We devise a specific construction of space-scale blocks to ensure that the residual dependencies of Leaders, both in space and in scale, are approximately reproduced by the bootstrap resamples. Also, these space-scale blocks ensure that bootstrapped wavelet Leaders preserve their key property for multifractal analysis, namely Eq. [START_REF] Schertzer | Multifractal cascade dynamics and turbulent intermittency[END_REF].

For an overview of bootstrap and bootstrap methods for dependent data, see e.g. [START_REF] Lahiri | Resampling Methods for Dependent Data[END_REF][START_REF] Zoubir | Bootstrap Techniques for Signal Processing[END_REF][START_REF] Zoubir | The bootstrap and its application in signal processing[END_REF].

Space-scale blocks of wavelet Leaders.

For ease of notation and without loss of generality, we assume square images (N 1 = N 2 = N ). The space-scale blocks of Leaders are defined as 3D boxes, extending over all scales, with a square base of a fixed number of pixels. All Leaders lying within such a box form one space-scale block. The blocks are constructed overlapping and on circularized Leaders, i.e. on L c X (j, k 1 , k 2 ) = L X (j, k 1 mod n j , k 2 mod n j ), such that each Leader has the same probability to be within a resample. More precisely, the collection of Leaders L X (j ′ , k ′ 1 , k ′ 2 ) that form a space-scale block I(k 1 , k 2 ), of 2l×2l pixels, located at position (k 1 , k 2 ), is given by (cf. Fig. 2):

I(k 1 , k 2 ) = {L X (j ′ , k ′ 1 mod n j ′ , k ′ 2 mod n j ′ ) : |k 1 -k ′ 1 2 j ′ | ≤ l, |k 2 -k ′ 2 2 j ′ | ≤ l, 1 ≤ j ′ ≤ j max }. ( 26
)
4.2 Bootstrap resampling and estimation.

Each single bootstrap resample of Leaders L * (r)

X , r = 1, • • • , R, is obtained as follows.
First, B = ⌈ N 2l ⌉ blocks of Leaders I * (k 1 , k 2 ) are drawn at random, independently and with replacement, from the available blocks

I(k 1 , k 2 ), k 1 , k 2 = 1, • • • , N .
Then, blocks are concatenated in space, such that each resampled Leader remains located at its original scale 2 j . The bootstrap estimation procedure re-applies the estimation procedure described above to each of this bootstrap resamples L * (r) X : First, the bootstrap structure functions S(2 j , q) * and Ĉ(2 j , q) * are obtained by applying Eqs. ( 3) and [START_REF] Decoster | A wavelet-based method for multifractal image analysis. II. Applications to synthetic multifractal rough surfaces[END_REF]. Then, bootstrap estimates ζ(q) * , D(q) * , ĥ(q) * and ĉ * p are calculated as in Eqs. ( 11), ( 13) and [START_REF] Bacry | Singularity spectrum of fractal signals from wavelet analysis: Exact results[END_REF]. These estimation and bootstrap estimation procedures are sketched in Fig. 3. Note that the bootstrap resampling and estimation procedures can as well be applied directly to wavelet coefficients instead of Leaders.

Bootstrap confidence limits

For parameters θ ∈ {ζ(q), D(q), h(q), c p }, we use the equi-tailed (1 -α) bootstrap percentile confidence interval,

real world {L X } (3) (8) s s s y y s s s s {L * (r) X } (8) (3) O O O ' ' O O O O O bootstrap world Ĉ(2 j , p) (12)
S(2 j , q) Overview of estimation procedure (left) and bootstrap estimation procedure (right). "resample" corresponds to the spacescale block bootstrap resampling procedure described in Section 4, "CI" to confidence limit calculation from empirical bootstrap distributions as in Section 4.3, and (•) on arrows indicates the equation involved in estimation. The scheme can be written equivalently for D(q) and ĥ(q).

defined as CI θ = θ * α 2 , θ * 1-α 2
, where θ * α is the α quantile of the empirical distribution of θ * obtained by the bootstrap estimation procedure described above. Double bootstrap estimates could be obtained and used for more sophisticated confidence limit and hypothesis test constructions, such as studentized or adjusted confidence limits and tests. Such issues have been considered in [START_REF] Wendt | Bootstrap for empirical multifractal analysis[END_REF][START_REF] Wendt | Multifractality tests using bootstrapped wavelet leaders[END_REF] for the 1D case. An extension to images is currently investigated.

5 Performance assessment on synthetic multifractal images

Methodology: Monte Carlo simulations

We assess the performance of the proposed estimation procedures by applying them to a large number N M C of realizations of synthetic stochastic 2D processes of size (N × N ) with a priori known and controlled multifractal properties. The aim of the numerical study is to address the following issues: Do the estimation procedures exhibit satisfactory statistical performance? Should one prefer wavelet coefficients or Leaders for the estimation of multifractal attributes of images? Are the bootstrap confidence limits reliable, i.e., do they reproduce targeted coverages?

Synthetic multifractal processes

We FBM. We use the 2D Gaussian self-similar Fractional Brownian Motion, as defined in [START_REF] Stein | Fast and exact simulation of fractional brownian surfaces[END_REF]. Its statistical properties are completely determined by a single parameter, the Hurst exponent H. FBM has scaling properties as in Eq. ( 5), with ζ(q) = qH and therefore c 1 = H, c p ≡ 0, p ≥ 2. Its multifractal spectrum reduces to one single point h = H where: D(H) = 2, hence, FBM is monofractal.

CMC-LP. The multiplicative cascades of Mandelbrot (cf. [START_REF] Mandelbrot | Intermittent turbulence in self-similar cascades; divergence of high moments and dimension of the carrier[END_REF]) have for long been the only example of multifractal processes practically available. Their construction is based on an iterative split-and-multiply procedure on an interval. In this work, we use a binary cascade and log-Poisson multipliers s = 2 γ exp (ln(β)π λ ), where π λ is a Poisson random variable with parameter λ = -γ ln(2) (β-1) . In order to ensure that the process has a minimum regularity, the cascade is integrated fractionally with order η > 0. The CMC-LP are multifractal with ζ(q) = (η -γ)q + γ(β q -1)

β-1 , c 1 = η + γ ln(β) β-1 -1 and c = -γ β-1 (-ln(β)) p , D(h) = 2 + γ β-1 + -η+γ+h ln β • ln (-η+γ+h)(β-1)
γ ln β -1 .

Simulation setup

Results are reported here with one specific selection of process parameter settings. Similar results have been obtained for a range of different process parameters, and are not reported here for space reasons. Parameters for numerical simulations are set to N M C = 500 and N = 1024. All results are obtained with Daubechies' wavelets with N ψ = 2 vanishing moments. It has been checked that using wavelets with larger N ψ yields identical results and conclusions. Linear regressions are performed over the scales 2 3 ≤ 2 j ≤ 2 7 with weights b j = n j , as proposed in [START_REF] Wendt | Bootstrap for empirical multifractal analysis[END_REF]. We use R = 99 bootstrap resamples and block lengths l = 85 (FBM) and l = 128 (CMC-LP).

The target significance for bootstrap confidence limits is set to 90%. The process parameters are fixed to H = 0.7 for FBM, and β = 0.84, η = 0.42, η = 0.5 for CMC-LP, such that [c 1 , c 2 , c 3 ] = [0.538, -0.080, 0.014]. CMC-LP

1 3 5 7 -1 0 1 q=-2 log 2 S(2 j
,q) -j⋅ζ(q) j 1 3

,q) -j⋅ζ(q)

1 3 5 7 -1 0 1 q=-2 log 2 S(2 j
,q) -j⋅ζ(q) j 1 3

,q) -j⋅ζ(q)

Figure 5: Structure Functions. Structure functions log 2 S(2 j , q) for q = 2 (top row) and q = -2 (second row), and deviation of structure functions from theoretical slope, log 2 S(2 j , q) -j • ζ(q), for q = 2 (third row) and q = -2 (bottom row), obtained by mean over Monte Carlo realizations for FBM (left) and CMC-LP (right). Red circles correspond to wavelet Leader, blue crosses to wavelet coefficient based estimations. The error bars correspond to 1.96 Monte Carlo standard deviation.

Results

Structure functions

Fig. 5 shows, for FBM (left) and CMC-LP (right), means over Monte Carlo realizations of structure functions log 2 S(2 j , q) for q = 2 (top row) and for q = -2 (second row), and the deviations of the structure functions from their theoretical slope, log 2 S(2 j , q) -jζ(q), for q = 2 (third row) and q = -2 (bottom row). The 95% asymptotic confidence limits are obtained by Monte Carlo simulation. As structure functions for wavelet coefficients diverge for negative qs, only positive qs are shown. Scaling range. A first investigation (top and second row) suggests that for both wavelet coefficients and Leaders, structure functions display scaling behavior as in Eq. ( 5) over the range of scales 2 2 ≤ 2 j ≤ 2 7 . A closer look at the deviations from the theoretical slope (third and last row) confirms this observation for coefficients. However, it reveals that for Leaders, log 2 S(2 j , q) becomes a linear function of j only for 2 j ≥ 2 3 . This can be explained by the fact that, theoretically, a Leader is defined as the sup of coefficients at all finer scales down to infinitely fine scales, whereas practically, the sup can be taken only down to the finest available, first scale. Hence, in practice, Leaders need one or two scales for initialization, whereas coefficients do not. Projection step. The non-scaling behavior of S(2 j , q) of coefficients at the first scale j = 1 is due to the fact that the prefiltering or projection step theoretically necessary for a clean wavelet analysis has been omitted (cf. [START_REF] Veitch | Meaningful MRA initialization for discrete time series[END_REF]). Regressions. We observe further that the error bars for structure functions for CMC-LP are substantially larger than those for FBM, suggesting a smaller variability of the S(2 j , q) and better subsequent estimation performance for the monofractal process. Moreover, it is interesting to note that whereas the error bars for FBM behave approximately as 1/ √ n j , this is not the case for CMC-LP, where the size of the confidence intervals varies only slightly with j. This confirms that the choice of weights for weighted linear fits, as proposed in [START_REF] Wendt | Bootstrap for empirical multifractal analysis[END_REF], is appropriate for FBM. Though not optimal for multiplicative cascades, it has been shown to perform better than non-weighted regression.

Performance of parameter estimation

We quantify the performance of the estimators θ ∈ { ζ(q), ĉp } by their root mean squared error,

mse θ N M C = E N M C θ -θ 2 + Var N M C θ, (27)
where E N M C and Var N M C stand for the sample mean and variance, respectively, over N M C independent realizations. Estimation performance results are summarized in Tab. 2 and illustrated in Fig. 6, both for wavelet coefficients and Leaders, for medium size images (N = 1024). In addition, Tab. 2 shows complementary results for large images of CMC-LP (N = 2048). On overall, we observe that Leader based estimations are more efficient than coefficient based ones.

Positive statistical moments q. For q > 0, both estimation procedures have approximately equal mse and std. Coefficient based estimations are slightly more efficient for monofractal FBM while Leader based ones are for multifractal CMC-LP.

Negative statistical moments q. For q ≤ -1, the wavelet coefficient based estimates of ζ(q) are not meaningful as they exhibit very large mse. Therefore, a wavelet coefficient based multifractal analysis of images allows to explore exclusively the range q > -1 and thus, in practice, only the increasing part of the multifractal spectrum D(h). In contrast, the proposed wavelet Leader based procedure permits a complete analysis of the multifractal properties of an image, and notably of the decreasing part of the multifractal spectrum (cf. Fig. 6). Log-cumulants. Tab. 2 shows that the wavelet Leader based estimations of c p exhibit consistently smaller std and mse than their wavelet coefficient based counterparts. Whereas the difference in performance is only small for the estimation of c 1 , it becomes more significant for c 2 and c 3 , with gains in mse of up to more than one order of magnitude. This is of crucial importance, since non-zero c 2 and c 3 discriminate self-similar from multiplicative cascade processes. Further, for N = 2048, Tab. 2 (last line) enables to deduce that an asymptotic 85% confidence interval for the Leader based estimation of c 3 excludes zero, hence that the real c 3 is different from zero with high probability. To the best of our knowledge, this had never been achieved on multifractal images. It also clearly shows that for images of smaller size, the estimation of parameter c 3 , a fortiori of higher order c p , should be used with care: CIs might be so large that they may not exclude 0 even if the parameters are non zero.

Self-similar vs. multiplicative cascade processes. Finally, for both the Leader and coefficient based procedures, estimation is more difficult for multifractal CMC-LP than for FBM, resulting in larger mse for estimations on the multifractal process.

Bootstrap estimation performance

We evaluate the reliability of the bootstrap percentile confidence limits CI θ by their recentered empirical coverage:

C θ = E N M C I θ + E N M C θ -θ ∈ CI θ , (28)
where I {•} is the indicator function of the event {•}. The recentered coverage C θ allows to assess the performance of the confidence limits without the contribution of a potential bias in the estimation of θ. The results for θ ∈ {ζ(q), c p } are summarized in Tab. actual coverage are very close to the target coverage, with error in coverage smaller than 5% both for ζ(q) and c p , and for both processes. The performance of the wavelet Leader confidence intervals is slightly inferior, with an average coverage error of approximately 9%. This can be interpreted as the result of the non linear operation underlying the construction of Leaders from coefficients, yielding an additional difficulty for bootstrap estimation. Still, the bootstrap confidence limits remain satisfactorily reliable. The use of more sophisticated and potentially more accurate double bootstrap confidence limits is currently under consideration.

Analysis of real-world image

The proposed estimation procedures have been used for the systematic analysis of large databases of real-world images (such as the texture image database in [START_REF] Xu | A projective invariant for texture[END_REF]). For illustrating its practical use, we present here the performance of the pro-

X (j, k 1 , k 2 )| with respect to scales j indicates the relevance of estimation of h min , which is found to be clearly negative. Therefore, the image has been fractionally integrated with an order η = 0.5 before analysis to ensure a minimum regularity (this order is assumed to be sufficient, since the estimated uniform regularity exponent for this image is h min = -0.24). Multifractal attribute estimates are obtained by weighted fits over the scales 2 3 ≤ 2 j ≤ 2 6 , for which Fig. 7 (center) indicates that log 2 S(2 j , q) are approximately linear functions of j. Results demonstrate that the proposed method can be readily applied to real-world images for their complete multifractal characterization. They demonstrate further that the bootstrap percentile 90% CIs are of significant practical usefulness: First, the confidence limits on structure functions help to verify that the chosen scaling range is appropriate. Second, confidence limits on multifractal attributes allow to decide that the image shows with high probability a multifractal signature, since confidence intervals for c 2 and c 3 exclude zero. The confidence limits thus provide information that might be of significant importance for e.g., detection or classification tasks, or for the understanding of physical or biological processes underlying the data. Finally, note that the entire estimation procedure with WMLF (without bootstrap) takes less then a minute on a standard PC, whereas the same estimation with MMWT increases computation time by a factor larger than 20. Table 4: Real-World Image Log-Cumulants. Log-cumulant estimates and bootstrap percentile 90% confidence limits (R = 99) for the real-world image in Fig. 7 (top).

Discussion

Numerical simulations, equivalent to those reported above, have been performed on other multifractal processes and lead to similar conclusions: Wavelet coefficients do not allow to meaningfully explore negative qs and thus to measure the decreasing part of the spectrum D(h) whereas wavelet Leaders do, estimation of c p is better when based on Leaders than on coefficients, and significantly so for c p , p ≥ 2.

Selection of scaling range. We have seen in Section 6.1 that the range of scales over which the linear fits are to be performed are likely to be narrower (requiring the use of a larger j 1 ) for Leaders than for wavelet coefficients. For practical multifractal analysis and real-world images, the choice of this regression range is a crucial, difficult and controversial issue, further complicated by the use of Leaders. In practice however, bootstrap confidence intervals for the structure functions, as obtained by the proposed procedure and illustrated in Fig. 5, constitute a precious support for solving this non-trivial issue and can be regarded as a first step towards statistical procedures for automatized or intelligently assisted detection of scales over which data are scale invariant.

Vanishing moments of wavelet. Another critical practical issue is the choice of the number of vanishing moments N ψ for the wavelet with which the data are analyzed. Condition N ψ > h, where h is the largest singularity exponent present in the data, is expected to be sufficient for a relevant multifractal analysis. However, in practice, the choice of N ψ results from a trade-off: A larger N ψ stabilizes the estimates of the negative q structure functions and enables to get rid of potentially superimposed smooth trends such as polynomial, hence improves estimation and brings robustness; A larger N ψ also implies a larger support for the wavelet and thus produces border effects of wider size, such that no coefficients may remain unaffected at large scales, hence in itself degrades estimation performance. Therefore, a reasonable practical rule of thumb is to choose the smallest N ψ for which the estimated multifractal attributes do not significantly change when N ψ is increased (i.e, remain within confidence intervals, emphasizing again the need for and the importance of such confidence intervals).

Real-world data and multifractal analysis. In Section 2.2, we recalled how the estimates of ζ(q), c p , D(q) and h(q) are theoretically intimately tied to the multifractal spectrum D(h) of the field X and hence to the analysis of its local regularity fluctuations. This relation is known to hold theoretically for all the synthetic images used here for Monte Carlo simulations and illustrations. For real-world images, this interpretation in terms of multifractal spectrum and singularity description might not always be completely relevant. However, this, in no way, prevents practitioners to make use of the measured ζ(q), c p , D(q) and h(q) to analyze the data in terms of a less mathematically stringent formulation of scale invariance property, or to perform standard image processing tasks such as classification or retrieval based on such quantities.

Computational costs. The proposed 2D wavelet Leader multifractal analysis is simple both conceptually (2D DWT, Leaders and linear regressions) and practically (very low computational cost, both with respect to time and memory, of the order of a 2D DWT). As an example, a 2048 × 2048 gray level image is processed in a couple of seconds on a standard PC under MATLAB. The bootstrap based CIs are obtained essentially by repeating the linear regression procedures. The larger the repetition number and therefore the computational load, the higher the precision of the CIs. This is the unavoidable price to pay to obtain reliable CIs from a single image. All the analysis and estimation procedures are designed by ourselves in MATLAB. The 2D DWT is performed using (a corrected version of) the Rice Wavelet Toolbox.

Conclusions and perspectives

The present contribution proposes a procedure to perform a fast, efficient and accurate analysis of the scaling and multifractal properties of images. It outperforms significantly the previous propositions based on 2D DWT or 2D CWT in estimation performance, and 2D-WTMM in computation time, memory cost and implementation complexity. It is fast and efficient since simply based on a 2D orthogonal DWT. Our procedure enables an accurate and complete characterization of the (ir)regularities of the texture of an image thanks to the use of original multiresolutions quantities called wavelet Leaders, and is backed up by a strong mathematical framework. Furthermore, the present contribution proposes and validates procedures for assessing function space model conditions for data -notably for the bounded function requirement -and handling them in practice.

The low memory and time costs together with the satisfactory estimation performance of the procedure opens for the first time the possibility to perform the multifractal analysis of voluminous databases of images with possibly large sizes. Therefore, multifractal analysis may be incorporated in procedures aiming at image retrieval, computer vision or robotic purposes. We are currently investigating such databases. Moreover, for applications such as surface roughness, where one wants to understand the fracture process, or biomedicine, where one wants to detect and classify pathologies, confidence intervals are crucial. The bootstrap approach, tailored here to match image processing via the proposition of a 2D space-scale block bootstrap procedure, provides reliable and straightforwardly usable confidence intervals. This is, to the best of our knowledge, the only procedure that has been proposed to provide accurate and operational confidence intervals for multifractal analysis of images. It can be further extended to devise statistical tests. This is currently investigated.

Figure 1 :

 1 Figure 1: Definition of Wavelet Leaders. The wavelet Leader L X (j, k 1 , k 2 ) at scale 2 j and position (k 1 , k 2 ) (black cross) is defined as the largest of the wavelet coefficients |d

Figure 2 :

 2 Figure 2: Space-scale block construction. The space-scale block of Leaders I(k 1 , k 2 ) at position (k 1 , k 2 ) consists of the collection of Leaders (fat black dots) that are within a box with square base of base length 2l, centered at (k 1 , k 2 ) and extending over all scales (red volume).

3 Figure 3 :

 33 Figure 3: Estimation and bootstrap estimation procedures.Overview of estimation procedure (left) and bootstrap estimation procedure (right). "resample" corresponds to the spacescale block bootstrap resampling procedure described in Section 4, "CI" to confidence limit calculation from empirical bootstrap distributions as in Section 4.3, and (•) on arrows indicates the equation involved in estimation. The scheme can be written equivalently for D(q) and ĥ(q).

  use two stochastic processes: Fractional Brownian Motion (FBM) and Canonical Mandelbrot's Multiplicative Cascade with log-Poisson multipliers (CMC-LP). They provide us with simple yet representative examples of Gaussian monofractal processes and multifractal processes, respectively. Example of such fields are shown in Fig. 4.

Figure 4 :

 4 Figure 4: Synthetic processes. One realization of FBM (left) and of CMC-LP (right). The bottom plane presents the images, the axes on top an amplitude view of the images (amplitudes normalized to [0, 1]), illustrating the (ir)regularity of the images.

Figure 6 :

 6 Figure 6: Estimates ζ(q) and D(h). Estimates ζ(q) (top) and D(h) (bottom) for FBM (left) and CMC-LP (right), obtained from Leaders (red circles, dashed lines) and wavelet coefficient (blue crosses, dashed lines) by mean over Monte Carlo realizations, and the corresponding theoretical values (black solid lines with dots). The error bars correspond to 1 Monte Carlo standard deviation. The error bars for coefficient based estimates for D(h) are not shown for better visibility of the spectra. The insert for D(h) of FBM corresponds to the area shaded around the theoretical spectra location.

Figure 7 :

 7 Figure 7: Analysis of Real-World Image. Real world image (top) and estimates: Structure functions (center) for q = -2 (left) and q = 2 (right), scaling exponents ζ(q) (bottom left) and spectrum D(h) (bottom right). The bootstrap percentile 90% confidence limits are obtained from R = 99 resamples.

Figure 8 :

 8 Figure 8: Real-World Image h min . Estimate of h min Eq. (18) for the real-world image in Fig. 7 (top).

Table 2 :

 2 Estimation performance.

	3 for coefficient (top) and Leader (bottom)
	based estimation and for a target coverage of 90%. We ob-
	serve that for wavelet coefficient based estimation, the perfor-
	mance of the bootstrap confidence limits is excellent. Their

Table 3 :

 3 Bootstrap estimation performance. Recentered empirical coverage (in %) of bootstrap percentile confidence limits for wavelet coefficient (top) and Leader (bottom) based estimation. The empirical coverage equals the percentage of realizations for which the theoretical value of an estimate lies within the estimated confidence limits (target coverage 90%).

	Recentered coverage of confidence limits -Target: 90%
	DWT	ζ(-2) ζ(-1) ζ(1) ζ(2)	1 c 2	c 3
	FBM	91.8 88.8 94.4 87.8 89.4
	CMC-LP	89.0 89.2 88.2 89.8 88.2
	LWT	ζ(-2) ζ(-1) ζ(1) ζ(2)	c 1 c 2	c 3
	FBM	79.4 79.8 79.0 78.6 78.8 82.0 85.0
	CMC-LP	80.8 82.0 82.6 81.2 81.0 81.2 83.0

posed estimation procedures for one real-world image of size 1024 × 768, taken by ourselves with a standard digital camera. It consists of trees in the Reichraminger Hintergebirge, Upper Austria (cf. Fig.

7

, top). Fig.

7

shows structure functions (center), scaling exponents (bottom left) and the spectrum (bottom right) estimated from this image, obtained with the Leaders based estimation procedure proposed in this work, with a N ψ = 3 Daubechies' wavelet. Results on log-cumulants are summarized in Tab. 4. Fig.

8

illustrates the uniform regularity exponent h min Eq. (

18

) estimate. The linear behavior of log 2 sup m,k1,k2 |d

A precise mathematical statement of this result requires further elaboration beyond the scope of this contribution, see e.g.[START_REF] Jaffard | Wavelet leaders in multifractal analysis[END_REF][START_REF] Jaffard | Wavelet techniques in multifractal analysis[END_REF].


Claire Meymandi-Nejad

Esteban Perrotin
email: eperrotin@laas.fr

Ariane Herbulot
email: herbulot@laas.fr

Michel Devy
email: devy@laas.fr

Aircraft navigation on taxiways: evaluation of line detection algorithms proposed for automotive applications

Keywords: line detection, feature extraction, particle filter, Hough transform, aircraft navigation

While working on aircraft navigation on taxiways, the line detection is one of the main challenging problems to be solved. This subject has been widely studied in the literature in the automotive field. In this paper, we propose a comparison of three line detection algorithms based on methods validated in the automotive field but transposed in aeronautics where this subject has not been widely addressed. Some problematics appear: the tarmac environment differs from the usual road model and the camera's position impacts the visibility on the image. The first method presented here uses a particle filter while the second one is based on the Hough transform. In the second method, we perform a color-based detection and introduce a method to compute the reference color, using technical specifications for airport markings. The last method is the LaneNet neural network. Criteria such as the precision or the max range of the detection are computed and exploited to discuss the algorithms relevance. The comparison is performed on both simulated images (from a product of the OKTAL-SE company) and real ones (from Airbus Operations S.A.S.).

INTRODUCTION

Lane detection plays a significant role in driver assistance systems. We can note several use cases: axis keeping, vehicle position estimation, lane departure warning, road modelling, and automatic lane following. There are two main types of geometric methods commonly used for lane detection: feature-based and model-based methods. Another possibility is to use neural network based methods.

In feature-based algorithms, the most used feature is the contrast variation between marking/road or road/grass. This assumes that the contrast between the object and its environment is obvious. In the case of the contrast between the line marking and the road, the ground lines are designed to have more reflectivity than the road. That is why most of the studies ( [START_REF] Sun | Hsi color model based lane-marking detection[END_REF], [START_REF] Lipski | A fast and robust approach to lane marking detection and lane tracking[END_REF], [START_REF] Hota | A simple and efficient lane detection using clustering and weighted regression[END_REF], [START_REF] Kazemi | L*a*b* color model based road lane detection in autonomous vehicles[END_REF]) try to detect line markings for lane detection applications, either by using color models, contrast or orientation of surfaces. The low computational complexity of those algorithms is an advantage but they are sensitive to shadows or occlusions. Model-based methods are more robust but require more assumptions on the road model (number of lines in the image, expected width and curvature, ...). Those algorithms are more consuming in terms of computational resources.

Neural network based methods have been developed for lane detection in the automotive field, based on machine learning and convolutive neural networks. In those studies, we can separate the networks in two specific fields: lane ( [START_REF] Neven | Towards end-to-end lane detection: an instance segmentation approach[END_REF]) and line ( [START_REF] Li | Line-cnn: End-to-end traffic line detection with line proposal unit[END_REF]) detections, where the first's objective is to define the usable area for the vehicle whereas the second's is specific to the detection of road markings. The neural network must learn how to extract important information from images in order to recognize lines. This method requires a larger dataset.

[10] reviews a large number of geometric line detection methods in the automotive fields and offers a good comparison in terms of advantages and drawbacks. Fewer studies concern the scene interpretation in aeronautics such as [START_REF] Theuma | An image processing algorithm for ground navigation of aircraft[END_REF] or [START_REF] Beda | Methods and systems for providing taxiway stop bar information to an aircrew[END_REF]. The distinction between the automotive and aeronautic fields is important because they differ in terms of line models and varying environments. Those specifications complicate the transfer of algorithms from one domain to the other.

One example is the region of interest (ROI) where lines must be detected in images. In the aeronautic field, images are provided by different cameras where the ROI is not always easy to define. Figure 1 presents the two cameras used for this study. The image from the fin camera is largely occupied by the sky and the plane. The area of detection is limited and far from the sensor (for a Airbus A380, which is our model for the simulation, lines to be detected are at least 70 meters from the camera). Hence, the number of pixels containing useful information is greatly reduced.

Another important difference is that, in the automotive field, the line detection is most of the time based on the assumption that the marking to be detected is white. This assumption enables simplifications such as using a grayscale or luminance image: it allows to focus on the edge detection. In the taxiway areas, the line to be detected is yellow so the line appearance is not as easy to exploit. As white lines also exist on the tarmac and should not be tracked, the chrominance information is important. The method also needs to be robust to the noise (as the plane fuselage is visible in some images and could be painted in yellow). This article aims to discuss line detection in the airport areas context. The discussion will focus on three main points: a comparison between three methods based respectively on the particle filter (PF), the Hough transform (HT) and the LaneNet neural network (LNN); the usefulness of the inverse perspective mapping (IPM) preprocessing and the performance of the three methods under degraded weather conditions. The algorithms will be compared on either simulated images or both simulated and real images.

This paper is organized as follows: in Section 2, we present first the IPM preprocessing, applied or not to every image before detecting lines, and then, the three selected methods. Section 3 details the criteria used for the comparison. The results are summed up and discussed in Section 4 before the conclusion in Section 5.

METHODS

Before performing one of the three methods to be compared, several preprocessing functions are automatically applied to our images such as distortion correction, white balance, calibration or homographic transformation. However, we decided to separate two cases of study; on one hand we use the resulting image as an input for the line detection algorithms, on the other hand, we apply an IPM transformation as another preprocessing.

The two geometric methods use a color selection function for the line detection. When working on line detection, the reference color can easily vary from one location to another, or from the scene lighting. Some authors ( [START_REF] Mammeri | A real-time lane marking localization, tracking and communication system[END_REF], [START_REF] Kazemi | L*a*b* color model based road lane detection in autonomous vehicles[END_REF]) propose to work in other color spaces, such as HSV, Lab or YCbCr, that have the particularity to differentiate lightness information from chromatic information. Other authors ( [START_REF] Saudi | Fast lane detection with Randomized Hough Transform[END_REF] and [START_REF] Kheyrollahi | Automatic real-time road marking recognition using a feature driven approach[END_REF]) try to use a local adaptative threshold in order to be robust to shadowing or lightness variation. However, all of these methods still require to define the desired color which is chosen arbitrary.

We decided to implement these functions as parts of the methods instead of considering them as a preprocessing, because methods do not share the same color assumption. The particle filter method is based on [START_REF] Meymandi-Nejad | Lane detection and scene interpretation by particle filter in airport areas[END_REF] which includes the color assumption as part of the whole process, while the color assumption used for our implementation of the Hough Transform method is specific to our application and based on technical specifications defined in the xyz color space. The Hough transform takes as an input a binary image obtained from features extraction. Since we want to use the color as a feature, we have to define a threshold (i.e. a subset from the color space) to obtain the binary image. In contrast, the particle filter requires a weight map representing the probability of pixels to be on a line. The color extraction is done by computing the distance between the pixels' color and a reference color, different for the two methods.

The output of the three algorithms is a list of clusters. A cluster is a part of the image considered as a line (or a marking). It can be composed of all the consecutive pixels selected as part of the line (when using the particle filter or the LaneNet algorithm) or the two line extremities (when using the Hough transform).

Preprocessing: The IPM transform

We do not detail the IPM transform but it can be found in [START_REF] Meymandi-Nejad | Lane detection and scene interpretation by particle filter in airport areas[END_REF]. The result of the IPM transform can be observed in Figure 2. We decided to study the result of the three line detection algorithms with and without the IPM transform preprocessing because one of our assumption is that the IPM transform corrects the perspective effects on the lines which could simplify the detection. Moreover, the IPM transform creates additional synthetic information. It can be considered as a noise addition and degrading for the image. However, distant parts of a line are represented by few pixels in the raw image. The IPM transform will increase the number of pixels representing the line. Our second assumption is that this phenomenon, while adding noise to the image, should ameliorate the line detection on distant areas.

First Method: Particle Filter

The first method we proposed for this comparison is an algorithm using the particle filter principle and presented in [START_REF] Meymandi-Nejad | Lane detection and scene interpretation by particle filter in airport areas[END_REF]. The philosophy of this method is to consider a pixel as a particle and a line of the image as an instant t. The particle filter is launched from the bottom line to the top line of the image. For each instant t, a number of particles are generated and a weight is assigned to them. The weight given by the weight map computed beforehand represents an observation whereas the pixel is a possible explanation of the hidden state of the system. The observation helps determine which pixel is prone to be a description of a line in the image. The weight map contains information from color and edges. A prediction of the particle positions on the next line (at time t + 1) is made following a transition probability density to determine. In this implementation, the state is assumed to spread up vertically in the image but other cases can be thought upon such as in [START_REF] Tran | A versatile object tracking algorithm combining Particle Filter and Generalised Hough Transform[END_REF].

Second Method: Hough Transform

The Hough transform is based on the extraction of the color and edges information. The color extraction is dependent on the definition of the reference color which is more difficult to express than the edges. In our approach, we proposed to use a more specific color definition.

In aeronautics, airport signs and markings are well defined by the ICAO recommendation [START_REF]AERODROMES: aerodromes design and operations[END_REF]. From these recommendations, the color of airport lines are defined by a subset of the xyz color space. With that definition, it is possible to select a threshold to extract most of the lines from the image. However, images used for this comparison are defined in the sRGB space (a sub-space of the xyz color space). The yellow definition provided by ICAO is defined for the full xyz space. For our implementation we applied an empirical modification of this definition.

The first step of this method is to apply a chromatic thresholding using this yellow definition. A binary image is produced and the edge extraction method applied is the morphological gradient which is more adapted in this case than the Sobel or Canny filters usually found in the literature ( [START_REF] Xu | Study on a robust approach of lane departure warning algorithm[END_REF], [START_REF] Srivastava | Improved lane detection using hybrid median filter and modified hough transform[END_REF]). The line detection is performed by a Hough transform of the resulting image and the clustering of the results are provided by a hierarchical clustering algorithm such as presented in [START_REF] Hota | A simple and efficient lane detection using clustering and weighted regression[END_REF].

Third Method: Neural network approach

CNNs (Convolutional Neural Network) provide the best performances for computer vision applications including lane detection. They obtain the best results on most of the automotive benchmarks. However, the nature of airport areas is quite different from automotive ones and could degrade CNNs performances.

We chose to use the LaneNet CNN ( [START_REF] Neven | Towards end-to-end lane detection: an instance segmentation approach[END_REF]). It is a CNN architecture composed of three networks, one encoder and two decoders. The first decoder is used to provide a binary image containing pixels with potential line information. The second decoder is used to construct a pixel embeddings map that will be used to cluster pixels from the binary image in different lanes. Since our comparison is based on pixel comparison, (see Section 3), we only require the binary image.

Because of the low number of images in our dataset, we could not train this network again but test the already pretrained version of this network from TuSimple. This comparison is to argue about the possibility of using algorithms from the automotive field directly for aeronautics.

CRITERIA

In order to compare those three algorithms, we decided to base our analysis on multiple criteria that can be divided in categories, depending on the added knowledge value:

1. Based on detection range

• Maximum range of detection

One of the objectives in aeronautics is to have a wide detection field, because the aircrafts speed implies that decisions have to be taken as soon as possible. This value is given in pixels, computed on the overall image and represents the highest line in the image where pixels are detected.

Based on detection precision

• Mean of the clusters' pixels' distances to the ground truth This criterion is the main one to quantify the algorithm precision. It computes, for each cluster, the distance between the proposed pixels to represent a marking and the ground truth.

• Recall

This criterion is used to analyze the pertinence of an algorithm, dividing the number of true positives by the total number of detections.

The aim of this article is to compare methods and analyze their relevance on real images. We are not concerned by the computation time because we are not working on the final implementation of the algorithm (in GPU or FPGA).

RESULTS AND DISCUSSION

Multiple datasets are available in the automotive field for the benchmark of algorithms (e.g. KITTI, Oxford Robotcar or TuSimple data sets). In aeronautics, images acquired by cameras embedded on aircrafts are more difficult to recover. Thanks to a joint project with AIRBUS, we can use the OKTAL-SE simulator to build synthetic images on airport areas. The three methods have been compared mostly on simulated images but also from few real images acquired by cameras placed at different spots in the aircraft, such as shown in Figure 1.

The aim of this comparison is to discuss three points. The first is to compare the results of the three methods on several images, varying the environment and image quality. The second point is to question the use of an IPM transform as a preprocessing and its contribution to the improvement of the algorithms performances. The third point of this study is to challenge the three algorithms by providing images with degraded meteorological conditions.

Performance of the three methods

Figure 3 presents the different precision criteria computed from the results obtained from both real and simulated images without IPM preprocessing nor degraded weather conditions. Figures 5(a) and 5(b) present the precision criteria with IPM processing and Figures 5(c) and 5(d) present the precision criteria with degraded weather conditions for images provided by the cockpit and fin cameras. Figure 4 presents the maximum range of detection for each use case.

The particle filter method recall is mostly higher than the recall of the other methods. This value should be compared with the mean distance between the detection and the ground truth. The particle filter has a higher mean distance because it returns a bigger set of points in each cluster. It increases the recall after the convergence of the clusters around the patterns to detect but also impact the precision. On contrary, the Hough Transform method returns few points, which are mostly well detected but a few false detections can lower significantly the recall of the method. LaneNet performances seem bad at first view. However, these performances have to be put in opposition with our annotation system. Because of its nature, the network tries to extrapolate more lines to describe a whole lane. Our aim is to detect only lines and the addition of lines from the LaneNet algorithm to recreate a lane strongly impacts the recall and precision.

Impact of the IPM transform on the performances

The IPM transformation requires a good estimation of the camera position on the aircraft, which we have not been provided with for the real images. We evaluate the interest of the IPM transform for synthetic images, acquired from simulated fin and cockpit cameras with a resolution of 1280pX960p and an horizontal field of view of 80 degrees. We also provide results for the real images, with an approximation of the position and field of view of the camera in real images but it should be noted that precision in this information can greatly impact the result of the IPM transform and the performances. Figure 5 gives an overview of the impact of the IPM and the effect of degraded weather conditions on the results. The use of the IPM transform does not increase massively the recall results from the three methods. For LaneNet it considerably degrades the results, probably because of the nature of its training (TuSimple). It is predictable as the IPM transform adds noise to the image when trying to reconstruct an information not present in the initial image.

While this preprocess increases the detection range, it also decreases the precision and recall of each methods. It is due to the interpolations used to commute from the initial view to the bird's eye view. A pixel chosen in the IPM view is computed from two interpolations. When compared to the ground truth in the initial view, it is assumed that the precision of the detection will be impacted. The results of the IPM show that this preprocess can benefit the detection when the position of the camera is well known. However, while applied to the real images, with an approximated position of the camera, we noticed that the IPM transform is more likely to degrade the detection.

We observed an impact of the IPM on simulated images on the long range detection. Results from the top part of the images have been greatly improved. On contrary, the IPM impact is not relevant on this part of the images from the cockpit camera. Thoses pixels are far from the airplane and they are not distinguishable for the two cameras. The IPM transform adds noise, increasing the detection range for those pixels.

Impact of the degraded weather conditions

For the first two methods, we can observe a drop in maximum detection range. This was predictable as these methods strongly rely on color assumptions that are impacted by variations of the visibility in the environment. However, LaneNet does not seem to be strongly affected, in terms of results only. Qualitatively, the algorithm detects less lines and it virtually increases the precision (see Figure 6). Figure 6(b) represents the detection for clear weather. We can see that the Particle Filter and Hough Transform methods detect lines homogeneously in the image. LaneNet encounters difficulties with this image as the pattern of the lines is specific to aeronautics and is not found in the automotive field. It performs at the top of the image where we can found two parallel lines, a mostly found configuration in automotive datasets. The detection in Figure 6(c) is sparse. The lines detections are reduced as the contrast between the line and the tarmac is more difficult to observe.

Qualitative comparison

In Figure 7, we selected a fin camera image and a cockpit camera to present a qualitative comparison of the algorithm where the results of the particle filter are represented in red, the results of the Hough transform are represented in blue and the results of LaneNet are represented in green.

We observe that the three methods offer a satisfactory performance on real images. However, the particle filter is subject to false detection on far range pixels due to the color assumption. The Hough transform method does not perform as well as the other methods on far range detection for fin camera images. LaneNet recreate lines based on its training. It enables the algorithm to be resistant to occlusions but degrades virtually the precision.

CONCLUSIONS

In this paper, we described three algorithms used in the literature for line detection. We compared those algorithms on simulated and real images, obtained by mounting cameras on several positions in an aircraft. We also studied the impact of performing an IPM transformation on the input image before running the algorithms and the impact of degraded weather conditions. The results presented above have been discussed based on average values of the detection on several images. For the particle filter method, we observe that results on simulated images are satisfactory but this evaluation must be mitigated on real images. A number of thresholds have to be fixed in this method which implies difficulties to adapt for real images or degraded weather conditions. The Hough Transform method is based on a strong color assumption. It is easy to fix on simulation but harder to fix in real images, as many other computer vision applications. LaneNet provides unsatisfactory results as expected. The change of field from automotive to aeronautics prevents the network from detecting a good number of lines because of the differences in line patterns. The IPM Transform requires a precise setup. It can improve detection range in simulated images but minor changes in parameters on real images can heavily degrade the results. Lastly, the methods still give good results at a near range but the degraded weather conditions have a great impact on middle and far range detections.

In future studies, we will work on the color assumption. Combining ICAO information and a method to perform an automatic detection of a reference color in the image should improve results on real images. We also plan to implement a new algorithm for the particle filter based on an idea presented in [START_REF] Tran | A versatile object tracking algorithm combining Particle Filter and Generalised Hough Transform[END_REF] where the particle can be expressed as a window of pixels instead of a simple pixel. Combined with a better appreciation of the line color, it should reduce the false detections and improve the precision. We plan to increase the number of test images as it will be beneficial for further CNN-based studies. It will enable us to retrain the LaneNet network to adapt to the aeronautic context. We are currently working on an estimation of the angular variation of the aircraft position in regards to the centerline of the tarmac. It can be added as another criterion to select the preferred line detection algorithm and it could be used to perform a lane departure estimation algorithm.

Figure 1 :

 1 Figure 1: Camera positions in the aircraft fin and in the cockpit

Figure 2 :

 2 Figure 2: (a) Original image, (b) IPM transform

Figure 3 :

 3 Figure 3: Mean distance to the ground truth and recall performances for real and simulated images from cockpit and fin cameras. Results in blue for Hough Transform, red for the Particle Filter and green for LaneNet.

Figure 6

 6 Figure 6 is based on Figure2(a) and offers a comparison of line detection of the same scene in two different weather conditions (either normal conditions or dusk and some fog). Figure6(b) represents the detection for clear weather. We can see that the Particle Filter and Hough Transform methods detect lines homogeneously in the image. LaneNet encounters difficulties with this image as the pattern of the lines is specific to aeronautics and is not found in the automotive field. It performs at the top of the image where we can found two parallel lines, a mostly found configuration in automotive datasets. The detection in Figure6(c) is sparse. The lines detections are reduced as the contrast between the line and the tarmac is more difficult to observe.In Figure7, we selected a fin camera image and a cockpit camera to present a qualitative comparison of the algorithm where the results of the particle filter are represented in red, the results of the Hough transform are represented in blue and the results of LaneNet are represented in green.We observe that the three methods offer a satisfactory performance on real images. However, the particle filter is subject to false detection on far range pixels due to the color assumption. The Hough transform method does not perform as well as the other methods on far range detection for fin camera images. LaneNet recreate lines based on its training. It enables the algorithm to be resistant to occlusions but degrades virtually the precision.

Figure 4 :Figure 5 :

 45 Figure 4: Maximum range detection for all situations (blue for Hough Transform, red for Particle Filter and green for LaneNet)

Figure 6 :Figure 7 :

 67 Figure 6: (a) Degraded conditions. Results of yellow area (blue for Hough Transform, red for the Particle Filter and green for LaneNet) for (b) clear weather and (c) degraded weather

ACKNOWLEDGEMENTS

We would like to thank the OKTAL-SE company for providing the simulation tool, used to build synthetic images. C. Meymandi Nejad is supported by the project COCOTIER funded by the Direction Générale de l'Aviation Civile depending on the French government and E. Perrotin is supported by a grant from AIRBUS OPERATIONS. We would also like to thank the AIRBUS OPERATIONS company for providing real images from test planes.