Mathieu De Craene

Martino Alessandrini

Pascal Allain

Stephanie Marchesseau

I Waechter-Stehle

J Weese

E Saloux

H G Morales

R Cuingnet

Hervé Delingette

Generation of ultra-realistic synthetic echocardiographic sequences

Keywords: Modeling of image formation, Ultrasound imaging, Motion detection and tracking

de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

GENERATION OF ULTRA-REALISTIC SYNTHETIC ECHOCARDIOGRAPHIC SEQUENCES

M. De Craene * , M. Alessandrini , P. Allain * , S. Marchesseau ∧ , I. Waechter-Stehle ∨ , J. Weese ∨ , E. Saloux • , H. G. Morales * , R. Cuingnet * , H. Delingette ∧ , M. Sermesant ∧ , O. Bernard , J. D'hooge

INTRODUCTION

Quantification of cardiac motion and deformation/strain from 3D cardiac Ultrasound (3D US) is an active area of research. Despite encouraging preliminary studies, a widespread dissemination of these tools in clinical diagnostics is held back because no solid quality assurance methodology exists [START_REF] Jasaityte | Current state of threedimensional myocardial strain estimation using echocardiography[END_REF].

The validation of strain imaging methods has been addressed in several ways [START_REF] Jasaityte | Current state of threedimensional myocardial strain estimation using echocardiography[END_REF]: in vivo in human subjects by using tagged magnetic resonance (tMRI) as reference method; in vivo in open-thorax animal preparations by using sonomicrometry as a reference; in vitro by using a physical cardiac phantoms and sonomicrometry as reference and in silico using fully synthetic sequences [START_REF] Duan | Simulation of 3D US with a realistic electro-mechanical model of the heart[END_REF][START_REF] Butakoff | Simulated 3D US LV cardiac images for active shape model training[END_REF][START_REF] De Craene | 3D strain assessment in US (Straus): A synthetic comparison of five tracking methodologies[END_REF].

In particular this latter solution presents the following advantages: it does not require inter-modal registration (unlike tMRI), and it provides the reference motion in every point of the heart (unlike sonomicrometry). Unfortunately, the level of realism of existing simulation techniques is highly unsatisfactory: surrounding structures such as papillary muscles and heart valves are typically neglected as well as typical artifacts such as reverberations, clutter noise, signal dropout and local intensity variations due to changing cardiac fiber orientation (see Fig. 1).

For this reason, recently two new simulations approaches were proposed in which a clinical ultrasound recording was used as an intensity template for the simulated data set [START_REF] Alessandrini | Simulation of realistic echocardiographic sequences for ground-truth validation of motion estimation[END_REF][START_REF] Prakosa | Generation of synthetic but visually realistic time series of cardiac images combining a biophysical model and clinical images[END_REF]. However, both solutions suffer from intrinsic limitations.

In [START_REF] Alessandrini | Simulation of realistic echocardiographic sequences for ground-truth validation of motion estimation[END_REF] the simulated cardiac motion is the result of a motion estimator applied to the template sequence. As a consequence no control on the simulated motion is possible. Also, errors of the adopted estimator reflect directly on the simulated motion. Finally, a bias can be expected in the comparison of different algorithms towards estimators of the same class as the one used to determine the synthetic motion.

On the other side, in [START_REF] Prakosa | Generation of synthetic but visually realistic time series of cardiac images combining a biophysical model and clinical images[END_REF] the synthetic sequence is obtained by warping directly the grayscale images of the template acquisition according to the desired motion field and this unavoidably introduces unrealistic warping artifacts in the resulting speckle texture.

In this paper we present a framework for the generation of 2D or 3D ultra-realistic synthetic echocardiographic sequences overcoming the limitations of existing approaches. an average shape model is deformed to match intensity profiles in the image with a training database [START_REF] Ecabert | Automatic model-based segmentation of the heart in CT images[END_REF]; (b) from the resulting surface mesh, a volumetric tetrahedral mesh is generated, (c) a sequence of meshes is generated through the SOFA framework [START_REF] Marchesseau | Preliminary specificity study of the Bestel-Clement-Sorine electromechanical model of the heart using parameter calibration from medical images[END_REF], here showing only the overlay between end-diastolic (gray) and end-systolic meshes (red).

cal model (E/M) [START_REF] Marchesseau | Preliminary specificity study of the Bestel-Clement-Sorine electromechanical model of the heart using parameter calibration from medical images[END_REF]. The E/M model provides the reference for heart motion and derived quantities and can generate a continuous range of pathological patterns by progressively altering the value of contraction and/or conduction parameters. As in [START_REF] Alessandrini | Simulation of realistic echocardiographic sequences for ground-truth validation of motion estimation[END_REF] a real clinical recording is adopted as template for local tissue echogenicity, which allows simulating surrounding structures as well as US artifacts. Finally, an US simulation environment is used to model the image formation [START_REF] Gao | A fast convolution-based methodology to simulate 2-D/3-D cardiac US images[END_REF], allowing full control over the transducer's parameters and ensuring realistic speckle patterns.

PROPOSED SIMULATION FRAMEWORK

In our framework, the generation of the synthetic ultrasound image sequence requires to define at every frame a set of discrete scattering sites and their respective echogenicity values. The E/M simulation affects the coordinates of scatterers inside the myocardium. Echogenicity values are obtained from the intensity of a real clinical recording. Integrating the mechanical simulation result with an ultrasound image sequence was performed in this paper by combining two elements. First, a patient-specific volumetric mesh is constructed by segmenting the first image frame (see Sect. 2.1). All E/M simulations are then performed on this geometry. Second, the echogenicity of scatterers moving in the simulation space is computed from the image intensities. This requires to compute a set of transformations that relate a coordinate from the simulation space to the image (see Sect. 2.2). The combination of this set of transformations with the mechanical simulation result is described in Sect. 2.3.

4D image segmentation, mesh generation and mechanical simulation

Left and right ventricular geometries are first extracted from the first frame of a template 3D US image sequence (i.e. a real clinical recording). Hereto, the four chambers were automatically extracted from the 3D image using a model-based segmentation technique initially developed for CT images [START_REF] Ecabert | Automatic model-based segmentation of the heart in CT images[END_REF] and recently extended to ultrasound. From the segmentation result on the first image frame, two outputs are further produced. The first output is a dynamic sequence of LV endo-and epicardial segmentations produced by tracking and manual editing. The editing is done in 4D to ensure a smooth propagation of manual updates. The result is a set of landmarks over time sampling the LV endoand epicardial surfaces regularly in the circumferential and longitudinal directions. This processing is part of the QLab software (3DQ Advanced plugin) [START_REF] Jacobs | Rapid online quantification of left ventricular volume from real-time threedimensional echocardiographic data[END_REF].

The second output is a volumetric tetrahedral mesh obtained by applying the NetGen (http://www.hpfem. jku.at/netgen/) software to the surface mesh. Prior to mesh generation, the right ventricle was first extruded along the mesh normals to produce a RV epicardial surface. We then applied the simulation technique described in [START_REF] Marchesseau | Preliminary specificity study of the Bestel-Clement-Sorine electromechanical model of the heart using parameter calibration from medical images[END_REF] and implemented in the SOFA framework (http: //www.sofa-framework.org/) for generating two simulations of healthy and ischemic cases. The ischemic case is produced by altering the contractility of segments perfused by the left circumflex artery [START_REF] De Craene | 3D strain assessment in US (Straus): A synthetic comparison of five tracking methodologies[END_REF].

Spatial and temporal matching between image and simulation sequences

The image and simulation sequences have different dynamics that have to be matched. A transformation between image and simulation times is computed by linear interpolation between temporal landmarks defined on both sequences. For each simulation frame indexed by t k simu , we compute t k imag being the interpolated time in the real sequence. The closest image frames are defined as t k+ imag = t k imag and t k- imag = t k imag . From t k+ imag and t k- imag , two motion fields can then be de-rived for matching coordinates from the first image frame to these two times. These two mappings will be denoted as T imag (•, t k- imag ) and T imag (•, t k+ imag ) where the first argument is a spatial coordinate in the domain of the first image frame. In our implementation, these mappings were computed from a Thin Plate Spline (TPS) transform parametrized by the 4D segmentation landmarks (Sect. 2.1).

In the simulation space, a set of transformations mapping each frame to the first frame are computed similarly as in the image space. They are noted as T simu (•, t k simu ) where the first argument is a spatial coordinate at frame t k simu . All mappings and coordinates notations are summarized in Fig. 3.

Scatter map generation

The first step is to generate a scatter map on the first frame. The scatter map contains the coordinates and echogenicity of each scatterer. As the E/M mesh matches the first image frame, the mesh can be used to differentiate foreground (within myocardium) from background scatterers. The scatterers coordinates on the first simulation frame are obtained by uniformly sampling N points through the image domain. Parametric E/M mesh coordinates of foreground scatterers are stored for further use. Their echogenechogenicity is obtained by interpolating the intensity of the first frame of the template image sequence.

For generating the scatter map at time t k simu > 0, we first move the N M scatterers within the myocardium by computing the absolute coordinates from the parametric coordinates computed on the first frame. We then redraw N B = N -N M samples using an uniform distribution through the image domain, excluding the myocardial volume.

For each scatterer at coordinate x simu and time t k simu , we then define the position x 0 of that scatterer at time 0. For foreground scatterers, x 0 is directly available since they are indexed consistently over time. For background scatterers, it can be computed as

x 0 = T simu (x simu , t k simu ) . (1)
The echogenicity of every scatterer can then be defined as

a =      (1 -w)I t k- imag • T imag (x 0 , t k- imag ) +w • I t k+ imag • T imag (x 0 , t k+ imag ) if x 0 / ∈ Ω M I 0 (x 0 ) if x 0 ∈ Ω M (2
) where t k+ imag and t k- imag are defined in Sect. 2.2, Ω M stands for the myocardial domain at time 0 and w is an interpolation weight defined as

w = t k imag -t k- imag . (3)
Amplitudes computed from Eq. 2 are further calibrated using where a M AX is the maximum a value, K = 70 is a constant and the power operation compensates the log compression traditionally performed in the ultrasound device prior to display. The normalized scatter map is set as the input of a fast ultrasound simulation [START_REF] Gao | A fast convolution-based methodology to simulate 2-D/3-D cardiac US images[END_REF] that computes radio-frequency signals on each scan line, and further performs envelope detection and scan conversion.

a = 10 K 20 ( a a M AX -1) (4)

RESULTS

We present two synthetic sequences obtained with the proposed framework corresponding to one healthy and one ischemic motion pattern. A standard 3D US sequence was acquired from a healthy volunteer using a iE33 (Philips Medical Systems, Bothell WA, USA) echocardiography system. The first frame of the image sequence was considered to be the end of diastole while the end of systole was visually determined based on the aortic valve closure. The synthetic motion was obtained from the mechanical simulation. For the ischemic simulation, peak contraction and stiffness values were modified in the territory of the left circumflex artery. To mimic the rigidity of the mitral annulus, basal segments were associated to high stiffness values. The E/M model provides the benchmark strain curves that can be used to evaluate 3D speckle tracking algorithms. The resulting radial strain curves obtained for the two cases are plotted in Fig. 4. Basal strain curves (in dashed) have lower values in both the ischemic and the healthy cases. Ischemic regions showed reduced radial strain values, they are underlined in the Bull's eye and plotted in bold. Taking as input the image sequence and the sequence of E/M meshes, a scatter map was generated as described in Sect. 2.3. After simulation of the ultrasound acquisition [START_REF] Gao | A fast convolution-based methodology to simulate 2-D/3-D cardiac US images[END_REF], the resulting synthetic image sequence is plotted in Fig. 5, with the E/M meshes overlaid on the image and two links to visualize the sequences of images and meshes. Generated images and meshes are also attached to this submission. The realism of this simulated data has to be compared with state of the art strategies, as plotted in Fig. 1. The inclusion of a real clinical recording in the pipeline adds realistic challenges: presence of pericardium, valves, etc. The speckle texture inside the myocardium is realistic and can be adapted to several probe designs [START_REF] Gao | A fast convolution-based methodology to simulate 2-D/3-D cardiac US images[END_REF]. Motion within the myocardium is unbiased to any image registration algorithm. Note that despite the fact that two cases have been illustrated here, several intermediate cases can be generated to study the sensitivity of any image-based strain quantification technique under a controlled and reproducible environment.

CONCLUSION

We proposed a novel framework for the generation of highly realistic 3D ultrasound images. The realism of the proposed data represents an important shift with respect to state of the art strategies. First exemplary datasets are available for download. They include all generated meshes, providing ground truth landmarks for validating tracking algorithms, and synthetic images. The generated sequences and meshes can be downloaded from http://bit.ly/synthusdata.

As future work, we will extend the database to span over a range of location and contraction values. Other pathological shapes (hypertrophic, dilated) will be included by selecting representative 3D ultrasound sequences. This extended data will serve as basis for the comparative and quantitative evaluation of different strain quantification algorithms.

Fig. 1 .

 1 Fig. 1. Synthetic echocardiographic images used in [2](a), [3](b), [4](c) and using the framework proposed in this paper (d).

Fig. 2 .

 2 Fig. 2. From a 3D US image to a mechanical simulation result: (a)

Fig. 3 .

 3 Fig. 3. Spatiotemporal mappings used to estimate the echogenicity of a background scatterer, from simulation to image spaces.

Fig. 4 .

 4 Fig. 4. Radial strain curves for the healthy and ischemic simulations. See the Bull's eye diagram for color legend. Dashed curves indicate basal segments. Non-basal ischemic segments are plotted in bold and underlined in the Bull's eye diagram.

Fig. 5 .

 5 Fig. 5. Simulated sequence (top row: long axis view and bottom row: short axis view) corresponding to a healthy case. The 'benchmark' E/M mesh is superimposed over the simulated image. A video showing the entire simulated sequence is available at http://bit.ly/sim_us. An animation comparing ischemic and healthy simulations can be downloaded from http://bit.ly/sim_us2.


Gaël Deest

Tomofumi Yuki

Olivier Sentieys

Steven Derrien

Toward Scalable Source Level Accuracy Analysis for Floating-point to Fixed-point Conversion

In embedded systems, many numerical algorithms are implemented with fixed-point arithmetic to meet area cost and power constraints. Fixed-point encoding decisions can significantly affect cost and performance. To evaluate their impact on accuracy, designers resort to simulations. Their high running-time prevents thorough exploration of the design-space. To address this issue, analytical modeling techniques have been proposed, but their applicability is limited by scalability issues. In this paper, we extend these techniques to a larger class of programs. We use polyhedral methods to extract a more compact, graph-based representation of the program. We validate our approach with a several image and signal processing algorithms.

I. INTRODUCTION

Custom embedded hardware platforms require the programmer (designer) to address many challenges that do not appear in general purpose setting. One reason is that such implementations must satisfy tight design constraints such as performance, cost, energy, reliability, and so on. The importance of tools for efficient exploration of the design space is further emphasized by the market pressure towards shorter design cycles.

Most custom embedded hardware platforms are designed to provide high computing power for a low area/energy budget. Hardware support for floating-point operations is rarely provided, necessitating the use of fixed-point arithmetic to achieve performance. However, most digital signal and image processing specifications (in Matlab or C/C++) use floating-point. The implementation stage requires a costly (experiments [START_REF] Clark | Accelerating Fixed-Point Design for MB-OFDM UWB Systems[END_REF] report up to 50% of design effort) floating-point to fixed-point conversion step, performed manually.

During the conversion process, designers try to derive implementations that maximize performance while retaining enough accuracy to preserve the functionality of the system. For most algorithms, the impact of encoding decisions on the overall system behavior is difficult to predict, and the performance/accuracy exploration therefore relies on simulations to evaluate encoding decisions. This simulation time is often a bottleneck, limiting the exploration of the solution space, leading to sub-optimal implementations.

To address this issue, analytical accuracy models have been proposed [START_REF] Constantinides | The Multiple Wordlength Paradigm[END_REF], [START_REF] Ménard | Analytical Fixed-Point Accuracy Evaluation in Linear Time-Invariant Systems[END_REF], [START_REF] Caffarena | SQNR Estimation of Fixed-Point DSP Algorithms[END_REF]. These methods focus on signal processing algorithms (e.g., FIR, IIR, FFT) that are modeled using This work is co-funded by the European Union under the 7th Framework Programme under Grant Agreement ICT-287733.

Signal Flow Graphs (SFGs). These approaches fit well with model-based design flows (e.g., graphical dataflow languages) but are difficult to adapt in a compiler context, where the system model is implicitly specified as program source code.

Existing tools solve the problem by building an SFG model out of the program [START_REF] Ménard | Analytical Fixed-Point Accuracy Evaluation in Linear Time-Invariant Systems[END_REF], [START_REF] Rocher | Analytical Approach for Numerical Accuracy Estimation of Fixed-Point Systems Based on Smooth Operations[END_REF], but are currently restricted to kernels operating on relatively small, one-dimensional inputs. The goal of this work is to extend the applicability of automatic analytical modeling techniques to a wider class of programs (in particular image processing algorithms). We focus on how to systematically build such models from source code in a more efficient and more scalable manner, rather than improving the quality (i.e., accuracy of the model) with respect to prior work.

The key idea in this work is the use of convolutions as a high-level operator to characterize input programs. Linear Time-Invariant (LTI) systems (and its generalization to higher dimensions) are convolutions by constants. Accuracy models for LTI systems rely on characterizing the system either with impulse response or transfer functions, which is significantly simplified by representing the program as combinations of convolutions. Specifically, our contributions are:

• a method to automatically build abstract models of the input program consisting of convolutions, employing a series of techniques from compilers, • use of the abstract view of the program to model noise propagation, as well as the effect of quantization decisions, and • prototype implementation and empirical validate our approach with several signal and image processing kernels.

The rest of this paper is organized as follows. In Section II, we present the necessary background on fixed-point accuracy analysis. We present the main limitations of existing approaches and highlight our contributions in Section III. Section IV introduces the Convolution System as a compact representation of the input program, which we use in Section V to build an analytical accuracy model of accuracy. We evaluate our approach with signal and image processing applications in Section VI. We discuss related work in Section VII, and conclude in Section VIII.

II. BACKGROUND AND MOTIVATION

In this section, we introduce the problem of floating-point to fixed-point conversion, its impact on system accuracy, and the optimization problem it raises. We then give the necessary background on Linear-Shift Invariant (LSI) systems, the kind of systems targeted by our approach.

A. Floating-Point to Fixed-Point Conversion

In embedded systems, most multimedia and signal applications are implemented using fixed-point arithmetic instead of floating-point. In fixed-point arithmetic, real numbers are represented using scaled integers with a fixed number of bits used for the integer part and for the fractional part. We write Q m,n for the fixed-point format with m bits for the integer part (exclusive of the sign bit) and n bits for the fractional part. Such a format has [-2 -m , 2 -m -2 -n ] as its dynamic range, with 2 -n as a quantization step.

When converting an initial program from a floating-point specification to a fixed-point implementation, one must choose a Q m,n format for every variable in the program. This choice must be done carefully as it can significantly alter the behavior of the system. Underestimating m may lead to overflows, i.e., operation output falls out of the dynamic range, which have high impact on system behavior. Underestimating n may lead to significant losses of precision due to quantization errors.

To handle these two types of errors, floating-to fixed-point tools operate in two stages:

1) Estimation of the dynamic range (one-time process).

2) Design space exploration (iterative process) that minimizes cost (energy, area) while satisfying execution time and accuracy constraints. The constraints on accuracy are expressed as Signal to Quantization Noise Ratio.

During the design space exploration, a large number of encoding combinations are considered before reaching a satisfying solution. Rapid evaluation of the effect of encoding decisions on system SQNR is therefore critical. For example, when targeting application-specific hardware, the designer has complete control over the operator and datapath wordlength. Using shorter fixed-point formats usually translates into both additional performance and lower area. The same holds (yet not as significantly) for embedded processors with sub-word SIMD execution units or non-uniform register sizes. However, using shorter encodings impacts the overall system accuracy leading to a complex trade-off between accuracy (i.e., application quality degradation) cost, and performance, which is part of the design space exploration. During this exploration, designers try to obtain the best implementation while constraining the output noise to stay below some threshold. One category of approaches for evaluating the SQNR is based on simulations (e.g., [START_REF] Belanovic | Automated Floating-Point to Fixed-point Conversion with the fixify Environment[END_REF], [START_REF] Kim | Finite wordlength effects analysis and wordlength optimization of a multiplier-adder based 8× 8 2D-IDCT architecture[END_REF]). The program is simulated a number of times with different inputs and the outputs from these runs are combined to evaluate the system accuracy. The main drawback is obviously simulation time: evaluating the impact on accuracy for one set of encoding decisions requires the program to be executed many times, slowing down the entire design space exploration. Another set of approaches analytically model the noise power as a mathematical expression depending on data wordlength to avoid long simulation time [START_REF] Constantinides | Wordlength optimization for linear digital signal processing[END_REF], [START_REF] Ménard | Analytical Fixed-Point Accuracy Evaluation in Linear Time-Invariant Systems[END_REF], [START_REF] Rocher | Analytical Approach for Numerical Accuracy Estimation of Fixed-Point Systems Based on Smooth Operations[END_REF]. These approaches are discussed in Section VII.

B. Error Propagation in Linear Shift-Invariant Systems

Our approach, like most others, relies critically on the system being Linear and Shift-Invariant (LSI). The main difference lies in our ability to efficiently handle complex multidimensional systems such as 2D recursive filters, while most approaches are, in effect, limited to one-dimensional Linear Time-Invariant systems.

LSI systems are a convenient framework to study error propagation. In such a system, signals and errors do not interfere and may be considered independently. Let T be a linear multidimensional system, X an input signal and X = X + ε be a perturbated input. Then:

T ( X) -T (X) = T (ε).
In other words, the propagated error is simply the output of the system when applied to the error itself.

A linear system is said Linear and Shift-Invariant if a shift at its input translates in a shift at its output. More formally:

∀ u, T (τ u (X)) = τ u (T (X))
where τ u represents a translation by vector u. LSI systems can be fully characterized by their impulse response h, the response of the system to a unit impulse at the origin:

T (X) = h * X,
where * denotes the convolution operation. This work relies on a compact representation of a system, based on convolutions, to efficiently estimate its impulse response. This response is then used to propagate the statistical properties of errors and compute output noise power.

III. CONTRIBUTIONS

In this paper, we present a technique for automatic evaluation of accuracy by program source code analysis. The key distinction from prior work in this area is that our method operates on a higher level of abstraction-convolutions-by explicitly modeling them as part of the program representation. The use of higher level abstraction give the following benefits:

• Improved scalability. The program representation becomes much more compact making the entire process much more scalable. Specifically, there is no need to flatten the control flow (unrolling of loops, etc.), which greatly reduces the number of nodes in our internal representation.

• Broadened applicability. Convolutions allow us to extend the applicability beyond LTI systems to its generalization to higher dimensional linear systems.

Combined with scalability improvement, our approach can handle image processing applications.

As a simple example, consider an implementation of block FIR filter:

for (n=0; n<NBLOCKS; n++) y[n] = a*x[n] + b*x[n-1] + c*x[n-2];
Existing approaches for automatic accuracy evaluation represent the program with a low-level graph-based representation where each node models an operation [START_REF] Ménard | Analytical Fixed-Point Accuracy Evaluation in Linear Time-Invariant Systems[END_REF], [START_REF] Rocher | Analytical Approach for Numerical Accuracy Estimation of Fixed-Point Systems Based on Smooth Operations[END_REF]. In order to construct such representation, the entire control structure in the source program must be flattened. For this example, the for loop must be completely unrolled, which requires that the number of blocks (NBLOCKS) is a known constant. The number of nodes in the graph corresponds to the number of operations in the program; NBLOCKS × 5 for this example.

In our work, the above program is represented as a single one-dimensional convolution over x using [a b c] as kernel coefficients. Such a high-level representation of the program is sufficient for our purposes, since the most important property of the system, i.e., impulse response, is captured.

The low-level representation quickly becomes infeasible for image processing applications as they have far more operations. For instance, Gaussian blur filters for N × M images have O(NM) number of operations. However, it can be represented as a single two-dimensional convolution when modeled with convolutions. While many computations may not be "convolutions" from a strict, mathematical, point of view, any weighted sum can be viewed as convolutions. Perceiving them as convolutions enables us to analyze many signal/image processing applications in a unified manner. This is the key insight behind this work.

Note that we do not address the closely related problem of dynamic range estimation in this work. It is important to emphasize that the dynamic range needs to be determined only once and is not part of the design space exploration. Thus, its execution time is not as critical as accuracy evaluation, and methods based on profiling (e.g., Chen and Singh [START_REF] Chen | Profile-Guided Floating-to Fixed-point Conversion for Hybrid FPGA-processor Applications[END_REF]) can be used. Nonetheless, estimation of dynamic range based on convolutions is an interesting direction of future work.

IV. CONVOLUTION SYSTEMS

In this section, we present how we model input programs as convolutions. We first use dependence analysis techniques from polyhedral compilation to represent programs as Systems of Recurrence Equations, and then identify convolutions from the equational view of the program.

A. Systems of Recurrence Equations

Systems of Recurrence Equations (SREs) is a formalism that has been extensively studied in the context of systolic array synthesis [START_REF] Quinton | Automatic Synthesis of Systolic Arrays from Uniform Recurrent Equations[END_REF], [START_REF] Rajopadhye | On Synthesizing Systolic Arrays from Recurrence Equations with Linear Dependencies[END_REF].

Given a set of functions a 0 , a 1 , a 2 , ..., indexed by integer vectors, a system of recurrence equations is a set of equations of the form:

a m (z) = expr (a i1 ( f i1 (z)) , a i2 ( f i2 (z)) , ...)
that define the functions a i . Functions f i j map integer vectors to other integer vectors, and are restricted to affine functions. In addition, the functions a x are associated with a validity domain, characterized by an integer polyhedron.

B. Extracting SREs from Loop Programs

The key link between SREs and programs is the Array Dataflow Analysis [START_REF] Feautrier | Dataflow Analysis of Array and Scalar references[END_REF]. ADA is a powerful dependence analysis that gives exact dependence information. It captures instance-wise and element-wise dependence, and hence is said to be exact. Instance-wise means that each iteration of the loop, or each instance of an operation, is distinguished from the others. Element-wise means that which element of the array is being accessed are captured in the analysis.

As a simple example, consider the following: for (i=0; i <N; i++) S0: x[i] = 0; for (j = 0; j <M; j++) S1:

x

[i] += a[j];
Without instance-wise information, the analysis result only says S1 depends on S0 and S1. Without element-wise information, an instance of S1 at i, j depends on all instances of S0, and S1, since they both write to the variable x.

With array dataflow analysis, the dependence information found for the read x[i] in S1 is the following:

Producer of x[i] at S1 i, j = j = 0 : S0 i j > 0 : S1 i, j -1
The result of array dataflow analysis can also be viewed as a system of recurrence equations [START_REF] Feautrier | Dataflow Analysis of Array and Scalar references[END_REF]. Intuitively, each statement can be considered to store its result in its own array, indexed by its instance vector. With exact dependence information, every read can be mapped to the appropriate element of this single assignment array, effectively constructing a functional specification. The above example can be re-written as the following SRE:

S0(i) = 0 S1(i, j) = j = 0 : S0(i) + a( j) j > 0 : S1(i, j -1) + a( j)
The domains of S0 and S1 are constructed from loop bounds. When the loop bounds are affine expressions, the domain of a statement is always a polyhedron. If statements with affine guards simply add additional constraints to the polyhedron. For example, S0 in the following code:

for (i=0; i <N; i++) if (i>c) S0 has domain D S0 = {i|0 ≤ i < N ∧ i > c} = {i|c < i < N}.

C. The Convolution Operator

We define the convolution operator ( * ) as the following:

y = x * c ≡ y( i) = ∑ k∈D c c( k)x( i + k)
where c and x are equations in the SRE that respectively define the convolution coefficients and the input dataset. For a d-dimensional convolution, the equations, c, x, and y are all defined over d-dimensional domain indexed by an integer vector i. The domain of c, D c defines the window where the convolution is applied, and c is expected to be constants. The equation x should refer to the input dataset or outputs from convolution, and hence transitively to input.

We call an SRE a convolution system if it consists entirely of convolutions, and its combinations by additions. Figure 1 illustrates its graph representation.

D. Identifying Convolutions

Identifying convolutions is an instance of the algorithm recognition problem, which is known to be undecidable in general. Hence, we do not claim to have a general solution, but our approach can find many variations of convolutions found in signal/image processing applications.

Starting from an SRE, convolutions are identified by the following three steps: 1) Preprocessing.

2) Alignment of domains.

3) Extraction of kernel coefficients.

We illustrate the identification process with the Deriche filter in the following. Figure 2 shows fragments of code implementing the Deriche filter, and corresponding SREs.

Preprocessing: This step consists of inlining simple copy equations, distributing multiplications, and merging sequence of additions as a single n-ary addition (summation). Since the convolutions we target are those by constant coefficients, they can be viewed as weighted sums. This is why we first extract summations as a preprocessing step. The summations are then reorganized into those that share the same data variable.

Constant coefficients are either given directly as constants or as references to scalars and/or arrays. We assume that the names of constant coefficient scalars/arrays (such as a1 through a8 in Figure 2) are known.

Applying the preprocessing to the equations in Figure 2 gives the following:

S1(i, j) = sum (a1 × in[i][ j], a2 × in[i][ j -1]) + sum(b1 × S1(i, j -1), b2 × S1(i, j -2)) (1) S5(i, j) = S1(i, j) + y2(i, j) S6( j, i) = sum (a7 × S5(i + 1, j), a8 × S5(i + 2, j)) + sum(b1 × S6(i + 1, j), b2 × S6(i + 2, j)) S11(i, j) = y1(i, j) + S6( j, i)
where the summations are convolution candidates.

Alignment: The aim of the alignment step is to "normalize" the domains of equations in the SRE. In this implementation of the Deriche filter, the horizontal passes visit columns of the image in the inner loop, where the vertical passes visit rows of the image (after horizontal filter) in the inner loop. This can be viewed as applying horizontal pass on the transpose of the image. In fact, the code for horizontal and vertical passes are almost identical due to this transposed view of the image. Thus, it is important to align the domains of equations in the SRE, or the "view" of the image, to appropriately identify the kernel coefficients, including its orientation, of the convolutions.

The alignment problem in our context is closely related to that in the context of uniformization/localization of dependences [START_REF] Chen | On Uniformization of Affine Dependence Algorithms[END_REF], [START_REF] Van Dongen | Uniformization of Linear Recurrence Equations: a Step Toward the Automatic Synthesis of Systolic Arrays[END_REF]). Uniformization is a process of transforming the equations such that all accesses in the definition of the equations are within some constant distance. The definitions of S6 and S11 does not satisfy this property (e.g., S11(i, j) uses S6( j, i), which is a sign of transposition). We use an adaptation of the alignment algorithms used for uniformization to align the equations. We obtain the following after alignment for S6 and S11 in our running example:

S6(i, j) = sum(a7 × S5(i + 1, j), a8 × S5(i + 2, j)) + sum(b1 × S6(i + 1, j), b2 × S6(i + 2, j)) S11(i, j) = y1(i, j) + S6(i, j)
Note that all references in the right hand side of the equations are i, j plus some constant offset. For the purpose of convolution detection, it is sufficient to have all dependences share the same linear part, and not necessarily uniform.

Kernel Extraction: Once the equations are aligned, we are ready to extract the kernel coefficients that characterize the convolution. The domain of the kernel, i.e., window of the convolution, is defined by the access functions to the data variables. As an example, a summation The equations for S1 and S6 are identified as a combination of four convolutions:

y(n) = sum (a × x(n), b × x(n -2), c × x(n -3)),
S1(i, j) = (in * c1) + (S1 * c2) S6(i, j) = (S5 * c3) + (S6 * c4)
where the coefficients c1 through c4 are as follows:

c1(i, j) = [a1 a2] , D c1 = {i, j| -1 ≤ i ≤ 0 ∧ j = 0} c2(i, j) = [b1 b2] , D c2 = {i, j| -2 ≤ i ≤ -1 ∧ j = 0} c3(i, j) = a8 a7 , D c3 = {i, j|i = 0 ∧ 1 ≤ j ≤ 2} c4(i, j) = b2 b1 , D c4 = {i, j|i = 0 ∧ 1 ≤ j ≤ 2}
Note that the domain of convolution is two-dimensional, even though it is effectively one-dimensional. This is necessary to distinguish horizontal and vertical passes, and the domains of coefficients always have the same dimensionality as the data being convolved in our representation of convolutions.

The constructed convolution systems have corresponding graph representations, as illustrated in Figure 3.

[i][j] = y1[i][j] + y2[i][j];
(b) Code for vertical passes Fig. 3: Graph representation of the convolution system for Deriche fragment using compact representation. Unlabeled edges can be viewed as "convolutions" by identity coefficient.

S1(i, j) = a1 × in[i][ j] + a2 × S2(i, j -1) + b1 × S4(i, j -1) + b2× S3(i, j -1) S2(i, j) = in[i][ j] S3(i, j) = S4(i, j -1) S4(i, j) = S1(i, j) S5(i, j) = S1(i, j) + y2(i, j) (c) SRE corresponding to horizontal passes S6( j, i) = a7 × S8( j, i + 1) + a8× S7( j, i + 1) + b1 × S10( j, i + 1) + b2× S9( j, i + 1) S7( j, i) = S8( j, i + 1) S8( j, i) = S5(i, j) S9( j, i) = S10( j, i + 1) S10( j, i) = S6( j, i) S11(i, j) = y1(i, j) + S6( j, i) (d) SRE corresponding to vertical passes
The nodes y1 and y2 would be identical to nodes S1 and S6 if the full code were analyzed.

E. Scope of the Analysis

The powerful dependence analysis is made possible by focusing on a subset of programs where the control flow is static. Informally, the scope of the analysis is defined to be loop nests where all array accesses, loop bounds, and if guards are (quasi-)affine expressions of the surrounding loop indices. Quasi-affine expressions are affine expressions with modulo by constants, which allow us to express conditionals such as: if (i%2==0). Moreover, the analysis handles sequences of imperfectly nested loops.

In addition, it is not necessarily for the full application to satisfy the above restriction. For instance, house keeping code that does not affect the main computation can be excluded from the analysis. Such code regions can be identified by less accurate but more general dependence analysis techniques (e.g., simple name-based dependences, and/or dependence vectors).

Recall that the boundary conditions have been ignored in this section. Although it is possible to construct convolution systems with all the boundary cases, an approximation of the system without the boundaries are used in the subsequent sections as well. The main branches of SREs can usually be identified by inspecting the conditions of each branch, since boundary conditions have less "effective" dimensions due to equalities. By constructing systems that consists of these branches, the resulting system becomes an approximation the underlying system that focus on the steady state.

V. ACCURACY MODEL CONSTRUCTION

This section describes how the convolution system introduced in the previous section is transformed to an accuracy model of the program. This model, parameterized by quantization decisions, is well suited to automatic wordlength optimization. In this section, we use a subset of the convolution system developed in Section IV to illustrate our approach. The sub-graph corresponding to the first horizontal pass is reconstructed in Figure 4.

in(i, j) S1 y [(a 1 , 0 0 ); (a 2 , 0 -1 )] [(b 1 , 0 -1 ); (b 2 , 0 -2 )]
id Fig. 4: Convolution system of the horizontal stage (right) of the Deriche filter. Edges are labeled with a compact representation of the convolution kernel. This graph corresponds to the computation y(i, j) = a 1 in(i, j) + a 2 in(i, j -1) + b 1 y(i, j -1) + b 1 y(i, j -2) in Equation 1.

e in S1 e y e + [(a 1 , 0 0 ); (a 2 , 0 -1 )]
id

[(b 1 , 0 -1 ); (b 2 , 0 -2 )]
id Fig. 5: Transformed system with noise sources. Input is replaced by e in , the quantization error of in. Node e y represents the sum of quantization errors in the computation of y.

A. Adding Noise Sources

The first step is to modify the system to model the computation at the quantization error level. This process is illustrated in Figure 5. Inputs are replaced by Quantization Noise Sources (QNSs) to model quantizations of inputs. A QNS is also added to every internal and output node to model the errors in the computation represented by this node. Note that these noise sources may correspond to several quantizations, as quantizations may occur after every arithmetic operation.

B. Modeling Contributions of Noise Sources to Outputs

In a linear system, the contribution of each noise source to an output may be considered independently. For every pair of noise source s and output o, a sub-system (sub-graph) of the system is constructed to independently model the contribution of s to o. The set of sub-systems modeling the contributions to o is denoted as Err o .

These subsystems are defined as follows. Let G be the graph representing the contribution of a noise source s to o. A node v belongs to G if there is a path from s to v and from v to o. An edge

v 1 → v 2 belongs to G if v 1 , v 2 ∈ G.
This process is illustrated in Figure 6.

C. Propagating Noise Characteristics

To estimate output accuracy, the statistical properties of the noises are propagated to the outputs. Let µ s and σ 2 s be, respectively, the mean and variance of the noise of some noise source s. Instances of this noise are supposed independent and identically distributed. Let h s,o be the (possibly multidimensional) impulse response of the system from s to o. The mean and variance of the propagated noise are:

µ s,o = µ s ∑ v h s,o ( v) σ 2 s,o = σ 2 s ∑ v h 2 s,o ( v), (2)
where vector v spans the space of coordinates. id

[(b 1 , 0 -1 ); (b 2 , 0 -2 )]
id Fig. 6: Decomposition of the horizontal Deriche filter graph. Each rounded box defines a subgraph which captures the propagation of a noise source to output y.

To enable rapid computation of µ s,o and σ 2 s,o , the above sums must be pre-computed by evaluating the impulse response on a large window around the origin. By definition, the impulse response is the response of the system to a unit impulse. We can thus obtain it by executing the system (SRE) with a unit impulse as input data-set, and then use the corresponding output as the system impulse response (this can be seen as a form of abstract simulation).

Since this impulse response can be infinite, some approximations may be required. If the system is numerically stable, it is assumed that the impulse response reaches zero beyond some reasonable range. We use this fact to bound the size of the impulse response obtained through our method.

If the graph does not contain any cycle, its impulse response is finite and can be computed exactly. The same is true for the parts of the graph that are not reachable from any cycle, such the as the first convolution in the Deriche filter in Figure 4. This observation can be used to further optimize the computation. Finally, this process can be easily parallelized.

D. Computing the Noise Power Function

Noise mean and variance at output o are given by summing those of individual noise contributions:

µ o = ∑ s∈Erro µ s,o σ 2 o = ∑ s∈Erro σ 2 s,o (3)
Noise power is given by:

P o = µ 2 o + σ 2 o (4)

E. Applying the Model

To compute noise power for a given fixed-point specification, (µ s , σ 2 s ) must be determined for every noise source s in the graph. It can then be computed with equations 2, 3 and 4.

As we observed in Section V-A, a noise source can actually correspond to several quantizations. Consider the Equation 1 from Deriche example (Figure 4). Quantizations may occur after every operation in the expression and at the end of the computation, when the result is constrained to the format of its final destination. In our example, the total error is thus the sum of six potential quantization errors. The properties (mean, variance) of these intermediate noises can be determined from the quantization mode (truncation, rounding, convergent rounding), the precision of the operands and the quantization step. As these noises are independent, µ s and σ 2 s can be computed by summation as in Equation 3.

VI. EXPERIMENTAL VALIDATION

The proposed analysis has been implemented within a compiler framework, and validated with signal and image processing applications. We contrast our scalability to another accuracy analysis tool based on the work by Ménard et al. [START_REF] Ménard | Analytical Fixed-Point Accuracy Evaluation in Linear Time-Invariant Systems[END_REF], called ID.fix 1 . We used a Linux machine with Intel core i7 (3.0GHz) and 8GB of memory for our experiments.

Table I shows the time required to construct the accuracy model for our benchmark applications. Although ID.Fix is not designed to handle 2D algorithms, it can still process them by seeing each pixel as a different input. The difference in scalability can be clearly observed with image processing applications. Our approach is unaffected by the input image size whereas ID.fix is unable to handle the Deriche filter on 16 × 16 images. 2D recursive filters, such as Deriche, have wide impulse response, and 16 × 16 image would not be sufficient to approximate the error.

Once the model is constructed, evaluating the model for a given quantization decision takes less than a second in all cases. The models constructed by ID.fix also take less than a second to evaluate. Table II demonstrates the accuracy of the constructed models validated against simulations.

VII. DISCUSSION AND RELATED WORK

Previous attempts to perform accuracy analysis from program sources have significant limitations that restrict their use to small kernels. Most tools [START_REF] Caffarena | SQNR Estimation of Fixed-Point DSP Algorithms[END_REF], [START_REF] Novo | Accuracy vs Speed Tradeoffs in the Estimation of Fixed-point Errors on Linear Time-invariant Systems[END_REF] ignore the program control-flow and remain at the basic-block level (i.e., DFG level). Hence, they fail at capturing inter basic-block interactions and it is unclear how such tools can handle sequence of loop nests accessing same arrays with non-trivial indexing functions, which are common in image/video processing.

The work by Ménard et al. [START_REF] Ménard | Analytical Fixed-Point Accuracy Evaluation in Linear Time-Invariant Systems[END_REF] addresses this problem by completely flattening the program control-flow to obtain a single (large) SFG for the whole program, leading to scalability issues as pointed out in Section III. In contrast, our approach is able to capture more complex propagation patterns (interblock, instance-wise) without resorting to flattening/unrolling.

A. Accuracy Analysis from a Compiler Point of View

From a strict accuracy modeling perspective, our approach may appear as a straightforward extension of earlier work on LTI systems to multi-dimensional signal processing systems (LSI). Indeed, we do use the same type of domain specific knowledge, which is understanding that the underlying system is LSI. Such knowledge can directly be exposed in the system specifications when appropriate model-based design approaches are used (e.g., Simulink), but are not available when working with arbitrary source programs.

Our work addresses two fundamental issues, (i) how to infer such domain specific knowledge from input source programs (i.e., from a compiler point of view), and (ii) how to represent the system in a compact way such that the analysis scales to programs beyond simple signal processing kernels. In this respect, our contribution also has some connections with algorithm recognition problems [START_REF] Alias | On the Recognition of Algorithm Templates[END_REF] and its scope goes beyond accuracy modeling techniques (at least as defined by the signal processing community).

One may argue that all these domain specific knowledge can be expressed in the source program as annotations/pragmas. However, it is extremely difficult to define a reasonably small set of annotations that can capture the domain specific knowledge, when the tool is expected to handle many different ways of implementing an algorithm. Our observation is that trying to specify sufficient information through annotations ends up asking the programmer to specify the from system à la model-based design, eliminating the benefits of automation. Moreover, the tool is incapable of verifying if the given annotations are "correct" if it relies all important information to be supplied by the user.

B. Extensions to non-LTI Systems

As already pointed out by previous authors, the limitation to LSI systems prevents from using the technique when non linear operators are involved in the algorithm or when the system impulse response changes over time, as is the case for adaptive filters. However, several extensions to the noise propagation model have been proposed to handle such situations. For example, Constantinides [START_REF] Constantinides | Perturbation Analysis for Word-length Optimization[END_REF] proposed to handle non-LTI, non-recursive filters by linearizing noise contributions. A variety of methods was developed to compute the linearization coefficients [START_REF] Ménard | Analytical Fixed-Point Accuracy Evaluation in Linear Time-Invariant Systems[END_REF], [START_REF] Caffarena | SQNR Estimation of Fixed-Point DSP Algorithms[END_REF], [START_REF] Shi | A Perturbation Theory on Statistical Quantization Effects in Fixed-Point DSP with Non-Stationary Input[END_REF]. Recently, Rocher et al. showed how to handle any recursive and non-recursive filter composed of smooth operations [START_REF] Rocher | Analytical Approach for Numerical Accuracy Estimation of Fixed-Point Systems Based on Smooth Operations[END_REF].

Most of these extensions follow the same principle: they use additional profiling information2 (statistic over signal values, or cross-correlation information) to built a LTI/LSI approximation of the real system. This approximation is then used to derive the error noise statistics. Since they ultimately recast the problem into a LTI/LSI formulation, we believe that such extensions can be easily integrated within our approach.

C. Working with non-polyhedral Programs

As explained in Section IV, our approach operates on programs with regular control structure that can be represented using the polyhedral model. Many/most image signal processing algorithms fall into this category, and those that do not are unlikely to correspond to LTI systems anyway. However, there exist a few exceptions, the most notable one being the Fast Fourier Transform. In such a case, the best solution to the accuracy analysis problem is probably to rely on a domain/application specific noise propagation model that would capture the recursive nature of the algorithm (borrowing from the philosophy of tools such as Spiral [START_REF] Koutsoyannis | Improving fixed-point accuracy of FFT cores in O-OFDM systems[END_REF]).

D. Arithmetic Accuracy Bounds in Numerical Computation

The problem of numerical accuracy analysis in program is a widely studied topic. Most work in the field aim at deriving safe error bounds between finite and infinite precision algorithms (or between two finite precision specifications) [START_REF] Goubault | Static Analysis of Finite Precision Computations[END_REF], or to help finding the cause of numerical inaccuracies [START_REF] Benz | A Dynamic Program Analysis to Find Floating-point Accuracy Problems[END_REF]. Boland and Constantinides [START_REF] Boland | A scalable approach for automated precision analysis[END_REF] showed how scalability issues could be addressed in that context by trading the quality of bounds for faster accuracy evaluation. The problem we address here is very different as we are interested in the statistical properties (such as mean and power) of the noise error rather than its bounds. Our work could be extended to derive bounds on error instead of statistical properties.

VIII. CONCLUSION

In this work, we have proposed a scalable approach to the problem of accuracy modeling for floating-point to fixedpoint conversion. Our approach works by inferring high-level convolution operations from the original source code, and explicitly modeling them as part of the program representation. The two main outcomes of the technique are (i) improved scalability with respect to prior work and (ii) the ability to handle higher dimensional linear systems. Future work directions include extending the approach to handle correlated quantification noises, and to support non-linear, but differentiable, operations.

Fig. 1 :

 1 Fig. 1: Graph representation of convolution system z(i) = (x1 * c1) + (x2 * c2). In the explicit representation, edges represent flow of data, and nodes represents collections of operations or data. Note that a single node represents more than one operation/data. In the compact representation, edges represent a convolution by the kernel associated with the edge, and the sum of the convolution results becomes the value of a node if it has multiple incoming edges. Indices of nodes are not shown except for input nodes in the compact representation.

  has a window of length 4; -3 ≤ n ≤ 0, where the coefficients are [c b 0 a].

Fig. 2 :

 2 Fig. 2: Code fragments from Deriche filter and corresponding SREs. Deriche filter applies total of 4 IIR filters in different directions (right, left, up, down). The SREs shown exclude boundary cases (e.g., when j = 0) to simplify our presentation.

TABLE I :

 I Model construction time for our tool and ID.fix.

	Algorithm	ID.Fix (s) Our tool (s)
	IIR8	23.1	20.5
	Sobel (32 × 32)	169.1	9.2
	Sobel (64 × 64)	2173.1	9.7
	Sobel (128 × 128)	-	9.4
	Gaussian blur (32 × 32)	160.1	10.2
	Gaussian blur (64 × 64)	2010.9	9.5
	Gaussian blur (128 × 128)	-	9.4
	Deriche blur (16 × 16)	-	6.5

TABLE II :

 II Validation of our model against simulations.

	Algorithm Simulation Our tool	Error
		(dB)	(dB)	(dB/%)
	IIR8	-17.80	-17.84	-0.04/-0.2%
	Sobel	11.62	12.04	0.42/3.6%
	Gauss	3.78	3.78	<0.01/0.1%
	Deriche	-18.01	-18.06 -0.05/-2.78%

http://idfix.gforge.inria.fr/

This profiling stage is performed only once on the initial floating point data-set, it should not be confused with simulation based approaches.


Dalla-Salda Guillermina
email: guillermina@inta.gob.ar

María Elena Fernández

Anne-Sophie Sergent

Philippe Rozenberg

Eric Badel

Alejandro Martinez-Meier

Dynamics of cavitation in a Douglas-fir tree-ring: transition-wood, the lord of the ring?

Keywords: conductivity loss, embolism, annual ring, wood density, earlywood, latewood

The objective of this work was to investigate the dynamics of embolism formation within a Douglasfir tree-ring. Four resistant and four vulnerable 10year-old trees were selected among 50 trees, based on their P 50 . Stem samples, taken next to those used to obtain the vulnerability to cavitation curves, were collected and submitted to increasing positive pressures, in order to simulate increasing tension caused by water stress in the xylem. Then the conductive surface of the samples was stained and scanned and the images were analyzed. X-ray microdensity profiles were obtained on the same samples. The microdensity profiles of the 2011 ring were analyzed in three parts, earlywood, transitionwood and latewood. The dynamics of embolism propagation was observed separately in these three parts. Our results showed that the initiation and the propagation of the cavitation follow a discrete trend, with at least two successive initiation events: first cavitation initiates and propagates rapidly in the latewood. Then, a second cavitation event begins and spreads in the earlywood and eventually propagates to the transition-wood, which remains the last conductive part in the ring before full embolism. We observed that resistant to cavitation trees showed lower transition-wood density than vulnerable to cavitation trees.

Introduction

Water transport in trees can be limited by cavitation. Cavitation is the breaking of the water column that follows atmospheric demand and transpiration increase and/or soil water content decrease, particularly in drought conditions. Cavitation initiates in the cells of the conductive xylem (tracheids in conifers) when the cell water copes with increasing tensions (Hacke and Sperry 2001) and turns from a liquid metastable state to a more stable gaseous state [START_REF] Cruiziat | The hydraulic architecture of trees: an introduction[END_REF][START_REF] Tyree | Xylem structure and the ascent of sap[END_REF], Mayr et al. 2014). According to the air-seeding hypothesis, cavitation propagates from cell to cell as air is aspirated from a non-functional air-filled tracheid into a conductive one through the inter-tracheid pits [START_REF] Delzon | Mechanism of water-stress induced cavitation in conifers: bordered pit structure and function support the hypothesis of seal capillary-seeding[END_REF]). Little information is available about where exactly cavitation springs up in a tree-ring and how it propagates inside it. [START_REF] Domec | How do water storage differ in coniferous earlywood and latewood?[END_REF] have studied vulnerability to cavitation of earlywood and latewood in Douglas-fir (Pseudotsuga menziesii (Mirb) Franco) finding that latewood was more vulnerable to cavitation at high trunk water potentials than earlywood. This suggests that cavitation could start in latewood, and then spreads to earlywood. The dynamics of cavitation initiation and propagation within the growth ring may have important functional implications since if water is released from cavitated latewood this could decrease earlywood tension (capacitance function, [START_REF] Domec | Safety factors for xylem failure by implosion and airseeding within roots, trunks and branches of young and old conifer trees[END_REF]) leading to ks maintenance during drought. However, if cavitation initiates on the largest tracheids of earlywood, this may result in a strong and early decrease in ks. Earlywood is produced in temperate climate during the first part of the growing season. In conifers, earlywood is formed of large-lumen, narrow-wall low-density tracheids. Latewood is produced later during the growing season, generally when atmospheric demand increases, soil water content depletes and photoperiod decreases. Latewood is formed of narrower-lumen, thicker-wall denser tracheids. The transition from earlywood to latewood is more or less gradual according to the species [START_REF] Rozenberg | Clonal variation of indirect cambium reaction to within-growing season temperature changes in Douglasfir[END_REF]) and the cambial age [START_REF] Zobel | Juvenile wood in forest Trees[END_REF]. Conventionally, tree-rings are cut into earlywood and latewood. This action arbitrarily splits the transition zone and mixes it with earlywood and latewood. Some authors have suggested that it may be useful to separate the transition zone from earlywood and latewood, cutting the ring into three parts [START_REF] Rozenberg | Clonal variation of indirect cambium reaction to within-growing season temperature changes in Douglasfir[END_REF][START_REF] Franceschini | Effect of ring width, cambial age, and climatic variables on the within-ring wood density profile of Norway spruce Picea abies (L.) Karst[END_REF], thus defining more homogeneous earlywood and latewood. As a matter of fact, [START_REF] Dalla-Salda | Genetic variation of xylem hydraulic properties shows that wood density is involved in adaptation to drought in Douglas-fir (Pseudotsuga Menziesii (Mirb.))[END_REF] found a stronger relationship between wood density and resistance to cavitation when they used the homogeneous first part of the earlywood rather than the conventionally-defined heterogeneous earlywood. Changes in ring structure may result in changes in the hydraulic architecture of plants with direct consequences on water supply to the leaves and tree growth and/or survival [START_REF] Froux | Xylem hydraulic efficiency versus vulnerability in seedlings of four contrasting Mediterranean tree species (Cedrus atlantica, Cupressus sempervirens, Pinus halepensis and Pinus nigra)[END_REF][START_REF] De Micco | Wood anatomy and hydraulic architecture of stems and twigs of some Mediterranean trees and shrubs along a mesic-xeric gradient[END_REF], Martinez-Meier et al. 2008[START_REF] Hoffmann | Hydraulic failure and tree dieback are associated with high wood density in a temperate forest under extreme drought[END_REF]. Tree-ring structure can be thoroughly described using X-ray microdensity profiles [START_REF] Polge | Etablissement des courbes de variations de la densité du bois par exploration densitométrique de radiographies d'échantillons prélevés à la tarière sur des arbres vivants[END_REF]. Wood microdensity is determined by cell dimension: it describes the proportion of cell wall in a given woody tissue. This proportion reflects wood mechanical and functional properties [START_REF] Zobel | Wood variation its causes and control[END_REF], Hacke and Sperry 2001, Hacke et al. 2001) that are related with the three wood functions: the mechanical support of the body of the plant [START_REF] Zobel | Juvenile wood in forest Trees[END_REF], the storage of water and biological chemicals [START_REF] Stratton | Stem water storage capacity and efficiency of water transport: their functional significance in a Hawaiian dry forest[END_REF] and the safe and efficient transport of water from the roots to the leaves [START_REF] Maherali | Adaptive variation in the vulnerability of woody plants to xylem cavitation[END_REF], Pitterman et al. 2006,). Deterministic relationship between wood density and vulnerability to cavitation is unclear. Several studies found positive correlations (Hacke and Sperry 2001[START_REF] Rosner | Hydraulic and mechanical properties of young Norway spruce clones related to growth and wood structure[END_REF], Dalla-Salda et al. 2009, 2011), while others reported a lack of relationship (Lamy et al. 2012). The present study had three objectives: 1) to investigate the dynamics of embolism formation within the tree-ring of Douglas-Fir; 2) to compare the dynamics of tracheids embolism within the tree-ring between resistant and vulnerable trees; and 3) to determine whether the variation found between resistant and vulnerable trees (if any) in terms of dynamic of embolism formation could be related with the microdensity structure of the ring. Initiation and propagation of embolism in the ring were studied for homogeneously-defined earlywood, latewood and a separate transitional zone by means of analyses of images of stained wood surface and X-ray microtomography of wood samples taken at variable embolism levels from exposure to different pressure intensities. X-ray microtomography appears now as a reference technique for embolism visualization in xylem [START_REF] Cochard | X-ray microtomography (micro-CT): a reference technology for high-resolution quantification of xylem embolism in trees[END_REF]. We applied the staining method analysis to eight Douglas-fir trees selected for their differential resistance to cavitation: the four more resistant and the four more vulnerable among a set of fifty trees. Microdensity profiles of wood samples of the same 8 trees were studied with x-ray methodology [START_REF] Polge | Etablissement des courbes de variations de la densité du bois par exploration densitométrique de radiographies d'échantillons prélevés à la tarière sur des arbres vivants[END_REF]). To our knowledge, there is no study documenting the intra-ring dynamic of embolism formation in relation with microdensity.

Materials and Methods

Plant Material

The studied plants were 11 years old Douglas-firs from a genetic trial belonging to INTA (Instituto Nacional de Tecnología Agropecuaria, Argentina) near Trevelin, Patagonia -Argentina (43° 06' S -71° 32' W) Firstly, stem vulnerability to cavitation (VC) curves were obtained with the air injection method [START_REF] Cochard | Vulnerability of several conifers to air embolism[END_REF]) on 50 trees from 10 different genetic entities. In order to characterize the trees for their resistance to cavitation, the wood samples were taken from the stem portion corresponding to the last ring (year 2011), following the methodology described in Dalla-Salda et al. (2009, 2011). Based on the P 50 (50% percent of conductivity loss) calculated on their VC curves (see Table 1 in Results section), eight trees were selected for a detailed study of embolism propagation within the growth-ring: the four most resistant and the four most vulnerable to cavitation.

Dynamics of embolism propagation studied through staining method:

From the eight selected trees, 30 cm long longitudinal samples from the 2011 ring were taken, next to the samples used to obtain the vulnerability curves. Each of these longitudinal samples was cut under water into six 5-cm-long segments. Each portion was subsequently subjected to a different pressure level (0, 1, 2, 2.5, 3 and 3.5 MPa) by means of a pressure chamber (Model 1000 Pressure Chamber Instrument, pms instruments, Co. USA) with its distal end protruding from the chamber. After being submitted to a given pressure level the sample was connected to a water column with deionized and distilled water colored with Phloxine 1% in order to stain the conductive area. The samples were then oven dried and cut in half. One of the middle surfaces was polished and digitalized. The obtained images were divided in earlywood (EW), transition-wood (TW) and latewood (LW) zones as described in the next section, and then the percentage of the stained (i.e. conductive) area was measured in each type of wood with ImageJ image analysis software [START_REF] Rasband | [END_REF]) (Fig. 1). Subtracting this area to the total conductive area, the variable named "% of non-conductive area" was estimated assuming that it represents a measure of the proportion of cavitated (embolized) cells in the sample.

Wood microdensity:

Microdensity profiles were obtained from wood samples of the 8 selected trees using indirect x-ray microdensitometry (Polge 1960). The resulting films were scanned and processed with WinDENDRO software [START_REF] Guay | A new automatic and interactive tree ring measurement system based on a line scan camera[END_REF]. A function written in Rlanguage (R Core Team 2014) allowed dividing each 2011 ring into three parts: earlywood, transition-wood and latewood (hereafter EW, TW and LW, respectively). The function firstly identifies the minimum (MID) and maximum ring density (MAD), calculates the average density of these two extreme values and locates it along the x axis (Mdp: red vertical line in Fig. 2) as in [START_REF] Mothe | Analyse microdensitométrique appliquée au bois: méthode de traitement des données utilisée à l'Inra-ERQB (programme Cerd)[END_REF]. At this point the ring is divided in two parts, corresponding to the conventional EW and LW portions. Subsequently the R function repeats the extreme average procedure described above in each part of the ring, defining Mdp1 (Mdp1: blue vertical line in Fig. 2) and Mdp2 (Mdp2: green vertical line in Fig. 2). The TW is finally defined as the portion of the ring located between Mdp1 and Mdp2. For each ring the following microdensity variables were calculated: ring width, ring density, EW width, EW density, TW width, TW density, the slope describing the increment of wood density within the transition-wood portion (slope), LW width and LW density. We obtained the length of each portion by multiplying the number of density points in each part of the ring per 24 µm (distance between each point of microdensity measurement). These measures were plotted on the stained wood surface images and were used to divide them into three parts matching the EW-TW-LW division in the corresponding ring microdensity profile (Fig. 1). A t-test (p<0.05) was assessed to compare percentage of non-conductive area at the different pressures between the parts of the ring and between the resistant and vulnerable trees. The same procedure was used to compare the differences between the resistant and vulnerable trees for the microdensity variables. A standardized tree-ring microdensity profile was computed to describe the intra-annual density variation along a common scale for both groups of trees. Each of the 2011 tree-rings was described by 33 density points by gathering and averaging microdensity values. For this purpose, significant difference between resistant and vulnerable trees was assessed by a t-test at the level of each density point (p< 0.05).

Embolism visualization by microtomography:

Xylem embolism was visualized by using x-ray microtomographic scans recorded in the PIAF laboratory (INRA, Clermont-Ferrand, France), under the same conditions as in [START_REF] Cochard | X-ray microtomography (micro-CT): a reference technology for high-resolution quantification of xylem embolism in trees[END_REF]. [START_REF] Brodersen | Visualizing wood anatomy in three with high-resolution X-ray micro-tomography (MCT)[END_REF] samples without cutting and surface preparation. One longitudinal sample was taken as described above from a single tree selected in a resistant to cavitation family of a Douglas-fir progeny trial located at Ecouves, Normandy, France (unpublished data). The sample was cut under water in three portions, 5 ×5 mm in the transversal section and about 20 mm in the longitudinal direction. Two of them were submitted to hydraulic tensions by centrifugation with the Cavitron technique [START_REF] Cochard | Evaluation of a new centrifuge technique for rapid generation of xylem vulnerability curves[END_REF].Two centrifugation speeds were chosen in order to generate -2.5 and -4 MPa that caused 50 and 88 % loss of hydraulic conductivity, respectively. After centrifugation, the samples were immediately dropped in liquid paraffin wax and cooled in order to prevent their drying during the x-ray scan. The samples were then placed in the Xray microtomograph (Nanotom 180 XS, GE, Wunstorf, Germany). X-ray source settings were 50 kV and 275 μA. After 3D reconstruction, the field of view was a 6.2 × 6.2 × 6.2 mm 3 with a spatial resolution of 2.8 × 2.8 × 2.8 μm per voxel. One radial-tangential 2D cross section was virtually extracted from the middle of the volume in order to visualize the part of the ring located in the middle of the sample, using VGStudioMax software (Volume Graphics, Heidelberg, Germany). Then, images were analyzed using the ImageJ software in order to compute the percentage of cross section area that was embolized.

Results

Vulnerability to cavitation:

Mean and standard deviation of P 12 , P 50 , P 88 and slope of the 50 trees for which the vulnerability to cavitation curves were obtained are shown in table 1. The resistant and vulnerable trees (n=4) were selected by their contrasting P 50 values. Mean ± standard deviation of P 50 of resistant and vulnerable trees were -3.32 ± 0.18 and -2.44 ± 0.18 MPa, respectively (Fig. 3).

Propagation of embolism within the growth ring:

Considering all analyzed trees, in the EW the % of non conductive area is very close to 0% at atmospheric pressure and at -1 MPa. It increases slowly to around 10% at -2 MPa, and then rapidly rises to more than 80% between -2 and -3.5 MPa (see Fig. 6 in next section, which shows the common general pattern of vulnerable and resistant trees). At all pressure levels the % of non conductive area is significantly lower in the EW than in the LW, except at -3.5 MPa where there is no significant difference between both parts of the ring (Table 2).

In the TW the % of non conductive area is significantly different from the one in the EW at -2 MPa, where it is significantly higher, and at -3.5 MPa, where it is significantly lower (Table 2).

Figure 3 Vulnerability to cavitation curves of 50 Douglasfir trees (grey lines) constructed with the air injection method. The four most vulnerable (red lines) and the four most resistant to cavitation (green lines) were selected for the present study. Black line with black full circles shows the % of non conductive area for the whole ring constructed by the staining method.

Partly in agreement with the observations made from staining wood samples subjected to different pressure levels, the microtomography image of a sample at atmospheric pressure shows that all tracheids were water filled (are conductive) in the three parts of the ring (Fig. 4A). The same conclusion was obtained for EW and TW from the staining method, but a high proportion of non-conductive area was observed in LW even at atmospheric pressure, which could be interpreted as embolized tracheids if microtomography images were not available. This novel method indicates that LW tracheids had water in their lumens at atmospheric pressure. Therefore, they were not stained at the very low pressures of ks determination due to their high resistance to water flow and not because they were cavitated.

The figure 4B shows the microtomography image of a sample that experienced a -2.5 MPa pressure that corresponds to P 50 (Fig 1 andTable 1). At this pressure, the x-ray observation enlightens that the LW region was greatly embolized. The EW portion showed low embolism and the TW was still fully water saturated.

- Finally, the figure 4C shows a sample submitted to a -4 MPa pressure that corresponds to P 88 (Fig. 1, Table 1). The x-ray observation indicated clearly the EW and LW regions were largely embolized and only the TW region showed water in the lumens of the tracheids.

Similarly, observation of the radial-longitudinal cross section of the sample submitted to -4 MPa allows to appreciate the dark grey strips corresponding to the empty lumens of the embolized tracheids both in the EW and in the LW parts of the ring (Fig. 5). Bright voxels (light grey strips) correspond to free water in the lumens, which are only located in the TW part of the ring (Fig. 5).

Resistant vs vulnerable trees:

Figure 6 shows the variation of the % area of non conductive tracheids as a function of the applied pressure in the EW, TW and LW for the resistant and vulnerable trees. In the LW, the % area of non conductive tracheids is nearly 80% at atmospheric pressure and increases to almost 100% at -3 and -3.5 MPa. As stated before from microtomography images, the high % of non conductive area at low pressures was not actually reflecting embolism, but only the high resistance to water flow in this portion of the growth ring. There was no significant difference between the resistant and the vulnerable trees in the proportion of non conductive area in the LW portion of the ring. In addition, there were also no significant differences in EW patterns of non conductive area between the resistant and vulnerable trees (Fig. 6).

In contrast, at -2 MPa there was a significant difference between the resistant and the vulnerable trees in the TW: the resistant trees showed a significantly lower % of non conductive area (Fig. 6). Regarding microdensity parameters, significant differences were found between the resistant and the vulnerable trees for ring width, EW width, TW width, and the slope and mean density of the TW (Fig. 7, only the variables with significant differences between the resistant and the vulnerable trees are shown). No significant difference was found for EW density nor for LW density. The resistant trees showed a significantly less dense TW with a less steep slope.

The resistant trees showed also a wider 2011 ring with significantly wider TW and EW (Fig. 7).

Figure 7:

The microdensity variables for which significant differences were found between the vulnerable and the resistant to cavitation trees. The significant differences between the resistant and the vulnerable trees are also viewed observing the standardized ring microdensity profiles (Fig. 8). The vulnerable trees have significantly higher microdensity values from the points 23 to 30 that effectively correspond to the TW part of the ring (Fig. 8).

Discussion

Up to now, most VC curves published for woody plants were performed on branches, seedlings and roots [START_REF] Cochard | Vulnerability of several conifers to air embolism[END_REF]Sperry and Ikeda, 1997;Kavanagh et al., 1999). As far as we know, little information is available from direct measures obtained on stem xylem of trees [START_REF] Domec | Relationship between growth rates and xylem hydraulic characteristics in young, mature and old-growth ponderosa pine trees[END_REF].

Our measurements of P 50 , achieved on trunk wood samples are consistent with those by other authors (Domec and Gartner, 2001) for the same species and are consistent with those of our previous studies (Dalla-Salda et al. 2009[START_REF] Dalla-Salda | Genetic variation of xylem hydraulic properties shows that wood density is involved in adaptation to drought in Douglas-fir (Pseudotsuga Menziesii (Mirb.))[END_REF]. The VC curves of 50 Douglas-fir trees in this study were performed by the use of the air injection method. With 8 trees (the 4 most vulnerable and the 4 more resistant to cavitation) we computed the % of non conductive area for the 2011 whole ring by the staining method.

Both methods seem to show the same types of curves and the % of non conductive area seem to be well represented by the % loss of conductivity (Fig. 3).

The dynamics of cavitation in the 2011 ring can be -in part-deduced from the plot in figure 6, which presents the curves describing the variation of the % of non conductive area with pressure. If we assume that the pressure increase along the x-axis of the plot simulates the condition of the start and development of a water deficit event, then these curves are likely to characterize the initiation and the propagation of cavitation in natural conditions. The observations of the stained samples images and of the microtomography images are generally consistent: the stained portions of the sample surface correspond to water filled lumen in the microtomography images, and the non-stained portions correspond mostly to embolized cells, with the exception of LW at atmospheric pressure. The results based on the stained samples images show that the LW is already mostly non-conductive at atmospheric pressure, while in the same conditions the EW and TW are either completely or nearly completely conductive. Although, the microtomography image (Fig. 4 A) shows that at atmospheric pressure the three parts of the ring are full of water. This may be due to the fact that the microtomography is an image of the inside of a sealed sample, while the image of the stained surface is that of an open sample which has been handled and more exposed to air. Other hypothesis is that the difference found in the LW portion between the two methods is due to LW great resistance to water flux preventing stained water from arriving to that middle part of the sample were the images were done. Both methods show the same trend in the initiation and propagation of the cavitation in the ring: the cavitation starts very early into the LW reaching more than 80% of non conductive area at -1 MPa and 100% between -2.5 and -3.5 MPa. This result confirms the findings by [START_REF] Domec | How do water storage differ in coniferous earlywood and latewood?[END_REF] who state that the LW is more vulnerable to cavitation than the EW at high water potential.

The trend is completely different in the EW and TW, where the cavitation appears later, then propagates slowly between atmospheric pressure and -2 MPa, except in the TW of the vulnerable trees, where it spreads significantly faster than in the EW between -1 and -2 MPa. The % of non conductive area eventually reaches a high level in the EW (over 80%), similar (not significantly different) to that reached in the LW, while the TW remains with a significantly lower % of non conductive area (between 60 and 70%). This trend is similar in the microtomography images, except that the cavitation seems to start and propagate in the EW before propagating to the transition zone located between the EW and the LW, which remains mostly conductive at -4 MPa. Hence it seems that cavitation occurs first in the LW, then "jumps" to the EW and TW (according to the stained samples images) or to the EW only (according to the microtomography images) and eventually to the TW. Thus, there are at least two successive events of cavitation initiation, the first one in the LW, then the second one in the EW. Whether there is a third initiation event in the TW or whether cavitation propagates from the EW to the contiguous TW cannot be inferred from our results. Assuming that the suggest that, in Douglas-fir, cavitation dynamic does not follow a trend linearly related to tracheid dimensions. According to [START_REF] Domec | Bordered pit structure and function determine spatial patterns of air seeding thresholds in xylem of Douglas-fir (Pseudotsuga menziesii, Pinaceae) trees[END_REF], pit number and dimension rapidly decrease from the EW to the LW. To our knowledge, no information has been published about the characteristics of the Douglas-fir TW pits. If we assume that the pits characteristics continuously vary from the EW to the LW, then cavitation propagation does not follow a trend linearly related to pits number and dimensions.

Published studies about the distribution of the pits and of their characteristics in radial-longitudinal (RL) and tangential-longitudinal (TL) cell-walls, show that in conifer species bordered pits are mostly located in the RL walls [START_REF] Leitch | Evidence for auxin regulation of bordered-pit positioning during tracheid differentiation in Larix Laricina[END_REF][START_REF] Domec | Bordered pit structure and function determine spatial patterns of air seeding thresholds in xylem of Douglas-fir (Pseudotsuga menziesii, Pinaceae) trees[END_REF][START_REF] Kitin | Anatomical features that facilitate radial flow across growth rings and from xylem to cambium in Cryptomeria japonica[END_REF]) rather than in the TL wall. In spite of we cannot extrapolate this pattern to all conifers, at least in Cryptomeria japonica, a temperate conifer, most of the EW tracheids lack of bordered pits in the TL walls, with the exception of the EW tracheids positioned next to growth-ring boundaries, whose tangential pits are connected with the previous ring's LW tracheids [START_REF] Kitin | Anatomical features that facilitate radial flow across growth rings and from xylem to cambium in Cryptomeria japonica[END_REF]. In contrast to the EW, the LW tracheids commonly have bordered pits in the TL walls not only at the growth-ring boundary but also in the tracheid layers apart from it [START_REF] Kitin | Anatomical features that facilitate radial flow across growth rings and from xylem to cambium in Cryptomeria japonica[END_REF]. The TW tracheids also show bordered pits in the TL walls. The pit membranes in the RL walls have torus and margo with loosely arranged fibrils, with spaces ('pores') in between the fibrils. In contrast, no torus-margo structure and no pores are visible in the pit membranes of the TL walls of either the EW or LW tracheids. The bordered pits in the RL walls of the EW tracheids are much larger than those in the TL walls of the boundary EW tracheids. In the LW and TW, however, pits in the RL and TL walls are similar in size. In another study, [START_REF] Utsumi | Seasonal and perennial changes in the distribution of water in the sapwood of conifers in a sub-frigid zone[END_REF] found in three other conifers, that most of the bordered pits of tracheids are located in the RL wall; large-scale cavitation would be facilitated in the tangential direction when such tracheids are exposed to water stress. Moreover, these cells were formed at the same time and experienced the same growing conditions, thus we suggest that they developed the same sensitivity to cavitation. These results support our hypothesis of, at least two independent cavitation initiation events in the LW and the EW, followed by within EW and LW tangential spreading of cavitation that could be promoted by the characteristics of the pits in the TL tracheids walls.

The difference between the stained samples data and the observations on the microtomography images for the cavitation propagation trend could be related to the fact that in the stained samples images, the definition of the TW is based on the ring microdensity profile.

Figure 1 shows that this so-called TW partly overlaps the stained portion of the ring and the non-conductive portion located on the LW side of the ring. In the tomography images, the zone located between the EW and the LW (Fig. 4C) seems on average much more conductive than the TW delineated on figure 1. This explains why on the figure 1 the % of non conductive area of the microdensity-defined TW is around 70% at -3.5 MPa while on the figure 4C the major part of the area of the zone located between the EW and the LW is formed of cells with water-filled lumens at 4 MPa.

There is no significant difference between the resistant and the vulnerable trees for the % of non conductive area in the EW and in the LW, whatever the pressure value. According to the tomography image, which was taken in a resistant tree, the TW remains the last conductive part during the progression of the cavitation in the ring. Both stained samples and microtomography results suggest that the TW would play a major role in the cavitation resistance.. [START_REF] Dalla-Salda | Genetic variation of xylem hydraulic properties shows that wood density is involved in adaptation to drought in Douglas-fir (Pseudotsuga Menziesii (Mirb.))[END_REF] found that more resistant to cavitation trees had a higher mean density of an alternative EW density variable. In this study we also found significant differences between the resistant and the vulnerable trees for some alternative microdensity variables, but, these differences were found only for the TW portion: width, density and slope. Furthermore we found that the TW of the vulnerable trees is significantly denser than that of the resistant trees. The comparison of the average standardized rings of the resistant and vulnerable trees show the same trend (Fig. 8). At the level of the 50 trees among which were selected the four resistant and the four vulnerable trees of this study, there is no relationship between P 50 and TW density (data not shown). Hence the significant difference observed in this study between the resistant and the vulnerable trees for TW density could be probably a selection artifact. Conversely, there is a significant negative relationship between P 50 and TW width (r=-0.48 and p-value = 0.0008, data not shown) and there is a significant positive relationship between P 50 and TW slope (r=0.50 and p-value = 0.0005, data not shown). These results are consistent with the significant difference found in this study between the resistant and the vulnerable trees for TW width and TW slope.

It is well-known that not only lumen dimension, but also pit number, dimension and characteristics play a major role in resistance to cavitation [START_REF] Pittermann | The relationships between xylem safety and hydraulic efficiency in the Cupressaceae: the evolution of pit membrane form and function[END_REF][START_REF] Delzon | Mechanism of water-stress induced cavitation in conifers: bordered pit structure and function support the hypothesis of seal capillary-seeding[END_REF][START_REF] Bouche | A broad survey of hydraulic and mechanical safety in the xylem of conifers[END_REF][START_REF] Tixier | Modelling the mechanical behaviour of pit membranes in bordered pits with respect to cavitation resistance in angiosperms[END_REF]. The air-seeding hypothesis predicts that xylem embolism resistance is linked directly to bordered pit functioning [START_REF] Choat | Structure and function of bordered pits: new discoveries and impacts on wholeplant hydraulic function[END_REF][START_REF] Lens | Embolism resistance as a key mechanism to understand adaptive plant strategies[END_REF][START_REF] Jansen | Plasmodesmatal pores in the torus of bordered pit membranes affect cavitation resistance of conifer xylem[END_REF]. In any case, cavitation dynamics is still a controversial issue and not fully elucidated [START_REF] Rockwell | Cavitation and its discontents: opportunities for resolving current controversies[END_REF]. According to [START_REF] Domec | Bordered pit structure and function determine spatial patterns of air seeding thresholds in xylem of Douglas-fir (Pseudotsuga menziesii, Pinaceae) trees[END_REF], the cell-walls of the Douglas-fir LW are punctured by a low number of small size bordered pits, while the EW cell-walls are punctured by a large number of big pits. The higher rigidity of the LW pits membrane would explain the LW low resistance to cavitation. Conversely, despite its high lumen-diameter and large number of big pits, the higher flexibility of the EW pits membranes would explain the EW higher resistance to cavitation. But what could explain the even higher resistance to cavitation of the TW? The TW is by definition a transition zone between the EW and the LW. Hence cell and pits characteristics could be intermediate between those of the EW and the LW: by its own nature, the TW cell-lumen diameter is intermediate between that of the EW and the LW. If we hypothesize that the TW pit dimension and pit number are also smaller than in the EW, but that the TW pit membrane is as flexible or nearly as flexible as that of the EW pits, then these variables would tend to confer to the TW a resistance to cavitation that would be higher than in the EW.

In addition we found that the resistant trees have a significantly wider 2011 ring, with significantly wider EW and TW and no change in the LW width: hence the increased resistance to cavitation of the resistant trees could be explained by this larger amount of more resistant to cavitation EW and especially TW. They also have a less steep microdensity slope, i.e. a more gradual lumen diameter decrease, which could also tend to increase the TW resistance to cavitation.

Conclusion

According to our results, the initiation and the propagation of the cavitation in a Douglas-fir tree-ring seems to follow a discrete trend, with at least two successive initiation events: at the very start of a water deficit episode, cavitation appears to initiate and propagate rapidly in the LW. Later, when the water deficit persists and increases, a second cavitation event begin and spread in the EW and eventually propagate to the TW, which remains the last conductive part in the ring before full cavitation. Whether cavitation actually propagates from the EW to the TW or whether it initiates for the third time in the TW is not known. Our results suggest that the TW may play a leading role in the between-tree variation of cavitation resistance and thus should be carefully considered when studying adaptation to drought. The sampling procedure carried out in this research benefited from a "smart sampling" approach [START_REF] Rosner | Wood density as a screening trait for drought sensitivity in Norway spruce[END_REF]): the samples used for measuring vulnerability to cavitation, conductivity and microdensity came all from the same annual ring (here the 2011 ring), increasing the probability to highlight significant relationships and differences. We used an original method to study the within-ring microdensity variation, defining a transition zone between the EW and LW, increasing the precision of the study and highlighting the major role of this intermediate part of the ring. Improving the definition method of the TW or dividing the ring into more than three parts will likely contribute to better explain the initiation and the propagation of the cavitation in the ring, as well as better explain other functional processes occurring in the ring. While little data is available, within-ring variation of pits characteristics seem to have a high explanatory potential for propagation of cavitation in the ring. The ratio between the pit aperture and the torus diameters is an important correlate with cavitation resistance in conifers and could perhaps help elucidating differences in the resistance of EW and TW tracheids. In the near future, studies of resistance to drought should concentrate efforts on large scale investigation of within-ring variation of pits characteristics and of the relationships between such pits characteristics and cavitation.

Figure 1 :

 1 Figure 1: Example of wood samples (photographs taken in a binocular loupe and corresponding scanned images) subjected previously to different pressure levels (MPa) showing conductive (red-stained) and non-conductive areas. Pressures were exerted using a Scholander bomb (PMS, instruments).

Figure 2 :

 2 Figure 2: Ring microdensity vs. radial distance in the stem (one ring was analyzed). The following variables are shown on the graph: the minimum and maximum wood density (MID and MAD, respectively), the mean density point of the ring (Mdp), the mean density point for the intial and final part of the ring (Mdp1 and Mdp 2, respectively).

  [START_REF] Torres-Ruiz | Vulnerability to cavitation in Olea europaea current-year shoots: more support to the openvessel artefact with centrifuge and air-injection techniques[END_REF] were the first to show xylem vulnerability to embolism using x-ray microtomographic scans. The great advantage of x-ray microtomography is the threedimensional observation of internal structure of

Figure 4 :

 4 Figure 4: X-ray microtomographic images of a cross section of Douglas-fir. A: image of the control sample submitted to atmospheric pressure. B: sample submitted to a pressure of -2.5 MPa. C: sample submitted to a pressure of -4 MPa. The grey levels refer to the local x-ray absorption property. The black spots (1) correspond to the empty lumens of the embolized tracheids; the grey spots correspond to the water filled lumens (2) and the brighter clear segments correspond to the water saturated cell walls (3). The spatial resolution is 2.8 µm and the scale bar represents 500 µm.

Figure 5 :

 5 Figure 5: x-ray microtomographic radial-longitudinal cross section of Douglas-fir. The sample was submitted to a pressure of -4MPa. Dark grey strips (1) correspond to the empty lumens of the embolized tracheids. The bright strips (light grey) correspond to free water columns in the tracheids (2). Bright vertical lines correspond to water saturated cell walls (3). Spatial resolution is 2.8 µm. Scale bar is 500 µm.

Figure 6 :

 6 Figure 6: Percentage of non-conductive area assessed by the staining method. Figure shows each type of wood (black, red and blue lines for the LW, TW and EW, respectively) in the resistant (diamonds with continuous lines) and vulnerable (squares with dotted lines) trees after being submitted to different pressure levels.

Figure 8 :

 8 Figure 8: Standardized ring profile for the vulnerable (red diamonds) and the resistant (green squares) trees for the 2011 ring. Asterisks show significant differences (p<0.05) between of the resistant and the vulnerable trees for the microdensity values.

  % of no conductive area corresponds to the propagation of the xylem cavitation, these results

Table 1

 1 Mean and standard deviation of the variables describing the VC of the fifty evaluated trees.

	Variable	Mean	Standard deviation
	P 12	-0.97	0.55
	P 50	-2.99	0.34
	P 88	-4.04	0.32
	slope	36.03	11.31

Table 2

 2 Probabilities associated to Student t-test for the comparison of the percentage of mean non-conductive area between the different parts of the ring (EW= earlywood; TW= transition-wood; LW=latewood)

				Pressures (MPa)
	Comparisons	0	1	2	2.5	3	3.5
	EW -TW	ns	ns	0.0045	Ns	Ns	0.0461
	EW -LW	0.0003 2.77 -7 2.94 -12 0.0005 0.0047	ns
	TW -LW	0.0005	4.81 -8 1.61 -6	0.0011 0.0009 0.0317

Acknowledgements

The authors would like to thank INTA-Programa Nacional-Forestales, PNFOR-1104073 and CONICET for founding part of this study. The x-ray tomography scans were partly founded by the PitBulles project (ANR no. 2010 Blanc 171001). Special thanks to Dr. Victor Mondino (INTA-Trevelin) who facilitated the access to the plant material and helped with the sample collection. Finally we would like to thank the two anonymous reviewers that contributed with their corrections and suggestions.