Abdourrahmane M Atto

Fractional Order Stochastic Image Integration Models and Deep Learning of Multi-Fractional Interactions

Keywords: { Fractional integration, Stochastic elds, Spectral analysis, Texture synthesis, Spectral learning. [27], [28]

Fractional Order Stochastic Image Integration Models and Deep Learning of Multi-Fractional Interactions

I. Introduction

V ISUAL information is a concept where the notion of variation plays a central role: for both static and dynamic vision systems, perceptual information is associated with edges and changes in optical ow in order to infer shapes of objects (in a single image or a static scene) or scene reconguration (in a video or a dynamic vision system). When edges are abundant in a given scene, the observation refers to the notion of texture and the scene under consideration is more complex for visual analysis/compression/interpretation due to a large amount of variations to handle. A theory of visual information has to take into consideration, parsimonious descriptions for textures in terms of their descriptions with as few parameters as possible.

The issue addressed in this paper is some contributions in visual information theory by the proposal of a parsimonious framework for synthesizing non-trivial textures: 1) from basic inputs such as white noise samples and 2) when considering fractional integration operating on very few parameters. White noise has, by denition, a at spectral response whereas natural images exhibit wide spectral variability, specically in presence of texture. Starting from a white noise, one can build integration procedures so as to constrain the output for presenting shapes, forms and textures that deserve interests in visual perception.

Most relevant integration solutions in terms of deriving rich spectral content and visual texture impact involve stochastic fractional order frameworks [START_REF] Atto | -stationary texture synthesis from random eld modeling[END_REF], [START_REF] Polisano | \Texture modeling by gaussian elds with prescribed local orientation[END_REF], [START_REF] Kroese | \Spatial process generation[END_REF], [START_REF] Pu | \Fractional dierential mask: A fractional dierential-based approach for multiscale texture enhancement[END_REF]. Early works in the 1980s have concerned feature extraction from a one-dimensional fractional Brownian motion model [START_REF] Lundahl | \Fractional brownian motion: A maximum likelihood estimator and its application to image texture[END_REF], [START_REF] Kube | the imaging of fractal surfaces[END_REF], [START_REF] Stewart | \Fractional brownian motion models for synthetic aperture radar imagery scene segmentation[END_REF]. The 1990s decade has established the foundations of the rst spatial extensions in terms of isotropic fractional Brownian elds [START_REF] Kaplan | improved method for 2d self-similar image synthesis[END_REF], [START_REF] Pesquet-Popescu | \Wavelet packet decompositions for the analysis of 2-d elds with stationary fractional increments[END_REF], [START_REF] Kaplan | fractal analysis for texture classication and segmentation[END_REF] or estimating a fractal dimension from object and scene observations [START_REF] Penn | \Estimating fractal dimension with fractal interpolation function models[END_REF], [START_REF] Balghonaim | maximum likelihood estimate for two-variable fractal surface[END_REF]. In the following decade (2000s), eorts have been concentrated on dierent issues raised by fractional integral image processing, for instance: synthesis schemes in cases of either isotopic guidelines [START_REF] Pesquet-Popescu | \Synthesis of bidimensional α-stable models with long-range dependence[END_REF] or anisotropic concerns (limited to a pre-specied elliptical form) [START_REF] Eom | -correlation image models for textures with circular and elliptical correlation structures[END_REF], but also the highlights and benets of fractional order modeling for biomedical applications [START_REF] Jennane | \Fractal analysis of bone x-ray tomographic microscopy projections[END_REF], [START_REF] Campisi | texture modelling and synthesis using fractal processes[END_REF], [START_REF] Hollingsworth | \The best bits in an iris code[END_REF]. The 2010s have been ages of:

• optimization of algorithms for simulating [START_REF] Kroese | \Spatial process generation[END_REF] and characterizing transformations [START_REF] Atto | -d wavelet packet spectrum for texture analysis[END_REF] of the oldest but \consistent" isotropic fractional Brownian eld model; • investigations of dierent methods for parameter estimation [START_REF] Pu | \Fractional-order euler-lagrange equation for fractional-order variational method: A necessary condition for fractional-order xed boundary optimization problems in signal processing and image processing[END_REF], [START_REF] Combrexelle | \Bayesian estimation of the multifractality parameter for image texture using a whittle approximation[END_REF], [START_REF] Ivanovici | \Fractal dimension of color fractal images[END_REF]; • spreading out fractional analysis to image processing applications such as feature enhancements [START_REF] Pu | \Fractional dierential mask: A fractional dierential-based approach for multiscale texture enhancement[END_REF], [START_REF] Yu | denoising algorithm based on entropy and adaptive fractional order calculus operator[END_REF], [START_REF] Zachevsky | \Statistics of natural stochastic textures and their application in image denoising[END_REF], [START_REF] Zachevsky | \Single-image superresolution of natural stochastic textures based on fractional brownian motion[END_REF] and feature extraction for clustering, segmentation or pattern recognition purposes [START_REF] Uss | \Ecient rotation-scaling-translation parameter estimation based on the fractal image model[END_REF], [START_REF] Khan | fractional moments of multilook polarimetric whitening lter for polarimetric sar data[END_REF],

A. Discrete domain fractional order spatio-stochastic integration 1) Basics on discrete spatial fractional order integration: An integrated eld X can be seen as combinations of spatial contributions of an elementary eld Z (white noise in general) measured via a spatial shift operator B.

Discrete domain integrated elds involve sums:

• that depend only on integer powers of B (integer order integrals) or • may depend on non-integer real values at the power of B (fractional order integrals). The shifts associated with operator B require computing a sub-pixel measure for the latter case.

In practice, not any sum forms an integrated process: for integration terminology, operator B must be chosen so that (1l -B)f(t, s) behaves approximately as a discrete partial spatial derivative of the function f(t, s). The following provides spatial operators B that respect this requirement.

Denition 1 (Gradient shift operators): Operators

B Grad-V : X(t, s) -→ X(t -1, s) (1) B Grad-H : X(t, s) -→ X(t, s -1) (2) B Grad-D : X(t, s) -→ X(t -1, s -1) (3)
are associated with Vertical (Grad-V), Horizontal (Grad-H) and Diagonal (Grad-D) gradient components.

These gradient based operators apply strict \backshifts" and one can check that since

B Grad-V X(t, s) = X(t -1, s), then (1l -B Grad-V )f(t, s) = f(t, s) -f(t -1, s)
is the standard discrete partial derivative of function f with respect to variable t. Note also that

B k Grad-V f(t, s) = B Grad-V B Grad-V . . . B Grad-V k terms f(t, s) = f(t -k, s)
so that shifting to position t -k involves k times composition of the elementary operator B Grad-V .

Because the gradient operators B Grad-V , B Grad-H , B Grad-D apply back-shifts, they are called causal (their shifting does not involve an f(t , s ) with t > t or s > s). Shift operator B can also be chosen anti-causal or \neither causal, nor anti-causal" with possible neighborhood consideration: for instance, since the standard Laplace operator satises:

∆f = f(t + 1, s) + f(t -1, s) + f(t, s + 1) + f(t, s -1) -4f(t, s)
then one can derive, by seeking B Lap such that (1 -

B Lap )f(t, s) ∝ ∆f(t, s):
Denition 2 (Laplace shift operator): The operator B Lap dened by:

B Lap X(t, s) = 1 4 X(t + 1, s) + X(t -1, s) +X(t, s + 1) + X(t, s -1) (4)
is called Laplace shift operator.

The Laplace shift operator B Lap is a combination that applies on the 4 closest variables pertaining to the neighborhood of X(t, s), without taking the latter into consideration. The following shift operators called Canny and Prewitt are compositions including the random variable X(t, s) and are neither causal, nor anti-causal.

Denition 3 (Canny shift operators): Canny Vertical (Can-V) and Horizontal (Can-H) shift operators are dened by:

B Can-V : X(t, s) -→ X(t, s) + X(t -1, s) -X(t + 1, s) (5) B Can-H : X(t, s) -→ X(t, s) + X(t, s -1) -X(t, s + 1) (6)
Denition 4 (Prewitt shift operators): Prewitt Vertical (Pre-V) and Horizontal (Pre-H) shifts consist in: [START_REF] Stewart | \Fractional brownian motion models for synthetic aperture radar imagery scene segmentation[END_REF] and

B Pre-V X(t, s) = X(t, s) + 1 k=-1 X(t -1, s -k) - 1 k=-1 X(t + 1, s -k)
B Pre-H X(t, s) = X(t, s) + 1 k=-1 X(t -k, s -1) - 1 k=-1 X(t -k, s + 1) (8)
Remark 1 (From derivation to integration): For B = B Grad-V dened by Eq. ( 1), we have seen that

(1 -B Grad-V )f(t, s) = f(t, s) -f(t -1, s) is a discrete partial derivative. When considering the inverse (thus integral) operator (1 -B Grad-V ) -1
, we obtain, provided that function g used below is chosen so that the following sum exists:

(1 -B Grad-V ) -1 g(t, s) = k 0 B k Grad-V g(t, s) = k 0 g(t -k, s)
Remark 1 above illustrates that one can switch from partial derivation induced by the dierence operator (1 -B) to partial integration by inverting the latter and this, for any of the shift operators given by Eqs. (1), ( 2), (3), (4), ( 5), ( 6), [START_REF] Stewart | \Fractional brownian motion models for synthetic aperture radar imagery scene segmentation[END_REF] and [START_REF] Kaplan | improved method for 2d self-similar image synthesis[END_REF].

More generally, we can consider a derivation scheme represented by a function of operator, say G(B), and equations having the form G(B)f(t, s) = h(t, s). The following properties then holds true: if G(B) is a derivation operator, then f is an integral of h and vice versa. Note that G(B) can operate integer or fractional order derivation, for instance:

• (1 -B Grad-V )f [from integer power \1" at the power of (1 -B Grad-V )] is the integer order derivation of f whereas • (1 -B Grad-V ) 1/2 f
is the fractional order \1/2" derivation of f with respect to B. Finally, when h is the solution of (1 -B Grad-V ) λ f = h where 0 < λ < 1, then f is said to be associated with a fractional order λ integral of h.

2) Fractional order spatio-stochastic K-factor Gegenbauer series: In the following, we consider equations with the form G(B)X(t, s) = Z(t, s) where X, Z : N × N × Ω -→ R are random elds and G(B) behaves as a fractional order derivation operator (see Section II-A1). Field X will be said fractional integral of Z.

Denition 5 (K-factor Gegenbauer eld): Field X(t, s) will be called spatial (2D) K-factor B-Gegenbauer (K-B-G) eld if it satises:

K k=1 1l -2ν k B + B 2 b k X(t, s) = Z(t, s) (9)
where |ν k | 1 for any k = 1, 2, . . . , K and Z(t, s) is assumed to be a zero-mean white Gaussian eld.

In Eq. ( 9), a term 1l -2νB + B 2 b behaves as (1l -B) 2b when ν is close to 1. Hence K-factor Gegenbauer elds include standard fractionally integrated elds especially when one among the ν k , k = 1, 2, . . . , K is close to one. Assuming a revertible power series expansion of any 1l -2ν k B + B 2 b k , synthesis of a K-factor Gegenbauer involves computing:

X(t, s) = K k=1 1l -2ν k B + B 2 -b k Z(t, s) (10)
where:

1l -2ν k B + B 2 -b k = +∞ =0 b k † (ν k ) B (11)
and the Gegenbauer function b k † (•) being:

b k † (x) = /2 m=0 (2x) -2m (-1) m Γ (b k + -m) m! ( -2m)! Γ (b k ) (12)
with denoting the oor function ( z is the largest integer less than or equal to z) and the Gamma function is dened by:

Γ (x) = R + u x-1 e -u du (13)
Samples of K-factor B-Gegenbauer elds are given by Figure 1 when B is one among the operators dened by Eqs. ( 4), ( 5), ( 6), [START_REF] Stewart | \Fractional brownian motion models for synthetic aperture radar imagery scene segmentation[END_REF] and [START_REF] Kaplan | improved method for 2d self-similar image synthesis[END_REF]. Depending on the operator used, one can obtain either a stochastic texture without any privileged orientation (case of B-Gegenbauer) or textures associated with a main Vertical/Horizontal/Diagonal edge orientation. Note that white noise shows no edge and the texture properties is driven by the operator B selected. More precisely, these texture feature orientations can be inferred from the spectral properties of K-factor B-Gegenbauer elds given below.

Proposition 1 (Power Spectral Densities, PSD, of K-B-G elds [START_REF] Atto | et processus stochastiques[END_REF]): The PSD associated with K-factor B-Gegenbauer elds are given below (PSD notation γ

K-B-G X Lap X Can-V X Can-H X Pre-V X Pre-H
Fig. 1. 5-factor Gegenbauer (see Eq. ( 9)) Laplacien (see Eqs. 4), Canny (see Eqs. ( 5) and ( 6)) and Prewitt (see Eqs. [START_REF] Stewart | \Fractional brownian motion models for synthetic aperture radar imagery scene segmentation[END_REF] and ( 8)).

for eld notation X K-B-G ) respectively for the gradient, Laplace, Canny and Prewitt operators, where

θ k = 1 + ν k and ξ k = 1 -ν k .
(i) PSD of X K-B-G for gradient shift operators (see Eqs.

(1), ( 2), (3), respectively):

γ K-B Grad-V -G (ω 1 , ω 2 ) = σ 2 K k=1 2 -2b k K k=1 |cos ω 1 -ν k | 2b k (14) γ K-B Grad-H -G (ω 1 , ω 2 ) = σ 2 K k=1 2 -2b k K k=1 |cos ω 2 -ν k | 2b k (15)
γ K-B Grad-D -G (ω 1 , ω 2 ) = σ 2 K k=1 2 -2b k K k=1 |cos(ω 1 +ω 2 ) -ν k | 2b k (16)
(ii) PSD of X K-B-G for Laplace shift operator (see Eq. ( 4)):

γ K-B Lap -G (ω 1 , ω 2 ) = σ 2 K k=1 cos 2 ω 1 2 +cos 2 ω 2 2 -θ k 2 + θ k ξ k 2b k (17)
(iii) PSD of X K-B-G for Canny shift operators (see Eqs. [START_REF] Atto | -d wavelet packet spectrum for texture analysis[END_REF], [START_REF] Pu | \Fractional-order euler-lagrange equation for fractional-order variational method: A necessary condition for fractional-order xed boundary optimization problems in signal processing and image processing[END_REF], respectively):

γ K-B Can-V -G (ω 1 , ω 2 ) = σ 2 K k=1 2 -2b k K k=1 4 sin 4 ω 1 -4 ν k ξ k sin 2 ω 1 + ξ 2 k b k (18) γ K-B Can-H -G (ω 1 , ω 2 ) = σ 2 K k=1 2 -2b k K k=1 4 sin 4 ω 2 -4 ν k ξ k sin 2 ω 2 + ξ 2 k b k (19)
(iv) PSD of X K-B-G for Prewitt shift operators (see Eqs.

γ K-B Pre-V -G (ω 1 , ω 2 ) = σ 2 K k=1 2 -2b k K k=1 4(1 + 2 cos ω 2 ) 2 sin 2 ω 1 -ξ 2 k 2b k (20) γ K-B Pre-H -G (ω 1 , ω 2 ) = σ 2 K k=1 2 -2b k K k=1 4(1 + 2 cos ω 1 ) 2 sin 2 ω 2 -ξ 2 k 2b k (21)
Directed vertical, horizontal or diagonal edges occur when a given PSD has its energy concentrated on certain frequency bands or varies only across one single variable (case for gradient and Canny operators). In particular, gradient elds exhibit an innite number of poles along many axes but in the same direction towards:

• frequencies (ω 1 , ω 2 ) ∈ R 2 /ω 1 = acos(ν k ), k = 1, 2, . . . , K for eld X K-B Grad-V -G , • frequencies (ω 1 , ω 2 ) ∈ R 2 /ω 2 = acos(ν k ), k = 1, 2, . . . , K for eld X K-B Grad-H -G and • frequencies (ω 1 , ω 2 ) ∈ R 2 /ω 1 + ω 2 = acos(ν k ), k = 1, . . . , K for eld X K-B Grad-D -G .
K-factor Canny-Gegenbauer elds present strong similarities with their analog gradient elds by construction, but they admit no pole when |ν k | < 1 for any k.

From Proposition 1, only Laplace and Prewitt elds show more complex PSD. However, K-factor Laplace-Gegenbauer elds X K-B Lap -G have no pole (dominant frequency) when the absolute value of every ν k is < 1 and are isotropic up to the mapping: ω i → cos(ω i /2) (same property whatever the direction). The motivation of the following section is the proposal of spatial eld models with richer PSD content in terms of spectral variability as well as limiting the number of poles.

B. Fractionally integrated elds from double indexed Gegenbauer polynomials

For being able to generate textures with parsimonious spectral information, we need to control PSD behavior in specic frequencies of interest. This cannot be achieved straightforwardly when using Denition II-A2. To overcome this limitation, the main contribution of this section is the construction of double indexed spatial multifractional integral models. Denition 6 (K-factor doubly indexed Gegenbauer eld):

Field X(t, s) will be called K-factor Doubly Indexed Gegenbauer (K-DI-G) eld if it satises one of the equations given below: [START_REF] Zachevsky | \Statistics of natural stochastic textures and their application in image denoising[END_REF] where |µ k | 1 and |ν k | 1 for any k = 1, 2, . . . , K; B 1 and B 2 are spatial shift operators chosen conveniently; Z(t, s) is a zero-mean white Gaussian noise and 1 , 2 are positive real valued parameters.

K k=1 1 B -2 1 1l -2µ k B 1 + B 2 1 2 + 2 B -2 2 1l -2ν k B 2 + B 2 2 2 b k X(t, s) = Z(t, s) (22) or K k=1 1 B 2 2 1l -2µ k B 1 + B 2 1 2 + 2 B 2 1 1l -2ν k B 2 + B 2 2 2 b k X(t, s) = Z(t, s)
Eqs. ( 22) and [START_REF] Zachevsky | \Statistics of natural stochastic textures and their application in image denoising[END_REF] show a balance between two Gegenbauer polynomials, the balance being governed by parameters 1 and 2 . When B 1 and B 2 are causal operators (gradient operators for instance), then Eq. ( 22) presents backward and forward shifts whereas Eq. ( 23) highlights only backward shifts.

Synthesis of sample X(t, s) satisfying Eq. ( 22) requires computing:

X(t, s) = K k=1 1 B -2 1 1l -2µ k B 1 + B 2 1 2 + 2 B -2 2 1l -2ν k B 2 + B 2 2 2 -b k Z(t, s) (24)
thus, decomposing the fractional double indexed polynomial as series involving integer powers of B 1 and B 2 (eld Z(t, s) is available only for t, s ∈ Z in the discrete domain framework). For this purpose, we write:

X(t, s) = K k=1 1 B -1 1 + B 1 2 -µ k 2 + 2 B -1 2 + B 2 2 -ν k 2 -b k Z(t, s) (25)
and set, for i = 1, 2,

B i = √ i (B -1 i + B i )/2 (26)
In this respect, we are concerned by:

X(t, s) = K k=1 B 1 - √ 1 µ k 2 + B 2 - √ 2 ν k 2 -b k Z(t, s) (27)
which amounts, by noting

µ k = √ 1 µ k , ν k = √ 2 ν k and c k = µ 2 k + ν 2 k , to: X(t, s) = K k=1 c k -2µ k B 1 -2ν k B 2 + B 2 1 + B 2 2 -b k Z(t, s) ( 28
)
The latter is an extension of the Gegenbauer polynomial introduced in [30, Eqs. ( 1) and ( 18)]. Denote, for k = 1, 2, . . . , K,

X k (t, s) = c k -2µ k B 1 -2ν k B 2 + B 2 1 + B 2 2 -b k X k-1 (t, s) (29)
with the convention X 0 = Z. The recursive Eq. ( 29) implies X K = X and we have:

Proposition 2: X k (t, s) = lim M,N→+∞ X k,M,N (t, s) with X k,M,N (t, s) = M m,n=0 C b k m,n (µ k , ν k )B m 1 B n 2 X k-1 (t, s) (30)
where function C b m,n is dened by:

C b m,n (µ, ν) = m/2 p=0 n/2 q=0 (-1) p+q Γ (b + m + n -p -q) Γ (b)p! q! (m -2p)! (n -2q)! × (2µ) m-2p (2ν) n-2q c m+n-p-q+b (31)
and the Gamma function is given by Eq. ( 13).

Proof: See Appendix A.

At this stage, we should specify B 1 and B 2 in Eq. ( 30) before continuing. When

B 1 = B r 1 Grad-V and B 2 = B r 1 Grad-H ,
we have from Eq. ( 26),

B 1 B 2 Y(t, s) = √ 2 B 1 Y(t, s -r 2 ) + Y(t, s + r 2 ) 2
thus:

B 1 B 2 Y(t, s) = √ 1 √ 2 4 Y(t -r 1 , s -r 2 ) + Y(t + r 1 , s -r 2 ) +Y(t -r 1 , s + r 2 ) + Y(t + r 1 , s + r 2 )
In practice, computing

B m 1 B n 2 Y(t, s)
consists in recursive applications of operator B 1 B 2 followed by either mcompositions with B 1 or n-applications of B 1 , where = min{m, n}. Textures given in Figure 2 have been obtained by using r 1 = r 2 = 1.

Figure 2 highlights that the doubly indexed polynomials of Eq. ( 22) yield a more powerful texture synthesis framework than the single one given by Eq. ( 9). In order to give a more formal proof of the latter assertion, let us derive the PSD associated with K-DI-G. We have:

Proposition 3: For B 1 = B r 1 Grad-V and B 2 = B r 1
Grad-H , the PSD associated with a K-DI-G given by either Eq. ( 22) or ( 23) is:

γ(ω 1 , ω 2 ) = σ 2 K k=1 2 -2b k K k=1 1 (cos(r 1 ω 1 ) -µ k ) 2 + 2 (cos(r 2 ω 2 ) -ν k ) 2 b k (32)
Proof: The proof follows from straightforward computations, by taking the moduli of Eqs. ( 22), [START_REF] Zachevsky | \Statistics of natural stochastic textures and their application in image denoising[END_REF] and by replacing in these quantities, B by e ir ω for = 1, 2.

The PSD of Eq. ( 32) shows poles at frequencies

(ω 1 , ω 2 ) = acosµ k r 1 , acosν k r 2
for k = 1, 2, . . . , K. Thus, we can specify specic frequency points associated with energy concentration, which makes K-DI-G a more concise framework for the synthesis of textures with multiple orientated features.

III. Continuous spatial domain fractional order integration models

In this section, t, s are continuous time variables, t, s ∈ R + , R, . . . and X : (t, s) → X(t, s) = X ω (t, s) = X t,s (ω) ∈ R or C is a spatial stochastic eld assumed to be centered in what follows. We recall that a stochastic eld X = (X(•, •)) is Gaussian if the random variable X(t, s) is Gaussian for every couple (t, s) pertaining to the domain of interest. Assuming that X is Gaussian, we can focus on its second-order statistics and since X has been assumed centered, then this eld is entirely characterized by its autocorrelation function. Thus, for synthesizing sample realizations of X, we only need to know the analytic form of this autocorrelation function.

As in previous section (dedicated to discrete domain models), the aim addressed in reference [START_REF] Atto | -stationary texture synthesis from random eld modeling[END_REF] was the derivation of a continuous domain fractional eld framework associated with parsimonious spectral contents. This aim has led [START_REF] Atto | -stationary texture synthesis from random eld modeling[END_REF] to propose a class of generalized fractional Brownian elds, the generalization concerning a model with 1) an arbitrary number of Hurst parameters and 2) an arbitrary number of spectral singularity points. However, the mathematical sense of these generalized fractional Brownian elds remains unspecied at present time. This section proposes a sense to these elds as a main contribution.

Section recalls the construction that has led to generalize fractional Brownian elds and Section shows that the convolution involved in their construction converges in quadratic mean sense.

A. Generalized Fractional Brownian Fields (GFBF) GFBF are associated to spatial domain convolution and modulation operators over a sequence of Fractional Brownian Fields (FBF). Process G H,(u 0 ,v 0 ) is called Modulated Fractional Brownian Field (MFBF) if it can be written in the form:

G H,(u 0 ,v 0 ) (x 1 , x 2 ) = e iu 0 x 1 e iv 0 x 2 F H (x 1 , x 2 ) (33)
where F H is a centered isotropic FBF with autocorrelation

R F H (x 1 , x 2 , y 1 , y 2 ) = σ 2 2 (x 2 1 + x 2 2 ) H + (y 2 1 + y 2 2 ) H -(x 1 -y 1 ) 2 + (x 2 -y 2 ) 2 H (34)
Field G H,(u 0 ,v 0 ) is complex valued and results as the interaction between a stochastic FBF (see rst row of Figure 3) and a deterministic phase term (see second row of Figure 3). Field G H,(u 0 ,v 0 ) is centered and its autocorrelation

R G H,(u,v) (x 1 , x 2 , y 1 , y 2 ) = E[G H,(u,v) (x 1 , x 2 ) G H,(u,v) (y 1 , y 2 )] (35) simplies as R G H,(u 0 ,v 0 ) (x 1 , x 2 , y 1 , y 2 ) = R F H (x 1 , x 2 , y 1 , y 2 )e iu 0 (x 1 -y 1 ) e iv 0 (x 2 -y 2 ) (36)
where R F H (x 1 , x 2 , y 1 , y 2 ) is given by Eq. [START_REF] Bahr | \Inequalities for the rth moment of a sum of random variables, 1 r 2[END_REF]. Thus,

R G H,(u,v) (x 1 , x 2 , y 1 , y 2 ) = σ 2 2 e iu(x 1 -y 1 ) e iv(x 2 -y 2 ) (x 2 1 + x 2 2 ) H +(y 2 1 + y 2 2 ) H -(x 1 -y 1 ) 2 + (x 2 -y 2 ) 2 H (37)
The GFBF is an integral representation involving interactions of several MIFBF: consider two sequences, the rst

H Q = {H 1 , H 2 , . . . , H Q }
composed of Hurst parameters and the second

{(u 1 , v 1 ), (u 2 , v 2 ), . . . , (u Q , v Q )} associated with modulation parameters (poles) of Q in- dependant MIFBF G Hs,(us,vs) s=1,2,...,Q having auto- correlations: R G Hs ,(us ,vs ) s=1,2,...,Q . A GFBF E H Q has X 1 X 2 X 3 X 4 X 5
Fig. 2. 1 to 5-DI-G elds generated from random selection of µ k , ν k by using Eq. ( 30) with B 1 = B Grad-V and B 2 = B Grad-H dened by Eq. ( 1) and ( 2) respectively. .

H = 0.2 H = 0.5 H = 0.8
||E H 2 || ||E H 3 || ||E H 4 || ||E H 5 || ||E H 6 ||
Fig. 4. GFBF samples E Hq for q = 2, 3, 4, 5, 6. Texture E H 6 includes all features issued from E Hq , q 5, whether these features are visible or not.

been dened in [START_REF] Atto | -stationary texture synthesis from random eld modeling[END_REF] as the stochastic Gaussian process with autocorrelation given by the convolution:

R H Q = Q s=1 R G Hs ,(us ,vs ) (38)
GFBF E H Q can be seen as the spatial convolution, with respect to spatial variables, of Q independent MIFBF:

E H Q = Q s=1 G Hs,(us,vs) (39)
and this makes synthesis possible from convolution of samples of MIFBF (see examples of GFBF given by Figure 4). However, this convolution cannot be understood as a standard function operator since otherwise, non-regularity (due to stochasticity) and innite energies (autocorrelation of MIFBF has fractional polynomial form) make straightforward calculus inecient. The following section provides a mathematical sense to the convolution of two MIFBF in terms of integrals associated to kernel based measures.

B. On convolution of MIFBF (to GFBF)

The MIFBF is G H,(u 0 ,v 0 ) (x, y) = e i(u 0 x+v 0 y) F H (x, y), a zero-mean Gaussian process with autocorrelation given by Eq. (37). We want to provide a sense to

G H 1 ,(u 1 ,v 1 ) G H 2 ,(u 2 ,v 2 ) (t, s) = R G H 1 ,(u 1 ,v 1 ) (x, y)G H 2 ,(u 2 ,v 2 ) (t -x, s -y)dµ(x, y)
This sense will be associated to quadratic mean convergence. For the monovariate stochastic process case the complete steps of proof are given in Appendix B. Hereafter, we focus on the specic extension of Appendix B to the multivariate spatial case. The extension will be performed by:

• considering stochastic integrals Y k (t, s) = R G H k ,(u k ,v k ) (x, y) dµ Φ (t,s) (x, y) (40)
• providing a quadratic mean sense to the following convolution of these integrals

Y 1 Y 2 (t, s) = R Y 1 (x, y)Y 2 (t -x, s -y)dxdy
• associating the quadratic mean convergence of

G H 1 ,(u 1 ,v 1 ) G H 2 ,(u 2 ,v 2 ) to that of Y 1 Y 2 from properties of Φ.
The convergence of Y 1 Y 2 will be obtained hereafter by considering in Denition of Eq. ( 40), a separable measure dµ Φ (α,β) (x, y) = φ α (x)φ β (y)dxdy where φ is the function described in Appendix B-A (having suciently fast decay in time and frequency). From this measure and by imposing that φ annihilates at frequencies ±(u k , v k ): φ(±u k ) = φ(±v k ) = 0, we obtain:

1) Field Y(t, s) dened by Eq. ( 40) is a second order stochastic process : Indeed,

R Y (x, y, t, s) = R 4 R G H,(u,v) (x 1 , x 2 , y 1 , y 2 ) dµ Φ (x,y) (x 1 , y 1 )dµ Φ (t,s) (x 2 , y 2 ) = R 4 R G H,(u,v) (x 1 , x 2 , y 1 , y 2 )φ x (x 1 )φ y (y 1 ) φ t (x 2 )φ s (y 2 )dx 1 dx 2 dy 1 dy 2 (41)
where R G H,(u,v) is given by Eq. ( 37) and (u, v) ∈ {(u 1 , v 1 ), (u 2 , v 2 )}. Let us analyze the 3 additive terms involved by expanding Eq. (37) in Eq. (41) separately. First

R 4 e iu(x 1 -y 1 ) e iv(x 2 -y 2 ) x 2 1 + x 2 2 H φ x (x 1 )φ y (y 1 ) φ t (x 2 )φ s (y 2 )dx 1 dx 2 dy 1 dy 2 = Constant × φ y (u) φ s (v) = Constant × e -iyu e -isv φ(u) φ(v) = 0 (42) Second R 4
e iu(x 1 -y 1 ) e iv(x 2 -y 2 ) y 2 1 + y 2 2 H φ x (x 1 )φ y (y 1 )

φ t (x 2 )φ s (y 2 )dx 1 dx 2 dy 1 dy 2 = Constant × φ x (-u) φ t (-v) = Constant × e -ixu e -itv φ(-u) φ(-v) = 0 (43)
Thus, R Y reduces to the contribution of the third term, that is:

R Y (x, y, t, s) = (44) - R 4
dx 1 dx 2 dy 1 dy 2 e iu(x 1 -y 1 ) e iv(x 2 -y 2 ) φ x (x 1 )φ y (y 1 )φ t (x 2 )φ s (y 2 )

(x 1 -y 1 ) 2 + (x 2 -y 2 ) 2 H (45)
and from a change of variable:

R Y (x, y, t, s)= (46) - R 4 dz 1 dz 2 dy 1 dy 2 e iuz 1 e ivz 2 z 2 1 + z 2 2 H φ x (z 1 + y 1 )φ y (y 1 )φ t (z 2 + y 2 )φ s (y 2 ) (47)
which exists and is nite provided that z 2 φ(z) ∈ L 1 (R) ∪ L 2 (R) (fast decay required for φ). In addition, after some steps of calculus similar to those of Appendix B-B, we derive:

R Y (x, y, t, s) = 1 4π 2 R 2 e i(x-t)ω 1 e i(y-s)ω 2 γ H,(u,v) (ω 1 , ω 2 ) (48) φ(ω 1 ) 2 φ(ω 2 ) 2 dω 1 dω 2 where γ H,(u,v) (ω 1 , ω 2 ) = σ 2 2 1+2H sin(πH)Γ 2 (1 + H) ((ω 1 -u) 2 + (ω 2 -v) 2 ) H+1 (49)
Thus, R Y (x, y, t, s) = R Y (x -t, y -s) so that stationarity holds true for Y and we can use the following twodimensional version for its autocorrelation:

R Y k (t, s) = 1 4π 2 R 2 e itω 1 e isω 2 γ H,(u,v) (ω 1 , ω 2 ) φ(ω 1 ) 2 φ(ω 2 ) 2 dω 1 dω 2 (50)
associated to the following PSD:

γ Y (ω 1 , ω 2 ) = γ H,(u,v) (ω 1 , ω 2 ) φ(ω 1 ) 2 φ(ω 2 ) 2 (51)
2) Convolution Y 1 Y 2 is a second order stochastic process: The convolution Y 1 Y 2 being with zero-mean, it is a well-dened second order process if

R Y 1 R Y 2 (x 1 , x 2 , y 1 , y 2 )
exists and is nite for every x 1 , x 2 , y 1 , y 2 ∈ R. Indeed, in this case:

E [Y 1 Y 2 (x 1 , x 2 ) Y 1 Y 2 (y 1 , y 2 )] = E R 2 Y 1 (x, y)Y 2 (x 1 -x, x 2 -y)dxdy R 2 Y 1 (t, s)Y 2 (y 1 -t, y 2 -s)dtds = R 4 dxdydtdsE[Y 1 (x, y)Y 1 (t, s)] E[Y 2 (x 1 -x, x 2 -y)Y 2 (y 1 -t, y 2 -s)] = R 4 dxdydtdsR Y 1 (x, y, t, s) R Y 2 (x 1 -x, x 2 -y, y 1 -t, y 2 -s) R Y 1 R Y 2 (x 1 , x 2 , y 1 , y 2 ) (52)
The existence of the convolution R Y 1 R Y 2 of autocorrelation functions is thus the main criteria for providing a sense to the existence of Y 1 Y 2 . Since Y 1 and Y 2 are stationary (see Section III-B1), we can focus on the two-dimensional versions of their autocorrelations. The problem then becomes the existence of:

R Y 1 R Y 2 (x, y) = R 2 R Y 1 (t, s) R Y 2 (x-t, s-y)dtds (53)
where R Y k is given by Eq. (50). To solve this problem, we thus impose, in addition to annihilating frequencies ±(u k , v k ), that φ is well localized in such a way that γ Y dened by Eq. (51) belongs to L 1 (R 2 ) ∩ L 2 (R 2 ). In this respect, both R Y 1 and R Y 2 pertain to L 1 (R 2 )∩L 2 (R 2 ) and we can conclude that Y 1 Y 2 is a second order stochastic eld that converges in quadratic mean.

Moreover, we have from Eq. ( 53):

γ Y 1 Y 2 (ω 1 , ω 2 ) = γ Y 1 (ω 1 , ω 2 ) × γ Y 2 (ω 1 , ω 2 ) = γ H 1 ,(u 1 ,v 1 ) (ω 1 , ω 2 ) × γ H 2 ,(u 2 ,v 2 ) (ω 1 , ω 2 ) × φ(ω 1 ) 4 × φ(ω 2 ) 4 (54)
The quadratic mean sense of

G H 1 ,(u 1 ,v 1 ) G H 2 ,(u 2 ,v 2 )
follows from the fact that we can provide a sequence

φ [n] t n
(of wavelets for instance) where every φ

[n] t ∈ L 1 (R)∩L 2 (R)
which converges to a Dirac distribution: the result follows as a limit case and makes PSD association to G H 1 ,(u 1 ,v 1 ) G H 2 ,(u 2 ,v 2 ) as (set φ = 1 in Eq. ( 54)):

γ G H 1 ,(u 1 ,v 1 ) G H 2 ,(u 2 ,v 2 ) (ω 1 , ω 2 ) = σ 2 1 2 1+2H 1 sin(πH 1 )Γ 2 (1 + H 1 ) ((ω 1 -u 1 ) 2 + (ω 2 -v 1 ) 2 ) H 1 +1 × σ 2 2 2 1+2H 2 sin(πH 2 )Γ 2 (1 + H 2 ) ((ω 1 -u 2 ) 2 + (ω 2 -v 2 ) 2 ) H 2 +1 (55)
Moreover, by iterating convolutions, we derive the PSD associated with a GFBF:

γ E H Q (ω 1 , ω 2 ) = Q k=1 σ 2 k 2 1+2H k sin(πH k )Γ 2 (1 + H k ) Q k=1 (ω 1 -u k ) 2 + (ω 2 -v k ) 2 H k +1 (56)
PSD which highlights Q spectral poles when all couples (u k , v k ) for k = 1, 2, . . . Q are dierent.

IV. Multi-fractional PSD pole identification

Fractionally integrated elds presented in Sections II and III share one common property: they have PSD admitting many frequency points associated with innite variances called poles. In practical applications requiring multi-fractional model selection, determining the number of interacting elds requires counting the number of frequency points with very high energy (exponential growth near the pole). Several solutions can be investigated for this problem. The easiest one involves selecting local maxima of the PSD and aecting them to poles. This has many drawbacks in terms of under-and over-determination (not any local maximum is eectively a pole and some poles can be squeezed when they are located near a pole with very high exponential growth).

We propose hereafter, after a long run search, a solution based on Convolutional Neural Network (CNN). The experimental setup concerns associating an observed multi-fractional GFBF 1 eld to a category indexed by Q ∈ {1, 2, . . .} corresponding to the number Q of poles used for sample generations.

The database generated for experimental tests contains 1200 images per specied value of Q: thus a total of 4800 GFBF images when the number of poles pertains to the category labels {1, 2, 3, 4}. For any category, poles and Hurst parameters are generated randomly, following the gamma distribution so as to impose a signicant energy concentration in ]0, π/2[×]0, π/2[ for poles and ]0, 1/2[ for Hurst exponents. An overview of this database 2 , denoted D hereafter, is given by Figure 5.

A. Direct learning of multi-fractional texture features

Because multi-fractional interactions are obtained by convolution operators (ltering white noise so as to impose a desired structure), deep CNN architectures seems 1 GFBF are chosen because their synthesis on the basis of fast Fourier transform requires less computational cost than the deluge of Γ computations for synthesis of K-DI-G elds. 2 The complete database is available for download at: http://am. atto.free.fr/AlbumTexturesGFBF.htm.

G H 1 -→ = E H 1 E H 1-• G H 1 G H 2 -→ = E H 2 E H 2-• E H 2 G H 3 -→ = E H 3 E H 3-• E H 3 G H 4 -→ = E H 4 E H 4-•
Fig. 5. Sample elements of GFBF database D, where every E H Q is associated with a number Q of distinct poles.

adapted per nature 3 to their discrimination. In this respect, we have rst investigated several existing architectures that are trained on millions of images to discriminate textures in D after retraining only the last convolutional layer. But none of them outruns more than 16% of good classication, the architecture achieving these 16% being \AlexNet" CNN [START_REF] Krizhevsky | classication with deep convolutional neural networks[END_REF] (trained on ImageNet [START_REF] Deng | A large-scale hierarchical image database[END_REF], a database with 1000 categories and 1.2 million training images) which is known to yield more than 90% accuracy in standard image classication tasks! These poor results are due to that multi-fractional textures are not objects with deterministic geometries: for the latter, feature recognition implies identifying similar structures up to elementary transformations whereas multifractional textures can present similar structures (Figure 5 shows almost more visual inter-category similarities than intra-category similarities), while having a dierent number of interactions: the issue is not thus a standard pattern recognition, but counting the number of visible pseudo-patterns.

We have thus been forced to retrain full CNNs by varying dierent features of the network. The most relevant architecture is still very close to [START_REF] Krizhevsky | classication with deep convolutional neural networks[END_REF], but with a signicant improvement of convolution lter lengths to ensure the results given in Table I. These results are unsatisfactory and the conclusion has been seeking multifractional feature interactions from an indirect learning (see next section).

B. Indirect learning from PSD as inputs

As highlighted in the previous section, we have to force learning in not focusing on a recognized fractional feature since the same feature can appear in textures pertaining to dierent categories. For this purpose, we propose to perform learning of multi-fractional interactions from PSD features. Table II presents the best relevant results derived from the CNN described in Appendix C.

It follows from Table II that counting the number of multi-fractional interactions can be performed with good performance from CNN based approach associated with spectral inputs. A study of the rst convolutional layer of the CNN shows specializations in 3 main convolution categories (see Figure 6):

• very selective lters in frequency and fast decay near the frequency selected (40% of rst level convolution lters), • selective lters in compact frequency blocks (rounded and wide support near the central frequency, 30% of rst level convolution lters), • multiple frequency selection lters (30% of rst level convolution lters). 3 Because they involve many layers of convolution lters and several non-linear transfer functions Note that from Table II, increasing the number of interactions leads to a signicant loss of performance. The reason is linked to the spectral content estimated from GFBF textures: 4 interactions means that 4 frequency points are associated with innite variance (energy for zero-mean eld) in ]0, π/2[×]0, π/2[. In addition, because of the exponential decays near these singular points, certain poles dominate others (the latter become almost nonvisible in spectra as illustrated in Figure 7). Moreover, Figure 7 highlights that when the number of interactions is large, then:

• only few poles (bright values dominating their surrounding neighborhoods) are `visible' in the PSD, • the overlays of poles yield bright intersections that can be confused with poles.

V. conclusion

The paper has shown that fractional order stochastic integration is a powerful tool for modeling with few parameters, image edge and texture information. Two categories of fractional eld models corresponding to discrete and continuous domain integrals have been studied.

For discrete domain modeling, the main contribution proposed by the paper is the derivation of a category of elds where any eld is associated with an arbitrary number of poles having their location controlled by a double indexed Gegenbauer polynomial.

For continuous domain modeling, the main contribution proposed by the paper is the proposal of a mathematical framework that provides a concise sense to convolutions of modulated fractional Brownian elds, any term involved in the convolution having the role of creating a spectral pole.

Because both discrete and continuous domain models can admit a pre-specied number of spectral poles, they are multi-fractional (a Holder exponent is associated with every pole) and this property raises the issue of counting the number of poles, given an image observation and for model selection purpose.

The paper has then proposed a third contribution consisting in providing a deep neural network architecture involving spectral information inputs and large size convolution lters at rst network layers so as to make pole identication possible. This network is shown ecient for discriminating up to 4 neighbor poles, which is a good performance in a tricky problem: every pole is a spectral point associated with an innite theoretical variance and its presence aects signicantly the observability of other poles located in a close neighborhood.

Appendix A Spatial Gegenbauer Fields

The problem expressed by Eq. ( 28) is expanding c -2µT -2νS + T 2 + S 2 -λ as a discrete series involv- Filters with multiple selection of frequencies Fig. 6. First convolutional layer of the CNN network described in C when learning has been addressed from GFBF PSD features. ing integer powers of T and S. This can be done by noting that the Gamma function given by Eq. ( 13) has equivalent expression (change of variable u → c -2µT -2νS + T 2 + S 2 x):

E H 1 E H 2 E H 3 E H 4 E H 5 E H 6 E H 7 E H 8 E H 9 E H 10 E H 1 E H 2 E H 3 E H 4 E H 5 E H 6 E H 7 E H 8 E H 9 E H 10 γ E H 1 γ E H 2 γ E H 3 γ E H 4 γ E H 5 γ E H 6 γ E H 7 γ E H 8 γ E H 9 γ E H 10
Γ (λ) = c -2µT -2νS + T 2 + S 2 λ × +∞ 0 e -(c-2µT -2νS+T 2 +S 2 )x x λ-1 dx (57)
so that we can write:

c -2µT -2νS + T 2 + S 2 -λ = 1 Γ (λ) +∞ 0 e -cx x λ-1 e (2µT-T 2 )x e (2νS-S 2 )x dx
The latter can be written in terms of the Hermite-Kamp e de F eriet polynomials [START_REF] Appell | \Fonctions hyperg eomtriques et hypersph eriques polynômes d'hermite[END_REF]:

c -2µT -2νS + T 2 + S 2 -λ = 1 Γ (λ) ∞ m=0 ∞ n=0 1 m! n! T m S n × +∞ 0 e -cx x λ-1 H m (2µx, -x)H n (2νx, -x)dx
where

H k (x, y) = k! k/2 j=0 x k-2j y j (k -2j)! j! (58) Thus c -2µT -2νS + T 2 + S 2 -λ = ∞ m=0 ∞ n=0 C λ m,n (µ, ν)T m S n (59)
with

C λ m,n (µ, ν) = 1 m! n! 1 c λ+m+n 1 Γ (λ) × +∞ 0 e -x x λ-1 H m (2µx, -cx)H n (2νx, -cx)dx (60)
Now, by using Eq. ( 58), we derive that:

+∞ 0 e -x x λ-1 H m (2µx, -cx)H n (2νx, -cx)dx = m! n! m/2 p=0 n/2 q=0 (2µ) m-2p (2ν) n-2q (-c) p+q (m -2p)! (n -2q)! p! q! × +∞ 0 e -x
x λ+m-p+n-q-1 dx which, by noting that

+∞ 0 e -x x λ+m-p+n-q-1 dx = Γ (λ + m + n -p -q)
leads us to conclude that Eq. ( 60) is equivalent to:

C λ m,n (µ, ν) = 1 c λ+m+n 1 Γ (λ) × m/2 p=0 n/2 q=0 Γ (λ + m + n -p -q)× (-1) p+q c p+q (2µ) m-2p (2ν) n-2q p! q! (m -2p)! (n -2q)! (61)
In particular, when b=1, we note C 1 m,n C m,n , which reduces to:

C m,n (µ, ν)= 1 c 1+m+n × m/2 p=0 n/2 q=0 (m + n -p -q)! × (-1) p+q c p+q (2µ) m-2p (2ν) n-2q p! q! (m -2p)! (n -2q)! (62)
since Γ (1) = 1 and Γ (k + 1) = k! when k is a non-negative integer.

Appendix B On convolutions of modulated fractional

Brownian motions

A. Context Consider a stochastic zero-mean process:

X : (t, λ) ∈ R×Ω or Z×Ω -→ X(t, ω) = X(t, λ) = X t (λ) = X λ (t) X(t) ∈ R.
Assuming that X is a second-order stochastic means, for all t:

X(t) 2 L 2 (Ω) = E |X(t)| 2 = Ω |X(t)| 2 dP(λ) < ∞ ( 63
)
where Ω = R in the following and P is a probability measure on elements of B(R). Under property given by Eq. (63), functions [mean]: t → EX(t) and [autocorrelation] (t, s) → E[X(t) X(s)] = R(t, s) are well-dened and real valued.

For a second-order stochastic process, continuity (that is X(t) → X(s) as t → s) is related to the so-called quadratic mean sense (use of the norm associated to L 2 (Ω)) and consists in:

lim t→s X(t) -X(s) 2 L 2 (Ω) = lim t→s E|X(t) -X(s)| 2 = 0
A consequence of using the quadratic norm is the following: X is continuous (in quadratic mean) if functions t → EX(t) and (t, s) → E[X(t) X(s)] = R(t, s) are continuous (standard function continuity). We assume that X is with zero-mean in what follows.

A second-order zero-mean stochastic Gaussian process is completely specied by its autocorrelation function. For the Fractional Brownian Motion (FBM) denoted X H with Hurst exponent/parameter H, 0 < H < 1, this autocorrelation function is given by:

R H (t, s) = σ 2 2 |t| 2H + |s| 2H -|t -s| 2H (64)
The FBM Modulation (FBMM) by a deterministic exponential wave, denoted X H,u 0 (x) = e iu 0 x X H (x) is a Gaussian process with zero-mean (as a consequence of zero-mean assumption on X H ) and autocorrelation:

R X H,u 0 (t, s) = E[X H,u 0 (t) X H,u 0 (s)] = e iu 0 (t-s) R H (t, s) (65)
where R H is given by Eq. (64).

Consider now two independent FBMM X H 1 ,u 1 and X H 2 ,u 2 . We want to provide a quadratic mean sense to the following convolution operator between FBMM:

X H 1 ,u 1 X H 2 ,u 2 (t) ? = R X H 1 ,u 1 (x)X H 2 ,u 2 (t -x)dµ(x) (66)
In practice, if we assume that µ is the Lebesgue's measure, it would have been suced that R X H 1 ,u 1 and R X H 2 ,u 2 pertain to L 2 (R 2 ). However, this is not true for FBMM as it can be observed by taking Eq. (64) and Eq. (65) into account. Nevertheless, we can go beyond these functional limitations by selecting an appropriate parametric measure that transforms

X H 1 ,u 1 and X H 2 ,u 2 into Y 1 and Y 2 : Y k (t) R X H k ,u k (x) dµ t (x) (67)
for k = 1, 2 and with convergence of the latter integrals in quadratic mean sense. Note that for the convergence of Eq. (67) to hold true, it suces to provide a kernel φ satisfying

• φ ∈ L 1 (R) ∪ L 2 (R), • (1 + x 2 )φ(x) ∈ L 1 (R) ∪ L 2 (R),
and let dµ t (x) = φ t (x)dx where φ t (x) = φ(x -λt). For such a measure,

R 2 e iu 0 (x-y) |x| 2H k φ t (x)φ s (y)dxdy and R 2 e iu 0 (x-y) |x -y| 2H k φ t (x)φ s (y)dxdy
exists and are nite, thus,

Y k (t) 2 L 2 (Ω) R X H k ,u k (x) dµ t (x) 2 L 2 (Ω) = R 2 R X H k ,u k (x, y)φ t (x)φ s (y)dxdy < ∞
in case where R X H k ,u k has the form given by Eq. ( 65). The quadratic mean sense of X H k ,u k follows from this last property as a limit case, given that we can provide sequence φ [n] t n such as wavelets where every φ

[n] t ∈ L 1 (R) ∩ L 2 (R)
: φ(u k ) = φ(-u k ) = 0
where we have used the following denition of the Fourier transform:

φ(ω) = R e -iωx φ(x)dx (68) Consider the autocorrelation function R Y k of Y k : R Y k (t, s) = R 2 R X H k ,u k (x, y)φ t (x)φ s (y)dxdy = σ 2 2 R 2 |x| 2H k + |y| 2H k -|x -y| 2H k
e iu k (x-y) φ t (x)φ s (y)dxdy (69) Under frequency ±u k blocking by φ, we obtain:

R Y k (t, s) = - σ 2 2 R 2 e iu k (x-y) |x -y| 2H k φ t (x)φ s (y)dxdy
This integral can be re-written, from a change of variable:

R Y k (t, s) = - σ 2 2 R 2 e iu k z |z| 2H k φ t (z + y)φ s (y)dzdy
By decomposing |z| 2H k with respect to [START_REF] Bahr | \Inequalities for the rth moment of a sum of random variables, 1 r 2[END_REF], we deduce:

R Y k (t, s) = - σ 2 2 Γ (2H + 1) sin(πH) π × R 2 e iu k z R 1 -cos(zξ) |ξ| 2H k +1 dξ φ t (z + y)φ s (
y)dzdy (70) and by noting:

γ H (ξ) = σ 2 Γ (2H + 1) sin(πH) |ξ| 2H+1 (71) we derive: R Y k (t, s) = - 1 2π R γ(ξ)dξ× R 2 e iu k z (1 -cos(zξ)) φ t (z + y)φ s (y)dzdy (72)
By taking again into account, the kernel φ frequency ±u k blocking property, we have:

R 2 e iu k z φ t (z+y)φ s (y)dzdy = e iu k (t-s) φ(-u k ) φ(u k ) = 0
and thus:

R Y k (t, s) = 1 2π × R dξγ(ξ) R 2 e iu k z cos(zξ)φ t (z + y)φ s (y)dzdy (73) which is equivalent to R Y k (t, s) = 1 2π × R dξγ(ξ) R 2 e iu k z e izξ + e -izξ 2 φ t (z + y)φ s (y)dzdy
The latter corresponds, after some Fourier based calculus:

R Y k (t, s) = 1 2π R dξγ(ξ) 1 2 × φ t (-ξ -u k ) φ s (-ξ -u k ) + φ t (ξ -u k ) φ s (ξ -u k ) thus, R Y k (t, s) = 1 2π R γ(ξ + u k ) φ t (ξ) φ s (ξ)dξ
Now, since we have assumed that t in notation φ t denotes a translation parameter, then: φ t (ξ) = e -ictξ φ(ξ) and φ s (ξ) = e -icsξ φ(ξ) so that:

R Y k (t, s) = 1 2π R γ(ξ + u k )e -ic
γ Y k (ξ) = γ(ξ -u k ) φ(ξ) 2
Assuming that ideal φ behaves approximately as a Dirac distribution (limit case as discussed at the end of the previous section), we can associate the following PSD to X k :

γ X k (ξ) = γ(ξ -u k ) = σ 2 Γ (2H + 1) sin(πH) |ξ -u k | 2H+1 C. Convolution between FBMM X H 1 ,u 1 , X H 2 ,u 2 via characterization of Y 1 Y 2
The sole condition required for the convolution

R Y 1 (x)Y 2 (t -x)dx Y 1 Y 2 (t)
4 For a zero-mean second order stochastic process X, stationarity (second-order or Wide Sense Stationarity, WSS) holds true if: R X (t, s) = R X (t -s, 0) ≡ R X (t -s = x) to make sense as the denition of a second order stochastic process Y 1 Y 2 is the existence of:

E[Y 1 (x)Y 2 (t -x)Y 1 (y)Y 2 (s -y)]dxdy = E[Y 1 (x)Y 1 (y)]E[Y 2 (t -x)Y 2 (s -y)]dxdy = R Y 1 (x, y)R Y 2 (t -x, s -y)dxdy = R Y 1 R Y 2 (t, s) (74) If R Y 1 , R Y 2 ∈ L 2 (R 2 )
, then this condition is satised. The latter is not very restrictive since φ can be chosen welllocalized in time frequency.

Furthermore, under the stationarity induced by φmeasure on Y 1 , Y 2 , then:

• from the change of variable t ← t -s and x ← x -y, we can write:

R Y 1 Y 2 (t) = R Y 1 (x)R Y 2 (t -x)dx
• and, provided that R Y 1 Y 2 obtained just above admits a Fourier transform, we can deduce the PSD

γ Y 1 Y 2 of Y 1 Y 2 as: γ Y 1 Y 2 (ξ) = R Y 1 R Y 2 (t)e -iξt dt = γ Y 1 (ξ)γ Y 2 (ξ) = γ(ξ -u 1 )γ(ξ -u 2 ) φ(ξ) 4 (75)
Finally, the PSD associated 5 to X H 1 ,u 1 X H 2 ,u 2 via kernel φ is:

γ X H 1 ,u 1 X H 2 ,u 2 (ξ) = γ(ξ -u 1 )γ(ξ -u 2 ) = γ X H 1 ,u 1 (ξ) × γ X H 2 ,u 2 (ξ) = σ 2 1 σ 2 2 Γ (2H 1 + 1)Γ (2H 2 + 1) sin(πH 1 ) sin(πH 2 ) |ξ -u 1 | 2H 1 +1 |ξ -u 2 | 2H 2 +1
Appendix C Convolutional Neural Network (CNN) architecture for multi-fractional pole identification

The CNN proposed for multi-fractional interaction count is presented in Table III. This CNN is characterized by a rst convolution layer with long-size impulse response lters so as to possibly allows for a concise frequency selection and surprisingly, up to 40% of the lters have been specialized (starting from a white noise initialization) in frequency selection (see Figure 6).

In the CNN of Table III, sizes of convolution lters have been decreased progressively, from layer to layer, in order 5 Deduction must be understood as a limit case involving an L 1 (R) ∩ L 2 (R) sequence φ [n] t n that converges to a Dirac distribution, thus φ [n] t n tends to constant function 1.

to limit the overall CNN complexity. From Table III, we have the following correspondences:

• Rectied Linear Unit (ReLU), function

x → (x) + =

x if x > 0 0 if x 0

• Normalization (cross-channel): where N is the number of kernels used in the layer and n is the normalization neighborhood size (same spatial position, n adjacent channels); • Softmax:

X = (X i ) n i=1 → e X i / n i=1 e X i n i=1
this loss function is used for the prediction of a single category between Q mutually exclusive Q-fractional categories.

1

  LISTIC, EA 3703, Universit e Savoie Mont Blanc -Universit e Grenoble Alpes, France Email: Abdourrahmane.Atto@univ-smb.fr -Phone: +334 50 09 65 27 -Fax: +334 50 09 65 59 ||||||||{ The work was supported by PHOENIX ANR-15-CE23-0012 grant of the French National Agency of Research.

( 20 )

 20 and (21), respectively):

Fig. 3 .

 3 Fig. 3. Examples of MFBF samples. Row 1: isotropic FBF F H associated with dierent Hurst parameters H.. Row 2: phase terms of G H,(u 0 ,v 0 ) .

  Selective filters in frequency (fast decay)Selective filters in frequency (slow decay)

Fig. 7 .

 7 Fig. 7. Curse of power for GFBF E H Q associated with PSD γ E H Q : PSD is displayed as log(1+log(1+WP Spectrum)) to allow observing some poles. But when the number of interactions is high, only the few dominant poles are visible in the PSD! [Sizes] / Input GFBF: 4092×4092 ; GFBF PSD: 1024×1024 ; Supports displayed: 64 × 64 -→ [0, π/16] × [0, π/16].

x m → y m = x

  . PSD association to FBMM X H k ,u k via that of Y k In the following, we assume that the following property holds true for k=1, 2:[Blocking frequency ±u k ] function φ introduced in section above satises at the specic real values ±u k

and parameterized by an upper index n which makes them converge to Dirac distributions. B

  (t-s)ξ φ(ξ) ∞ for some η > 0. In this case, Y k is a stationary4 stochastic process with autocorrelation in the reduced onedimensional form:

			2	dξ
	and nally:
	R Y k (t) =	1 2π R	γ(ξ -u k ) φ(ξ)

R Y k (t, s) = 1 2π R γ(ξ -u k )e ic(t-s)ξ φ(ξ) 2 which converges provided that sup |ξ-u k | η φ(ξ)/ξ < 2 e ictξ dξ

and associated with PSD (derived from inverse Fourier identication in the latter form):

Test confusion matrix from a direct learning: average retrieval per texture category and average inter-category confusion (in %). The D 3 database is composed by 3 GFBF categories (1800 textures for training, 1800 for learning). The D 4 database is composed by 4 GFBF categories (2400 textures for training, 2400 for learning). The CNN used is described in Appendix C.


Aurélie Martin
email: aurelie.remy@inrs.fr

Pascal Wild

Bivariate left-censored measurements in biomonitoring: a Bayesian model for the determination of Biological Limit Values based on Occupational Exposure Limits

Keywords: Biomonitoring, Biological Limit Values, Limit of Detection, Occupational Exposure Limit, Bayesian model, left-censored data

Biological Limit Values (BLV) are often determined from the occupational exposure limits (OEL) in modelling biological data obtained on a number of exposed subjects based on measurements of air exposure. In order to obtain such BLVs, biomonitoring studies are conducted collecting simultaneously biological and airborne measurements to these substances in exposed workers.

One obstacle in the modelling of such data is the often large number of values below the Limit of Detection (LOD) for both biological and airborne measurements (left censored measurements). A second difficulty, which is also a strength, is that multiple measurements are obtained for the same workers, leading to non-independence of the data.

In this paper, we propose a statistical method based on Bayesian theory making use of measurements below the LOD for both dependent (biological) and independent (air exposure) data, and taking into account multiple measurements on the same worker.

This method relies on the modelling of the airborne exposure measurements using standard random effect models adapted for values below LOD and the simultaneous modelling of the biological measurements assumed to be linearly (on the log scale) related to the airborne exposure while accounting for between worker variability. This method is validated by a simulation study in which up to 50% of the measurements are censored for both variables in realistic settings. This simulation study shows that the proposed method is uniformly more efficient than the candidate alternative we considered (MLE method) that did not make use of a data with airborne measurements below the LOD.

When the method is applied on a real biomonitoring data set among electroplating workers exposed to chromium with 54% censored airborne measurements and 20% censored urinary measurements, the slope is steeper when incorporating these data using the proposed Bayesian method leading to different BLV estimations depending on the OEL used.

INTRODUCTION

According to the CDC, biomonitoring is defined as a method for assessing human exposure to chemicals by measuring the chemicals and/or their metabolites in human tissues or specimens, such as blood or urine. In occupational settings, its ultimate objective is to be able to propose suitable measures (improvement of technical, organisational and personal prevention) in order to reduce the exposure and its potential impact on health. For many hazardous substances, in particular when entry paths are multiple, the individual exposure sustained can only be quantified, and therefore assessed, by means of biomonitoring. One means of ensuring that workers are adequately protected is by setting Biological Limit Values (BLVs) which are not to be exceeded. Ideally, as for instance for Cadmium and Lead, BLVs are based on the relation between the internal exposure and its health effects observed in human studies. This is however only rarely the case.

For substances for which the main entry is inhalation, in particular when the exposure is to fumes or aerosols, the Biological Limit Values (BLV) can be derived from Occupational Exposure Limit values (OELs). The French Agency for Food, Environmental and Occupational Health & Safety (ANSES) states for instance that (our translation) "in absence of such (health-based) data, for substances with a no-effect level, the BLV will be computed from the expected biological concentration, when the worker is exposed at the 8 hour OEL".

Estimating such an expected biological concentration relies on the statistical modelling on the biological concentration among exposed workers as a function of simultaneously measured airborne exposure concentrations [START_REF] Lison | BIOLOGICAL MONITORING OF WORKERS EXPOSED TO COBALT METAL, SALT, OXIDES[END_REF] [START_REF] Pierre | EFFECT OF DIFFERENT EXPOSURE COMPOUNDS ON URINARY KINETICS OF ALUMINUM AND FLUORIDE IN INDUSTRIALLY EXPOSED WORKERS[END_REF] [START_REF] Pierre | Biomonitoring of two types of chromium exposure in an electroplating shop[END_REF]. Such studies typically include a limited number of workers for which repeated measurements are obtained. A relationship is usually documented [START_REF] Acgih | ACGIH 2015 . TLV's and BEI's Based on the Documentation of the Threshold Limit Values for Chemical Substances and Physical Agents and Biological Exposure Indices[END_REF] by establishing a linear relationship between log transformed atmospheric exposure and biological exposure using simple linear regression with air concentrations as independent variable and biological concentrations as dependent variable. Logtransformation is used because the statistical distributions of both air and biological data are highly skewed. The logarithmic transformation has the consequence to make the distribution more symmetric.

There are a series of issues with such linear regression of log-transformed data, which need to be accounted for by specific statistical methods.

First, a majority of these data are repeated measurements on the same workers, leading to the non-independence of the data. Indeed, the elimination rates in different workers may differ significantly (between-worker variability in biological measurements) even with the same airborne exposure. On the other hand, as has been shown some time ago (Kromhout, worker variability in excess of the day-to-day variability.

Second, although analytical methods are increasingly sensitive, the improvement of worker protection and collective protection leads to lower concentration levels and increasing numbers of exposure measurements are below the Limit Of Detection (LOD) that in the statistical jargon are referred to as left-censored data. Note that sometimes exposure measurements are reported as below the Limit of Quantification (LOQ) rather than LOD.

Up to now, in the industrial hygiene context, methods of data censored treatment have mainly focused on the univariate analysis of the airborne exposure measurements [START_REF] Helsel | Nondetects And Data Analysis: Statistics for censored environmental data New York[END_REF]) [START_REF] Hewett | A comparison of several methods for analyzing censored data[END_REF] [START_REF] Ganser | An Accurate Substitution Method for Analyzing Censored Data[END_REF] [START_REF] Huynh | Comparison of Methods for Analyzing Left-Censored Occupational Exposure Data[END_REF]) [START_REF] Huynh | A Comparison of the beta-Substitution Method and a Bayesian Method for Analyzing Left-Censored Data[END_REF]. In this context, [START_REF] Jin | Analysis of Lognormally Distributed Exposure Data with Repeated Measures and Values below the Limit of Detection Using SAS[END_REF]) performed a simulation study comparing different approaches to take account of up to 80% left-censored measurements while at the same time acknowledging a between-worker variability. Their study compared the widely used substitution method, replacing data below the LOD for the independent variable by LOD/2 or LOD/sqrt(2), with mixed-effects models (based on Maximum Likelihood Estimation, MLE) to log-normally distributed occupational exposure data. They found that the substitution methods led to large biases (underestimation or overestimation) of geometric means and standard deviations, especially when the proportion of censored measurements was large, but that this was not the case with the MLE methods. However, we know of no paper in the context of industrial hygiene that takes into account of censoring for both biological and airborne concentrations. In the context of modelling the relation between biological concentrations as a function of airborne concentrations, censored data for the independent variable, i.e. the airborne concentrations, can only be ignored, leading to a loss of information, loss of statistical power, and loss of precision of the relation assessment or substituted by LOD/2 or LOD/sqrt(2) leading to bias (see for instance [START_REF] Ong | Biomarkers of exposure to low concentrations of benzene: a field assessment[END_REF] or [START_REF] Ghittori | Evaluation of occupational exposure to benzene by urinalysis[END_REF]).

The purpose of the present paper is to present a statistical method making use of leftcensored measurements for both dependent (biological) and independent (airborne) concentrations while accounting for between-worker variability both for airborne and biological measurements in order to improve the precision of their relation. Such an improved method would yield a BLV that better corresponds to the OEL.

METHODS

In this paper, we propose a statistical method based on Bayesian theory to use bivariate leftcensored data, that is data with measurements below the LOD for both dependent (biological) and independent (air exposure) data. This is done by modelling simultaneously the two variables. This proposed method is validated in a simulation study contrasting its performance to that of standard methods (based on MLE) that ignore all data corresponding to atmospheric data below the LOD.

Finally the proposed Bayesian method is applied to a real data set of chromium measurements among electrolytic platers.

Description of the data structure

We assume that we measured repeatedly the occupational airborne exposure as well as the biological exposure (e.g. in urine) for a series of workers. We assumed that these workers belong to the same Homogeneous Exposure Group i.e. a group of workers with one geometric mean airborne exposurethe difference between workers being characterized by a between-worker geometric standard deviation.

Both types of measurements are assumed to follow log-normal distributions and for both, there are within and between variance components. While the between worker variance component for air exposure is mainly due to differences in work practices, the betweenworker variance component for biologic exposure is due to intrinsic between-worker differences in metabolism.

Finally, both measurements are measured with a certain precision and have therefore LODs and the biologic exposure measurement is assumed to depend, for any given worker, only on its atmospheric exposure.

The log-transformed airborne measurements are denoted by X and the log-transformed biological measurements are denoted by Y. The corresponding respective limits of detection (on the log scale) are denoted by and .

In a first data pattern, we simulated 20 subjects each with 11 measurements for both airborne and biological measurements (balanced data). In a second data pattern, we considered 80 subjects, with 1 measurement for 50 subjects, 2 measurements for 20 subjects and 5 measurements for 10 subjects (unbalanced data).

Bayesian modelling

Let us first present the full Bayesian model ignoring censoring Following [START_REF] Kromhout | A COMPREHENSIVE EVALUATION OF WITHIN-WORKER AND BETWEEN-WORKER COMPONENTS OF OCCUPATIONAL EXPOSURE TO CHEMICAL-AGENTS[END_REF], the airborne exposure data X are modelled by a random effect model with a random effect representing between-worker variability and a residual error representing within-worker variability. A similar model is used for the biological data Y, with random effects for between-worker variability and a residual within-worker variability. However we assume that, given these random effects, the log biological exposure measurement depends linearly on the (log-transformed) airborne exposure measurement.

For the jth observation of the ith worker, this can be written as:

with (1) with (2)
where and are respectively the between-worker standard deviation and the within-worker standard deviation for the airborne log exposure measurements and the population mean log airborne exposure; and are respectively the between-worker standard deviation and the within-worker standard deviation for log biological measurements, is the population log biological exposure for a zero log airborne exposure (intercept), and is the slope of the relation.

Censoring

In order to accommodate measurements below the LOD, we introduce following notation for a truncated normal distribution , a normal distribution truncated to the interval . Thus, in our model corresponding to (1), a left censored value observed for jth observation of the ith worker would be specified as ( ) . In our model corresponding to (2), a left censored value observed for jth observation of the ith worker would be specified as .

Prior distribution of the parameters

Weakly informative prior distributions of the variability parameters are based on inverse gamma distribution with parameters (0.001, 0.001). It was chosen as it is the conditionally conjugate distribution for the inverse-chi2 distribution and leads to less computation time.

However, in order to examine the sensitivity to these choices of priors, we also used noninformative Half-Cauchy distributions for these parameters instead of Gamma distributions in a limited set of simulations.

Estimation of model parameters

The Bayesian model was fitted using Gibbs Sampling, a special case of Markov chain Monte Carlo (MCMC) methods (see [START_REF] Wild | Confidence intervals for probabilities of exceeding threshold limits with censored log-normal data[END_REF] for Gibbs Sampling in the context of censored data, and [START_REF] Lunn | The BUGS Book -A practical Introduction to Bayesian Analysis[END_REF]) for more recent general discussion of MCMC methods) that by default take the censored data into account. We used the freely available RJags software (JAGS 4.2.0 and R version 3.3.2) we compared to other freely available softwares. In the supplementary material available online, Appendix 1 contains the code of the model used, respectively in Jags, OpenBugs and WinBugs.

We ran the Gibbs Sampling algorithm for N=50000 samples, a thinning of 5, and discarded the 10000 to account for "burn-in". Convergence was checked by running two different chains. The starting points of the parameters were based on their true value for the simulation study. When analysing the actual data, we ran three chains starting from overdispersed starting points.

Simulation Study

In order to assess the performance of the above proposed Bayesian model in comparison with previously published methods, we did a simulation study by repeatedly simulating airborne exposure measurement and biological exposure measurement in settings as realistic as possible. We thus fixed the parameters of the data generated by our simulations to values that were abstracted from the scientific literature or, when not available, from the analysis of our own biomonitoring data given rise to a series of scenarios detailed below Data generation scenarios Specifically, we defined the following main scenario by: exposure measurements and the biological exposure measurement on the original (non log-transformed) scale.



= 1.25 corresponding to a =3.48, that is the median of between-worker GSD of historic urinary data of our own biomonitoring laboratory.  For each combination of these parameters, we fix the theoretical probability of censoring in X and Y of 30 and 50%.

 = 0.
Fixing these parameters, determines (for given values of the LODs) the intercept and the mean log-exposure . Without any loss of generality we arbitrarily fixed the (logtransformed) LOD for atmospheric measurements to 0. The intercept of the regression can thus be interpreted as the predicted mean urinary measurement when the atmospheric measurement is at the LOD.

This main scenario was simulated for balanced and unbalanced data pattern.

Subsequently, eight scenarios were then explored only for the balanced data pattern, modifying from the main scenario one of the value of the 4 standard deviations at a time:

 is then set to 0.71 and 1.6 ( =2 and =4.95 ) , that is the minimum and the maximum of between-worker GSD of historic urinary data of our own biomonitoring laboratory. Each selected combination of parameters was simulated N=1000 times

Data analysis

For each simulated dataset, we applied -the full proposed model based on Bayesian approach as described above, using RJags (JAGS 4.2.0 and R version 3.3.2),

-the reference MLE model (mixed tobit model), using STATA version 14.0 (program xttobit), similar to the method described in [START_REF] Jin | Analysis of Lognormally Distributed Exposure Data with Repeated Measures and Values below the Limit of Detection Using SAS[END_REF]). However, this last method leads to the exclusion to all data below the LOD for X.

Retained variables of interest:

 : the intercept of the relation that is the predicted mean urinary measurement when the atmospheric measurement is set to the LOD  : the slope of the relation.

 The derived BLVs corresponding to 3 theoretical values of OEL, defined by the theoretical percentage of exceeding them (0.1%, 2.5% and 10%).

Evaluation metrics

In order to compare the two methods, we computed the method-specific biases and RMSEs (Root Mean Square Error) summarizing both bias and variance over all 1000 simulated data sets. For each scenario, for and as well as for the derived BLVs we computed the bias and the RMSE according to following formulas

̂ √ ̂

For each parameter, we also computed the ratio between the RMSEs :

Chromium exposure in electrolytic plating workers

The proposed Bayesian method was applied to a real dataset of electrolytic plating workers exposed to chromium and contrasted with the reference MLE method.

Between 2007 and 2012, the INRS (French Institute for research and Safety), conducted a biomonitoring study to evaluate urinary chromium of workers exposed to electrolytic plating.

The study was carried out on 47 male workers exposed to Cr VI during electrolytic plating process.

Workers were from 5 companies. Full shift air samples were collected, over one to 5 days, and urinary samples were collected at the end of each shift.

The measure of airborne exposure to soluble CrVI (inhalable fraction) (CrA) was obtained by personal sampling. The samples were taken by aspiration of the ambient air through a device consisting of a 37 mm diameter filter (Whatman QMA quartz fiber) placed in a closed cassette (Millipore) at a flow rate of 2 L / min provided by a controlled flow pump [START_REF] Inrs | Echantillonnage des aérosols. Echantillonnage individuel d'un aérosol par cassette fermée -Méthode METROPOL fiche H2[END_REF]. The LOQ of CrA ranged between 0.16 and 1.16 µg/m 3 depending on the sampled volume. Analysis of CrA was performed by optical emission spectroscopy analysis (ICP-OES) according to the method described in the French METROPOL database (INRS 2008).

Chromium urinary analyses were performed in a controlled atmosphere room. The method consisted of the graphite tube atomization followed by detection by an atomic absorption spectrometer with Zeeman correction of the type AA220 (Varian, Melbourne Australia) [START_REF] Dube | Determination of chromium in human urine by graphite furnace atomic absorption spectrometry with Zeeman-effect background correction[END_REF]). The method developed and validated in the laboratory has a LOQ ranging between 0.3 and 0.83 µg/ L.

Creatinine in urine samples was determined by the Jaffe method (JAFFE 1886) using a clinical biology apparatus of the type Daytona (Randox, Crumlin, UK).

The proposed Bayesian approach was first applied as described above with the variant specifying 5 different Homogeneous Exposure Groups (HEG), one for each company. Each HEG had a different geometric mean for the airborne exposure but we assumed that the 5

HEGs had the same between and within worker GSDs.

On the basis of the current French OEL (1 µg/m 3 ) of CrVI (INRS 2006), and on the current American TLV® (50 µg/m 3 ) of soluble CrVI (ACGIH 2015), the corresponding BLVs were computed from the regression parameters estimated with the two methods.

RESULTS

Simulation study

Figure 1 shows the distribution of the slope estimations in the simulated data sets of the main scenario with balanced and unbalanced data patterns. As expected, the two methods were unbiased. However, the dispersion of the estimates increases with X-censoring but is always lower using our proposed Bayesian method, whatever censoring and data pattern. Note that for the unbalanced data pattern the dispersion is greater and the difference in dispersion between the two methods more important. See also the corresponding scatter plots in Appendix 2 of the supplementary material available online.

Table 1 shows the results of the simulations. Overall these simulations show that the RMSE is always lower when using the proposed Bayesian approach as documented by the ratios of RMSEs which are always greater than 1. This is particularly apparent for the unbalanced data.

Considering the balanced data, the mean ratio over the scenarios for the intercept is only 1.04 thus the gain in using the proposed Bayesian procedure is moderate and does not vary much according to the censoring pattern nor according the scenario.

For the slope β, the gain in using the proposed method is much more obvious (mean ratio over all scenarios 1.27), and is particularly important when the censoring is 50% in X (mean 1.37 with 30% censoring in Y and 1.35 with 50% censoring in Y). This ratio of RMSEs does not depend to a great extent on , it decreases with and with and increases with The maximal ratio 1.66 was obtained for scenario 4 ( =1.6, =1.34, =1.45, = 3.49) and 50% censoring for both X and Y.

With respect to the BLVs, the gain is almost negligible when the OEL is exceeded in 10% of the airborne measurements (BLV1). The gain is more important for (2.5% excess of OEL) and is most important for BLV3 (0.1% excess of OEL). In the latter case, the mean ratio of RMSEs is most important for 50% censoring for both X and Y (1.10). Again, scenario 4 shows the highest gain in efficiency.

We didn't observe any noteworthy differences on the posteriors when using Half-Cauchy non-informative priors compared with non-informative gamma distribution (observed on the main scenariodata not shown). Moreover, influence of the informative priors was very slight even when the prior was very different from the true values.

Chromium exposure in electrolytic plating workers

The studied dataset comprised 166 pairs of urinary and airborne chromium measurements, CrU (in µg/g creatinine) and CrA (in µg/m 3 ) respectively. This dataset contained a large number of measurements below the respective LOQs.

Figure 2 provides a schematic view of the distribution of these measurements. At the bottom left, the figure shows that measurements below the LOQ for both CrU and CrA represent 20% of the dataset. At the upper left, we have data below the LOQ for CrA and above LOQ for CrU, which represent 34% of the dataset. As mentioned above, 54% of airborne measurements were below LOQ and would thus have been excluded using standard statistical methods, although for many of them, urinary chromium measurements were available. Note that we considered LOQs and not LODs as the data were only available in this format.

Table 2 shows the descriptive results of our data set. One can notice that the airborne exposure is quite high, with 69 measurements (42%) above the present French OEL=1 µg/m 3 but also 4 (2%) above the American TLV ® = 50 µg/m 3 , which was the French OEL at the time of the study.

The two statistical methods give similar estimates for the urinary GSDs and for within-worker airborne exposure GSD but not for between-worker GSD (Table 3). It must however be noted that the analysis of the airborne exposure data were not part of the reference MLE model but were modelled independently and were therefore not influenced by the urinary data.

DISCUSSION

In this paper, we proposed a statistical method based on Bayesian theory making use of measurements below the LOD for both dependent (biological) and independent (airborne exposure) data, in the context of repeated measurements per worker. The simulation of realistic settings showed that the proposed method is uniformly more efficient than the candidate alternative we considered (reference MLE method) that did not make use of a data with airborne measurements below the LOD. This efficiency gain was most important for unbalanced data, it increased with the percentage of airborne measurements below the LOD.

With 50% censoring on airborne measurements, the RMSE of the estimated BLV decreased by up to 31% using the proposed Bayesian method. When the method was applied on a real biomonitoring data set with 54% airborne measurements below the LOQ, the slope was steeper and had a smaller standard error when incorporating these data using the proposed Bayesian method leading to different BLV estimations depending on the OEL used.

For the slope parameter, the simulation study shows a decrease of the RMSE with the proposed Bayesian method compared to reference MLE method, whatever the level of the variability parameters. However, the results are less clear-cut for the intercept of the relation;

although RMSE is always lower with proposed Bayesian method, the gap with reference MLE method is less than for the slope estimation. This might be because, with 30 or 50% censoring, the LOD is close to the middle of the airborne measurements and the intercept is thus well known in any estimation method.

Simulations were also run with lower percentage of censoring (10% -data not shown), on X and Y. In this case, there is still a gain using the proposed Bayesian method, but is not large enough to justify the added complexity of its use.

For derived BLV, the gain is lower, especially for BLV1 and BLV2, corresponding to OEL defined by the theoretical percentages 10% and 2.5% of exceeding them, respectively. The lower gain could be explained by the fact that, as for the intercept, these values of OEL correspond to observed data which lead to a better stability of response, whatever the method used. BLV3 corresponds to 0.1% excess and thus involves some extrapolation that is more sensitive to slope estimation.

We checked (data not shown) that when discarding all censored data for airborne measurements, Bayesian model and MLE method gave very similar results. We can therefore safely conclude that the observed gain in RMSE is due to the inclusion of these censored data. On the other hand, when discarding all censored data (airborne and biological), we observed a rather important bias in all estimated parameters.

Our proposed Bayesian method is a rather complex model, integrating random effects in order to quantify the variability both between and within subject. This type of model has been proposed already 20 years ago [START_REF] Kromhout | A COMPREHENSIVE EVALUATION OF WITHIN-WORKER AND BETWEEN-WORKER COMPONENTS OF OCCUPATIONAL EXPOSURE TO CHEMICAL-AGENTS[END_REF] but we are not aware that the proposed extension to biomonitoring data has ever been used even without consideration of data censoring. Indeed, in actual biological data, while the between-subject variability is often lower than the within-subject variability, we think that this is due to differences in day-to-day within-subject exposure. For any given airborne exposure level, we think that the between-subject variability in data is greater than the within-subject variability if inhalation is the predominant route of entry into the organism. When analysing historical data for metals (mostly unpublished), collected by our laboratory in a mixed model including the airborne exposure, the variance component for subject was indeed greater than the withinworker variance component adjusted on airborne exposure. The reason might be that the subject-specific physiological parameters like height, weight, BMI and metabolism vary widely between subjects, whereas within a given subject, the urinary excretion, adjusted on creatinine, depends mostly on the airborne exposure. Thus this between-subject variability cannot be ignored when analysing biomonitoring data, involving repeated measurements by subject. On the opposite, for the airborne exposure model, within-subject variability, due to day-to-day differences is usually greater (e.g. 77% of the 165 situations considered in [START_REF] Kromhout | A COMPREHENSIVE EVALUATION OF WITHIN-WORKER AND BETWEEN-WORKER COMPONENTS OF OCCUPATIONAL EXPOSURE TO CHEMICAL-AGENTS[END_REF]) than the between-subject variability which depends mostly on subject-specific work practices.

In our model, we did not consider subject-specific random slopes although this might be the reality. It would have added complexity to our already complex model and might not be applicable to real data because of the low number of available data per subject.

However, the proposed Bayesian method has some limitations.

First, to run this method, we need a large number of measurements, to be able to estimate the different parameters of variability.

Second, our proposed Bayesian method has been only compared with the mixed-effect tobit model fitted using MLE [START_REF] Jin | Analysis of Lognormally Distributed Exposure Data with Repeated Measures and Values below the Limit of Detection Using SAS[END_REF]. Several analyses of left-censored occupational exposure data have been developed in literature, like substitution methods, β-substitution methods [START_REF] Ganser | An Accurate Substitution Method for Analyzing Censored Data[END_REF], Kaplan-Meier methods, multiple imputation based model, even Bayesian method [START_REF] Huynh | Comparison of Methods for Analyzing Left-Censored Occupational Exposure Data[END_REF][START_REF] Huynh | A Comparison of the beta-Substitution Method and a Bayesian Method for Analyzing Left-Censored Data[END_REF]. However, while βsubstitution methods have been found optimal for parameter estimation in [START_REF] Huynh | Comparison of Methods for Analyzing Left-Censored Occupational Exposure Data[END_REF]) when considering only airborne exposure measurements, none of these methods (except the Bayesian) can be directly extended for the purpose of establishing the relationship between two left-censored variables. Moreover, none of these methods (except [START_REF] Jin | Analysis of Lognormally Distributed Exposure Data with Repeated Measures and Values below the Limit of Detection Using SAS[END_REF])) take into account non-independence in repeated measurements by subject.

The only method that can be directly extended to this context is the mixed-effect tobit model fitted using MLE we refer to as the reference MLE model.

Third, measurements are often made on a single week, and we have to be sure that the halflife of the substance is short, (less than a few hours), to possibly interpret biological measurements to reflect the airborne exposure of the same day. In this case only, does a regression approach make sense. If the half-life of the substance were greater, the biological measurements would be autocorrelated.

Fourth the random effects model for the log transformed urinary concentrations is probably not in accordance with underlying physiological mechanisms. It is precisely at low air/urinary levels that the assumed model might provide a poor representation. In particular, the urinary concentration for many, if not most, substances will be an additive combination of occupational and non-occupational (background) sources (such as diet and environmental exposures). Thus as airborne exposures approach zero, the urinary concentrations do not, but instead approach a 'background' or non-occupationally exposed level. The assumed model on a log-scale does not share this behaviour. Additionally, in the proposed model, the between-worker random effects for the (log) biological monitoring values translate to multiplicative random effects on the natural scale. Whilst this makes sense for the portion of the between-worker effect that is biological, e.g. weight and metabolism, another component in the between-worker variation are variations in the non-occupational sources of exposure, which (on the natural scale) are additive. It was beyond the scope of this paper to include these aspects in our model or in our simulations. In our opinion, what our model most importantly adds, is that when the airborne measurements are below LOD, but the urinary measurements are quantified, we can make use of these data. Data points where both measurements are below LOD will have less influence on the results although limited simulations showed that suppressing them would lead to an attenuation of the slope between air and biological monitoring. Thus the issues mentioned in this paragraph may not predominate especially when we consider LOQs instead of LODs.

Fifth, the assumption of linearity on the log scale, between the two variables, is strong, it cannot be checked on the data below the LOD and it is possibly not realistic for all substances. This must be judged for each practical case, when trying to derive a BLV from existing OELs. An added complexity might be to acknowledge, that the slope for a single substance may differ according to the industrial process that might generate aerosols with different physical and chemical properties, which might influence the biological exposition.

Finally, our simulation study, like all simulation studies, was conducted in a simplified setting.

For example, we did not consider the above-mentioned additive components and only one Homogenous Exposure Group was considered; in practice, as in our real example, this is rarely the case.

As a last point, we consider our practical example. This example was mainly illustrative, but we can discuss some specifics of this data set. The main hypotheses on which this model is based are reasonable. Indeed, the main half-life time is 4.5 hours and the linearity assumption on the log scale fits our data in the range of data above LOQ. However, one limitation could be the assumption of linearity when the airborne exposure is more than an order of magnitude lower than the LOQ. It is probably not reasonable to assume that so small airborne Cr VI exposures would imply level of urinary chromium between 0.05 and 1 µg/g creatinine. Indeed, as mentioned previously, when airborne exposures approach zero, the urinary concentrations approach a non-occupationally exposed level. This is why we commented on the reanalysis excluding the two least exposed companies. Figure 3 shows that the model fit is more reasonable in this case. So that the very similar estimation results

were obtained, validate our method.

Finally, we can compare our estimated BLVs of 5.1 µg/g creatinine to the value (1.8 µg/g creatinine) published by the French Agency for Food, Environmental and Occupational

Health & Safety (ANSES) on the basis of the French recommended OEL 1 µg/m 3 . On the other hand, the ACGIH defined a BEI of 25 µg/L (about 18 µg/g creatinine) based on the 50 µg/m 3 TLV, which is close to our estimate of 23.6 µg/g creatinine. One point to be mentioned is that our censored data are reported as below LOQ rather than below LOD. From a statistical point of view, this does not change anything, if all we know is that some measurements are below a given value. (11.4-25.9) 32 (19%) 23.6 (15.1 -33.9) 15 (9%)

is the between-worker geometric standard deviation of urinary chromium data is the within-worker geometric standard deviation of urinary chromium data is the between-worker geometric standard deviation of airborne exposure chromium data is the within-worker geometric standard deviation of airborne exposure chromium data Table 3: Estimated parameters and BLVs using French OEL or TLV® according to the two methods



  is then set to 0.26 and 0.53 ( =1.3 and =1.7 ) , that is the minimum and the maximum of within-worker GSD of historic urinary data of our own biomonitoring laboratory  is set to 0 and 0.5 ( =1 and =1.65)  is set to 0.47 and 1.53 ( =1.6 and =4.65), defined from the minimum and the maximum of the ratio between and for vapours and non-vapourscomponents (McNally et al. 2014).

Figure 1 :Figure 3 :

 13 Figure 1: Distribution of the slope estimation in the simulated data sets of the main scenario, with balanced and unbalanced data patterns

  The prior distribution of the variability parameter is informative as information exists in the literature on this quantity. Following[START_REF] Mcnally | Advanced REACH Tool: A Bayesian Model for Occupational Exposure Assessment[END_REF], we chose a lognormal with Geometric Mean GM=0.29 and Geometric Standard Deviation GSD=2.82 expressing that lies between 0.036 and 2.22 in 95% of the cases.

	The prior distribution of the variability parameter	is also informative as information exists
	in the literature on this quantity. Following (McNally et al. 2014), we chose a lognormal with
	Geometric Mean GM=0.92 and Geometric Standard Deviation GSD=1.64 expressing that
	lies between 0.34 and 2.38 in 95% of the cases.
	Prior distributions of the parameters	,	and	are non-informative normal
	distributions.

Table 1 :

 1 Pattern with 20 subjects each with 11 measurements for both airborne and biological measurements (**) Pattern with 80 subjects, with 1 measurement for 50 subjects, 2 measurements for 20 subjects and 5 measurements for 10 subjects Description of explored scenarios, (defined by censoring on X and Y, and by variability parameters), and results of the comparison between the two methods for each scenario, on variables of interest:, β , BLV1, BLV2 and BLV3.

	SCENARIO	% Censure X	% Censure Y
	1 Balanced Data (*)	30 50 30 50	30 50	3.49	1.45	1.34	2.66	1.015 1.027 1.028 1.049	1.213 1.399 1.170 1.344	1.005 1.005 1.018 1.032	1.010 1.014 1.023 1.041	1.026 1.053 1.037 1.078
	1 Unbalanced Data (**)	30 50 30 50	30 50	3.49	1.45	1.34	2.66	1.177 1.384 1.193 1.376	1.305 1.614 1.249 1.549	1.057 1.125 1.089 1.132	1.090 1.178 1.105 1.170	1.153 1.310 1.145 1.286
	2	30 50 30 50	30 50	3.49	1.45	1.00	2.66	1.011 1.021 1.031 1.048	1.273 1.379 1.122 1.371	1.009 1.004 1.029 1.031	1.019 1.014 1.033 1.040	1.042 1.052 1.044 1.078
	3	30 50 30 50	30 50	3.49	1.45	1.65	2.66	1.023 1.002 1.027 1.040	1.107 1.361 1.162 1.293	1.017 0.987 1.021 1.022	1.019 1.000 1.027 1.032	1.027 1.044 1.043 1.068
	4	30 50 30 50	30 50	3.49	1.45	1.34	1.60	1.052 1.051 1.137 1.130	1.364 1.547 1.314 1.665	1.039 1.035 1.137 1.121	1.046 1.049 1.144 1.139	1.068 1.097 1.163 1.194
	5	30 50	30	3.49	1.45	1.34	4.65	1.005 1.014	1.124 1.235	1.000 0.998	1.004 1.004	1.016 1.030
		30	50					1.003	1.050	1.004	1.007	1.015

Table 2 :

 2 Description of urinary and atmospheric chromium measurements

		Reference MLE model		Proposed Bayesian Model
		Estimate	95 % CI		Estimate	95 % CI
	Intercept [µg/g creatinine]	1.95	1.63	2.26	1.62	1.41	1.85
	Slope	0.23	0.10	0.36	0.39	0.29	0.49
		1.88	1.46	2.43	1.68	1.36	2.19
		1.30	1.23	1.38	1.43	1.35	1.54
		1.57	1.19	2.06	1.71	1.20	2.50
		2.03	1.77	2.34	2.1	1.85	2.50
		Estimate (95 % CI)	Urinary measurements above BLV N (%)	Estimate (95 % CI)	Urinary measurements above BLV N (%)
	BLV estimation
	[

µg/g creatinine] French OEL = 1 µg/m 3

	7.0 (5.1-9.8)	69 (41%)	5.1 (4.1-6.4)	78 (47%)

BLV estimation [µg/g creatinine] TLV® = 50 µg/m 3 17

.1

ACKNOWLEDGMENTS

We thank Nadège Jacoby for providing the Chromium biomonitoring data that were used in our paper. We thank also Dr Nancy B Hopf for her critical review of an earlier draft of this paper.

DECLARATION

Funding for this project was entirely based on the core funding of the INRS which is the sole employer of the authors. The authors declare no conflict of interest relating to the material presented in this Article. Its contents, including any opinions and/or conclusions expressed, are solely those of the authors.

The main differences are however in the estimate of the slope and in the intercept (corresponding to an exposure of 1 µg/m 3 which happens to be the current OEL).

The most striking result is that incorporating the data with atmospheric measurement below the LOQ results in a significantly steeper slope. This difference in slope is illustrated in Figure 3, which includes virtual measurements below the LOQs sampled randomly from the model.

Note that the estimate of the slope is not only steeper but has a smaller standard error (in relative terms) using the proposed Bayesian method while the intercept has a similar precision in both methods.

Note that if the OEL used is 1 µg/m 3 (current French OEL), the corresponding BLV is 7.0 µg/g creatinine for the reference MLE model and 5.1 µg/g creatinine for the proposed Bayesian method when using informative prior on standard deviation of airborne exposure.

With non-informative prior on standard deviation of airborne exposure, the proposed Bayesian method gave a similar BLV. If we use 50 µg/m 3 TLV ® , the corresponding BLV is 17.1 µg/g creatinine for the reference MLE model and 23.6 µg/g creatinine for the proposed Bayesian method with informative prior on airborne exposure standard deviation, and 25 µg/g creatinine with non-informative priors (data not shown). Very similar values were obtained although with a slightly larger confidence interval when discarding the two least exposed companies (circles and diamonds in Figure 3).