**Abstract** : 1. Motivation Scale invariance has been observed in numerous applications involving data of various and very different natures. It can be operationally defined as the power law behavior with respect to scale of the structure functions, which are given by the empirical moments of the absolute value of the multiresolution coefficients of the data at a given scale (cf. Eq. (1)). The estimation of the exponents characterizing these power laws - termed scaling exponents - constitutes the ultimate goal of the practical analysis of scale invariance (also called scaling analysis). These scaling exponents are then commonly involved in standard signal processing tasks, such as detection, identification, or classification. In practice, scaling analysis is often conducted within the theoretical framework of multifractal analysis. In a nutshell, multifractal analysis aims at characterizing the fluctuation (in time or space) of the local regularity of the process under analysis through analysis of the (power law) behavior of the structure functions in the limit of fine scales, Though multifractal analysis can theoretically be extended to dimensions higher than 1 without technical difficulties, most practical implementations remain restricted to one dimensional signals. This is mainly due to the fact that multifractal analysis requires the use of a range of both positive and negative empirical moments, hence demanding for multiresolution quantities with adequate properties. To date, the only practically available procedure for the multifractal analysis of 2D signals, hence images, is the so called Wavelet Transform Modulus Maxima (WTMM) procedure (based on the skeleton of a continuous wavelet transform (CWT)). Yet, the WTMM procedure suffers from a number of theoretical and practical drawbacks: It has a high computational cost; The calculation of the 2D CWT skeleton requires involved theoretical definitions as well as a cumbersome practical procedure; It is still lacking a theoretical support. Therefore, in numerous applications where the data are naturally images, multifractal analysis remains restricted to 1D slices of the data and hence incomplete. In the present contribution, elaborating on previous results obtained for (1D) signals, we propose a practical multifractal analysis method for (2D) images based on two key features: The use of a 2D Discrete Wavelet Transform (DWT) (instead of a 2D CWT); The replacement of wavelet coefficients with wavelet Leaders. This yields two major benefits: The computation cost is very low; Wavelet Leaders have been shown to yield a complete and rigorous analysis of the multifractal analysis of bounded functions. This is because wavelet Leaders consist of monotonous increasing quantities that finely account for the irregularities of the analyzed function. The aims of the present contribution are twofold: First, studying the necessary theoretical elements, validity and limitations of a wavelet Leader based multifractal analysis of images, and second, the evaluation of its practical statistical performance. 2. Multifractal Analysis: Theory Multifractal analysis of a bounded function consists of describing the fluctuation of its local regularity with respect to space. The local regularity is measured by the Holder exponent. The variability of the Holder exponents with respect to space is characterized by the multifractal spectrum. It is defined as the Hausdorff dimension of the ensemble of points at which the function has the same Holder exponent. Practical multifractal analysis essentially aims at obtaining the multifractal spectrum from a single finite size image. The corresponding theoretical elements are given in Section 2.1. Since the multifractal spectrum is numerically not directly accessible, so-called multifractal formalisms have to be used: They recast the estimation to quantities that can actually be numerically calculated in practice and allow - under certain mathematical hypotheses - the estimation of the multifractal spectrum, For functions that are uniform Holder (i.e., functions whose minimal regularity exponent, as defined in Eqs. (3) and (4), is positive), it has been shown that the wavelet Leaders locally reproduce the Holder exponent through a power law behavior in the limit of fine scales. They should hence constitute the key ingredient for the wavelet Leader multifractal formalism for images. Wavelet Leaders are defined as the supremum of the coefficients of a discrete wavelet transform in a small (time or space) neighborhood over all finer scales. First, the wavelet Leader structure functions are calculated from the image, both for positive and negative moments, and the corresponding scaling exponents are estimated. By taking their Legendre transform, one obtains the Legendre spectrum. The multifractal formalism consists now of stating that the Legendre spectrum is identical with the multifractal spectrum. This is described in Section 2.2. A practical procedure to measure the minimal regularity exponent is given in Eq. (5). This interpretation of the Legendre spectrum as the multifractal spectrum of the image is however not of general validity. Yet, it has been proven to be valid for self-similar finite variance stationary increments processes, and is conjectured for multiplicative martingale based processes. Also, it has be shown that the Legendre spectrum is independent of the choice of the mother wavelet, and hence provides an intrinsic characterization of the data, regardless of its interpretation. It is shown in Section 2.3 that the multifractal formalism holds without further hypothesis only for monofractal processes, i.e., when both the Legendre and multifractal spectra reduce to a point. Also, it is shown that one can always construct a process with any admissible Legendre spectrum, whose multifractal and Legendre spectra do not coincide. 3. Positive measures and tractional integration Numerical images can naturally be modeled with positive measures. This violates the uniform Holder assumptions, and the wavelet Leader multifractal formalism can not be directly applied: Theoretically, wavelet Leaders become infinite; In practice, they are always finite, yet their use would yield meaningless results. This difficulty can be related to the existence of negative Holder exponents in the data (explicit definition detailed in Section 3.1). It is shown that performing a fractional integration of the data, with an order larger than the opposite of the minimal regularity exponent, guarantees uniform Holder regularity, and hence that the wavelet Leader based multifractal formalism can be applied, It is also explained in Section 3.1 how fractional integration and wavelet Leaders can be merged into a single computationally inexpensive and numerically unproblematic operation, referred to as pseudo-fractional integration. The question that arises naturally, though, is in how far the Legendre spectrum obtained by this practical solution can be related to the original image. It is argued here that even when the multifractal formalism holds for the integrated image, the common heuristic association with a translated (by the order of integration) version of the Legendre spectrum can not, in general be interpreted as the spectrum of the image, since the image is not uniform Holder, Yet, if the image contains only cusp-type singularities, the Legendre spectrum is nevertheless intrinsically related to the original image, thus characterizes it precisely, and can be used as such. However, it is shown and discussed in Sections 3.2 and 3.3 that if the image contains oscillating singularities, this reasoning is not valid: fractionally integrating with different orders and translating the Legendre spectra does not define one unique function. It is argued here that this fact can potentially be exploited for detecting the presence of oscillating singularities in data. 4. Empirical Multifractal Analysis In practice, scaling analysis amounts to estimating scaling exponents. However, their estimation is often advantageously replaced by the estimation of a smaller set of specific parameters, termed the log-cumulants. Section 4.1 states how these log-cumulants can be used to approximate the Legendre spectrum with only a small number of parameters, and without having to calculate the Legendre transform numerically. In Sections 4.2 and 4.3, explicit forms of the wavelet coefficient and Leader based estimation procedures for multifractal analysis of images are given. In Section 4.4, the use of wavelet coefficients and wavelet Leaders for multifractal analysis are compared. 5. Numerical Simulations Finally, the performance of the practical wavelet Leader based multifractal attribute estimation procedures for images are studied numerically, by means of Monte Carlo simulations conducted on two reference stochastic processes, whose multifractal properties are theoretically known: 2D fractional Brownian motions, Gaussian self-similar hence monofractal processes; Mandelbrot's multiplicative cascades (or martingales), celebrated examples of multifractal processes (with both log-normal and log Poisson multipliers). These processes are defined in Section 5.1. Simulation parameters are chosen to match realistic values observed in a number of real world data and situations (cf. Section 5.2), 6. Estimation Performance From a large number (1000) of 2D realizations of fractional Brownian motion, and of Mandelbrot's multiplicative cascades, the estimation performance of the wavelet coefficient and wavelet Leader based multifractal formalisms are measured in terms of biases, variances and mean squared error. The results demonstrate that the proposed wavelet Leader based procedures outperform those based on wavelet coefficients, while obtained at the same computational cost, in estimation performance: Whereas the former enable the relevant estimation of scaling exponents for both positive and negative statistical orders, and hence of the entire spectrum, that latter fail for negative statistical orders and hence the right part of the spectrum (cf. Section 6.1). Also, the former have significantly better performance for the estimation of the log-cumulants, in particular for cumulants of order larger than 1, i.e., for cumulants that enable practitioners to discriminate between finite variance self similar processes and multiplicative cascades (cf, Section 6.2). Yet, there is a price to pay for such improvements: First, the scales that can be involved in the estimations are slightly restricted when using wavelet Leaders (cf. Section 6.3), and second, the variance of Leader based estimates of the self-similarity parameter of Gaussian self similar processes is not independent of the parameter, as is the case for coefficient based estimates (cf. Section 6.4). Section 6.5 illustrates the behavior of the wavelet Leader based multifractal analysis procedure on a real-world image, consisting of natural fern, chosen in reference to the classical simulated fern example used as a paradigm fractal object. 7. Conclusion and discussion This contribution proposes the definition of a wavelet Leader multifractal analysis procedure for 2D signals, i.e., for Images, It is shown to have estimation performance that significantly outperform those obtained form wavelet coefficients and to have both a low computational cost and a low level in implementation difficulty. No other procedure combining this high performance level and low computational and conceptual difficulty levels has been proposed so far. Multifractal analysis consists of a conceptually involved characterization of an image: It relates, via the so-called multifractal formalism, the scale invariance properties of the data to their sample path regularity properties. The exact use of this theoretical framework may turn practically difficult. However, this does not prevent practitioners from using a set of estimated multifractal attributes as relevant image features. This contribution hence opens the route towards a systematic use of multifractal analysis for image detection and classification tasks, as is tentatively experimented in various applications (e.g., computer vision).