Bayesian Estimation of the Multifractality Parameter for Image Texture Using a Whittle Approximation

Texture characterization is a central element in many image processing applications. Multifractal analysis is a useful signal and image processing tool, yet, the accurate estimation of multifractal parameters for image texture remains a challenge. This is due in the main to the fact that current estimation procedures consist of performing linear regressions across frequency scales of the 2D dyadic wavelet transform, for which only a few such scales are computable for images. The strongly non-Gaussian nature of multifractal processes, combined with their complicated dependence structure, makes it difficult to develop suitable models for parameter estimation. Here, we propose a Bayesian procedure that addresses the difficulties in the estimation of the multifractality parameter. The originality of the procedure is threefold. The construction of a generic semiparametric statistical model for the logarithm of wavelet leaders; the formulation of Bayesian estimators that are associated with this model and the set of parameter values admitted by multifractal theory; the exploitation of a suitable Whittle approximation within the Bayesian model which enables the otherwise infeasible evaluation of the posterior distribution associated with the model. Performance is assessed numerically for several 2D multifractal processes, for several image sizes and a large range of process parameters. The procedure yields significant benefits over current benchmark estimators in terms of estimation performance and ability to discriminate between the two most commonly used classes of multifractal process models. The gains in performance are particularly pronounced for small image sizes, notably enabling for the first time the analysis of image patches as small as 64 × 64 pixels.

otherwise infeasible evaluation of the posterior distribution associated with the model.Performance is assessed numerically for several 2D multifractal processes, for several image sizes and a large range of process parameters.The procedure yields significant benefits over current benchmark estimators in terms of estimation performance and ability to discriminate between the two most commonly used classes of multifractal process models.The gains in performance are particularly pronounced for small image sizes, notably enabling for the first time the analysis of image patches as small as 64 × 64 pixels.

Index Terms
Texture characterization, Multifractal analysis, Wavelet leaders, Bayesian estimation, Whittle approximation, Multiplicative cascades, Fractional Brownian motion I. INTRODUCTION

A. Context and motivation
Since the early origins of digital image processing, texture has been recognized as one of the central characteristic features in images.There is no common definition for texture, and different paradigms have been introduced in the literature [1].Several authors have proposed to model texture using random fractals, scale invariance or self-similarity [2], [3].Indeed, it has been reported in the literature that scale invariant processes are relevant and effective models for textures associated with a large class of natural images, see, e.g., [4]- [6].
The concepts of scale invariance and self-similarity are tied deeply to the degree of pointwise singular behavior or local regularity of the image amplitudes [7], [8].It has long been recognized that multiscale and wavelet analyzes constitute ideal tools to study data regularity [8]- [12].It is therefore not surprising that these tools play a central role not only for the study of image contours (edges), but also for texture characterization [13]- [15].Yet, while contours are essentially isolated singularities, the texture models consist of densely interwoven sets of singularities of different regularity strength.
Multifractal analysis provides a mathematical framework for the study of such spatial fluctuations of local regularity and texture characterization is therefore nowadays often conducted using this tool [16], [17].
Multifractal analysis.The local regularity of an image X is commonly measured using the socalled Hölder exponent h(t) [8], [11].Qualitatively, the smaller h(t 0 ), the rougher X is at spatial location t 0 and the larger h(t 0 ), the smoother X is at t 0 .The goal of multifractal analysis is the estimation of the multifractal spectrum D(h), which provides a global description of the spatial fluctuations of h(t).It is defined as the collection of the fractal dimensions of the sets of points for which the Hölder exponent takes the same value [8], [11], cf., Section II-A.Multifractal analysis has when c 2 < 0, while D(h) = δ(h − c 1 ) when c 2 ≡ 0 (c 2 cannot be positive theoretically [8], [11], [27]).
Estimation of c 2 .The leading order coefficients c p provide a relevant summary of the multifractal properties of X in applications where it would often not be convenient to handle an entire function D(h) [12], [26]- [28].The first log-cumulant c 1 , for instance, is the mode of D(h) and can be read as a measure for the "average" smoothness of X.More importantly, the coefficient c 2 , referred to as the multifractality or intermittency parameter, is directly related to the width of D(h) and captures the multifractal signature (i.e., the fluctuations of the local regularity) of the image X.Its primary importance stems from the fact that it enables the identification of the two major classes of multifractal stochastic processes: self-similar processes for which c 2 = 0 and multifractal multiplicative cascade (MMC) based processes for which c 2 is strictly negative [29].While the former class is tied deeply to additive constructions, the latter is based on multiplicative constructions and is hence linked to fundamentally different physical principles [7], [8], [30].Moreover, the magnitude of c 2 quantifies the degree of multifractality of an image for the latter class.For an overview and details on scale invariant and multifractal processes, the reader is referred to, e.g., [8], [31] and references therein.
In the seminal contribution [27], it has been shown that the log-cumulants c p are tied to the quantities (j, k) through the key relation Cum p [log (j, k)] = c 0 p + c p log 2 j , where Cum p [•] is the p-th order cumulant.In particular Relation (2) leads to the definition of the current standard and benchmark estimator for the parameter c 2 , based on linear regression of the sample variance, denoted by Var, of log (j, k) over a range of where w j are suitably defined regression weights [12], [28].
Limitations.The use of multifractal analysis remains restricted to images of relatively large size (of order 512 2 pixels) because a sufficient number of scales j must be available to perform the linear regression (3).While a similar issue is encountered for the analysis of 1D signals, it is significantly more severe for images: indeed, modulo border effects of the wavelet transform, the number of available scales is proportional to the logarithm of the number of samples for 1D signals and to the logarithm of the square root of the number of pixels for an image.For instance, for a 1D signal with 256×256 = 65536 samples, j 2 = 13 or 14 scales can be computed, while j 2 = 4 or 5 for an image of N ×N = 256×256 pixels.In addition, the finest scale, j = 1, should not be used in (3), see, e.g., [32].
The practical consequences for the multifractal analysis of images are severe: First, images of small size and thus image patches cannot be analyzed in practice.Second, (3) yields modest performance for images when compared with 1D signals of equivalent sample size [12], making it difficult to discriminate between c 2 ≡ 0 and values c 2 < 0 that are encountered in applications (typically, c 2 lies between −0.01 and −0.08).The goal of this work is to propose and validate a novel procedure for the estimation of c 2 for images that addresses these difficulties.

B. Related works
There are only few works reported in the literature that attempt to overcome the limitations of multifractal analysis for images described above.The generalized method of moments has been proposed and studied in, e.g., [33]- [35] and formulates parameter inference as the solution (in the least squares sense) of an over-determined system of equations that are derived from the moments of the data.The method depends strongly on fully parametric models and yields, to the best of our knowledge, only limited benefits in practical applications.
Although classical in parameter inference, maximum likelihood (ML) and Bayesian estimation methods have mostly been formulated for a few specific self-similar and multifractal processes [36], [37].The main reason for this lies in the complex statistical properties of most of these processes, which exhibit marginal distributions that are strongly non-Gaussian distribution as well as intricate algebraically decaying dependence structures that remain poorly studied to date.The same remark is true for their wavelet coefficients and wavelet leaders, see, e.g., [38], [39].
One exception is given by the fractional Brownian motion (in 1D) and fractional Brownian fields (in 2D) (fBm), that are jointly Gaussian self-similar (i.e., c 2 ≡ 0) processes with fully parametric covariance structure appropriate for ML and Bayesian estimation.Examples of ML and Bayesian estimators for 1D fBm formulated in the spectral or wavelet domains can be found in [36], [37], [40], [41].For images, an ML estimator has been proposed in [42] (note, however, that the estimation problem is reduced to a univariate formulation for the rows / columns of the image there).
As far as MMC processes are concerned, [43] proposes an ML approach in the time domain for one specific process.However, the method relies strongly on the particular construction of this process and cannot easily accommodate more general model classes.Moreover, the method is formulated for 1D signals only.Finally, a Bayesian estimation procedure for the parameter c 2 of multifractal time series has recently been proposed in [44].Unlike the methods mentioned above, it does not rely on specific assumptions but instead employs a heuristic semi-parametric model for the statistics of the logarithm of wavelet leaders associated with univariate MMC processes.Yet, it is designed for and can only be applied to univariate time series of small sample size.

C. Goals and contributions
The use of fully parametric models for the data can be very restrictive in many real-world applications.Therefore, the goal and the main contribution of this work is to study a Bayesian estimation procedure for the multifractality parameter c 2 with as few as possible assumptions on the data (essentially, the relation ( 2)) that can actually be applied to real-world images of small as well as large sizes.To this end, we adopt a strategy that is inspired by [44] and develop the key elements that are required for its formulation for images.
First, we show by means of numerical simulations that the distribution of the logarithm of wavelet leaders log (j, k) of 2D MMC processes can, at each scale j, be well approximated by a multivariate Gaussian distribution.Inspired by the covariance properties induced by the multiplicative nature of cascade constructions, we propose a new generic radial symmetric model for the variance-covariance of this distribution.This second-order statistical model is parametrized only by the two parameters c 2 and c 0 2 in (2) and enables us to formulate estimation in a Bayesian framework.Second, we formulate a Bayesian estimation procedure for the parameter c 2 of images that permits to take into account the constraints that are associated with the proposed statistical model.To this end, an appropriate prior distribution is assigned to the parameter vector (c 2 , c 0 2 ) which essentially ensures that the variance (2) is positive.Additional prior information, if available, can easily be incorporated.The Bayesian estimators of c 2 associated with the posterior distribution of interest cannot be evaluated directly because of the constraints that the parameter vector (c 2 , c 0 2 ) has to satisfy.Therefore, we design a suitable Markov chain Monte Carlo (MCMC) algorithm that generates samples that are asymptotically distributed according to the posterior distribution of interest.These samples are in turn used to approximate the Bayesian estimators.More precisely, we propose a random-walk Metropolis-Hastings scheme to explore efficiently the posterior distribution according to the admissible set of values for c 2 and c 0 2 .Finally, the exact evaluation of the likelihood associated with the proposed model for the logwavelet leaders requires the computation of the inverse and the determinant of large dense matrices, which is numerically and computationally too demanding for practical applications.To obtain a stable and efficient algorithm that can actually be applied to images, following intuitions developed in the univariate case, cf.e.g., [36], we approximate the exact likelihood with a Whittle-type expansion that is adapted to the proposed model and can be efficiently evaluated in the spectral domain.
The proposed algorithm for the estimation of the multifractality parameter c 2 is effective both for small and large image sizes.Its performance is assessed numerically by means of Monte Carlo simulations for two classical and representative 2D MMC constructions, canonical Mandelbrot cascades (CMC) [7] and compound Poisson cascades (CPC) [45], using the most common multipliers and a large range of process parameters.Complementary results are provided for 2D fBms (that are self-similar but not MMC).Our results indicate that the proposed estimation procedure is robust with respect to different choices of process constructions and greatly outperforms (2), in particular for small images and for identifying a value c 2 ≡ 0. It enables, for the first time, a multifractal analysis of images whose sizes are as small as 64 × 64 pixels.
The remainder of this work is organized as follows.Section II summarizes the main concepts of multifractal analysis and the wavelet leader multifractal formalism.Section III introduces the statistical model and the Bayesian framework underlying the estimation procedure for the parameter c 2 of images, which is formulated in Section IV.Numerical results are given in Section V.In Section VI, the proposed procedure is applied to the patch-wise analysis of a real-world image, illustrating its potential benefits for practical applications.Finally, Section VII concludes this paper and presents some future work.

II. MULTIFRACTAL ANALYSIS OF IMAGES
Let X : R 2 → R denote the 2D function (image) to be analyzed.The image X is assumed to be locally bounded in what follows (see Section II-B for a practical solution to circumvent this prerequisite).

A. Multifractal analysis
Hölder exponent.Multifractal analysis aims at characterizing the image X in terms of the fluctuations of its local regularity, characterized by the so-called Hölder exponent, which is defined as follows [8], [11].The image X is said to belong to C α (t 0 ) if there exists α > 0 and a polynomial P t0 of degree smaller than α such that The Hölder exponent at position t 0 is the largest value of α such that this inequality holds, i.e., Multifractal spectrum.For large classes of stochastic processes, the Hölder exponents h(t) can be theoretically shown to behave in an extremely erratic way [11], [26].Therefore, multifractal analysis provides a global description of the spatial fluctuations of h(t) in terms of the multifractal spectrum D(h).It is defined as the Hausdorff dimension (denoted dim H ) of the sets of points at which the Hölder exponent takes the same value, i.e., For more details on multifractal analysis and a precise definition of the Hausdorff dimension, see, e.g., [11], [26].

B. Wavelet leader multifractal formalism
Historically, multifractal formalisms have been proposed based on increments or wavelet coefficients.These choices of multiresolution quantities lead to both theoretical and practical limitations, see [12], [28] for a discussion.Recently, it has been shown that a relevant multifractal formalism can be constructed from the wavelet leaders [11], [26], [28], which are specifically tailored for this purpose.
Wavelet coefficients.We assume that the image is given in form of discrete sample values X(k), yields the approximation coefficients D X (j, k), whereas the high-pass filters defined by X (j, k), m = 0, . . ., 3 are obtained by convolving the image X with G (m) , m = 0, . . ., 3, and decimation; for the coarser scales j ≥ 2 they are obtained iteratively by convolving G (m) , m = 0, . . ., 3, with D (0) X (j − 1, •) and decimation.For scaling and multifractal analysis purposes, the approximation coefficients D (0) X are discarded and it is common to normalize the wavelet coefficients according to the L 1 -norm so that they reproduce the self-similarity exponent for self-similar processes [17].For a formal definition and details on (2D) wavelet transforms, the reader is referred to [13], [46].
Wavelet leaders.Denote as the dyadic cube of side length 2 j centered at k2 j and the union of this cube with its eight neighbors.The wavelet leaders are defined as the largest wavelet coefficient magnitude within this neighborhood over all finer scales [11] (j, k) ≡ (λ j,k ) sup Wavelet leaders reproduce the Hölder exponent as follows where λ j,k (t 0 ) denotes the cube at scale j including the spatial location t 0 [11].It has been shown that ( 8) is the theoretical key property required for constructing a multifractal formalism, see [11] for details.In particular, it can be shown that the wavelet leader multifractal formalism (WLMF), i.e., the use of ( 1) with coefficients c p estimated using wavelet leaders, is valid for large classes of multifractal model processes, see [12], [28] for details and discussions.The WLMF has been extensively studied both theoretically and in terms of estimation performance and constitutes the benchmark tool for performing multifractal analysis, cf.e.g., [12], [28].
Negative regularity.The WLMF can be applied to locally bounded images (equivalently, to images with strictly positive uniform regularity) only, see [12], [28] for precise definitions and for procedures for assessing this condition in practice.However, it has been reported that a large number of realworld images do not satisfy this prerequisite [5], [12].In these cases, a practical solution consists of constructing the WLMF using the modified wavelet coefficients in (7).When α is chosen sufficiently large, the WLMF holds (see [12] for details about the theoretical and practical consequences implied by this modification).

III. BAYESIAN FRAMEWORK
In this section, a novel empirical second-order statistical model for the logarithm of wavelet leaders for 2D MMC processes is proposed.This model is the key tool for estimating the multifractality parameter c 2 in a Bayesian framework.

A. Modeling the statistics of log-wavelet leaders
Marginal distribution model.It has recently been observed that for 1D signals the distribution of the log-wavelet leaders can be reasonably well approximated by a Gaussian distribution [44].Here, we numerically investigate the marginal distributions of l(j, •) for 2D images.To this end, a representative selection of scaling processes (the MMC processes CMC-LN, CMC-LP, CPC-LN and CPC-LP as well as fBm) have been analyzed for a wide range of process parameters (see Section V-A for a description of these processes).Representative examples of quantile-quantile plots of the standard Normal distribution against empirical distributions of log-wavelet leaders (scale j = 2) associated with CPC-LN, CPC-LP and fBm are plotted in Fig. 1 (upper row).
Clearly, the normal distribution provides, within ±3 standard deviations, a reasonable approximation for the marginal distribution of log-wavelet leaders of images for both members of the MMC class.
It is also the case for the fBm, a Gaussian self-similar process that is not a member of MMC.Note that the fact that the marginal distributions of the log-wavelet leaders are approximately Gaussian for scale invariant processes confirms the intuitions formulated by Mandelbrot [47].However, it is not a trivial finding: There is no a priori reason for this property even if the analyzed stochastic process has log-normal marginals (as is the case for CMC-LN, for instance).Indeed, it is not the case for the logarithm of the absolute value of wavelet coefficients whose marginal distributions are significantly more complicated and strongly depart from Gaussian, cf., Fig. 1 (bottom row).
Variance-covariance model.We introduce a model for the covariance of the logarithm of 2D wavelet leaders for MMC processes at fixed scale j denotated as Cov[l(j, k), l(j, k + ∆k)].It is motivated by the asymptotic covariance of the logarithm of multiscale quantities generically associated with multiplicative construction (c.f.[7]), studied in detail for wavelet coefficients of 1D random wavelet cascades in [48], and also by recent numerical results obtained for the covariance of the logarithm of 1D wavelet leaders for MMC processes [44].These results suggest a linear decay of Cov[l(j, k), l(j, k + ∆k)] in log coordinates log ∆k, with slope given by the parameter c 2 .Numerical simulations with 2D MMC processes for a wide range of process parameters (detailed in Section V-A) indicate that the empirical intra-scale covariance is radially symmetric and decays as c 2 log ∆r with ∆r ||∆k|| for an intermediary range of values ∆r given by 3 < ∆r ≤ ∆r max j Cov[l(j, k), l(j, k + ∆k)] ≈ (1) where ∆r max j = √ 2( √ n j − 1) and n j ≈ N 2 /2 2j denotes the number of wavelet leaders at scale 2 j of an N × N image.The constant γ is found to be well approximated by using the heuristic condition  The theoretical variance of the log-wavelet leaders is given by C 2 (j) = C 2 (j; c 2 , c 0 2 ) defined in (2).Finally, the short-term covariance is modeled as a line connecting C 2 (j; c 2 , c 0 2 ) at ∆r = 0 and (1) Combining ( 2), ( 11) and ( 12) yields the following full model for the covariance, parametrized by two Here, only the positive portions of (1) j are considered for numerical reasons (conditioning of the covariance matrix).The proposed covariance model is illustrated in Fig. 2 (left) for CMC-LN.
The joint Gaussian model with covariance model ( 13) assumes limited information on the dependence between different scales; essentially the variance (2), and the corresponding covariance matrix model for log-wavelet leaders at several scales j ∈ [j 1 , j 2 ] has the block-diagonal structure.For convenience and without loss of generality, the formulations given below and in Section IV will be stated in block-diagonal form yet could be extended without difficulty to any other valid covariance matrix model.

B. Likelihood, prior and posterior distributions
We focus on the estimation of the parameter c 2 and therefore work with centered log-wavelet leaders below.Let l j denote the vector of the n j centered coefficients l(j, k) − E[l X (j, .)]at scale j ∈ [j 1 , j 2 ], organized in lexicographic order, where E[•] stands for the sample mean.Let Σ j (θ) denote the corresponding n j × n j covariance matrix whose entries are given by the 2D parametric covariance function model (13).For convenience of notation, all coefficients are stacked in a unique zero-mean vector L = [l T j1 , ..., l T j2 ] T .Likelihood.With the above notation and assumptions, the likelihood of l j is given by Using the independence between l(j, k) and l(j , k ) for j = j , the likelihood of L is given by Prior distribution.The parameter vector θ must be chosen such that the variances of l(j, k) are positive, C 2 (j) ≥ 0. We define the admissible set where and c m 2 , c 0,m 2 quantify the largest admissible values for c 2 and c 0 2 , parameters that need to be tuned by practitioners and may depend on the application considered.When no additional prior information is available regarding θ, a uniform prior distribution on the set A is assigned to θ Posterior distribution and Bayesian estimators.The posterior distribution of θ is obtained from the Bayes rule and can be used to define the Bayesian maximum a posteriori (MAP) and minimum mean squared error (MMSE) estimators given in ( 20) and ( 21) below.

IV. ESTIMATION PROCEDURE
The computation of the Bayesian estimators is not straight-forward because of the complicated dependence of the posterior distribution (18) on the parameters θ.Specifically, the inverse and determinant of Σ j in the expression of the likelihood ( 14) do not have a parametric form and hence (18) can not be optimized with respect to the parameters θ.In such situations, it is common to use a Markov Chain Monte Carlo (MCMC) algorithm generating samples that are distributed according to p(θ|L).These samples are used in turn to approximate the Bayesian estimators.

A. Gibbs sampler
The following Gibbs sampler enables the generation of samples {θ (t) } Nmc 1 that are distributed according to the posterior distribution (18).This sampler consists of successively sampling according to the conditional distributions p(c 2 |c 0 2 , L) and p(c 0 2 |c 2 , L) associated with p(θ|L).To generate the samples according to the conditional distributions, a Metropolis-within-Gibbs procedure is used with a random walk Gaussian proposal.For details on MCMC methods, the reader is referred to, e.g., [49].

B. Whittle approximation
The Gibbs sampler defined in subsection IV-A requires the inversion of the n j × n j matrices Σ j (θ) in (19) for each sampling step in order to obtain r c2 and r c 0 2 .These inversion steps are computationally prohibitive even for very modest image sizes (for instance, a 64×64 image would require the inversion of a dense matrix of size ∼ 1000 × 1000 at scale j = 1 at each sampling step).In addition, it is numerically instable for larger images (due to growing condition number).To alleviate this difficulty, we propose to replace the exact likelihood (15) with an asymptotic approximation due to Whittle [50], [51].With the above assumptions, the collection of log-leaders {l(j, •)} are realizations of a Gaussian random field on a regular lattice P j = {1, .., m j } 2 , where m j = √ n j .Up to an additive constant, the Whittle approximation for the negative logarithm of the Gaussian likelihood ( 14) reads [36], [51]- [53] − log p(l j |θ) ≈ p W (l j |θ) = 1 2 ω∈Dj log φ j (ω; θ) + I j (ω) n j φ j (ω; θ) (22) where the summation is taken over the spectral grid and φ j (ω; θ) is the spectral density associated with the covariance function j (∆r; θ), respectively.
Since no closed-form expression is available for φ j (ω; θ), it is obtained numerically by discrete Fourier transform (DFT) Frequency range.It is commonly reported in the literature that the range of frequencies used in ( 22) can be restricted.This is notably the case for the periodogram-based estimation of the memory coefficient of long-range dependent time series for which only the low frequencies can be used, see, e.g., [36], [54], [55].Similarly, in the present context, the proposed spectral density model φ j (ω; θ) yields an excellent fit at low frequencies and degrades at higher frequencies.This is illustrated in Fig. 3 where the average periodograms of l(j, k) for CPC-LN are plotted together with the model φ j (ω; θ).We therefore restrict the summation in (22) to low frequencies where the fixed parameter η approximately corresponds to the fraction of the spectral grid D j that is actually used.We denote the approximation of p W (l j |θ) obtained by replacing D j with D † j in ( 22) by p † W (l j |θ).The final Whittle approximation of the likelihood ( 15) is then given by the following equation up to a multiplicative constant.

V. NUMERICAL EXPERIMENTS
The proposed algorithm is numerically validated for several types of scale invariant and multifractal stochastic processes for different sample sizes and a large range of values for c 2 .

A. Stochastic multifractal model processes
Canonical Mandelbrot Cascade (CMC).CMCs [7] are the historical archetypes of multifractal measures.Their construction is based on an iterative split-and-multiply procedure on an interval; we use a 2D binary cascade for two different multipliers: First, log-normal multipliers W = 2 −U , where , where π λ is a Poisson random variable with parameter λ = − γ log 2 (β−1) (CMC-LP).For CMC-LN, the log-cumulants are given by c 1 = m + α, c 2 = −2m and c p = 0 for all β−1 − 1 and all higher-order log-cumulants are non-zero with Below, γ = 1.05 and β is varied according to the value of c 2 .Compound Poisson Cascade (CPC).CPCs were introduced to overcome certain limitations of the CMCs that are caused by their discrete split-and-multiply construction [4], [45].In the construction of CPCs, the localization of the multipliers in the space-scale volume follows a Poisson random process with specific prescribed density.We use CPCs with log-normal multipliers W = exp(Y ), where Y ∼ N (µ, σ) is a Gaussian random variable (CPC-LN), or log-Poisson CPCs for which multipliers W are reduced to a constant w (CPC-LP).The first log-cumulants of CPC-LN are given by c 1 = − µ + 1 − exp µ + σ 2 /2 +α, c 2 = −(µ 2 + σ 2 ), and c p = 0 for p ≥ 3. Here, we fix Fractional Brownian motion (fBm).We use 2D fBms as defined in [56].FBm is not a CMC process and is based on an additive construction instead.Its multifractal and statistical properties are entirely determined by a single parameter H such that c 1 = H, c 2 = 0 and c p = 0 for all p > 2, and below we set H = 0.7.

B. Numerical simulations
Wavelet transform.A Daubechies's mother wavelet with N ψ = 2 vanishing moments is used, and α = 1 in (9), which is sufficient to ensure positive uniform regularity for all processes considered.
Estimation.The linear regression weights w j in the standard esimator (3), denoted LF, have to satisfy the usual constraints j2 j1 jw j = 1 and j2 j1 w j = 0 and can be chosen to reflect the confidence granted to each Var nj [log (j, •)], see [12], [28].Here, they are chosen proportional to n j as suggested in [28].The Gibbs sampler is run with N mc = 7000 iterations and a burn-in period of N bi = 3000 samples.The bandwidth parameter η in (25) has been set to η = 0.3 following preliminary numerical simulations; these are illustrated in Fig. 4 where estimation performance is plotted as a function of η for LN-CMC (N = 2 8 bottom, N = 2 9 up) with two different values of c 2 (−0.02 in red, −0.08 in blue).As expected, η tunes a classical bias-variance tradeoff: a large value of η leads to a large bias and small standard deviation and vice versa.The choice η = 0.3 yields a robust practical compromise.
Performance assessment.We apply the LF estimator (3) and the proposed MAP and MMSE estimators (20) and (21) to R = 100 independent realizations of size N × N each for the above described multifractal processes.A range of weak to strong multifractality parameter values c 2 ∈ {−0.01, −0.02, −0.04, −0.06, −0.08} and sample sizes N ∈ {2 6 , 2 7 , 2 8 , 2 9 } are used.The coarsest scale j 2 used for estimation is set such that n j2 ≥ 100 (i.e., the coarsest available scale is discarded), yielding j 2 = {2, 3, 4, 5}, respectively, for the considered sample sizes.The finest scale j 1 is commonly set to j 1 = 2 in order to avoid pollution from improper initialization of the wavelet transform, see [32] for details.Performance is evaluated using the sample mean, the sample standard deviation and the root mean squared error (RMSE) of the estimates averaged across realizations  N = {2 8 , 2 9 }.The performance of the MAP estimator has been found to be similar to the MMSE estimator and therefore is not reproduced here for space reasons.

Estimation
First, it is observed that the proposed algorithm slightly but systematically outperforms LF in terms of bias.This reduction of bias does neither depend on a specific choice of the multifractal process or its parameters, or on the sample size.Second, and most strikingly, the proposed Bayesian estimators yield significantly reduced standard deviations, with a reduction of up to a factor of 3 as compared to linear regressions.The variance reduction is more important for small values of |c 2 | yet remains above a factor of 1.5 for large values of |c 2 |.
These performance gains are directly reflected in the overall RMSE values, which remain up to a factor of 2.5 below those of linear fits.Finally, note that the estimation performance for CMCs and CPCs with log-Poisson multipliers are found to be slightly inferior to those with log-normal multipliers.This may be due to an arguably slightly stronger departure from Gaussian for the former, cf.Fig. 2.
Performance for small sample size.For small sample sizes N ≤ 2 7 , the limited number of available scales forces the choice j 1 = 1.Results for N = {2 6 , 2 7 } are reported in Tab.II.They indicate that the performance gains of the proposed Bayesian estimators with respect to LF estimators are even more pronounced for small sample size, both in terms of bias and standard deviations, yielding a reduction of RMSE values of up to a factor of 4. In particular, note that LF yields biases that are prohibitively large to be useful in real-world applications due to the use of the finest scale j = 1, cf., [12].Notably, values c 2 = 0 cannot be reliably detected with LF.In contrast, the proposed Bayesian procedure yields sufficiently small bias and standard deviations to enable the estimation of the multifractality parameter c 2 even for very small images (or image patches) of size 64 × 64.The reported performance gains come at the price of an increased computational cost, with computation times of the order of 8s (N = 64) to 50s (N = 512) per image, respectively, on a standard desktop computer.
Performance for fractional Brownian motion.Self-similar fBms with c 2 = 0 do not belong to the class of MMC processes for which the proposed estimation procedure was designed.The correlation structure of the wavelet coefficients of fBms has been studied in, e.g., [57].This correlation is weak, i.e., it goes to zero fast with the distance between wavelet coefficients in the time-scale plane.

VII. CONCLUSION AND PERSPECTIVES
This paper proposed a Bayesian estimation procedure for the multifractality parameter of images.
The procedure relies on the use of novel multiresolution quantities that have recently been introduced for regularity characterization and multifractal analysis, i.e., wavelet leaders.A Bayesian inference scheme was enabled through the formulation of an empirical yet generic semi-parametric statistical model for the logarithm of wavelet leaders.This model accounts for the constraints imposed by multifractal theory and is designed for a large class of multifractal model processes.The Bayesian estimators associated with the posterior distribution of this model were approximated by means of samples generated by a Metropolis-within-Gibbs sampling procedure, wherein the practically infeasible evaluation of the exact likelihood was replaced by a suitable Whittle approximation.The proposed procedure constitutes, to the best of our knowledge, the first operational Bayesian estimator for the multifractality parameter that is applicable to real-world images and effective both for small and large sample sizes.Its performance was assessed numerically using a large number of multifractal processes for several sample sizes.The procedure yields improvements in RMSE of up to a factor of 4 for multifractal processes, and up to a factor of 10 for fBms when compared to the current benchmark estimator.The procedure therefore enables, for the first time, the reliable estimation of the multifractality parameter for images or image patches of size equal to 64 × 64 pixels.It is interesting to note that the Bayesian framework introduced in this paper could be generalized to hierarchical models, for instance, using spatial regularization for patch-wise estimates.In a similar vein, future work will include the study of appropriate models for the analysis of multivariate data, notably for hyperspectral imaging applications.
2D) orthonormal discrete wavelet transform (DWT) can be obtained as the tensor product of one-dimensional (1D) DWT as follows.Let G 0 (k) and G 1 (k) denote the lowpass and high-pass filters defining a 1D DWT.These filters are associated with a mother wavelet ψ, characterized by its number of vanishing moments N ψ > 0. Four 2D filters G (m) (k), m = 0, . . ., 3 are defined by tensor products of n j /4 ; c 2 )) = 0, where the operator truncates to integer values.

Fig. 4 .
Fig. 4. Influence of the bandwidth parameter, η, on estimation performance for two different sizes of LN-CMC and two different values of c2 values (−0.02 in red and −0.08 in blue).
performance.Tab.I summarizes the estimation performance of LF and MMSE estimators for CMC-LN, CPC-LN, CMC-LP, CPC-LP (subtables (a)-(d), respectively) and for sample sizes

Fig. 5 .
Fig. 5. Band #20 of hyperspectral datacube (a); estimates of c2 for overlapping 64 × 64 pixel patches obtained by MMSE (c) and LF (d), respectively, together with their histograms (e); zoom on the patches indicated by a red frame (b).The center of the image patches are indicated by white dots in the original image, the distance between two of the dots corresponds to one half of the patch size, axes labels indicate patch numbers.

FBm results are summarizedFig. 5 (
Fig.5(a), (c) and (d)).Specifically, considering the zoom in Fig.5(b) (equivalently, the corresponding textures in Fig.5(a) and (c)), the Bayesian estimates are spatially strongly homogeneous for the forested regions with visually homogeneous texture (e.g., upper right portions in Fig.5(b)) despite the modest patch size used for estimation and indicate a weak yet non-zero multifractality for these regions.Similar observations are obtained for the vegetation patches (e.g., bottom left corners in Fig.5(b)).Furthermore, the zones of mixed vegetation (e.g., upper left corner in Fig.5(b)) also yield spatially coherent and consistent estimates of c 2 , with more negative values (stronger multifractality).In contrast, the linear regression based estimates do not succeed in reproducing the spatial structure of the image and display strong variability throughout the image.Indeed, even for the homogeneous texture in the forested regions, LF yields strongly varying and spatially incoherent estimates.Finally, note that the strongly negative values of c 2 observed in the bottom left corner of Fig.5(c) correspond to regions consisting of both (textured) vegetations and of roofs of buildings (with close to zero amplitudes and no texture).