Herwig Wendt
email: herwig.wendt@ens-lyon.fr

Stephane G Roux
email: stephane.roux@ens-lyon.fr

Patrice Abry
email: patrice.abry@ens-lyon.fr

BOOTSTRAP FOR LOG WAVELET LEADERS CUMULANT BASED MULTIFRACTAL ANALYSIS

Multifractal analysis, which mostly consists of estimating scaling exponents related to the power law behaviors of the moments of wavelet coefficients, is becoming a popular tool for empirical data analysis. However, little is known about the statistical performance of such procedures. Notably, despite their being of major practical importance, no confidence intervals are available. Here, we choose to replace wavelet coefficients with wavelet Leaders and to use a log-cumulant based multifractal analysis. We investigate the potential use of bootstrap to derive confidence intervals for wavelet Leaders log-cumulant multifractal estimation procedures. From numerical simulations involving well-known and well-controlled synthetic multifractal processes, we obtain two results of major importance for practical multifractal analysis : we demonstrate that the use of Leaders instead of wavelet coefficients brings significant improvements in log-cumulant based multifractal estimation, we show that accurate bootstrap designed confidence intervals can be obtained for a single finite length time series.

MOTIVATION

Scaling or Multifractal analysis [START_REF] Jaffard | Multifractal formalism for functions, part 2: Selfsimilar functions[END_REF][START_REF] Abry | Wavelets for the analysis, estimation and synthesis of scaling data[END_REF][START_REF] Bacry | Singularity spectrum of fractal signals from wavelet analysis: Exact results[END_REF] is becoming a standard analysis procedure commonly available in empirical data analysis toolboxes. Scaling, or scale invariance, is indeed a property that has been extensively observed in empirical data produced from numerous applications of very different nature. Multifractal analysis mostly consists of measuring scaling exponents, whose values are then commonly involved in various detection, identification or classification tasks. Despite becoming increasingly popular in data analysis, multifractal estimation procedures remain poorly studied. Questions, that may appear natural or simple, such as should one prefer increments or wavelet coefficients?, should one perform weighted or non weighted regressions? or what are the typical sizes of the confidence intervals? still remain insufficiently addressed. However, for practical uses and purposes, elements of answers to such issues are crucial. Indeed, in many real life applications, the sizes of confidence intervals are as important as the values of the scaling exponents themselves, as no classification, discrimination or hypothesis testing are possible without them. In the present contribution, we elaborate on multifractal analysis in combining together three key improvements: wavelet Leaders, log-cumulants and bootstrap. First, it is now considered as classical and powerful to chose wavelet coefficients as the key multiresolution quantities multifractal analysis should be based on [START_REF] Abry | Wavelets for the analysis, estimation and synthesis of scaling data[END_REF][START_REF] Bacry | Singularity spectrum of fractal signals from wavelet analysis: Exact results[END_REF]. Very recent findings reported in [START_REF] Jaffard | Wavelet techniques in multifractal analysis[END_REF][START_REF] Jaffard | Wavelet leaders in multifractal analysis[END_REF] indicate that an accurate multifractal analysis should be based on wavelet Leaders rather than on wavelet coefficients. Indeed, the former enable to estimate exactly the entire multifractal spectrum and to analyze accurately processes containing oscillating singularities when the later do not. For further details, the reader is referred to [START_REF] Jaffard | Wavelet techniques in multifractal analysis[END_REF][START_REF] Jaffard | Wavelet leaders in multifractal analysis[END_REF]. Wavelet Leaders are defined in Section 4.

Second, multifractal estimation procedures are commonly based on structure functions (i.e., power law behaviors of the moments of multiresolution quantities) as in Eqs. [START_REF] Jaffard | Multifractal formalism for functions, part 2: Selfsimilar functions[END_REF] or [START_REF] Abry | Wavelets for the analysis, estimation and synthesis of scaling data[END_REF] below. However, it has been proposed to use instead the cumulants of the logarithm of the multiresolution quantities. This was originally introduced in the early nineties in [START_REF] Castaing | Log-similarity for turbulent flows[END_REF] and largely developed in [START_REF] Delour | Intermittency of 1d velocity spatial profiles in turbulence: A magnitude cumulant analysis[END_REF]. We follow here these promising developments. Cumulant based estimation procedures are described in Section 2. Third, we investigate on potential benefits resulting from the use of non parametric bootstrap for multifractal estimation. In bootstrapping, the distribution of an estimator is approximated through repeated resampling with replacement from the available data. The technique was introduced in the eighties [START_REF] Efron | The Jackknife, the Bootstrap, and Other Resampling Plans[END_REF] and has recently regained interest due to continuously growing computer facilities [START_REF] Zoubir | Bootstrap Techniques for Signal Processing[END_REF][START_REF] Zoubir | On confidence intervals for the coherence function[END_REF][START_REF] Davison | Bootstrap methods and their application[END_REF]. Bootstrapping has been used in the wavelet domain after the pioneering work reported in [START_REF] Percival | Wavestrapping time series: Adaptive wavelet-based bootstrapping[END_REF]. It has also been considered for the estimation of the Hurst parameter of self-similar processes [START_REF] Jach | Wavelet based confidence intervals for the self-similarity parameter[END_REF][START_REF] Sabatini | Wavelet-based estimation of 1/f -type signal parameters: Confidence intervals using the bootstrap[END_REF]. In the present work, we intend to explore the use of bootstrap in two respects: estimation procedure enhancement and confidence intervals derivation. Bootstrap procedures are detailed in Section 3. Therefore, the aims of the present article are to contribute to the answers of the two following questions: Does the use of (log-cumulant) wavelet Leaders improve multifractal estimation procedures ? Can bootstrap provide us with reliable confidence intervals ? To address these questions, we first compare the statistical performance of estimation procedures based on wavelet coefficients and on wavelet Leaders. Second, we use a coverage procedure to compare confidence intervals obtained from a simple bootstrap approach. Results are derived by applying our procedures to a large number of realizations of synthetic scaling processes with a priori known and controlled multifractal properties (cf. Section 5). In Section 6, we show that the use of wavelet Leaders instead of wavelet coefficients brings substantial improvements in multifractal estimation performance. We also clearly demonstrate that the use of bootstrap procedures enables us to obtain highly reliable confidence intervals. We end up with a practical procedure that provides us with both accurate multifractal estimates and confidence intervals and that can actually be used for analyzing a single run of empirical data with finite observation duration.

MULTIFRACTAL, CUMULANTS AND WAVELETS

Definitions Wavelet Coefficients.

Let X(t), t ∈ [0, n) denote the process under analysis and n its observation duration. ψ0(t) is a reference pattern with fast exponential decay, called the mother-wavelet. It is characterized by its number of vanishing moments, a strictly positive integer N ψ ≥ 1 defined as: ∀k = 0, 1, . . . , N ψ -1, R R t k ψ0(t)dt ≡ 0 and R R t N ψ ψ0(t)dt = 0. Let us further denote by {ψ j,k (t) = 2 -j ψ0(2 -j t -k), j ∈ Z, k ∈ Z} templates of ψ0 dilated to scales 2 j , and translated to time positions 2 j k. The wavelet coefficients of X are defined as dX (j, k) = ψ j,k |X . Scaling and Multifractal. A process X is said to possess scale invariance or scaling properties if, for some statistical orders q ∈ [q - * , q + * ] (cf. [START_REF] Lashermes | New insight in the estimation of scaling exponents[END_REF]), the time averages of |dX (j, k)| q taken at fixed scales display power law behaviors with respect to scales a = 2 j ,

|dX (j, •)| q = Fq|a| ζ(q) , (1)
over a wide range of scales a ∈ [am, aM ], aM /am ≫ 1. The ζ(q) are referred to as the scaling exponents of X and are closely related to its multifractal spectrum [START_REF] Jaffard | Wavelet leaders in multifractal analysis[END_REF]. When ζ(q) is linear in q, i.e. ζ(q) = qH, the process X is said to be monofractal. This is, for instance, the case for finite variance self similar processes such as fractional Brownian motion. When ζ(q) = qH, X is said to be multifractal. This is clearly only a poor and operational definition of multifractality. However, for the purposes of this article it is sufficient. We refer the reader to, e.g., [START_REF] Jaffard | Multifractal formalism for functions, part 2: Selfsimilar functions[END_REF], for a thorough introduction to multifractal analysis.

Cumulants. [START_REF] Castaing | Log-similarity for turbulent flows[END_REF][START_REF] Delour | Intermittency of 1d velocity spatial profiles in turbulence: A magnitude cumulant analysis[END_REF] For some processes, Eq. ( 1) is equivalent to

E|dX (j, •)| q = Fq|2 j | ζ(q) . (2)
Using a second characteristic function type expansion, one can rewrite Eq. ( 2) as :

ln Ee q ln |d X (j,•)| = ∞ X p=1 C j p q p p! = ln Fq + ζ(q) ln 2 j , (3)
where C j p stand for the cumulant of order p ≥ 1 of the random variable ln |dX (j, •)|. Combining Eqs (2) and (3) yields that the cumulants of ln |dX (j, •)| have to be of the form:

∀p ≥ 1 : C j p = c 0 p + cp ln 2 j (4)
and therefore that

ln Ee q ln |d X (j,•)| = ∞ X p=1 c 0 p q p p! | {z } ln Fq + ∞ X p=1 cp q p p! | {z } ζ(q) ln 2 j , (5)
where c 0 p and cp do not depend on the scale 2 j . Thus, the measurements of the scaling exponents ζ(q) can be interestingly replaced by those of the log-cumulants cp. This is mainly motivated by the fact that it emphasizes the difference between monofractal (∀p ≥ 2 : cp ≡ 0) and multifractal processes [START_REF] Castaing | Log-similarity for turbulent flows[END_REF][START_REF] Delour | Intermittency of 1d velocity spatial profiles in turbulence: A magnitude cumulant analysis[END_REF]. The next section describes estimation procedures for the cps.

Estimation Procedures

Commonly, the scaling exponents ζ(q) are estimated by linear fits performed in log 2 |dX (j, •)| q vs. log 2 |2 j | plots (see e.g., [START_REF] Abry | Wavelets for the analysis, estimation and synthesis of scaling data[END_REF]). In the present work, we explore the alternative estimation of the equivalent quantities cp.

Cumulant estimations.

Given nj coefficients d(j, k) and thus samples Yj(k) = ln |dX (j, k)|, the asymptotically unbiased and consistent standard estimators (see e.g., [START_REF] Kendall | The Advanced Theory of Statistics, Distribution Theory[END_REF]) are employed to obtain estimates Ĉj p for the cumulants of ln |dX (j, •)|. Linear regressions.

From these Ĉj p s, the cp can then be estimated by linear regression (cf. equation ( 4)),

ĉp = log 2 e j 2 X j=j 1 wj Ĉj p . (6)
Theoretical performance.

Since the Ĉj p s are asymptotically unbiased and consistent, the ĉps are asymptotically unbiased. As detailed in, e.g., [START_REF] Abry | Wavelets for the analysis, estimation and synthesis of scaling data[END_REF], the dX (j, k) of scaling processes are weakly correlated. Hence, one can approximate the variance of ĉp as: Var ĉp ≃ (log 2 e) 2 P j 2 j=j 1 w 2 j Var Ĉj p . Thus, the ĉps are as well consistent. Weights. The weights wj have to satisfy the constraints P j 2 j 1 jwj ≡ 1 and P j 2 j 1 wj ≡ 0 and can be expressed as wj

= 1 b j S 0 j-S 1 S 0 S 2 -S 2 1
, with Si = P j 2 j 1 j i /bj, i = 0, 1, 2. The positive numbers bj are freely selectable and reflect the confidence granted to each Ĉj p . We have chosen to compare three cases, corresponding respectively to i) non-weighted regression, ii) the dX (j, k) can be idealized to independent random variables (cf. [START_REF] Abry | Wavelets for the analysis, estimation and synthesis of scaling data[END_REF]), and iii) the confidence level is set proportional to the inverse of the estimated variances σ2 p (j) of Ĉj p (in the present case, these variances will be estimated by bootstrap): w0,j: bj = 1.

(non-weighted regression) w1,j: bj = 1/nj (assuming Ĉjm p , Ĉjn p uncorrelated) w2,j: bj = σp(j) 2 (σp(j) 2 : estimate of variance of κp,(Ya) ).

BOOTSTRAP

We use bootstrap generated nonparametric empirical distributions (see e.g. [START_REF] Efron | The Jackknife, the Bootstrap, and Other Resampling Plans[END_REF][START_REF] Zoubir | Bootstrap Techniques for Signal Processing[END_REF][START_REF] Zoubir | On confidence intervals for the coherence function[END_REF][START_REF] Davison | Bootstrap methods and their application[END_REF]) for first estimating the variance σ2 p (j) of Ĉj p , and second for constructing confidence intervals for ĉp. As the wavelet coefficients of scaling processes at a given scale are weakly correlated, we adopt a moving blocks bootstrap with overlapping blocks of length L. At each scale a = 2 j , the R bootstrap resamples D * (1) j

, • • • , D * (R) j are generated from the original sample Dj = {dX (j, 1), • • • , dX (j, nj)}. Each resample D * j = n d * (j,1) X (•), • • • , d * (j,n j ) X (
σ2 p (j) = d Var C j * (r) p , (7)
ĉ * (r)

p = log 2 e j 2 X j 1 wj Ĉj * (r) p . (8)
The σ2 p (j) are involved in the calculation of the weights w2,j. The ĉ * (r) p are used to construct 100(1-α)% confidence intervals for the ĉps, according to:

CIp = " Qp " α 2 " , Qp " 1 - α 2 "" = " ĉ * (r 1 ) p , ĉ * (r 2 ) p " (9)
Here Qp(α) is the α-th empirical quantile of the empirical distribution of the R estimates of ĉp, i.e., r1 = ⌊ Rα 2 ⌋ and r2 = R -r1 + 1. Alternatively to this simple bootstrap procedure, the use of pivotal statistics or variance stabilizing transformations can be considered (see e.g., [START_REF] Zoubir | Bootstrap Techniques for Signal Processing[END_REF]). This is currently being investigated.

WAVELET LEADERS

As indicated in Section 1, Wavelet Leaders consists of multiresolution quantities that present significant theoretical and practical qualities to perform Multifractal analysis. Notably, they enable the use of positive and negative q in Eq. 1 as well as a relevant analysis of chirp-type oscillating singularities and hence the correct analysis of the entire spectrum of multifractal properties of X. This has recently been proven in [START_REF] Jaffard | Wavelet techniques in multifractal analysis[END_REF][START_REF] Jaffard | Wavelet leaders in multifractal analysis[END_REF]. Therefore, in the estimation procedures described in Section 2, wavelet Leaders are used instead of wavelet coefficients. Let us introduce the indexing λ j,k = [k2 j , (k +1)2 j ) and the union 3λ j,k = λ j,k-1 ∪ λ j,k ∪ λ j,k+1 . The wavelet Leaders LX (j, k) are defined as

LX (j, k) = sup λ ′ ⊂3λ j,k |d λ ′ |, (10)
where the supremum is taken on the discrete wavelet coefficients dX (•, •) in the time neighborhood 3λ j,k over all finer scales 2 j ′ < 2 j . All relations in Subsection 2.1, in particular Eqs. (1-5) can be rewritten replacing the dX (j, k) by the LX (j, k). Thus, the estimation procedures detailed in Subsection 2.2 and in Section 3 can be rewritten, mutatis mutandis. For convenience, we introduce the superscript d and L to distinguish between estimates involving wavelet coefficients and Leaders, i.e., Ĉj, d p , ĉd p , CI d p involve the dX (j, k) while Ĉj, L p , ĉL p , CI L p involve the LX (j, k).

NUMERICAL SIMULATIONS

Monte Carlo Simulation.

We evaluate the performance of the proposed estimation procedures by applying them to a large number NMC of realizations of synthetic stochastic multifractal processes with known and controlled multifractal properties. From averages • over Monte Carlo realizations, we compute the standard deviations ŝp = p ĉ2 p (i) -ĉp(i) 2 , the biases βp = cp(i) -ĉp(i) and mean-square errors MSEp = q ŝ2 p + β2 p of the proposed estimators. To evaluate the reliability of the confidence intervals obtained from bootstrap, we investigate the coverages produced by re-centered confidence intervals,

CIp,R(i) = CIp(i) -βp, i = 1, • • • , NMC
i.e., confidence intervals that are corrected by the Monte Carlo estimates of the bias of the estimators. This allows us to determine the quality of the confidence intervals independently of the influence of a possible bias of the estimators. The empirical coverages of re-centered confidence intervals are then calculated as:

C emp p = ε (cp, CIp,R(i)) .
Here, ε (cp, CIp,R(i)) = 1 if cp ∈ CIp,R(i) and 0 otherwise: i.e., the empirical coverages C emp p equals the percentage of MC realizations for which the true cp fall within the corresponding re-centered confidence intervals. Scaling Processes.

We use two well known scaling processes, Fractional Brownian motion (FBM) and Multifractal random walk (MRW), chosen because they provide us with simple yet representative examples of Gaussian monofractal processes and non Gaussian multifractal processes respectively. FBM is defined as the only Gaussian exactly self-similar process with stationary increments. Its full definition as well as that of self-similarity can be found in e.g., [START_REF] Samorodnitsky | Stable non-Gaussian random processes[END_REF]. The statistical properties of FBM are entirely determined by the parameter H. FBM possesses scaling properties as in Eq. ( 2), with ζ(q) = qH, for q ∈ (-1, ∞). Thus, c1 = H and cp ≡ 0 for all p ≥ 2. MRW has been introduced in [START_REF] Bacry | Multifractal random walk[END_REF] as a simple multifractal (hence non Gaussian) process with stationary increments: k) , where GH (k) consists of the increments of FBM with parameter H. The process ω is independent of GH , Gaussian, with non trivial covariance:

X(k) = P n k=1 GH (k)e ω(
cov(ω(k1), ω(k2)) = λ ln " L |k 1 -k 2 |+1
" when |k1 -k2| < L and 0 otherwise. MRW has interesting scaling properties as in Eqs. ( 1) or ( 2

) for q ∈ h - p 2/λ, p 2/λ i (cf. [15]
), with ζ(q) = (H + λ)q -λ 2 q 2 /2, hence c1 = H + λ, c2 = -λ 2 and cp ≡ 0 for all p ≥ 3. One sees that the departure from a linear behavior in q is fully controlled by λ (or c2). Simulation Setup.

The results presented here are obtained using Daubechies wavelets with N ψ = 3. Parameters were set to NMC = 1000, n = 2 15 , R = 200, L = 6, H = c1 = 0.8 for FBM and (H, λ) = (0.72, √ 0.08), i.e. c1 = 0.8 and c2 = -0.08, for MRW.

RESULTS

Statistical Performance.

Tables 1 and2 compare the biases and MSEs (respectively) of ĉd p and ĉL p for p = 1 -5, obtained for 1000 realizations of FBM and MRW.

Bias.

Table 1 shows that, while ĉd 1 and ĉL 1 have comparable biases, for p ≥ 2, ĉL p systematically exhibits smaller biases. Note that this discrepancy increases with p and that for p = 5, the difference counts at least 3 orders of magnitude ! Clearly, ĉd p become useless in practise for p = 4, 5, whereas ĉL p continue to give estimates of high accuracy, similar to those produced for p = 1, 2. Also, ĉd 1 has a bias smaller than that of ĉL 1 for (monofractal) FBM, but a larger bias for (multifractal) MRW. Together, these arguments clearly indicate that when the deviations from linearity of ζ(q) are of interest, Leaders must be preferred to wavelet coefficients. From the weight choice point of view, biases are equivalent. Mean Square Error.

Standard deviations are an order of magnitude larger than biases, so that they mostly contribute to MSEs. Hence MSEs only are reported. Table 2 shows that the MSEs of the ĉL p are systematically much smaller than those of the ĉd p for both processes, all weights and all orders. Again, this difference grows with p (roughly as 10 p-1 ), showing that estimates for the cumulant of order 3 can no longer be obtained from wavelet coefficients and require the use of Leaders. For the weight choice issue, we note that the MSEs of ĉL p based on w1,j and w2,j are of comparable order of magnitude in all cases and much smaller than those obtained with w0,j. Four conclusions can be drawn: i) linear regressions for the cumulants must be weighted, ii) the independence assumption of the dX (j, k) underlying the choice w1,j remains valid for Leaders, iii) bootstrap does a good job in estimating the variances of the Ĉj p s, iv) as the choice w2,j involves a much higher computational cost for results equivalent to those obtained with the choice w1,j, this latter is preferred. Fig. 1 displays histograms of ĉd p and ĉL p for p = 1, 2, 3 (weights w1,j) obtained from 1000 realizations of MRW. Whereas the empirical distributions of ĉd p and ĉL p have similar shape and center, those of ĉL p have smaller spread, in particular for p = 2, 3. Again, this suggests the use of ĉL p rather than ĉd p .

Conclusion.

These results lead us to conclude that weighted wavelet Leader based estimators produce the most accurate estimates for log-cumulants. In particular, highly relevant estimates of log-cumulants of higher order are obtained. Hence, we recommend the use of ĉL p with weights w1,j for practical log-cumulant based multifractal analysis. Block Bootstrap.

From the definition of the Leaders LX (j, k) in Eq. [START_REF] Zoubir | On confidence intervals for the coherence function[END_REF], it is clear that they may display complicated intra-and inter-scale (i.e., joint time and scale) correlation structures. Therefore, one may expect poor results for leader-based simple bootstrap estimates. Our results show that this is not the case. Moreover, in Fig. 2 indication in favor of a relevant use of leader-based bootstrap estimation procedures, as described in Section 3. Joint time-scale block bootstrap is however under current investigation.

Confidence Intervals.

Table 3 summarizes the empirical coverages of the re-centered confidence intervals [START_REF] Zoubir | Bootstrap Techniques for Signal Processing[END_REF] for the estimators ĉd p and ĉL p of order p = 1-5. The targeted coverage is 95%. We observe that the bootstrap based procedure produces satisfactory confidence intervals in all cases. Moreover, we note that the coverages obtained with the choice w1,j are highly relevant and quasi systematically the best (or very close to the best). This may be because it allies relevant weighting without having recourse to estimated quantities. This yields two main conclusions. First, bootstrap approaches yield highly relevant confidence intervals for multifractal estimation. Therefore, we highly recommend their use in practical multifractal analyses. Second, the choice w1,j is to be favored. Improvements resulting from the use of pivotal statistics (see e.g., [START_REF] Zoubir | Bootstrap Techniques for Signal Processing[END_REF]) are under current investigation.

Practical procedure and regression range.

From these analyzes, we have designed a MATLAB routine that implement wavelet coefficient and wavelet Leader based log-cumulant multifractal analysis together with bootstrap based confidence intervals. Therefore, it enables us to obtain from a single observed times series with finite length, both estimates for the cps and error bars. FBM MRW Estimate w0,j w1,j w2,j w0,j w1,j w2,j 3. Empirical coverage of re-centered 95% confidence interval for log-cumulant estimates from wavelet coefficients (top) and wavelet Leaders (bottom) for FBM and MRW and p = 1 -5. Results closest to target coverage are marked in bold. This significantly improves already available practical multifractal estimation procedures as for most applications error bars are as important as estimates themselves. We see this as a major result of the present contribution. Fig. 3 illustrates this procedure at work and shows logscale diagrams (for p = 1, 2, 3): Ĉj, L p (or Ĉj, d p ) as a function of j for a single realization of MRW, together with corresponding regression lines and ±1.96σp(j) error bar estimates obtained from bootstrap. Whereas all Ĉj, L p s display a highly linear behaviour over a large range of scales j, even for large ps, a zone of linearity is more difficult to find for Ĉj, d 2 and Ĉj, d 3 , suggesting estimates of poorer quality in these cases. Clearly, the selection of the regression range of octaves is a key practical issue that adds an extra complication to the analyzes reported here. Again, one sees that for Leader based multifractal, selecting the regression range should be easier.

CI d 1,

CONCLUSION AND PERSPECTIVES

In the present work, we compared various log-cumulants multifractal estimation procedures. First, we demonstrated that the use of wavelet Leaders instead of wavelet coefficients brings substantial improvement in estimation performance. In particular, highly accurate estimates for log-cumulants of order p ≥ 2 can be obtained. To the best of our knowledge, this had never been illustrated clearly before. Second, we showed that the simple bootstrap approach provides us with highly relevant confidence intervals for the estimates ĉd p and ĉL p . This is another major improvement as, to the best of our knowledge, this is the first time that a non parametric confidence interval estimation procedure, with excellent performance, is obtained for multifractal analysis. Wavelet Leaders and bootstrap confidence intervals together lead to the design of a multifractal analysis procedure (available in MATLAB upon request) of primary interest for the exploration of empirical data with possibly multifractal properties.

A key feature of the results obtained in this work lies in the fact that they hold for both Gaussian monofractal and non Gaussian multifractal processes. This is very promising as it opens the track for the design of hypothesis tests aiming at discriminating between monoand multi-fractal processes, and between different multi-fractal processes, two major practical issues. These ideas are currently under investigation.

  •) o represents an unsorted collection of nj sample points, drawn blockwise and with replacement from the original sample. These collections D * (r) j are used to compute R bootstrap cumulant estimates r = 1, . . . , R, Ĉj * (r) p . In turns, these C j * (r) p are used for obtaining i) variance estimates for Ĉj p , and ii) R bootstrap ĉp:

3 Fig. 1 .

 31 Fig. 1. Histograms of ĉ1 (left), ĉ2 (center) and ĉ3 (right) for wavelet coefficients (top) and wavelet Leaders (bottom) for 1000 realizations of MRW and weights w1,j (note the different scale of the xaxis for ĉd 3 and ĉL 3 .

3 Fig. 2 .

 32 Fig. 2. Histograms of ĉL p obtained from 200 realizations of MRW (top row) and bootstrap resamples ĉL * p obtained from one single realization chosen at random (bottom row), for p = 1, 2, 3 (from left to right).

Fig. 3 .

 3 Fig. 3. Logscale diagram of Ĉj, d p (left) and Ĉj, L p (right) for a single realization of MRW, with weighted regression line with weights w1,j over scales j1 = 3 to j2 = 10. The error bars ±1.96σp(j) are estimated by bootstrap.

Table 1 .

 1 , histograms for ĉL p , based on 200 realizations of MRW (cf. below), and for 200 corresponding bootstrap estimates ĉL, * Bias (× 10 3 ) of estimators ĉp for FBM (left) and MRW (right) and p = 1 -5. Best results are marked in bold.

	BIAS		FBM			MRW
	Estimator	w0,j w1,j	w2,j	w0,j	w1,j	w2,j
	ĉd 1	0.4	0.8	-3.3	32.8	31.2	28.6
	ĉL 1	2.2	5.8	6.2	28.3	30.2	28.4
	ĉd 2	10.7	1.3	48.3	-0.1	-10.6	33.1
	ĉL 2	1.5	-2.8	0.6	-1.8	-10.0	1.0
	ĉd 3	-51.3	-3.5 -267.9	-36.7	-6.7 -258.9
	ĉL 3	0.9	1.2	1.4	2.2	-0.9	2.0
	ĉd 4	394.9 54.3 1059.6	352.8	86.0 1078.0
	ĉL 4	0.8	0.6	0.2	5.9	3.2	1.7
	ĉd 5	-3513.8 -783.3 -3235.6 -3344.8 -1036.9 -3232.7
	ĉL 5	0.0	-0.2	-0.4	7.8	4.2	3.0
	MSE		FBM			MRW
	Estimator	w0,j	w1,j	w2,j	w0,j	w1,j	w2,j
	ĉd 1	32.8	15.5	16.6	47.0	35.3	34.3
	ĉL 1	19.8	10.8	12.0	38.1	32.8	32.9
	ĉd 2	77.6	37.3	65.6	87.8	42.8	60.1
	ĉL 2	6.5	4.1	4.6	24.1	17.5	20.6
	ĉd 3	340.0 187.7 302.4	399.8	200.7 300.7
	ĉL 3	3.1	1.8	1.9	23.4	18.4	19.9
	ĉd 4	1822.4 1251.2 1125.4 2505.8 1366.4 1164.0
	ĉL 5	2.0	1.0	0.8	31.0	30.0	25.0
	ĉd 5	9622.5 9803.0 3621.3 15741.9 11068.3 3673.2
	ĉL 5	2.6	0.7	0.6	72.9	50.1	52.3

p , based on a single realization, are shown. The closeness of Monte Carlo and bootstrap empirical distributions for ĉL p provides us with a clear

Table 2 .

 2 MSE (× 10

3

) of estimators ĉp for FBM (left) and MRW (right) and p = 1 -5. Best results are marked in bold.

  Table

	R	85.7 92.1 89.5 90.1 98.6 95.5
	CI d 2,R	79.4 92.3 75.3 79.1 95.5 79.6
	CI d 3,R	72.3 91.2 66.9 71.2 92.6 72.9
	CI d 4,R	64.3 87.8 85.0 62.2 90.3 86.8
	CI d 5,R	85.7 87.8 76.6 83.8 90.8 82.5
	CI L 1,R	71.8 83.4 76.6 84.5 98.8 89.0
	CI L 2,R	76.1 90.3 73.3 75.7 97.0 71.3
	CI L 3,R	92.6 94.7 94.7 80.4 96.4 75.6
	CI L 4,R	96.3 97.2 89.7 85.1 95.4 79.2
	CI L 5,R	96.1 98.0 98.2 85.0 96.6 82.4


Jean Stawiaski

Etienne Decencière

François Bidault

Interactive Liver Tumor Segmentation Using Graph-cuts and Watershed

We present in this paper an application of minimal surfaces and Markov random fields to the segmentation of liver tumors. The originality of the work consists in applying these models to the region adjacency graph of a watershed transform. We detail the assumptions and the approximations introduced in these models by using a region graph instead of a pixel graph. This strategy leads to an interactive method used to delineate tumors in 3D CT images. We detail our strategy to achieve relevant segmentations of these structures and compare our results to hand made segmentations done by experienced radiologists. This paper summarizes our participation to the MICCAI 2008 3 workshop called: "3D segmentation in the clinic : A Grand Challenge II".

Introduction

Possible applications related to liver tumor segmentation are mainly radiotherapy and surgery planning. In both scenarios, the knowledge of the exact location and volume of the tumors is a key problem. Liver tumors present low contrasted boundaries and exhibit a large variability of shapes, sizes and locations in the liver. Due to these multiple difficulties, automatic or model based approaches seem to be inadequate. A few semi-automatic approaches have already been proposed in the literature. Recent approaches are mainly based on active contours [START_REF] Lu | Liver tumor volume estimation by semiautomatic segmentation method[END_REF], level-sets [START_REF] Popa | Tumor volume measurement and volume measurement comparison plug-ins for volview using itk[END_REF], as well as machine learning [START_REF] Li | A machine learning approach for locating boundaries of liver tumors in ct images[END_REF]. These methods have the interesting property to allow an interaction with its user through land marks positioning or interactive refinement of the segmentation. The interactive approach of medical images segmentation is, to our mind, the most reliable method to provide robust results. We propose in this paper an interactive segmentation method based on watershed and graph cuts to extract liver tumors boundaries in 3D CT images.

The studied CT images 4 were acquired on one 64-slice and two 40-slice CT scanners using a standard four-phase contrast enhanced imaging protocol. The resulting images have a slice thickness of 1mm or 1.5mm and an in-plane resolution of 0.6-0.9mm. Depending on several parameters such as the patient size and cardiac rhythm, the images present different contrasts in the liver. Several cases are illustrated in figure 1. The liver segmentation presents some difficulties that have to be taken into account to design a relevant segmentation protocol. Depending on the acquisition time, tumors can be in positive or negative contrast with the liver. Secondly, the tumors boundaries are not well defined and perceptual properties have to be used to define the exact contours of the tumors.

Images segmentation strategy

The first step of our methodology is the manual definition of a sub-volume containing one or more tumors that need to be segmented in order to reduce the computation time of the segmentation algorithm. Our strategy is then entirely based on the use of the region adjacency graph of an unsupervised watershed segmentation [START_REF] Beucher | The morphological approach to segmentation: the watershed transformation[END_REF], as originally proposed by Li et al. [START_REF] Li | Lazy snapping[END_REF]. From this first low-level segmentation, a region adjacency graph is extracted and used for the following optimization steps. Our protocol is motivated by the following simple observation: the liver presents two kinds of tissues: tumoral and healthy tissues. The classification of liver pixels in one of these two classes provides the tumor boundaries. To achieve this classification we model the liver pixels as a Markov random field and the classification is performed through the maximum a posteriori estimation. This classification step is supervised by user defined markers that specify both tumoral and healthy tissues. The markers are used to locate the tumors and to estimate the grey levels characteristics of these structures. However, the liver pixels classification needs also that the liver boundaries are extracted. This task is realized by computing a minimal surface based on user defined markers. The user has finally to specify normal liver tissues, tumoral tissues, and external tissues surrounding the liver. Note that the liver boundaries extraction is only necessary if the sub-volume considered for the segmentation contains also nonliver tissues. In any other cases, the liver boundaries are not computed, and only the classification step is performed. The different energy minimization strategies (minimal surfaces and Markov random field) are based on the computation of a minimal graph cut using Boykov and Kolmogorov algorithm described in [START_REF] Boykov | An experimental comparison of min-cut/maxflow algorithms for energy minimization in vision[END_REF] and implemented in the Boost Graph Library5 .

The watershed transform

An unsupervised watershed transform of the morphological gradient of the original CT image is used in our work to produce a region adjacency graph. Both minima of the gradient and the watershed transform are computed using the 6-neighborhood adjacency system. The watershed transform [START_REF] Beucher | The morphological approach to segmentation: the watershed transformation[END_REF] allows to obtain a partition of an image composed of small and numerous homogeneous regions. Moreover important contours of the image are preserved during the segmentation and regions of the partition are mostly composed of homogenous pixels (pixels of similar grey values). The quality of this first unsupervised segmentation is important to guarantee a minimal loss of information, an ideal situation would be that important information (contours and/or homogenous regions) about the original image is accessible from this first segmentation. These observations about the watershed transform are not theoretically guaranteed but are verified when working on natural images. Another important point is that the watershed transform algorithm based on hierarchical queues exhibits a linear complexity [START_REF] Meyer | Morphological segmentation[END_REF]. The time needed by the watershed transform is in practice negligible compared to a standard graph cuts algorithm [START_REF] Boykov | An experimental comparison of min-cut/maxflow algorithms for energy minimization in vision[END_REF].

Approximate minimal surfaces

We detail now how to extract the liver boundaries by an approximate minimal surface using a region adjacency graph [START_REF] Stawiaski | Computing approximate geodesics and minimal surfaces using watershed and graph-cuts[END_REF]. The combination of graph-cuts with a watershed low-level segmentation provides us with an explicit and efficient way to compute approximate minimal surfaces. Our basic assumption is that the minimal surface to be computed is embedded in the watershed low-level segmentation contours. We propose thus to solve the following combinatorial problem: finding a surface composed of a finite union of watershed contours such that the surface minimizes a given geometric functional. We solve this problem by using graph-cuts optimization on a region adjacency graph.

Following the formulation of Caselles et al. [START_REF] Caselles | Geodesic active contours[END_REF], we want to find a surface S defined by a finite union of watershed contours that minimizes the following energy function:

E(S) = S g( ∇I(x, y) )dxdy ( 1
)
where g is a positive and strictly decreasing function and ∇I(x, y) is the modulus of the gradient of the image I (image contrast) along the surface S.

Note that Cauchy-Crofton formulaes can be used to minimize the energy function E(S) by computing a minimal graph cut as described by Boykov et al. in [START_REF] Boykov | Computing geodesics and minimal surfaces via graph cuts[END_REF].

Let us consider G = (V, E, W ) as the pixel graph of an image I. Classically V is the set of nodes and represents the pixels of I, E is the set of edges representing neighborhood relations between pixels and W is a positive weight assigned to each edge of E. In our terminology, an edge linking two nodes i and j is written e i,j and the corresponding edge weight is denoted by w i,j . From the pixel graph, we define the region adjacency graph G R = (V R , E R , W R ) of the watershed transform where V R is the set of nodes (i.e the regions of the watershed transform). E R is the set of edges (i.e the neighborhood relation between regions) and W R is the weights of the edges. Let us define F (ri,rj ) as the set of edges of the pixel graph connecting two regions r i and r j of the low-level watershed segmentation:

F (ri,rj ) = {e m,n ∈ E | m ∈ r i , n ∈ r j } . (2)
Note that the set F (ri,rj ) depends on the adjacency system of the pixel graph G. The set of edges of the pixel graph describes also implicitly a set of surfaces between the regions r i and r j as illustrated in figure 2. Let S (ri,rj ) denote the set of surfaces that could cross the edges of F (ri,rj ) . Following Cauchy-Crofton formulaes with the V 6 adjacency system, the energy function E(S (ri,rj ) ) can be approximated by:

E(S (ri,rj ) ) ≈ (em,n∈F (r i ,r j ) ) g(max ( ∇I(m) , ∇I(n) )) , (3)
where ∇I(m) and ∇I(n) are the gradient magnitudes of the end points of e m,n . In the following, we consider the strictly positive and decreasing function g:

g( ∇I(p) ) = 1 1 + ∇I(p) k . (4)
The parameter k ∈ R + is a free parameter that can be used as a smoothing term as shown by Allène et al. in [START_REF] Allène | Some links between min-cuts, optimal spanning forests and watersheds[END_REF]. In our application this parameter was set to k = 2. The function g works as an edge indicator of the image I and takes a small value if neighbors pixels m and n take different grey values p m and p n .

The energy E(S (ri,rj ) ) of the boundary between two regions is simply obtained by summing the local contrasts along the boundaries between two regions. The edge weights of the region adjacency graph are then set such that the weight of a graph cut equals the energy function of the surface it implicitly defines:

w ri,rj = E(S (ri,rj ) ) . (5)
The liver boundaries are finally extracted by computing a minimal graph cut of the region adjacency graph with weights given by equation 5. The minimal cut is computed on the region adjacency graph with two additional nodes s and t, respectively connected to the markers of the liver and the markers of the external tissues. In the following, we denote the markers that specify the liver tissues as the set of regions M t and M h , respectively for tumoral and healthy tissues. The markers specifying the tissues surrounding the liver are denoted by M ext . The edge weights of the graph are summarized in table 1. 1. Edge weights for approximate minimal surfaces.

Edge Weight for

ws,r i +∞ ri ∈ Mt ws,r i +∞ ri ∈ M h wr i ,t +∞ ri ∈ Mext wr i ,r j E(S (r i ,r j ) ) ri ∈ VR, rj ∈ Nr i Table

Approximate maximum a posteriori estimation of a Markov random field (MRF)

This section details the segmentation method used to detect the tumors in the liver. We are now going to take into account a second assumption about the watershed transform: the unsupervised watershed transform of a natural image is composed of regions of homogenous grey level intensities. This assumption permits us to model the image to be segmented as a Markov random field [START_REF] Boykov | Markov random fields with efficient approximations[END_REF][START_REF] Boykov | Interactive graph cuts for optimal boundary and region segmentation of objects in n-d images[END_REF], where each random variable corresponds to the mean value inside a region of the watershed transform [START_REF] Stawiaski | Combining graph-cuts and morphological segmentation[END_REF].

Let us consider the pixel graph G = (V, E, W ) of an image, as well as the corresponding region adjacency graph G R = (V R , E R , W R ) of its watershed transform. The binary image restoration problem [START_REF] Geman | Stochastic relaxation, gibbs distributions, and the bayesian restoration of images[END_REF][START_REF] Boykov | Markov random fields with efficient approximations[END_REF] is classically defined as the labeling X of the nodes V :

X : V → {0, 1} i → X(i) = x i
such that X minimizes:

E(X) = i∈V -ln(P r(p i |x i )) + i∈V j∈Ni u i,j .δ(x i = x j ) , (6)
where p i ∈ [0, 1] is the grey level of the pixel i, and x i ∈ {0, 1} is a label that has to be assigned to the pixel i. N i is the set of neighbors of the pixel i and u i,j is a positive function. δ(x i = x j ) is the indicator function: δ(x i = x j ) = 1 if and only if x i = x j and δ(x i = x j ) = 0 otherwise.

This equation describes the classical formulation of the maximum a posteriori estimation of a Markov random field [START_REF] Boykov | Markov random fields with efficient approximations[END_REF]. The data term, also called the likelihood function P r(p i |x i ), ensures that dark pixels, p i ≈ 0, will be assigned the label x i = 0; and that bright pixels, p i ≈ 1, will be assigned the label x i = 1. This model is regularized with the prior function u i,j which is typically used to guarantee that the resulting segmentation has smooth boundaries. We detail now how each term of this energy function can be approximated by using the region adjacency graph of the watershed transform.

Since we assume that each region of the watershed transform is composed of region of homogenous intensities, we can approximate the likelihood term by:

i∈V -ln(P r(p i |x i )) ≈ i∈V R -|r i |.ln(P r(µ ri |x i )) , (7)
µ ri = 1 |r i | pi∈ri p i , (8)
where |r i | is the number of pixels inside the region r i .

Assuming that the data are corrupted with a white gaussian noise, the likelihood function can finally be written as:

P r(µ ri |(x i = 0)) = exp(- (µ ri -µ (xi=0) ) 2 2σ 2 0 ) , (9)
P r(µ ri |(x i = 1)) = exp(-

(µ ri -µ (xi=1) ) 2 2σ 2 1
) ,

where µ ri is the mean gray level of the region r i , and µ (xi=0) and µ (xi=1) are the mean values of the pixels expected to take the values x i = 0 and x i = 1.

In our experiments the values of σ 0 and σ 1 were empirically set to σ 0 = 0.25 and σ 1 = 0.2. In our model, σ 0 represents the grey level variance of the tumoral tissues and σ 1 represents the grey level variance of the healthy tissues. We have experimentally found that the variance of tumoral tissues is slightly higher than the grey level variance of the healthy tissues. On the other side, the values of µ xi are estimated from the user defined markers as:

µ (xi=0) = 1 |M t | ri∈Mt µ ri , (11)
µ (xi=1) = 1 |M h | ri∈M h µ ri , (12)
where |M h | and |M t | are respectively the number of regions marked as healthy or tumoral. We recall that the markers that specify the liver tissues are denoted as the set of regions M t and M h , respectively for tumoral and healthy tissues.

The markers specifying the tissues surrounding the liver are denoted by M ext .

On the other side, the prior function depends on the boundaries properties of the labeling x. We assume now that all pixels inside a region of the watershed transform are assigned the same label x i , the prior function does thus only depend of the pixels lying in the boundaries between two regions:

i∈V j∈Ni u i,j .δ(x i = x j ) ≈ i∈V R j∈Nr i |F (ri,rj ) |.u i,j .δ(x i = x j ) , (13)
where N ri is the set of neighbor regions of the region r i and |F (ri,rj ) | is the number of edges of F (ri,rj ) .

The prior function u i,j that we use in our application is a contrast sensitive function:

u i,j = β -β * (µ ri -µ rj ) n . ( 14
)
where n is a free parameter describing the strength of the term β * (µ ri -µ rj ) and takes also into account the local contrast. In the following we set the parameter to n = 4. Note that if n is very large, our contrast sensitive model is equal to the classical Ising model of ferromagnetism [START_REF] Ising | Beitrag zur theorie des ferromagnetismus[END_REF]. The prior function takes into account the contrast between two regions and is equal to a constant in the areas where the contrast is low. This function allows to detect correctly the boundaries of high contrasted tumors, whereas low contrasted boundaries are smoothed such that the surface of the object is minimized.

We have finally to minimize the following energy function:

E(x) = i∈V R |r i | (µ ri -µ xi ) 2 2σ 2 + i∈V R j∈Nr i |F (ri,rj ) |.u i,j .δ(x i = x j ) . (15)
As shown in [START_REF] Boykov | Markov random fields with efficient approximations[END_REF], we can minimize this energy function by computing a minimal graph cut. The minimal cut is computed on the region adjacency graph with two additional nodes s and t, respectively connected to the markers of the tumoral tissues and the markers of the healthy liver tissues. The edge weights of the graph are summarized in table 2.

Edge

Weight for 2. Edge weights for approximate maximum a posteriori estimation of a MRF.

ws,r i +∞ ri ∈ Mt wr i ,t +∞ ri ∈ M h ws,r i (µr i -µ (x i =0) ) 2 2σ 2 0 ri ∈ VR \ {Mt ∪ M h } wr i ,t (µr i -µ (x i =1) ) 2 2σ 2 1 ri ∈ VR \ {Mt ∪ M h } wr i ,r j |F (r i ,r j ) |.ui,j ri ∈ VR, rj ∈ Nr i Table

Post-processing

An additional post processing step is also proposed to the user after the extraction of the segmented tumors. This step consists in smoothing the segmentation by using a morphological opening of the object representing the tumor. The opening is computed with a structuring element of size 1 using the V6 adjacency system. This additional step permits to obtain slightly smoother tumors boundaries.

Example

Figure 3 illustrates our segmentation strategy on a single slice of a 3D CT image. The obtained results are in good concordance with the expected results obtained by a hand made segmentation. Alternatively, the user can add or delete markers if he is not satisfied with the computed segmentation.

Graphical user interface

We have developed a graphical user interface dedicated to 3D medical image segmentation and visualization. The software is entirely developed in python and C++ based on the VTK6 and our own image analysis library Morph-M7 .

Our software allows the user to explore a highly detailed view of the data-set for easy interpretation. Visualization is particularly important for segmentation validation purposes. Data sets can be explored through 2D orthogonal cuts of the 3D volume and 3D rendering of the whole image. The user can interactively provide the markers needed for the segmentation algorithms by drawing on slices of the 3D image, as illustrated in figure 4. The user can visualize the segmented image as a set of surfaces, one surface for each object and surimpose it on the original image. Combination of all this visualization methods allows an easy and fast interpretation of the segmentation result.

Results

Table 3 summarizes the evaluation scores of our method on a set of 5 CT images presenting 10 tumors with unknown hand made segmentations. The evaluation scores compare our results with the radiologists segmentations. The important point is that these results have been obtained without the knowledge of the hand made segmentations. The mean surface distance between our segmentations and the references is approximatively one and a half millimeter, which represents 2 to 3 voxels of the studied 3D CT images. The volumetric overlap error shows that approximatively 71 % of our segmentation volume is in perfect match with a hand made segmentation. Some evaluation results show that we have misunderstood some structures that had to be extracted. This problem leads to very low scores (see for instance tumor number 5). However our results are promising, considering the low quality of some images of the dataset. In comparison, the mean total score obtained on the training data set with known segmentations was equal to 88. The computation time of our method depends on the image size and the number of refinements of the segmentation. The first step, which consists in extracting a sub-volume, is typically done in one minute. The markers placement requires basically one to two minutes. In our experiments we have drawn markers in three orthogonal slices centered on the tumor in about two minutes. The computation time of the segmentation algorithms is then relatively fast: about 1 second for the watershed transform, 2 seconds for the liver extraction and 2 other seconds for the tumor extraction (for a typical sub-volume of size 100×100×100, computed on a common personal computer). The time needed for the refinement steps is then mainly due to the markers placement. In our experiments, two to five refinement steps were needed to provide the presented results. The typical time spent for these refinements varied from two to five minutes. The total time needed for the segmentation of a tumor is then approximatively equal to five, up to height minutes.

Conclusion

Our method exhibits promising results for the aimed application. However some open problems still remain. First, the segmentation of multiple tumors in the same liver often requires additional user markers to correctly separate the tumors. The developed method merges the tumors in a single object when different tumors are too close. This problem requires that the user adds markers between the merged tumors. This additional interaction speeds down the segmentation protocol. However the used methods (minimal surfaces and Markov random fields on a region adjacency graph) are fast and can be used interactively. Secondly we did not develop any preprocessing step such as filtering of the images. There is thus still some possible improvements of our methodology. Future work will be concentrated on the development of adapted filters to simplify the segmentation and the classification step of our methodology.

The use of a region adjacency graph offers a good trade-off between speed and precision for the computation of minimal surfaces and maximum a posteriori estimation of a MRF. We want also to point out that a region based approach is potentially richer than a pixel approach since a wide class of geometrical and statistical functionals can be computed on each region of the watershed transform and additional constraints, such as curvature smoothing term, could be added to the energy function to minimize. These approaches were already successfully used for various image segmentation problems for medical and material sciences applications [START_REF] Stawiaski | Computing approximate geodesics and minimal surfaces using watershed and graph-cuts[END_REF][START_REF] Sørensen | Virtual open heart surgery: Obtaining models suitable for surgical simulation[END_REF][START_REF] Stawiaski | Proceedings of Medicine Meets Virtual Reality 15[END_REF].

Fig. 1 .

 1 Fig. 1. Liver 3D CT images. Tumors are indicated by red arrows and the aorta is indicated by a blue arrow. The brightness of the aorta indicates the quantity of contrast agent present in it.

Fig. 2 .

 2 Fig. 2. (a) A region adjacency graph. (b) The set of nodes of the pixel graph considered to compute boundary properties between two regions, with a V4 adjacency system. (c) A curve crossing the edges of the boundary between two regions r1 and r2.

Fig. 3 .Fig. 4 .

 34 Fig. 3. Liver tumor segmentation. (a) User specified markers. In blue, the tumor markers, in red, the liver markers and in green, the external markers. (b) Results of our segmentation strategy. (c) Radiologist hand made segmentation and the liver contours extracted by our method.

Table 3 .

 3 Results of comparison metrics and scores for all ten test tumors.

		Overlap	Volume		Average	RMS	Maximum	Total
		Error	Difference	Surface	Surface	Surface	Score
					Distance	Distance	Distance
	Tumor	(%) Score (%) Score (mm) Score (mm) Score (mm) Score
	IMG05 L1 27,16 79 21,68 77 2,01 49 2,90 60 10,50 74	68
	IMG05 L2 36,41 72 24,98 74 1,36 66 1,83 74 6,15 85	74
	IMG05 L3 31,99 75 16,93 82 1,18 70 1,59 78 5,64 86	78
	IMG06 L1 33,19 74	1,86	98 0,84 79 1,09 85 4,00 90	85
	IMG06 L2 61,24 53 119,82 0	2,29 42 2,95 59 9,30 77	46
	IMG07 L1 21,83 83 10,68 89 2,65 33 3,58 50 18,31 54	62
	IMG07 L2 21,44 83	4,85	95 0,93 77 1,38 81 6,89 83	84
	IMG08 L1 16,02 88 14,47 85 1,84 54 2,70 62 9,73 76	73
	IMG09 L1 22,87 82	2,14	98 0,65 83 0,97 87 5,10 87	87
	IMG10 L1 22,78 82 21,27 78 1,21 69 1,73 76 7,33 82	77
	Average 29,49 77 23,87 78 1,50 62 2,07 71 8,29 79	73

MICCAI 2008 is the 11th International Conference on Medical Image Computing and Computer Assisted Intervention: http://miccai2008.rutgers.edu/

The images are provided by the organizers of the workshop "3D segmentation in the clinic : A Grand Challenge II": Wiro Niessen, Martin Styner, Simon K. Warfield and Xiang Deng. http://grand-challenge2008.bigr.nl/index.php.

http://www.boost.org/

www.vtk.org, The Visualization ToolKit (VTK) is an open source, freely available software system for 3D computer graphics, image processing, and visualization.

http://cmm.ensmp.fr/Morph-M/, Morph-M is the result of the work of several researchers at the Centre for Mathematical Morphology.

The author has requested enhancement of the downloaded file. All in-text references underlined in blue are linked to publications on ResearchGate. The author has requested enhancement of the downloaded file. All in-text references underlined in blue are linked to publications on ResearchGate.

Acknowledgements

We would like to thank Région Ile-de-France and the Cancéropôle Ile-de-France for funding our research on medical image segmentation.


Jean Salzemann

Vinod Kasam

Nicolas Jacq

Astrid Maass

Horst Schwichtenberg

Vincent Breton

Grid Enabled High Throughput Virtual Screening Against Four Different Targets Implicated in Malaria

Keywords: large scale deployment, Computational grids, Malaria, In silico docking, Virtual Screening, WISDOM

After having deployed a first data challenge on malaria and a second one on avian flu, respectively in summer 2005 and spring 2006, we are demonstrating here again how efficiently the computational grids can be used to produce massive docking data at a high-throughput. During more than 2 months and a half, we have achieved at least 140 million dockings, representing an average throughput of almost 80,000 dockings per hour. This was made possible by the availability of thousands of CPUs through different infrastructures worldwide. Through the acquired experience, the WISDOM production environment is evolving to enable an easy and fault-tolerant deployment of biological tools; in this case it is the FlexX commercial docking software which is used to dock the whole ZINC database against 4 different targets.

Introduction

WISDOM stands for World-wide In Silico Docking On Malaria. Malaria together with many other tropical and protozoan diseases is one of the most neglected diseases by the developed countries as well as by the pharmaceutical industries. Plasmodium is the protozoan genus causing malaria. Due to very high costs associated to the drug discovery process as well as due to late stage attrition rates, novel and cost effective strategies are absolutely needed for combating the neglected diseases, especially malaria [START_REF] Breman | Conquering the intolerable burden of malaria: what's new, what's needed: a summary[END_REF].

In silico screening of chemical compounds against a particular target is termed as Virtual Screening. The costs associated to the virtual screening of chemical compounds are significantly reduced when compared to screening of compounds in experimental laboratory. Beside the costs, virtual screening is fast and reliable [START_REF] Bleicher | Hit And Lead Generation: Beyond High-Throughput Screening[END_REF][START_REF] Boehm | Virtual Screening For Bioactive Molecules. Chapter 1, High Throughput Screening and Virtual Screening Entry Points To Drug Discovery[END_REF]. However, it is computationally intensive: docking a single compound within the active site of a given receptor requires about 1 minute CPU. With the development of combinatorial chemistry technology, millions of different chemical compounds are now available in digital databases [START_REF] Spencer | High throughput virtual screening of historic collections on the file size, biological targets, and file diversity[END_REF]. To screen all these compounds and store the results is a real data challenge. To address this problem computational grid infrastructures are used. WISDOM-I [START_REF] Jacq | Grid enabled virtual screening against malaria[END_REF] is the first large scale deployment of molecular docking application on EGEE grid infrastructure. It took place from August 2005 to September 2005 and achieved 41 million dockings which is equivalent to 80 CPU years. The docking was performed on Plasmepsins, a aspartic protease involved in haemoglobin degradation. On the biological front three scaffolds were identified, of them one is guanidino scaffold which is likely to be novel as it was not known as a plasmepsin inhibitor before [START_REF] Kasam | Design of Plasmepsin Inhibitors: A Virtual High Throughput Screening Approach on the EGEE Grid[END_REF].

With the success achieved by the WISDOM-I project both on the computation and biological sides, several scientific groups around the world proposed targets implicated in malaria, which led to the second assault on malaria, WISDOM-II.

Materials and methods

Virtual Screening by molecular docking requires a target structure, a chemical compound database and docking software. The targets used in the current project are and Glutatione-S-trasferase (GST, pdbid: 1Q4J) [START_REF] Perbandt | Native and inhibited structure of a Mu class-related glutathione S-transferase from Plasmodium falciparum[END_REF], Plasmodium falciparum Dihydrofolate reductase (DHFR) wild type (pdbid: 1J3I), quadrupule mutant (pdbid: 1J3K) [START_REF] Kongsaeree | Crystal Structure of Dihydrofolate Reductase from Plasmodium Vivax: Pyrimethamine Displacement Linked with Mutation-Induced Resistance[END_REF], Plasmodium vivax Dihydrofolate reductase wild type (pdbid: 2BL9), double mutant (pdbid: 2BLC) [START_REF] Yuvaniyama | Insights into antifolate resistance from malarial DHFR-TS structures[END_REF]. In another experiment the same structures of Plasmodium vivax Dihydrofolate reductase wild type (pdbid: 2BL9), double mutant (pdbid: 2BLC) but after minimization are used. The chemical compound database used is ZINC database [START_REF] Zinc Database | ZINC -A Free Database of Commercially Available Compounds for Virtual Screening[END_REF][START_REF]A Free database for virtual screening[END_REF] and the docking software used is FlexX. To store the results of docking, MySQL databases are used, but it is still in process. FlexX [START_REF] Rarey | Predicting Receptor-Ligand interactions by an incremental construction algorithm[END_REF]13] is an extremely fast, robust and highly configurable computer program for predicting protein-ligand interactions. During our experiment, after several control tests, standard parameter settings are used except for two cases: "Place particles" and "Maximum overlap volume".

Procedure

The goal of WISDOM II is two fold, the biological goal is to find the best hits against the targets implicated in malaria and the computational goal is to keep improving the relevance of computational grids in drug discovery applications. Here in this paper we are going to discuss in details the grid architecture and deployment.

Virtual screening experimental setup

The complete virtual screening experiment is segmented into five different phases.

i. Target preparation ii. Compound database iii. Validation of the docking experiment iv. Screening v. Result analysis

Target preparation

A standard protocol is used while preparing the target structures. The initial coordinates for all the target structures are obtained from Brookhaven protein database (www.pdb.org). Depending upon the inclusion of the significant residues, cofactors and the binding pocket, active site is defined as 8.0 Å -10.0 Å around the co-crystallized ligand.

Compound database

The Compound library used for WISDOM was obtained from the ZINC database [START_REF] Auvergrid | auvergrid[END_REF][START_REF] Eela | [END_REF]. The ZINC database is a collection of 4.3 million chemical compounds ready for virtual screening from different vendors. We have chosen to use the ZINC library because ZINC is open source database and the structures have already been filtered according to the Lipinski rules. Moreover the data are available in different file formats (Sybyl mol2 format, sdf and smiles). A total of 4.3 million compounds were downloaded from the ZINC database and screened against four targets.

Validation of the docking experiment

Re-docking against the co-crystallized compound is performed to check and tune the docking experiment requirements. Re-docking serves as a control for finally selecting the parameters for target structure, before subjecting it to large scale docking. Docking pose is validated at two levels; RMSD value (the lower, the better) and binding pose of ligand (the more similar the docking pose to the co-crystallized ligand, the better).

Grid infrastructure and Deployment

Grid Infrastructures

The deployments were achieved on several grid infrastructures: Auvergrid [START_REF] Auvergrid | auvergrid[END_REF], EELA [START_REF] Eela | [END_REF], EGEE [START_REF] Gagliardi | Building an infrastructure for scientific Grid computing: status and goals of the EGEE project[END_REF], EUChinaGrid [17] and EUMedGrid [18]. All these infrastructures are actually using the same middleware, gLite. EGEE is the main infrastructure offering the largest resources; they are all interconnected with EGEE, in the sense that all of these Grids share some resources with EGEE. In the case of Auvergrid, it is even more evident as all the resources available through the Auvergrid Virtual Organization (VO) are also shared with several EGEE VOs. The EUChinaGrid project for instance made available all the grid sites belonging to its infrastructure; seven Computing Elements in total and two Storage Elements were used to store the databases and result files on the EuChinaGrid.

WISDOM production Environment

WISDOM environment has been used two times in previous large-scale experiments, WISDOM-I in the summer 2005 [START_REF] Jacq | Virtual Screening on Large Scale Grids[END_REF] and a second deployment against avian flu in the spring 2006 [START_REF] Lee | Grid-Enabled High-Throughput In Silico Screening Against Influenza A Neuraminidase[END_REF]. WISDOM environment keeps evolving in order to make it more user friendly and easier to use by non grid expert. The main objective was also to improve the fault-tolerance of the system, in implementing, for instance, a persistent environment, that can be stopped and restarted at any time without risk of loosing significant information, which proved to be also very useful as it enables the whole maintenance of the scripts and code and improve the interactivity with the user, as the user can also manage jobs finely, for instance force the cancellation and resubmission of a scheduled job. Along with this, we tried to minimize the cost of the environment in terms of disk space and CPU consumption for the user interface. Most of the job files are now generated dynamically: this allows the user as well to modify on the fly the configuration of the resource brokers and the jobs requirements. This way, the user is sure that the next submissions will take these modifications into account. The Figure 1, shows the overall architecture of the environment. The user is interacting with the system through the two main scripts (widom_submit and wisdom_status) deployed on the user interface. These scripts will take care of job files generation, submission, status follow-up and eventually resubmission automatically. The jobs are submitted directly to the grid Workload Management System (WMS), and are executed on the grid computing elements and worker nodes (CEs and WNs). As soon as it is running, a job transfers all the files stored on the Storage Elements (SEs) via the Data Management System of the grid (DMS) with the gridFTP protocol. Once the job is finished, the outputs are stored back on the grid Storage Elements via the Data Management System and the useful docking results are inserted directly from the grid to a relational database where they can later be more easily queried and analyzed.

Data Challenge Deployment

The deployment was performed on the infrastructures listed under section grid infrastructures, and involved at least one manager to oversee the process on each of them. The three groups of targets (GST, Plasmodium vivax and falciparum DHFR) were docked against the whole ZINC database (4.3 millions of compounds). The database was actually cut into 2,422 chunks of 1,800 compounds each. This splitting was chosen because we wanted to have an approximated processing time ranging from 20 to 30 hours for each job (one docking process takes from 40s to 1min depending on the CPU power). The subsets were stored on the grid and replicated on several locations whenever possible to improve fault-tolerance. We define a WISDOM instance as one target structure docked against the whole ZINC database, with a given parameter set.

A total number of 32 instances were deployed, corresponding to an overall workload of 77,504 jobs, and up to 140,000,000 docking operations. On these total 32 instances, 29 instances ran on EGEE, and 3 were on Auvergrid, EELA and EuChinaGrid.

As shown in Figure 1, the environment included a FlexLm server that was providing the floating licenses for the FlexX commercial software. The FlexX software was already used during the first WISDOM deployment in 2005, and the license server was identified as a potential bottleneck and point of failure because we had just one server available at this time. For WISDOM-II, up to 3 servers were made available at the SCAI Fraunhofer institute (http://www.scai.fraunhofer.de), with 3,000 licenses available on each server.

As the average duration of a job was around 30 hours, we submitted 1 instance per day, with a delay of 30 seconds between each submission. As one instance was submitted in about 20 hours, the submission process was quite continuous during the first month of deployment. The jobs were submitted to 15 Resource Brokers (the components of the Workload Management System) in a round-robin order. At the end of a job, the results were stored on the grid Storage Elements and directly into a relational database. The job repartition was quite similar to the previous deployments, but here the UKI federation played an even bigger part. For instance, one of the British sites offered for quite a long period of time more than 1,000 free CPUs, which is half of the average used CPUs. Auvergrid, EELA, Euchinagrid and Eumedgrid contributed by running each 3% of the jobs.

Results and Statistics

Table 1, shows the overall statistics of the deployment. The number of jobs here are the number of awaited results, but far more jobs were actually submitted on the grid. When a job was done on the grid, the environment checked a status file specifying the final result of the job: a job can be done in the point of view of the worker node, without having produced the result files, in this specific case, the status of the job, which was stored on the grid as well, was labeled as failed, and the environment had to resubmit the job. In some cases, the environment failed at retrieving the status from the grid, and thus considered implicitly the job has failed, even if the job has succeeded. It explains why some jobs ran several times, and why the final completed docking number is bigger than the useful awaited dockings. Anyhow the average docking throughput is coherent with the crunching factor. The crunching factor is the ratio of the total CPU time over the duration of the experiment. It represents the average number of CPUs used simultaneously all along the data challenge and is a metric of the parallelization gain. If we consider 80,000 dockings per hour for 2,000, it means 40 dockings for one CPU per hour, which is coherent with the empiric observation of one docking process lasting approximately 1 minute on a 3.06 Intel Xeon processor. In the Table 1, the estimated grid success rate is the ratio of successful grid jobs on the total of submitted jobs. The success rate after output checking will consider just the jobs that succeeded in producing the result files, that's why this score is lower. One can notice that these values are very small, but there are several explanations for this. At the beginning of the data challenge, the observed grid success rate was about 80 to 90%, but it decreased constantly because of sites overload. Sometimes the available disk space was decreasing on some Resources Brokers, up to a point where some of the job data could not reach the Computing Element. In other cases, the sites were simply producing a lot of aborted job for an undetermined reason. The Resource Brokers failed again to balance reasonably the jobs on the Computing Elements, and some of them ended up with more than 500 jobs in queue, the site administrator had no other choice than kill all these jobs, producing in a single row more than 500 aborted jobs. Actually, because of the automatic resubmission, this information should not be taken as an overall significant way to evaluate the efficiency of the grid, because the automatic resubmission guaranteed a successful job, and the aborted jobs are not staying a long time on the grid consuming useful resources. The grid is a very dynamic system, and errors can occur at the last minute.

Observations and Issues

The scheduling efficiency of the grid is still a major issue. The resource broker is still the main bottleneck, and even if used in high number (>15), is always a source of trouble. Moreover things get worse as load is increasing on the grid. The « sink-hole » effects can result in sites overloading in a very short amount of time, and if not taken care quickly can lead to an impressive overhead caused by the long lasting waiting state of the jobs. Added to that the sometimes unreliable and incomplete information provided by the information system, which does not publish the available slots and VO limitations that would be mandatory to perform an efficient scheduling.

Another issue was that to be able to store and treat the data in a relational database, the machine hosting the database must have good performances or the number of queries coming from the grid may also overload the machine significantly. In this deployment we used a MySQL database and planned to put all the produced result in the same table, but finally we had to spit this database in several ones (one per target), because MySQL would not have been able to withstand the total number of records, It was generating CPU overloads on the machine, which lead to serious slowdowns.

All these elements demonstrate clearly that even if the grid can show very good result in comparison to very simple architecture it is still missing robustness and reliability, and can indeed be improved performance-wise.

Conclusions

We have demonstrated the role and significance of computational grids in the life science applications like structure based drug design. Large scale virtual screening on four different targets of malaria was performed in search for potential hits on several grid infrastructures: Auvergrid, EELA, EGEE, EUChinaGrid and EUMedGrid. One of our goals was to further demonstrate the impact of computational grids in life science applications like virtual screening where large amounts of computing power is required; we have achieved it by successfully screening the whole ZINC database for three malaria targets in 76 days instead of 413 CPU years. We have reached during this ten-week period an average docking throughput of 78,400 dockings. MySQL databases are used for the analysis of the docking results, which will ease the final analysis of the virtual screening data. On the biological front 1,738 TB of valuable data has been produced. Analysis of the results (identification of the best hits) is under way: the best hits will be post processed by molecular dynamic simulations and tested in the experimental laboratories.

Figure 1 .

 1 Figure 1. Schema of the WISDOM production environment.

Table 1 .

 1 Overall statistics concerning the deployment.

	Number of Jobs	77,504
	Total Number of completed dockings	156,407,400
	Estimated duration on 1 CPU	413 years
	Duration of the experiment	76 days
	Average throughput	78,400 dockings/hour
	Maximum number of loaded licences (concurrent	5,000
	running jobs)
	Number of used computing elements	98
	Average duration of a job	41 hours
	Average crunching factor	1,986
	Volume of output results	1,738 TB
	Estimated distribution efficiency	39%
	Estimated grid success rate	49%
	Estimated success rate after output checking	37%

Corresponding author: Vinod Kasam: Laboratoire de Physique Corpusculaire, Campus des Cézeaux -24 av des Landais, 63177 Aubière cedex, France; E-mail: kasam@clermont.in2p3.fr.

Acknowledgements

Information on WISDOM collaboration can be found at http://wisdom.healthgrid.org.

The authors acknowledge the contribution of Emmanuel Medernach (AuverGrid), Ignacio Blanquer (EELA), Gabri Aparico (EELA), Enrico Fattibene (EUChinaGrid), Kostas Koumantaros (EUMedGrid) and Domenico Vicinanza (EUMedGrid).

The Enabling Grids for E-sciencE (EGEE) project is co-funded by the European Commission under contract INFSO-RI-031688. Auvergrid is a project funded by the Conseil Regional d'Auvergne. The BioinfoGRID project is co-funded by the European Commission under contract INFSO-RI-026808. The EMBRACE project is funded by the European Commission within its FP6 Programme, under the thematic area "Life sciences, genomics and biotechnology for health", contract number LHSG-CT-2004-512092.

We gratefully acknowledge the support of BioSolveIT who provided more than 6,000 free licenses for their commercial docking program FlexX.


Christos Pappas
email: pappasch@inf.ethz.ch

Raphael M Reischuk
email: reischuk@inf.ethz.ch

Adrian Perrig
email: adrian.perrig@inf.ethz.ch

Forwarding Accountability: A Challenging Necessity of the Future Data Plane

Forwarding accountability mechanisms pinpoint the sending/forwarding properties of traffic to the entities that send and forward the traffic along a path. In this paper, we take flooding attacks as a use case and describe a proposal to hold senders accountable for the sending rates of their flows. Furthermore, we describe the corresponding challenges, potential solutions, and briefly present the literature in the area of forwarding accountability.

Introduction

The Internet started out as a small-scale network among scientists, and turned into a global-scale network for business and private communication alike. At the heart of this success story lies the best-effort delivery service of the network layer -a fundamental property of the Internet architecture. The network does neither provide guaranteed delivery of traffic, nor guaranteed quality of service. The simplicity of this design enabled multiple services and protocols to evolve on top of the minimalistic network layer.

Another property of the Internet architecture is the lack of feedback about the fate of packets. Functionality that detects whether and when packets were delivered is pushed to the end points without aid from the network. However, contrary to best-effort delivery, this design principle has raised trouble with respect to network security. To counter the problems of this design principle, we consider forwarding accountability a necessity for the future data plane.

In general, accountability mechanisms associate state and actions to entities, rendering misbehavior detectable, provable, and non-repudiable. Forwarding accountability in particular, associates the sending/forwarding properties of traffic (e.g., latency, bandwidth) to the entity that sends/forwards the traffic (e.g., a host, a router, a switch, or even an Autonomous System), constructing verifiable evidence about how traffic is sent/forwarded. This verifiable information can then be used by users or legal authorities to make informed decisions in cases of misbehavior or poorly performing networks entities.

In this paper, we highlight the importance of forwarding accountability by describing how it would aid in solving two burning issues for the networking community: network neutrality violations and flooding attacks. Network-neutrality violations. Network neutrality has become an increasingly hot subject in the networking community. Internet service providers (ISPs) have been accused of blocking [START_REF]AT&T Faces Formal FCC Complaint for Blocking Cellular FaceTime Use[END_REF]1] and slowing down traffic from specific content providers [START_REF]Netflix Performance on Verizon and Comcast Has Been Dropping for Months[END_REF].

Consider the dispute between Netflix and Comcast [START_REF] Comcast Vs | Netflix: Is This Really About Net Neutrality?[END_REF]. Netflix -backed by the media and Internet activists -accused Comcast of deliberately slowing down its video traffic, causing an unacceptable quality of experience for the customers. Comcast denied the blame and attributed the problem to the inability of Netflix's direct ISPs to handle the amount of traffic.

An accountable data plane could alleviate many concerns raised by the neutrality debate [START_REF] Pappas | Transparency Instead of Neutrality[END_REF]. As the exposed forwarding information would be trustworthy, the Internet community would obtain much richer feedback on how "neutral" an ISP really is. Verifiable information could be combined with other higher-level information (e.g., Service Level Agreements) to make an informed judgement about ISP's practices.

Flooding attacks. In recent months, we have observed an increase in the frequency and intensity of flooding attacks rooted in misconfigured or vulnerable Internet services: in February 2014, attackers used misconfigured time synchronization servers to attack Cloudflare with a peak of 400 Gbps. For 2015, Akamai reports a 116.5% increase in total DDoS attacks and a 42.8% increase in the average attack duration compared to the previous year [START_REF]Akamai: Q1 State of the Internet -Security Report[END_REF].

In an ideal Internet, users could enjoy the benefits of an accountable forwarding plane. Receivers could specify a traffic profile that sources need to adhere to, drawing a clear line for benign traffic and enabling misbehavior detection. In case of traffic profile violations, receivers could provide proofs of misbehavior to the origin and transit ISPs, and ISPs could ensure compliance for misbehaving hosts through traffic shaping.

We provide an example to demonstrate the virtues of forwarding accountability. Consider the topology depicted in Figure 1 and assume that a web server is located in AS n . We assume that an attacker launches a reflection attack against the server by exploiting the NTP protocol running on vulnerable servers in AS 0 . Specifically, the attacker fakes the victim's source IP address and sends NTP commands to the servers within AS 0 . Due to traffic amplification, the NTP servers generate traffic that overpowers the victim's resources.

With forwarding accountability in place, each packet is associated with a proof that can later remind every AS on the path that it forwarded the traffic. When the web server reports the attack to the ASes on the path by providing the per-packet proofs, the ASes can acknowledge or deny that they forwarded the malicious traffic. Based on the feedback from the ASes, it is possible to detect misbehavior and narrow down its location on the path. In our example, it becomes clear that the NTP servers in AS 0 sourced the malicious traffic. AS 0 can then drop or deprioritize AS 0 's traffic and thus protect the victim web server and the other networks on the path.

Contributions. This paper's focus is on forwarding accountability with respect to flooding attacks. We extend our recent work, FAIR [START_REF] Pappas | FAIR: Forwarding Accountability for Internet Reputability[END_REF], by outlining an ACcountability-based Ddos Protection framework -ACDP -to hold the sending hosts accountable for the sending rates of their originated flows. Furthermore, we describe the corresponding challenges and present the literature in the area of forwarding accountability.

Overview of ACDP

We provide a high-level overview of ACDP. In ACDP, communication proceeds in the following three stages (Figure 1).

-Stage 1 (Setup): Source and destination hosts set up a sending policy that dictates the sending rate for a specific flow between them.

-Stage 2 (Transmission):

The source sends its traffic to the destination.

Each AS on the path (including the source's and destination's ISP) inscribes information in the packet headers, which serves as a reminder to itself that it has forwarded the packets. -Stage 3 (Protest): If the destination host detects a sending-rate violation, it proceeds to the protest phase and hands the sending policy together with the packet headers to its own ISP. The ISP then contacts other ASes on the path by providing the aggregated proof; the proof eventually identifies the adversary.

Setup (Stage 1)

Before sending the actual data, communicating hosts set up a sending policy. The sending policy specifies the sending properties that a host should apply to its outgoing traffic towards the communicating peer. The sending properties can be formally expressed through the Token Bucket [START_REF]Cisco: Cisco Policing and Shaping Overview[END_REF] parameters (the average sending rate, the maximum burst size, and the measurement interval). We consider bidirectional communication channels and thus, both communicating hosts indicate their preferred sending properties. Specifically, the sending policy is constructed as follows:

1. The source1 initiates the policy setup and constructs a policy packet. It inscribes the sending properties that the destination should adhere to when sending traffic to the source. 2. Each AS on the path (including the host's ISP) indicates its presence on the path by inscribing its identifier in the policy packet. However, it does not interfere with the policy details. 3. The destination completes the policy by filling in its own desired sending properties that the source should adhere to. Then, it sends the policy packet back to the source. 4. Similar to step 2., each AS on the path back to the source -possibly different from the outbound path -indicates its presence.

We assume that communicating hosts sign their information with their private key, in order to make the policy non-repudiable. Furthermore, each AS uses a secret key to protect the integrity of its own information, so that it can later remind itself -without keeping state locally -that it witnessed the corresponding policy.

Transmission (Stage 2)

With the sending policy in place, hosts start exchanging traffic under the restrictions of the sending policy. We describe the data-plane operations performed by the source, transit, and destination ASes. These operations are applied to each packet.

Source AS. The border routers of the source AS inscribe the Autonomous System Number (ASN), a timestamp, and a sequence number. The timestamp is included to calculate the sending rate of the source in the next stage. Furthermore, in conjunction with the sequence number, it serves as a protection against packet replay from transit ASes. The ASN points to the AS that forwards a packet and constructs a trace of the AS path together with the ASNs of the next on-path ASes. Finally, the inscribed information is protected with a short MAC in the packet, computed with the secret key of the source AS.

Transit and Destination ASes. Each border router of the transit and destination ASes, performs the following operations:

1. The router verifies that the source's timestamp does not deviate from the local time. If the check fails, the packet gets dropped. This check ensures that the source AS does not collude with a customer host in order to conceal an attack by reporting false timestamps. 2. The router inscribes its own information in the packet: its ASN, a short nonce, and a MAC computed over the inscribed information. The nonce serves as protection against replaying the MAC of the AS.

The destination AS forwards the packet to the eventual recipient, who monitors each flow in order to detect policy violations. In case of a violation, the destination sends the received packets and the policy packet to its own ISP. The ISP proceeds to the protest phase, representing its customer.

Protest (Stage 3)

In the third stage, destination ASes provide proofs of misbehavior to the source AS and the other transit ASes. It is an offline procedure of at most two rounds.

In the first round, the destination AS sends the policy packet and the packet headers it received from its customer to the source AS. The policy packet contains the transmission properties, and the packet headers contain evidence about the actual transmission properties by the source. The source AS examines the evidence (i.e., verifies its own MAC inscribed in the packets) and approves or rejects the complaint. If the source AS approves the complaint, it can take measures against its misbehaving customer. However, a non-cooperating AS or a replay attack from a transit AS can lead to a rejected complaint in the first round.

A rejection in the first round leads to the second round. The destination AS sends the same information to all ASes on the path. They examine the evidence in the same way as the source AS and approve or reject the complaint. Based on the approvals and rejections, the ASes can determine the root of the problem, because each complaint is accepted at least by the benign cooperating ASes adjacent to the destination, as shown in FAIR [START_REF] Pappas | FAIR: Forwarding Accountability for Internet Reputability[END_REF].

Challenges

In this section, we discuss the major challenges related to our proposal. We sketch potential solutions and list the open problems to encourage future research with respect to deployment and performance.

Deployment Challenges

The first challenge is the required modification of hosts. With ACDP, end hosts have to perform additional functionality compared to the legacy communication paradigm (e.g., under TCP or UDP). Namely, end hosts have to perform a policy setup before a connection starts transmitting data. A change in the network stack of the host's operating system is an unrealistic requirement and we believe that this task can be delegated to a gateway between the host and its ISP. Typically, hosts connect to the Internet through their ISP-provided routers, which act as middleboxes and usually perform additional tasks (e.g., acting as a Network Address Translators or a firewall). Requiring middleboxes to interfere and perform the additional functionality keeps the hosts unmodified and provides a smoother deployment path.

A second challenge is the upgrade of the AS infrastructure required in order to inscribe the additional information in the packet headers; specifically, the MAC computation requires a hardware implementation of a cryptographic engine. Although it is impossible to circumvent this requirement, hardware cryptographic engines are readily available for commodity processors [START_REF] Gueron | Intel Advanced Encryption Standard (AES) New Instruction Set[END_REF]. Since the required technology exists at a low price, the required upgrades would not incur a high procurement cost for ISPs.

The third challenge is a viable business model that provides incentives to ISPs to adopt such a mechanism. We anticipate that security-concerned customers (especially enterprise networks) will be interested in buying service from an ISP that handles and forwards its customers' complaints to the sources of misbehavior. Hence, competition would be the key to promote Accountabilityas-a-Service [START_REF] Bender | Accountability as a service[END_REF]. However, a thorough economic analysis is required to explore the viability of such a security service.

Performance Challenges

The additional functionality required for ACDP introduces overhead with respect to processing, latency, and bandwidth.

A border router of an AS has to inscribe additional information, which includes the computation of a MAC. We conducted an experiment on one 10 GbE NIC port of a commodity server machine, simulating the required processing, and found that there was not a substantial drop in throughput. Specifically, for 64 byte packets (the minimum packet size, i.e., the maximum packet rate) the switch forwards at 95% of the line-rate. For 128-byte packets and larger, the switch achieves line-rate performance. The initial results indicate that forwarding performance would not suffer from such an accountability framework.

Another performance issue is the increase in latency for communicating hosts. Before the actual communication starts, end hosts must establish the sending policy, which translates to one Round-Trip-Time (RTT). This overhead can be significant for latency-sensitive applications (e.g., video streaming). More extensive research is required to optimize this aspect, but potential solutions include piggybacking the policy packet on the first data packets and embedding a default policy in the DNS records.

Our proposal comes with an increased packet size that leads to bandwidth overhead. The length increase is inevitable, but certain measures can limit the introduced overhead. For instance, a short MAC (4 bits) per ISP is enough to enable misbehavior detection in the context of flooding attacks [START_REF] Pappas | FAIR: Forwarding Accountability for Internet Reputability[END_REF].

Related Work

We present the main proposals in the area of forwarding accountability to date.

Goldberg et al. [START_REF] Goldberg | Path-quality Monitoring in the Presence of Adversaries[END_REF] propose end-to-end path quality monitoring in the presence of adversaries. Specifically, an alarm is raised when packet loss and delay exceed a given threshold. The proposal leverages secure sampling, which allows end points to coordinate their measurements of loss and delay when an onpath adversary delays or drops packets. An alternative protocol uses a sketch to exchange loss measurements securely and efficiently in adversarial scenarios, accompanied by a theoretical analysis about their accuracy vs. overhead tradeoffs. In addition, these protocols make sensible assumptions for networking environments: no symmetric paths, no processing at forwarding devices, and configurable storage overhead based on accuracy target. However, these protocols do not provide granular performance reports for smaller path segments and do not localize misbehavior.

FAIR [START_REF] Pappas | FAIR: Forwarding Accountability for Internet Reputability[END_REF] is a forwarding accountability mechanism that pushes stricter security policies to ISPs. The source and destination ASes set up a communication channel with a corresponding sending policy, which can specify sending properties (e.g., average sending rate) or forbid abnormal packet headers used for malicious activity (e.g., Christmas tree packets). Transit ASes on the path mark packets and in case of policy violations, the packets are used as a proof of misbehavior. FAIR comes with an implementation that introduces low bandwidth overhead and can switch packets at a line-rate of 120 Gbps. However, the proposal does not allow proving misbehavior at the granularity of flows, and thus cannot be used to identify individual misbehaving flows or hosts.

AudIt [START_REF] Argyraki | Loss and Delay Accountability for the Internet[END_REF] proposes an accountability interface, provided by ISPs, that gives loss and delay feedback to the traffic sources. The framework relies on statistics reports from ISPs, without requiring complicated key establishment. However, the proposal is based on aggregation of flow information, and thus ISPs can hide their lies since they report mean values.

Fig. 1 .

 1 Fig. 1. ACDP Operation.

We refer to the host that initiates the connection as the source; and to its communicating peer as the destination.