Olivier Rozenbaum
email: rozenbaum@cnrs-orleans.fr

3-D characterization of weathered building limestones by high resolution synchrotron X-ray microtomography

Keywords: Stone weathering, 3-D image analysis, X-ray microtomography, porous media, segmentation, topology. 1. Research aims

 by two-dimensional image analysis was selected. The 3-D images were analysed by a set of mathematical tools that enable the description of the pore and solid phase distribution. Results show that 3-D image analysis is a powerful technique to characterise the morphological, structural and topological differences due to weathering. The paper also discusses criteria for mathematically determining whether a stone is weathered or not.

Introduction

In a previous article [START_REF] Rozenbaum | 2-D image analysis: A complementary tool for characterizing quarry and weathered building limestone[END_REF], an experimental approach to characterize weathered stones by using some classical techniques was reported, and certain image analysis tools were presented. Image analysis provides a huge quantity of data, offering important insights to characterize the porous medium. Within recent decades, 2-D image analysis (images generally acquired by SEM) has been commonly used to characterize different porous materials such as porous silica, soils, concrete and stones (e.g. [START_REF] Cousin | Three-dimensional analysis of a loamyclay soil using pore and solid chord distributions[END_REF][START_REF] D'acqui | Study of soil porosity with mercury porosimetry and image analysis on backscattered electron scanning images (BESI). Application to tilled ''crusting soils'' in Zimbabwe. Soil micromorphology: studies in management and genesis[END_REF][START_REF] Levitz | Disordered porous solids: from chord distribution to small angle scattering[END_REF][START_REF] Török | Morphology and mineralogy of weathering crusts on highly porous oolitic limestone, a case study from Budapest[END_REF].

However, the transfer properties of porous media are related to 3-D characteristics of the porous network which cannot be appropriately described by 2-D analysis since the latter cannot provide complete topological information on the porous network. Since the middle of the '90s, 3-D X-ray microtomography has significantly improved and has been widely used in many fields of investigation such as materials science and geology [START_REF] Bentz | Microstructure and transport properties of porous building materials II. Three dimensional X-ray tomography studies[END_REF][START_REF] Biswal | Quantitative analysis of experimental and synthetic microstructures for sedimentary rock[END_REF][START_REF] Burlion | X-ray microtomography: Application to microstructure analysis of a cementitious material during leaching process[END_REF]Cnudde and Jacobs, 2004a;[START_REF] Coles | Pore level imaging of fluid transport using synchrotron X-ray microtomography[END_REF][START_REF] Ketcham | Acquisition, optimization and interpretation of X-ray computed tomographic imagery: applications to the geosciences[END_REF][START_REF] Lindquist | Investigating 3-D geometry of porous media from high resolution images[END_REF][START_REF] Lindquist | Quantitative analysis of three-dimensional x-ray tomographic images, Developments in X-Ray Tomography III[END_REF][START_REF] Pierret | 3-D reconstruction and quantification of macropores using X-ray computed tomography and image analysis[END_REF][START_REF] Spanne | Synchrotron computed microtomography of porous media: topology and transports[END_REF]).

Sandstones, for example, have been extensively characterized and the subject of many publications [START_REF] Biswal | Quantitative analysis of experimental and synthetic microstructures for sedimentary rock[END_REF][START_REF] Coles | Pore level imaging of fluid transport using synchrotron X-ray microtomography[END_REF][START_REF] Spanne | Synchrotron computed microtomography of porous media: topology and transports[END_REF]. This technique is underused, however, in the field of building conservation, despite the fact that it can provide this branch of research with new information, as has been shown recently (Cnudde et al., 2004b;[START_REF] Graef | A sensitivity study for the visualisation of bacterial weathering of concrete and stone with computerised X-ray microtomography[END_REF]. In these two studies, X-ray computed tomography (CT) analysis was used to visualize bacterial weathering of concretes and stones or to locate water repellents and consolidants in building stones. The technique was applied with a resolution of about 10 µm, as this resolution was considered appropriate for building stones. However, part of the porous network was not analyzed because the corresponding pores were much smaller than the resolution used. For some building stones such as the one selected for this study, a high proportion of the porous network requires a higher resolution than 10 µm if it is to be analyzed properly (i.e. less than 1 µm voxel size) and a high-contrast image. Currently, these specifications require the use of synchrotron radiation [START_REF] Gallucci | 3-D experimental investigation of the microstructure of cement pastes using synchrotron X-ray microtomography (μCT)[END_REF][START_REF] Colombo | Characterization of the morphology of cellular ceramics by 3D image processing of X-ray tomography[END_REF][START_REF] Denis | 3D complex shape characterization by statistical analysis: Application to aluminium alloys[END_REF][START_REF] Rolland Du Roscoat | Estimation of microstructural properties from synchrotron X-ray microtomography and determination of the REV in paper materials[END_REF][START_REF] Tafforeau | Nature of laminations and mineralization in rhinoceros enamel using histology and X-ray synchrotron microtomography: Potential implications for palaeoenvironmental isotopic studies[END_REF]) coupled with reliable image analysis.

The present work focuses on the analysis of 3-D images using mathematical tools and illustrates the interest of this kind of analysis in building stone conservation. The stone selected was the Tuffeau limestone which had been studied earlier (e.g. [START_REF] Rozenbaum | 2-D image analysis: A complementary tool for characterizing quarry and weathered building limestone[END_REF] and references herein) by classical techniques.

Material and methods

Materials

Tuffeau limestone is a yellowish-white porous sedimentary limestone, mainly composed of calcite, silica in the form of opal cristobalite-tridymite and quartz, and some secondary minerals such as clays and micas. The unweathered Tuffeau stone from a quarry (Saint-Cyren-Bourg, France) studied in this paper was mainly composed of calcite (0.503 g g -1 ), and silica (0.452 g g -1 ) and the total porosity was 0.481. It was a multi-scale porous medium since the equivalent pore size distribution ranged from 0.01 to 50 µm [START_REF] Rozenbaum | 2-D image analysis: A complementary tool for characterizing quarry and weathered building limestone[END_REF].

The weathered Tuffeau stone selected came from the Saint-Donatien church located in Orléans (France). Note that the origin of weathering was undetermined and that no building or environmental monitoring was done before sampling. Blocks were extracted from the northeast wall of the bell tower at an elevation of 20 m. Sampling occurred during restoration of the church, which consisted in replacing entire blocks by new ones. The exposed surface was a grey crust harder than the stone core, and the first 3 cm below the surface were crumbly with micro-cracks parallel to the surface. Gypsum was present in the first 30 mm below the surface but mainly concentrated in the first 20 µm, corresponding to the grey crust [START_REF] Rozenbaum | 2-D image analysis: A complementary tool for characterizing quarry and weathered building limestone[END_REF]. A morphologically representative block was selected and two sets of samples were then constituted. The first set was taken from the 0-3 mm layer below the exposed surface and corresponds to a weathered Tuffeau stone. It consisted of three samples W0-1, W1-2 and W2-3 from the 0-1, 1-2 and 2-3 mm layers below the exposed surface (Fig. 1). The second set was taken from the 80-81 mm layer underneath the exposed surface (UW sample) and was considered as a sample of unweathered Tuffeau stone based on the criteria established in a previous study [START_REF] Rozenbaum | 2-D image analysis: A complementary tool for characterizing quarry and weathered building limestone[END_REF], where "unweathered" was considered to mean that part of the stone (in depth) had undergone very little (or no) evolution and was in almost the same state as at the building of the wall. This was the case far from the surface (80-81 mm in the present study, following the previous article) where the stone had rarely (or never) undergone wetting/drying cycles. Unweathered stones were therefore considered as reference stones.

Obtaining the 3-D images

X-ray microtomography is a technique which enables high-resolution evaluation of 3D structural information. The microtomography images were collected at the MS-X04SA-Tomo beam line of SLS (Swiss Line Source, [START_REF] Stampanoni | High resolution X-ray detector for synchrotron-based microtomography[END_REF] and the energy used was 12 keV. Successive rotations of the sample, 1500 projections, corresponding to 1500 angular positions ranging between 0° and 180°, were acquired by a camera with 1024×1024 pixels.

The good quality of the beam made it possible to image samples at a 0.7 µm pixel resolution.

In order to remain within the field of view of the detector and to avoid additional artifacts, the samples had to measure less than 700 µm. The samples were prepared as cylindrical cores, 700 µm in diameter, and were mounted on a vertical rotator on a goniometric cradle.

Preparing limestone samples with these dimensions requires particular care. The samples were resin impregnated in order to make a 700 µm thin section [START_REF] Bruand | Effect of water content on the fabric of a soil material: an experimental approach[END_REF]. The thin section was then cut into 700 µm matches with a diamond saw, and the edges of the matches were filed down to obtain quasi cylinders.

Imaging time was approximately 120 minutes. A filtered back-projection algorithm was used to reconstruct horizontal slices 0.7 µm thick from the projections. Stacking the 1024 slices reconstructed produced a 3-D volume (Fig. 2). Each pixel from a tomographic slice corresponds to a specific voxel in the 3-D image. The whole data set corresponding to one 3-D image has an approximate size of one Gbyte (1024 3 voxels). Each voxel was characterized by its grey level (ranging from 0 to 255 for an 8-byte image) which depended on the X-ray attenuation coefficient of the element [START_REF] Baruchel | X-Ray tomomography in material science[END_REF].

Image analysis and segmentation

Image quantification involves correctly segmenting the solid and porous phases. However, it was not possible to obtain a direct threshold of the raw images from their histogram (Fig. 3, dot curve), as the expected three phases (pore, silica, and calcite) could not be identified. In an ideal image without any noise, two grey levels (ll and lh) would be selected. Hence, all the voxels belonging (i) to the [0, ll] range would represent the pore phase, (ii) to the [ll, lh] range would represent the silica phase and (iii) to the [lh, 255] range would represent the calcite phase. In the present case of a real image, the presence of noise [START_REF] Trong | A simple methodology to segment X-ray tomographic images of a multiphasic building stone[END_REF] meant that these two grey levels, which distinguish each phase from the other two, could not be selected. Unfortunately, the algorithms generally used for 2-D image segmentation were not efficient for most of the 3-D images due to finite resolution effects of noise, leading to similar grey levels for different phases. Many authors used their own algorithm, depending on the type of materials studied (e.g. [START_REF] Trong | A simple methodology to segment X-ray tomographic images of a multiphasic building stone[END_REF][START_REF] Lindquist | Quantitative analysis of three-dimensional x-ray tomographic images, Developments in X-Ray Tomography III[END_REF][START_REF] Vignoles | Image segmentation for phase-contrast hard X-ray CMT of C/C composites[END_REF]. Initial studies of our complex 3-D images showed that these algorithms were inefficient and new procedures to binarize the images were therefore created [START_REF] Trong | A simple methodology to segment X-ray tomographic images of a multiphasic building stone[END_REF].

Briefly, a two-step filtering approach was chosen. Alternate sequential filters were first used to remove the noise but, as they destroy the structural components smaller than the structuring element used, they cannot be pursued far enough. As a second step in the filtering process, a mosaic operator, based on a pragmatic but efficient marker determination, was therefore implemented to simplify the images further. This yielded the mosaic (Fig. 4b) of flat zones, which was much simpler than the starting image (Fig. 4a) and straightforward to segment.

Using the histogram (Fig. 3, solid curve) the threshold levels were easily identified (as the minimum values between the peaks). The values of ll and lh for each sample are reported in table 1.

Estimators

Different estimators were then calculated on the segmented images. Most of these estimators were described in [START_REF] Rozenbaum | 2-D image analysis: A complementary tool for characterizing quarry and weathered building limestone[END_REF] for 2-D images but their extension to 3-D images was straightforward. The estimators are briefly explained below; the reader is referred to references [START_REF] Levitz | Disordered porous solids: from chord distribution to small angle scattering[END_REF] and [START_REF] Rozenbaum | 2-D image analysis: A complementary tool for characterizing quarry and weathered building limestone[END_REF] for more detail.

Porosity

As the thresholds were distinguished, the porosity, calcite and silica volume fraction could be calculated (Table 1). The fraction of phase i within the volume was defined as:

T i i V V = φ , (1)
where V i and V T are the volume of phase i and the total volume of the sample, respectively.

For a 3-D digitized medium, V i and V T are the number of voxels corresponding to phase i and to all the phases, respectively. Hence, total porosity φ is calculated if phase i is the pore phase.

Chord distributions

Chord distributions are stereological tools used to describe the interface between pore and solid phases [START_REF] Levitz | Disordered porous solids: from chord distribution to small angle scattering[END_REF][START_REF] Rozenbaum | 2-D image analysis: A complementary tool for characterizing quarry and weathered building limestone[END_REF]. The chord distribution gives the probability of having a chord length between r and r+dr. Chord distribution can be calculated either for the pores (f p (r), or pore chord distribution) or for the solid (f s (r), or solid chord distribution). Furthermore, the first momentum of f p (r) (called l p ) and the first momentum of f s (r) (called l s ) are defined as:

∫ = dr ) r ( f r ) r ( l p p (2) and, ∫ = dr ) r ( f r ) r ( l s s . ( 3
)
The values of l p and l s can be considered to be the estimators of the mean size of the pore and solid phases.

Due to digitizing, the smallest pores (or solids) were badly defined and did not correspond to real objects: as the chord distributions for r < 3 pixels presented a linear increase, they were therefore not taken into account [START_REF] Cousin | Three-dimensional analysis of a loamyclay soil using pore and solid chord distributions[END_REF][START_REF] Rozenbaum | 2-D image analysis: A complementary tool for characterizing quarry and weathered building limestone[END_REF].

Autocorrelation function

The autocorrelation function

) (r R Z
r is a statistical description of a porous medium [START_REF] Levitz | Disordered porous solids: from chord distribution to small angle scattering[END_REF][START_REF] Rozenbaum | 2-D image analysis: A complementary tool for characterizing quarry and weathered building limestone[END_REF]. This function can be interpreted as the probability of finding two randomly selected points that are both in the same phase. For an isotropic

medium with r r v = , ) r ( R ) r ( R Z Z = r
, the autocorrelation function depends on distance alone.

Otherwise, the one dimensional autocorrelation functions for r r parallel to the Ox, Oy or Oz axis (where Ox, Oy, Oz define an orthogonal reference), respectively noted Zx R , Zy R and Zz R , give information about isotropy or anisotropy, since pore space anisotropy is revealed as a disparity between the one-dimensional autocorrelation functions along different directions [START_REF] Frykman | Anisotropy in pore networks analyzed with 2-D autocorrelation (variomaps)[END_REF]. Another characteristic length scale is the correlation length L, defined as the integral of the autocorrelation function. It represents a length scale characteristic of the spatial structure [START_REF] Ioannidis | Statistical analysis of the porous microstructure as a method for estimating reservoir permeability[END_REF]:

∫ ∞ = 0 ) ( dr r R L Z . ( 4
)
For anisotropic materials, L x , L y and L z are defined accordingly from Zx R , Zy R and Zz R .

Topological analysis

In addition to these estimators, which are also valid in 2-D, the topology of the pore space is a fundamental property as regards flow and transport in porous media (and can only be estimated on a 3-D volume), since pore connectivity strongly influences these properties.

Fortunately the main advantage of X-ray microtomography resides in 3-D structure analysis and topology analysis. The Hoshen-Kopelman algorithm (6-connectivity) was used to identify (or label) each pore with a number [START_REF] Hoshen | Percolation and cluster distribution I. Cluster multiple labelling technique and critical concentration algorithm[END_REF]. This is also a simple and efficient cluster identification algorithm for identifying connected clusters on a lattice (percolation theory). The size is determined by the number of voxels in each cluster. This algorithm was also used to label each calcite cluster and each silica cluster.

The basic topological properties are the number N of isolated pores, the number of redundant connections C within the pore space, and the number of completely enclosed cavities H [START_REF] Vogel | Quantitative morphology and network representation of soil pore structure[END_REF]. In other words, C represents the number of bonds that can be removed without creating an additional isolated cluster. The Euler-Poincaré characteristic ε V was then computed as ε V = N -C + H. For the topology of porous media such as stones, the quantity H is of minor importance since solid grains or particles completely surrounded by pores are scarce. H was therefore set to zero within the previous equation [START_REF] Vogel | Quantitative morphology and network representation of soil pore structure[END_REF]. As this characteristic is related to the total volume (V, in voxel units) and in order to compare images of different sizes, the specific Euler-Poincaré number was defined as

χ v = ε v /
V. This number can be considered as a connectivity measure: the greater the number is, the less connectivity can be expected.

Analysis and Discussion

Pore and solid phase morphology

As the grey level value of a voxel was linked to the X-ray absorption of the sample at the voxel position [START_REF] Baruchel | X-Ray tomomography in material science[END_REF], the pores appeared in dark grey, the silica compounds in medium grey and the calcite compounds in light grey (Fig. 2). The different phases were easily distinguishable to the naked eye: calcite was present in the form of large irregular grains (sparitic calcite) or small grains that look like crumble (micritic calcite), whereas silica had the form of large quartz crystals or small spheres of opal. The fine structure of the latter

was not visible at this pixel size, but their global shape made it possible to distinguish their spatial distribution. For example, on the upper right corner of Fig. 2f, some opals embedded in a large calcite grain can be seen.

Simple observation of the 3-D images such as those presented on Fig. 2 showed that the W0-1 sample was strongly altered due to the presence of numerous large cracks. For the W1-2 and W2-3 samples, some cracks or large pores were visible but, at first glance, the weathering was far less severe than for W0-1. Some pores were visible on the UW sample, as expected for an unweathered Tuffeau stone (no cracks or large pores). The four samples could not be distinguished by their solid compounds morphologies, only by the presence or absence of cracks or pores i.e. the UW sample was compact, the W1-2 sample slightly less compact (due to the presence of some large pores), the W2-3 sample was a little less compact again, and the W0-1 sample was not compact at all.

Porous phase analysis.

Average porosities calculated for the four 3D images are reported in Table 1 (volume fraction of pores). The porosity determined by image analysis for the unweathered stone (0.228) was smaller than the porosity determined from the density measurement conducted in [START_REF] Rozenbaum | 2-D image analysis: A complementary tool for characterizing quarry and weathered building limestone[END_REF], where a porosity of 0.481 was recorded for this sample. This difference may be related to the small pores which were not taken into account in the 3D image analysis because they were smaller than the pixel size (0.7 µm). These small pores correspond partly to those resulting from the assemblage of the silica and calcite grains, the structure and morphology of these grains (e.g opal spherolites) as well as fine particles such as clays. Nevertheless, a comparative study of these 3D stone images enabled weathering features to be detected.

First of all, as shown in Fig. 5, the pore and solid chord distributions for all samples presented exponential decreases. This type of porous medium, in which both the porous and solid chord distributions are exponential, is called a long-range random medium [START_REF] Levitz | Disordered porous solids: from chord distribution to small angle scattering[END_REF][START_REF] Rozenbaum | 2-D image analysis: A complementary tool for characterizing quarry and weathered building limestone[END_REF]. In other words, the pore and solid phases in the 3D image were randomly distributed. The absence of correlation peaks demonstrated that the pore and solid phases were heterogeneously distributed in size.

As shown on Table 1, the porosity of the surface sample (W0-1) was high (0.467) in comparison with the other porosities. In particular, from the exposed surface, porosities were much lower in the 1-2 mm layer (0.289) and in the 2-3 mm layer (0.296) (that were both very close) than in the 0-1 mm layer (0.467). The porosity of the unweathered Tuffeau stone (i.e.

between 80 and 81 mm below the exposed surface) was lower than that of the weathered Tuffeau stones (0.228 compared to the recorded porosity in the 0-1, 1-2 and 2-3 mm layers).

Similarly, the pore chord distribution (Fig. 5a) of the surface sample (0-1 mm layer) extended to the longest chord range (about 310 µm for a probability of 2×10 -6 ) while the pore chord distributions of the 1-2 mm and 2-3 mm layer samples extended to about 200 µm for a probability of 2×10 -6 . Finally, the pore chord distribution of the 80-81 mm layer extended to the shortest chord range (about 130 µm for a probability of 2×10 -6 ). This behaviour was in accordance with that of the total porosities and suggests that, statistically, there were strong modifications in the pore distribution within the W0-1 sample and non-negligible porous modifications in the W1-2 and W2-3 samples. Using the classification used in [START_REF] Rozenbaum | 2-D image analysis: A complementary tool for characterizing quarry and weathered building limestone[END_REF] and the results presented above and in the following (visual analysis of 3D volumes, porosity, chord length distributions, anisotropy etc…), the W0-1 sample could be considered as strongly weathered and the W1-2 and W2-3 samples could be considered as moderately weathered.

Furthermore, the Hoshen-Kopelmann algorithm showed that for each sample one percolating pore cluster existed and accounted for the quasi-totality (about 99%) of the total pore volume, demonstrating that the porous media can be considered as fully connected. This was confirmed by the negative values of the specific Euler-Poincaré characteristics χ V (Table 1).

The percolating pore cluster (Table 1) was the highest in the 0-1 mm layer (≈2.5×10 8 voxels)

and was the lowest in the 80-81 mm layer, i.e. in the unweathered Tuffeau stone (≈1.2×10 8 voxels). The remaining unconnected pores (4% to 0.4% of the total pore volume, depending on the sample) consisted of small clusters with sizes ranging from a few voxels up to about 80×10 3 voxels (i.e. considerably smaller than the percolating pore volume). The quantity of these remaining unconnected pores was the lowest for the weathered sample W0-1 (≈3400)

and was the highest for the unweathered UW sample (≈6200). In fact, because of weathering, initially unconnected pores were connected together, resulting in larger pore clusters. This also led to a decrease in the number of unconnected pores and an increase in the main percolating pore cluster volume.

Isotropy/anisotropy

The autocorrelation functions ( Zx R , Zy R and Zz R ) and the average value Z R for the UW sample are presented in Fig. 6a. The Oxy plane was taken parallel to the stone surface, the axis Oz being perpendicular to the rod sample (Fig. 1). All these functions presented a decreasing behaviour without any particular correlation. For the UW sample, no particular disparity was observed for the one-dimensional autocorrelation functions along the orthogonal

x-, y-, z-directions and Z R (Fig. 6a) and between the correlation lengths (Table 2). This means that the UW sample presents no, or only very slight, anisotropies. In contrast, the W0-1 sample autocorrelation functions (Fig. 6b) and the correlation lengths (Table 2) presented a non-negligible anisotropy. As shown on Fig. 2, this anisotropy was caused by cracks. The anisotropy decreased for the W1-2 sample and decreased yet further for W2-3 (curves not shown here), becoming nil for the UW core sample, as described previously.

Solid phase analysis

Calcite volume fractions (Table 1) were quasi constant (about 0.29) whatever the sample, except for the W0-1 sample which had a much lower volume fraction (0.19). In other words, the calcite volume fractions of the moderately weathered samples (W1-2 and W2-3) were similar but, surprisingly, were also similar to the calcite volume fraction of the unweathered stone. The silica volume fractions were the highest for the unweathered sample (0.48) and decreased for the moderately weathered stones (about 0.42). Finally the silica volume fraction of the W0-1 sample was the smallest (0.35).

The solid chord distribution (Fig. 5b) in the 0-1 mm layer extended to the shortest chord range (about 260 µm for a probability of 2.10 -6 ), whereas the solid chord distributions of the W1-2 and UW samples were both similar and extended to the highest chord length (about 425µm

for a probability of 2×10 -6 ). Moreover, the solid chord distribution of the W2-3 sample extended to an intermediate chord length (about 375 µm for a probability of 2×10 -6 ). This result meant that the distances between two interfaces within the solid phase were statistically equivalent in the W1-2 and UW samples and were both larger than (i) the W2-3 sample above 250 µm and (ii) the W0-1 sample above 75 µm.

The Hoshen-Kopelmann algorithm showed that for each sample, one percolating calcite cluster and one silica cluster existed. In other words, the silica grains (resp. calcite) adhered to each other in order to build a cohesive backbone that guarantees, a priori, the integrity of the stone. As for the calcite volume fraction analysis, the volumes of the W1-2, W2-3 and UW calcite clusters were quasi constant, about 1.5×10 8 voxels, accounting for more than 99% of the total calcite fraction (Table 1). On the other hand, the W0-1 sample had a much lower calcite cluster (about 0.4×10 8 voxels) representing "only" 37% of the total calcite fraction.

The unconnected grain parts were similar (0.4%) for the W1-2, W2-3 and UW samples, and some large unconnected large calcite grains (about 2×10 6 to 5×10 5 voxels i.e. much less than the percolating calcite volume) were found. Although the number of calcite clusters was very similar (about 2900) whatever the sample, the unconnected grains in the W0-1 sample were 10 times larger than in the other samples. This may result from the mechanical constraints suffered and/or partial dissolution of the initial calcite cohesive backbone of the unweathered Tuffeau stone as it appears in the 80-81 mm stone, partially cracking this backbone into relatively large unconnected calcite grains.

The percolating silica cluster (Table 1) was the smallest in the 0-1 mm layer (≈1.9×10 8 voxels) and was the largest in the 80-81 mm layer (≈2.6×10 8 voxels). The percolating silica clusters represent more than 99% of the total silica fraction for all the samples. The remaining unconnected silica grains consisted of about 4000 small grains with sizes ranging from a few voxels up to about 50×10 3 voxels for all the samples (i.e. much smaller than the percolating silica volume).

At this stage of analysis it seemed that some silica grains had been removed from their settings during the stone history. This phenomenon occurred, however, without any calcite dissolution for the moderately weathered stones.

2D/3D comparisons and weathering criteria

It was possible to make some qualitative comparisons between the 3-D image data and the earlier 2-D image data [START_REF] Rozenbaum | 2-D image analysis: A complementary tool for characterizing quarry and weathered building limestone[END_REF], bearing in mind however that (i) the pixel and voxel sizes were different (2.8 and 0.7 µm, respectively) and (ii) the studied stone thickness from the surface exposed to environmental agents was also different: for the 2-D images the samples were 1.2 to 1.4 mm thick, while for the 3-D images the rod samples were "only" 0.7 mm thick.

The present results (variation of porosity according to depth, anisotropy etc…) are consistent with those recorded earlier and enable a more detailed analysis of the differences resulting from weathering (topological analysis). Moreover, as expected for stereological tools, the general trends of the chord length distributions were similar to the 2-D chord length analysis.

One criterion (among others) to characterise the degree of weathering in [START_REF] Rozenbaum | 2-D image analysis: A complementary tool for characterizing quarry and weathered building limestone[END_REF] was obtained from the solid chord distribution analysis. The classification used in [START_REF] Rozenbaum | 2-D image analysis: A complementary tool for characterizing quarry and weathered building limestone[END_REF] established that the solid chord distribution of moderately weathered or unweathered samples tended to extend to the longest chord length, while the solid chord distribution of strongly weathered samples tended to extend to the shortest chord length (due mainly to cracks and large pores caused by partial grain dissolution and or mechanical constraints such as freezing/thawing cycles). This classification seemed to remain valid as the W1-2, W2-3 and UW solid chord distributions (Fig. 2) extended to longer chord lengths than those of the W0-1 sample. It was observed that the W2-3 solid chord distribution was modified above 230 µm. In other words, the solid chord distributions of the W1-2 and W2-3 samples show a slightly different solid phase distribution. This is rather surprising, as the volume fraction of pore, calcite and silica (Table 1) and the pore chord distributions of these two samples were identical.

One hypothesis is that during weathering, cracks and pores separated mineral grains, thus breaking the continuity of the solid. The result was a statistical decrease in the solid phase (silica+calcite) on the long solid chord range for the W2-3 sample (Fig 5b). In contrast (the case of the W1-2 sample), cracks "pushed" the solid phase, resulting in a higher concentration of the solid phase within some localised zones. Within these zones the solid phase was more compact (Fig. 2d) than within the W2-3 sample. Nevertheless, this analysis did not enable any conclusions to be drawn about the resistance of these layers, since although the grains were in contact as shown here, this gave no information about the bond strengths between these grains.

This hypothesis explains the difference between the solid chord distributions of the W1-2 and W2-3 samples and also explains why this behaviour was not observed in the previous study but was revealed now. In the present investigation (3-D), the probed sample thickness was 0.7 mm while the probed depth in the previous article (2-D) was twice as high. The latter depth was sufficiently large to contain both large cracks and compact zones (Fig 1e in [START_REF] Rozenbaum | 2-D image analysis: A complementary tool for characterizing quarry and weathered building limestone[END_REF]) in all the moderately weathered samples considered. The W1-2 sample studied here represented the same moderately weathered zone but, due to the smaller depth sampling, W1-2 sample was situated in a zone without large or long cracks. Note that in spite of this lack of long cracks, the porosity of W1-2 sample was higher than the porosity of the UW sample, due to weathering processes. The W1-2 sample contained some compact solid zones which created large solid chords, but also some zones with large pores which also created large pore chords. The same observations can be made for the W2-3 sample except that the solid phase was slightly more dispersed by the porosity. This could explain the small number of cohesive multilayer structures macroscopically observed in the first few millimetres below the surface.

In the previous paper [START_REF] Rozenbaum | 2-D image analysis: A complementary tool for characterizing quarry and weathered building limestone[END_REF], the ratio of the first momentum of chord length distributions was pointed out as a possible indicator of weathering. The lower the ratio was, the stronger the weathering. This is in accordance with the measurements presented here: l s /l p (W0-1)=1.1 and l s /l p (UW)=3.1 (Table 2), confirming the hypothesis that the W1-2 and W2-3 samples were moderately weathered as l s /l p (W1-2)=2.2 and l s /l p (W2-3)=2.1.

Nevertheless, if these ratios are indeed weathering estimators, what is the threshold limit to determine whether a stone is strongly, moderately or not at all weathered?

This requires a statistical analysis on more images. Only new images on weathered and unweathered stones will provide an answer to this question. This will be the next step and the subject of forthcoming microtomography measurements on the synchrotron.

Conclusion

The results recorded in this study confirm and complete the previous work [START_REF] Rozenbaum | 2-D image analysis: A complementary tool for characterizing quarry and weathered building limestone[END_REF] and show that 3-D image analysis can quantify calcite, silica and pore modifications in stones. These modifications are related to the morphological differences caused by weathering on a building limestone: weathering created cracks on the exposed surface and increased the porosity within the first millimetres. Compared to the unweathered Tuffeau stone (which appeared to be isotropic), the weathered Tuffeau stones appeared to be rather anisotropic, which probably led to changes in water transfer properties. One of the main advantages of a 3-D image is that a topological analysis of the different phases present can be performed.

Thus, it has been shown here that percolating pore clusters existed for weathered and unweathered samples, accounting for 99% of the total porosity. This topological analysis showed that weathering connected initially closed porosities, increasing the pore cluster volumes. This increase in porosity (notably by cracks) removed silica grains from their settings, induced below the surface some constraints on the solid phases, and locally compacted the grains. Hence, 3-D image analysis along the weathering depth of a building stone enhances our understanding of weathering mechanisms and, in the longer term, of transfer properties.

Finally weathering criteria have been proposed, but only a statistical analysis on more samples will be able to confirm their validity. This statistical analysis will also make it possible to determine the REV of each sample (moderately or strongly weathered and unweathered) and hence prove whether the samples are representative. Indeed, the main drawback of this technique and the use of high resolution is that the higher the resolution, the smaller the sample size. This will be the aim of the next proposal addressed to a synchrotron X-ray microtomography beam line.

Figure captions

Fig. 1: Schematic representation of the origin of the four cylindrical samples studied in this article. Three came from the weathered surface (between 0 and 3 mm) while the fourth came from the unweathered core (between 80 and 81 mm). where the silica phase is in white (the others are in black). e) thresholding of image (b) where the calcite phase is in white (the others are in black).

Table captions

Table 1: Grey level threshold, volume fraction, volume of percolating cluster (voxels) and specific Euler number for the four samples studied. The latter was calculated by using the library QuantIm [START_REF] Vogel | QuantIm -Some useful C/C++ functions for scientific image processing[END_REF].

Table 2: Ratios of chord length momentum, correlation length and anisotropy for the four samples studied.

Fig. 2

 2 Fig. 2: 3-D images of the weathered tuffeau originating from a church. The images were

Fig. 3 :

 3 Fig. 3: Greyscale histograms recorded for the original image and after the mosaic operator of

Fig. 4

 4 Fig. 4: a) 2-D zoom of a non segmented 3-D image (720×570 pixels). b) Image (a) after

Fig. 5 :

 5 Fig. 5: Comparison between a) the pore chord and b) the solid chord distributions.

Fig. 6 :

 6 Fig. 6: Autocorrelation functions of a) the UW sample and b) the W0-1 sample.

Acknowledgements

The author would like to thank Marco Stampanoni and Amela Grosso (PSI, SLS, Villingen) for scientific support concerning the tomography experiments and Ary Bruand (ISTO) for his comments and suggestions on this paper.

Sample

Grey level

Volume fraction Volume of percolating cluster χ V (µm


Gaëtan Guérin

Dominique Ambard

Pascal Swider

Implant Mechanobiology

Introduction

Conditions influencing bone growth in the early post-operative period include the surgical technique, mechanical [Prendergast, 1997] and biochemical factors [Bailón-Plaza, 2001]. Low performances of implant fixation were generally associated with a low mineralization or a strong heterogeneous distribution of bony structure in the new-formed surrounding tissue and the physicochemical properties of the implant surface might pay a significant role. We previously developed a mechanobiological model of healing coupling porous media mechanics to biomathematics [Ambard, 2006]. To go further, we hypothesized that such mathematical model could be completed to investigate the role of implant surface in cell proliferation, migration, and adhesion. The application concerned our stable canine implant [Vestermark, 2004].

Methods

The coupling of porous media mechanics and biomathematic allowed the diffusive-convectivereactive governing equations (1) to be derived; L, C, D, respectively were the local variation, the convection, the diffusion, and the source terms. (1 )

g g s o g C C (2) 2 [ ( ) (1 ) -(1-) ] o o s o o s o g 0 C N C C C (3)
The PMMA implant [Vestermark, 2004] was the reference and we compared with two other surface treatments: acid-etched and coarse grit blasted acidetched with RGDS peptide. Material properties are given in Table 1 [Dee, 1999], [Rausch, 2007]. Acid C-RGDS PMMA g (e -9 cell -0.5 /s) 5.22 10.1 2.82 D o (e -7 mm 2 /s) 1.75 1.45 1.75 Table 1: Material properties

Results

The implant healing showed a polar symmetry and we plotted in Figure 1 the radial distribution of mineralized fraction from the implant surface toward the surrounding bone after 4 weeks. We observed the influence of implant surface properties since mineralized fraction increased from 62% for the PMMA implant to 85% for the C-RGD surface.

Discussion

The TGFsynthesis coefficient and cell diffusion coefficient D o , were two main parameters that allowed the mechanobiological role of the implant surface to be predicted in time and space. The decrease of cell diffusion (D o ) and the increase of growth factor synthesis ( g ) improved the bone formation. We also noticed that it was decreasing for the highest value of the litterature, probably because of a too rapid accumulation of cells in the vicinity of the implant and an early haptotactic migration towards the surrounding bone where the porosity gradient stayed important.

  the structural (or mineralized fraction) s , the fluid fraction f , the growth factor concentration C g (TGF-) and the osteoblast concentration C o . Structural porosity, fluid flow and growth factors conditionned the cellular behavior (proliferation, chemotactic & haptotactic migrations, mineral fraction aposition). The cell adhesion influenced the motility through the cell diffusion coefficient D o dependant upon the substrate (bone or implant). The growth factor retention into the initial gap was modelled by a local diffusion coefficient. The source of growth factors involved the osteoblast concentration C o , the growth factor concentration C g to take into account the autocrine and paracrine modes of TGF-, and g dependant upon the osteoblast localization (bone or implant). The model of cell proliferation was similar (equ.3); N o being the proliferation threshold, C 0 the initial growth factor concentration.

Figure 1 :

 1 Figure 1: Distribution of mineralized fraction

Journal of Biomechanics 41(S1) 16th ESB Congress, Oral Presentations, Wednesday 9 July 2008


Patrick Cañadas

Bernard Maurin

Haimad Baudriller

Nadir Philippe Montcourrier

Nadir Bettache

Numerical model of the cytoskeleton structuration during cell spreading

teaching and research institutions in France or abroad, or from public or private research centers.

Introduction

The objective of this work is to propose a mechanical model and a numerical method devoted to the cytoskeleton (CSK) form-finding resulting from its structuration during cell spreading. It is now well assumed that the cell mechanics closely depends on its CSK architecture and on the tension and compression forces carried by its filaments. Several structural models have been therefore developed, especially those based on the tensegrity analogy. However, the observed topological and geometrical complexities of the CSK networks lead to difficulties for determining and using a realistic structure [Baudriller et al., 2006]. The presented method allows calculating such a complex architecture and the associated forces in the different CSK filaments of an adherent cell.

CSK modelling and calculation

A cell is modelled by using a granular medium with rigid grains representing respectively: (i) the nucleus, (ii) the centrosome, (iii) the CSK crosslinkers and (iv) some adhesion proteins like integrins. The grains are connected by at-distance mechanical interaction forces associated to the CSK filaments. These forces are initially tacit and the spreading process is then simulated by considering the growing of microtubules from the centrosome towards the integrins. These connect to the substratum by considering adhesive reactions. In the process, the microtubule growing pushes the other grains, thus changing the whole cell shape and generating interaction forces (Fig. 1). This result in a 3D network of compression forces and in an associated network of traction forces, which ensure the mechanical equilibrium of the overall structure connected to the substratum.

Results

The simulations lead to realistic compressive microtubule (star-like organization; Fig. 2) and associated tensile actin networks (Fig. 3). The results highlight basal higher tensions while the less ones appear on the apical face, consistently with in vivo distribution of strong stress fibres and of thin actin cortex.

Discussion and conclusion

The presented model allows determining a complex and realistic 3D structure of the CSK filaments. Their internal tension and compression forces are mechanically equilibrated, such as considered in the tensegrity analogy. This approach may provide significant information for developing a mechanical model of the CSK structure and its spatial distribution of internal forces, which could be useful in further mechanical studies on adherent cell.

Figure 1 :

 1 Figure 1: Steps of the simulation of cell spreading.

Figure 2 :

 2 Figure 2: Obtained microtubule-like compressive network.

Figure 3 :

 3 Figure 3: Resulting actin-like tensile network.


F Cherblanc

A.-L Lozano

F Ouedraogo

J.-C Bénet

Non-equilibrium liquid-gas phase change in hygroscopic porous media

Keywords: phase change, mass non-equilibrium, hygroscopic porous media, chemical potential

Introduction

When focusing on the drying of hygroscopic materials (soil, wood, gel …), liquid-gas phase change is the key phenomenon that enables the complete dewatering. Evaporation is often considered instantaneous when compared to other transport phenomena, what is expressed by the local equilibrium assumption. However, a volatilisation retardation time has been observed at the macroscopic scale under certain conditions [START_REF] Armstrong | Nonequilibrium mass transfer between the vapor, aqueous, and solid phases in unsaturated soils during vapor extraction[END_REF][START_REF] Bénet | Phenomenological relation of phase change of water in a porous medium: experimental verification and measurement of phenomenological coefficient[END_REF]. At the microscopic scale, the water phase is partially adsorbed on grains surface, and it cannot be considered in a "free"' state, where classical thermodynamics applied. In such situations, water displacement induced by pressure gradient is not conceivable and evaporation requires higher energy exchange. At this scale, the phase change phenomenon is generally represented by the Hertz-Knudsen law [START_REF] Eames | The evaporation coefficient of water: a review[END_REF], derived from the kinetic theory for ideal gases. Trying to upscale this relation confronts itself to various interrogations such as the way the phase change flux depends on hygroscopic effects, on liquid-gas interface curvature, on material composition and texture. In that complex framework, macroscopic phase change descriptions have been proposed, generally based on a firstorder kinetic associated with a mass exchange coefficient. This linearization should be valid close to equilibrium conditions, whereas the validity domain has not yet been investigated. To revisit these issues, a macroscopic formulation of liquidgas phase change has been developed. It is derived from the TIP ideas (Thermodynamics of Irreversible Processes), by referring to the chemical potential difference taken as the thermodynamic "driving force". This approach introduced a phenomenological coefficient that must be experimentally determined. An original experimental device has been designed allowing us to explore the complete range of water content, temperature, gas pressure. This leads us to differentiate two domains, a linear behaviour close to equilibrium and a non-linear one far from equilibrium. Eventually, this modelling is confronted to experimental drying kinetics of soil columns.

Phenomenological phase change formulation

As long as the water phase is continuous, the capillary pressure is defined and can be used. However, in the pendular or hygroscopic domain, water is partially adsorbed on the solid surface and the concepts of pressure field and capillary potential are meaningless [START_REF] Nitao | Potentials and their role in transport in porous media[END_REF]. In these situations, the chemical potential becomes the key parameter that governs the water transport. Developing the mass, momentum, energy and entropy balances in a tri-phasic porous medium leads to express the internal production rate of entropy as [START_REF] Bedeaux | Irreversible thermodynamics: a tool to describe phase transitions far from global equilibrium[END_REF][START_REF] Bénet | Phenomenological relation of phase change of water in a porous medium: experimental verification and measurement of phenomenological coefficient[END_REF]:

0 ≥       - - = T J l v µ µ σ
where J is the volumetric phase change flux, l and v are respectively the chemical potentials of water and its vapour. To be used in practical applications, the chemical potential must be expressed as a function of an experimentally measurable quantity, like the vapour pressure. In the framework of linear TIP, the thermodynamic flux is expressed as a linear function of the thermodynamic force by [START_REF] Bénet | Phenomenological relation of phase change of water in a porous medium: experimental verification and measurement of phenomenological coefficient[END_REF][START_REF] Ruiz | Phase change in a heterogeneous medium: comparison between the vaporisation of water and heptane in an unsaturated soil at two temperatures[END_REF]:

        - = veq v P P M R L J ln
where P v is the partial pressure of vapour, and P veq is the partial pressure of vapour in equilibrium conditions. For hygroscopic materials, it depends on the activity of the liquid phase and on the saturating vapour pressure through:

( ) ( )

T P w a P vs veq =
The activity can be determined from experimental desorption isotherms. This relation introduces a phenomenological coefficient L that depends on the state variables, such as the water content w, the temperature T ... The next section focuses on the experimental determination of this coefficient.

Materials and Methods

The porous media used is a natural soil (clayey silty sand), in which the clayey fraction (10%) causes its hygroscopic nature (the hygroscopic domain spreads until a water content of 7%). The principle of experiments consists in creating a thermodynamic non-equilibrium in a soil sample and, then, analysing the return to equilibrium. For that purpose, the gas pressure is lowered by connecting the sample pore network to a vacuum pump during 3 seconds. This operation aims to decrease the vapour pressure below its equilibrium value. Under these conditions, the phase change occurs until the equilibrium vapour pressure is attained. A sensor inside the cell provides a continuous recording of the variations in the total gas pressure and the temperature of the sample. For indication, the time required to get back to equilibrium ranges from 20 to 100 minutes, emphasizing the slowness of phase change in these situations. From these measurements combined with mass balance equations, the phenomenological coefficient L can be determined. For a typical experiment at fixed temperature and water content, this coefficient is plotted as a function of a non-equilibrium criterion: the ratio of the vapour pressure divided by its equilibrium value (Fig. 1). This figure clearly shows two domains. Close to equilibrium, where the ratio of vapour pressure is close to 1, the phase change behaves as a linear phenomenon with a constant mass exchange coefficient L eq . This defines the limit of validity in the framework of linear TIP. Beyond this limit as the non-equilibrium increases, a non-linear behaviour is observed, what is not surprising. The matter is that this situation is regularly encountered in practical drying application since the vapour pressure can easily drop far below its equilibrium value, especially when venting techniques are involved. As a first approximation, the phase change coefficient variation can be represented through the slope k.

For a set of environmental parameters (temperature, pressure, porosity), the phase change coefficient close to equilibrium L eq is determined as a function of the water content (Fig. 2). This is done for various total gas pressures by regulating the vacuum pump used to extract the gas phase. The third curve corresponding to atmospheric pressure is carried out by "completely" removing the gas phase and connecting the pore network with dry air. For each gas pressure, the variation of L eq reaches a maximum around 7%, which is roughly the upper limit of the hygroscopic domain. Above this maximum, the phase change velocity decreases since the liquid-gas interface is getting lower. For water content greater than 12%, the gas phase is occluded and phase change cannot be activated.

Below the maximum, when hygroscopic effects become predominant, the intensity of solid-liquid interactions increases in the adsorbed layers. The supplementary energy requires for water desorption decreases the phase change velocity, leading to a lower coefficient L eq . The influence of these micro-scale interactions have also been shown by adding a mass fraction of calibrated sand to the original material. The presence of sand, weakly hygroscopic, manifests itself through shifting the maximum to the left and increasing the phase change coefficient. The perfect gases theory is unable to predict the influence of gas pressure experimentally observed. At the interface scale, evaporation is predominantly stimulated by shocks with gas phase molecules [START_REF] Matsumoto | Molecular dynamics of fluid phase change[END_REF]. This confirms the greater phase change coefficient resulting from higher gas pressure as shown in figure 2.

In order to illustrate the interests of such investigation, macroscopic drying experiments have been conducted. 30 cm soil columns prepared at 8% water content were placed in a temperature-and moisture-regulated atmosphere. Drying kinetics were established through differential weighing of the column while water content profiles were measured by destructive uptakes. Numerical simulations of water transport were performed in which the phenomena taken into account are liquid-gas phase change and vapour diffusion. In this simple case without gas convection, the vapour pressure computed drops far below its equilibrium value, justifying a non-equilibrium description. When compared to experimental kinetics, a fair agreement is obtained, while some discrepancies are observed on the upper side of the column. At the interface with the atmosphere, phase change is highly activated and do not behaves as inside the sample. Further investigations are needed to develop accurate boundary conditions.

Conclusion

A macroscopic modelling of liquid-gas phase change in porous media has been developed from thermodynamics principles. In non-equilibrium conditions, the chemical potential appears to be the key parameter that governs the mass transport. This phenomenon has been experimentally characterized, emphasizing that hygroscopic effects have a significant influence by lowering the phase change velocity. Further microscopic thermodynamics developments should bring some essential physical descriptions of phase change.

Fig 1 :

 1 Fig 1: Phase change coefficient L as a function of the ratio of the vapour pressure divided by its equilibrium value.

Fig 2 :

 2 Fig 2: Phase change equilibrium coefficient L eq as a function of the water content w.


Astrium Services

Astrium Satellites

Astrium Space Transportation

The European prime contractor for civil and military space transportation and manned space activities A world leader in the design and manufacture of satellite systems

At the forefront of satellite services in the secure communications, geo-information and navigation fields

Feb. 2011 -7

Astrium: an integrated company

Man in space

The International Space Station (ISS) is an uplifting example of international co-operation and ambition not only for advances in space technology and engineering, but also for all kinds of scientific research It is organised around five major flagship pluri-annual projects, so called the "Top Five", receiving the biggest budgets, and several smaller projects (also plurianual), called "STEP".

The time frame for projects is typically 3 to 5 years.

ASTRIUM ST R&T PROGRAM This network is flexible, and our will is to extend it, according to the evolution of our needs.

NEEDS IN MATHEMATICS

Research & Partnerships in Mathematics

In the domain of mathematics, we have already some partnerships, covering at least partly the topics presented in the following pages

We are actively working on these topics as they are key for most of our major future programmes For some very innovative topics we are using the support of EADS Innovation Works research centre (IW)

  The yearly home funded ASTRIUM-ST R&T/D Program gets a budget higher than 40 millions €, shared between France and Germany

Feb

  up a large network of scientific partners, from the academic community, in many scientific domains.

Space Transportation Astrium Chief Financial Officer Institutional Relations Human Resources Strategy & Business Development Communications & Public Relations Satellites Services Central Functions Chief Technical Officer General Counsel / Corporate Secretary International Sales & Marketing

  Space systems have become key factors in defence and security systems. Astrium provides Space Defence systems, technologies and services. ASTRIUM ST works on :

	Ballistic Missiles
	Future Space Defence Systems
	Launchers Ballistic Missile Defence
	Space Police	Defence
	Orbital System & Space Exploration Space debris removal
		Propulsion & Equipment

ASTRIUM ORGANISATION

Feb. 2011 -8

Astrium Space Transportation

The European prime contractor for space transportation, orbital infrastructure and manned space ASTRIUM ST ORGANISATION CEO:

Alain Charmeau

Business Divisions:

Feb. 2011 -9

Access to space

Without launchers, no space activity would be possible Astrium provides Europe with a full, flexible launch capability Ariane 5 Soyuz Rockot ASTRIUM AT WORK Feb. 2011 -10 A safer world ASTRIUM AT WORK Feb. 2011 -11

Q & A

Questions & Answers

Any questions?

We would be happy to answer them! More details are available on: www.astrium.eads.net You may also contact directly : gerald.pignie@astrium.eads.net


J Benedicto

A Moreau

R Pollès

E Centeno

Comment on "Negative refraction in 1D

Keywords: A. Photonic crystal, D. Negative refraction, D. Group velocity

published or not. The documents may come

In a recent publication [START_REF] Srivastava | Negative refraction in 1d photonic crystals[END_REF], Srivastava et al. explain why they expect negative refraction to occur in the gap of a one-dimensional photonic crystal with a periodicity along the optical axis (see figure 1). Negative refraction occurs when the velocity group in the x direction is negative. Following [START_REF] Srivastava | Negative refraction in 1d photonic crystals[END_REF] we have established the analytical expression of that velocity, using the same notations :

V gx = c 2 β w ε 1 a k 1 + 1 2 k 1 k 2 + k 2 k 1 ε 2 b k 2 sin(k 1 a) cos(k 2 b) + ε 2 b k 2 + 1 2
k 1 k 2 + k 2 k 1 ε 1 a k 1 cos(k 1 a) sin(k 2 b) + 1 2 k 2 ε 1 k 1 - k 1 ε 2 k 2 1 k 2 2 - 1 k 2 1 sin(k 1 a) sin(k 2 b) -1 × a k 1 + 1 2 k 1 k 2 + k 2 k 1 b k 2 sin(k 1 a) cos(k 2 b) + b k 2 + 1 2
k 1 k 2 + k 2 k 1 a k 1 cos(k 1 a) sin(k 2 b) + 1 2 k 2 k 1 - k 1 k 2 1 k 2 2 - 1 k 2 1 sin(k 1 a) sin(k 2 b) . (1)
The latter expression is defined even in the gap of the photonic crystal, where it may present negative values. Since the transmission coefficient [START_REF] Uitham | Transmission coefficient of a one-dimensional layered medium from a light-path sum[END_REF] is not null in some gaps for a finite structure, the authors of [START_REF] Srivastava | Negative refraction in 1d photonic crystals[END_REF] conclude that negative refraction can occur when expression 1 is negative.

Using a numerical tool some of us have developed and which is freely available [START_REF] Krayzel | Simulation and analysis of exotic nonspecular phenomena[END_REF], we have computed the propagation of a gaussian beam in the structure when negative refraction is expected. The result is presented figure 2. It is obvious that the lateral shift of the transmitted beam is positive, so that there is definitely no negative refraction in this situation, although we have checked that expression 1 is actually negative. This underlines that in the gap the group velocity is not well-defined and its expression should then be considered meaningless.

We have finally two arguments which strengthen the idea that negative refraction cannot occur for such a structure.

First, when a finite structure with a negative refractive index is illuminated by a beam, this beam undergoes a negative lateral shift [START_REF] Moreau | Leaky modes of a lefthanded slab[END_REF]. This shift is linked to the opposite of the derivative of the transmission coefficient's phase (Artmann's formula). Since the phase of the transmission coefficient (see figure 3) is always decreasing, then the transmitted beam is always positively shifted and not negatively, as would happen in the case of a negative refraction. That result is in agreement with previous works on the lateral shift of beams reflected by one-dimensional photonic crystals [START_REF] Felbacq | Goos-haenchen effect in the gaps of photonic crystals[END_REF].

Second, when the group velocity is well-defined (outside a gap), it is related to the Poynting vector along the x axis by

V gx = P x ǫ d dz (2)
where ǫ d is the energy density and P x is the averaged Poynting vector [START_REF] Scalora | Negative refraction and sub-wavelength focusing in the visible range using transparent metallo-dielectric stacks[END_REF][START_REF] Gralak | Anomalous refractive properties of photonic crystals[END_REF]. After a simple calculation, the averaged Poynting vector can be written in TE polarization

P x = 1 2 k x ω µ 0 a 0 | E y 1 | 2 dz (3) + 1 2 k x ω µ 0 a+b a | E y 2 | 2 dz. (4)
Since, (

1 2 kx ω µ 0 | E y 1 | 2 ) 0 et ( 1 2 kx ω µ 0 | E y 2 |
2 ) 0, P x is always positive and so is the velocity group V gx , as long as it is well-defined. That is the case in TM polarization as well. Negative refraction can thus not occur outside the gap. Inside a gap, the Poynting vector cannot be averaged, but it is linked to the lateral shift. A positive Poynting vector whatever the layer means a positive lateral shift, as mentioned above.

  Comment on "Negative refraction in 1D photonic crystals" [Solid State Communications 147 (2008) 157-160] Benedicto Jessica, Rémi Pollès, Antoine Moreau, Emmanuel Centeno To cite this version: Benedicto Jessica, Rémi Pollès, Antoine Moreau, Emmanuel Centeno. Comment on "Negative refraction in 1D photonic crystals" [Solid State Communications 147 (2008) 157-160]. Solid State Communications, 2011, 151, pp.157. 10.1016/j.ssc.2010.12.026. hal-00582962

Figure 1 :

 1 Figure 1: As in [1], we study a stack of alternating layers. Layer A has an optical index of 1 (ε 1 = 1) and has a thickness of a=0.75 in arbitrary units. Layers of type B present an optical index of 2.3 (ε 2 = 5.29) and are b=0.25 thick. The finite structure contains 10 periods. The wavelength is of 1.0744 to get a negative group velocity.

Figure 2 :

 2 Figure 2: (Color online) Modulus of the field in a structure alternating a layer of optical index 1 and thickness 0.75 and a layer of optical index 2.3 with a thickness 0.25. The structure contains 10 periods. The wavelength is 1.0744, and the waist is 25 λ.

Figure 3 :

 3 Figure 3: Phase of the transmission coefficient in the same conditions.

Preprint submitted to Solid State Communications

We finally came to the conclusion that the purely dielectric structure described in [1] cannot produce negative refraction, contrary to metallodielectric structures [6] or to Bragg mirrors oriented in another direction [8,9].


Till Isambert

Jean-Paul Berroir

Isabelle Herlin

A multi-scale vector spline method for estimating the fluids motion on satellite images

Satellite image sequences visualize important patterns of the atmospheric and oceanographic circulation. Assessing motion from these data thus has a strong potential for improving the performances of the forecast models. Representing a vector field by a vector spline has been proven efficient for fluid motion assessment: the vector spline formulation makes it possible to initially select the locations where the conservation equation has to be taken into account; it efficiently implements the 2nd order div-curl regularity, advocated for turbulent fluids. The scientific contribution of this article is to formulate vector splines in a multiscale scheme, with the double objective of assessing motion even in the case of large displacements and capturing the spectrum of spatial scales associated to turbulent flows. The proposed method only requires the inversion of a band matrix, which is performed by an efficient numerical scheme making the method tractable for large satellite image sequences.

Introduction

Meteorological and oceanographic forecasts are based on the numerical solving of 3D circulation models that require an initial condition of their state variables, among which velocity, i.e. winds and currents. The estimation of the initial condition relies on real measurements ingested into models by means of data assimilation [START_REF] Korotaev | Retrieving Ocean Surface Current by 4D Variational Assimilation of Sea Surface Temperature Images[END_REF]. Operational forecast models make use of in situ observations for this purpose, provided by ground stations, balloons, drifting buoys, etc. One challenge of environmental modelling is the ability to derive motion estimates from dynamic satellite data in order to complement in situ data. Atmospheric and oceanographic satellite images constitute a 2D visualization of a 3D fluid flow. The link between the 2D motion and the 3D flow is non trivial: for instance, on water vapor meteorological images, the pixel value results from contributions of water particles in the vertical column, hence the 2D image motion is a complex combination of real 3D motions. Correlation-based methods constitute the current operational standard for assessing image motion from satellite imagery (the so-called Atmospheric Motion Winds). Their applicability to environmental forecast is however limited owing to the lack of reliable tracers identifiable on images and to uncertainties affecting their height assignment. Despite these limitations, image motion fields are considered to be an essential data source for environmental forecast.

Assessing image motion requires defining an image processing model, constituted of a 2D conservation law and regularity properties. State-of-the-art models for fluid flow estimation are characterized by the three following properties. [START_REF] Korotaev | Retrieving Ocean Surface Current by 4D Variational Assimilation of Sea Surface Temperature Images[END_REF] The conservation equation applies either to the luminance [START_REF] Horn | Determining optical flow[END_REF] or to the mass [START_REF] Béréziat | A generalized optical flow constraint and its physical interpretation[END_REF][START_REF] Wildes | Physically based fluid flow recovery from image sequences[END_REF]. [START_REF] Horn | Determining optical flow[END_REF] The regularity is modelled by the 2nd order div-curl constraint, which provides direct control of the Helmholtz decomposition of the motion field in terms of divergence and vorticity [START_REF] Gupta | Stochastic models for div-curl optical flow methods[END_REF]. A recent study [START_REF] Ruhnau | Optical Stokes Flow Estimation: An Imaging-based Control Approach[END_REF] formulates the regularity with effective physical constraints: the retrieved flow is supposed steady and satifying the incompressible Stokes equations. The applicability to non steady flows satisfying the Navier-Stokes equations remains however open. (3) Multiscale coarse-to-fine numerical schemes [START_REF] Anandan | A computational framework and an algorithm for the measurement of visual motion[END_REF][START_REF] Bergen | Hierarchical model-based motion estimation[END_REF][START_REF] Enkelmann | Investigation of multigrid algorithms for the estimation of optical flow fields in image sequences[END_REF][START_REF] Moulin | Multiscale modeling and estimation of motion fields for video coding[END_REF] are implemented in the case of large displacements, that prevent the linearization of the conservation equation. The multiscale scheme is further needed in the case of highly turbulent flows in which motion patterns co-exist at different spatial scales.

Solving the image processing model is usually addressed by a gradient-based minimization of an energy functional made up of two components: the confidence in data and the regularity of the result. The data confidence is often the integral over the image domain of the residual of the conservation equation. Applying such an approach to satellite atmospheric or oceanographic sequences faces two main difficulties. First, the conservation equation can be locally incorrect -3D motions such as upwelling or convection constitute typical examples in which the 2D conservation equation is not valid-or unworkable -missing data, occlusions by clouds, absence of contrast, motion parallel to contours. Accounting for the conservation equation in the whole image domain therefore leads to erroneous flow estimates. Second, the iterative minimization of the 2nd order div-curl regularity constraint leads to 4th order PDEs, hence to complex iterative minimization not guaranteeing the convergence to a minimum.

An elegant solution to these two problems is provided by spline-based methods. First, they make it possible to easily formulate models in which data confidence applies only at selected locations: the so-called control points. Second, they do not require an iterative energy minimization. Vector-valued thin-plate splines minimize the 2nd order div-curl regularity constraint and have been proved efficient to interpolate or approximate wind measurements provided by ground stations [START_REF] Amodei | A vector spline approximation[END_REF]. The use of thin-plate vector splines for image motion estimation has been further proposed by [START_REF] Suter | Motion estimation and vector splines[END_REF], who formulated a vector spline model in which the 2D conservation law is satisfied at control points. Thin-plate splines are, however, unsuitable for multiscale estimation: they are defined from a harmonic radial basis function that increases with the distance to the control point and hence does not provide a local information.

This paper presents an innovative multiscale formulation of the vector spline approach for motion estimation. It keeps the nice properties of vector splines (use of control points, div-curl regularity) but makes use of a parametric spline model based on a compactly supported and rapidly decaying radial basis function, thus adapted to multiscale representation. The solution is obtained by solving a sparse and well-conditioned linear system. The motion is computed on a pyramidal representation of images, as the sum of a coarse scale motion and increments from one scale to the immediately finer one. Results are presented to demonstrate the effectiveness of the characteristics of the multiscale vector spline: use of control points, div-curl regularity and multiscale coarse-to-fine motion estimation.

This paper is organized as follows: the section 2 recalls the vector spline theory applied to fluid motion estimation; the proposed multiscale vector spline is presented in section 3. Results are analyzed in section 4, and conclusions and prospects for future work are given in section 5.

Apparent motion estimation using vector splines

Vector splines have been initially introduced [START_REF] Amodei | A vector spline approximation[END_REF] for the interpolation and approximation of vector observations. In this context, the vector spline model is defined from: (1) a set of n control points x i in a spatial domain Ω; (2) a vector observation w i at each control point. The vector spline is solution of the following minimization problem:

   min Ω w 2 d w(x i ) = w i ∀i or: min i (w(x i ) -w i ) 2 + λ Ω w 2 d Interpolation Approximation (1)
In equation ( 1), the parameter λ of the approximating spline controls the compromise between regularity and confidence in data, w d denotes the 2nd order div-curl semi-norm, defined as:

w 2 d = α ∇div w 2 + β ∇curl w 2 (2)
It is a semi-norm which zero-space is the set of affine vector fields. It has been proven [START_REF] Amodei | A vector spline approximation[END_REF] that this minimization problem admits a unique solution: a thin-plate spline based on the harmonic radial basis function φ:

φ(x) = (128π) -1 x 4 log x (3)
The spline is defined by:

u = p(x) + n i=1 a i ( 1 α φ xx (x -x i ) + 1 β φ yy (x -x i )) + b i ( 1 α -1 β )φ xy (x -x i ) v = q(x) + n i=1 a i ( 1 α -1 β )φ xy (x -x i ) + b i ( 1 α φ yy (x -x i ) + 1 β φ xx (x -x i )) (4
) with p and q degree 1 polynomials. The coefficients of the spline (a i , b i and the 6 coefficients in p and q) are obtained by solving a linear system.

The use of vector splines for motion estimation from image data has been proposed in [START_REF] Suter | Motion estimation and vector splines[END_REF] for the luminance conservation equation, and in [START_REF] Isambert | Apparent motion estimation for turbulent flows with vector spline interpolation[END_REF] for mass conservation. In contrast to the previous formulation, only an indirect observation is provided by the conservation equation. The latter is formulated as:

Lw + I t = 0, L being a linear operator. If assuming luminance conservation, we have Lw = ∇I.w, and for mass conservation: Lw = (∇I + I∇).w. For both cases, at each control point x i , we have L i .w = -I t , L i denoting the observation operator at control point x i .

The vector spline model is then rewritten as:

   min Ω w 2 d L i w + I t = 0 or: min i (L i w + I t ) 2 + λ Ω w 2 d ( 5
)
for the interpolation and approximation cases. It has been proven [START_REF] Suter | Motion estimation and vector splines[END_REF] that the solution of ( 5) exists and is unique if the observation operators L i are linear and non zero, and if the control points are non aligned. The solution is a thin-plate spline, with the same basis function φ as in equation ( 3):

w = n i=1 c i L i φ( x -x i ) + 6 i=1 d i p i (x) (6)
p = (p i ) being the canonical basis of degree 1 polynomials. The vectors of coefficients c = (c i ) and d = (d i ) are solution of:

Q + λI d S S T 0 c d = -It 0 (7)
I d being the identity matrix. The general term of

Q (size n × n) is L i L j φ( x i - x j ) and S (size n × 6) is defined as: S = L i p.
The vector spline model used for image motion estimation shares the same interesting properties as the model used for interpolating or approximating vector data: its solution is obtained without iterative minimization and it has few parameters (i.e. α, β and λ). One will set α ≫ β when processing images of highly turbulent flows, thus favoring an estimation with large curl variability. Reversely, β ≫ α is appropriate for situations such as intense vertical activity causing locally divergent 2D flows.

One main difference with vector data is that there are no predefined locations where observations are available. The control points selection process must be defined prior to the estimation. The theoretical conditions guaranteeing the existence and uniqueness of the solution can easily be met by selecting these control points with a double thresholding: first on the magnitude of the spatial gradient, discarding low contrast areas; second on the magnitude of the motion index (ratio of temporal and spatial gradient), thus discarding visually still areas. Depending on the applicative domain, further locations in specific structures must be discarded. This is, for instance, the case of ejection filaments in oceanography, as in these structures the motion is parallel to the contours and the conservation equation is degenerated. The selection of control points is therefore highly application-dependent; they should be as evenly distributed as possible, but there is, to our knowledge, no criterion for defining an optimal distribution of control points.

Multiscale vector splines

Parametric spline model

Thin-plate vector splines minimize the 2nd order div-curl regularity, but are inappropriate for multiscale estimation as they are defined from a harmonic basis function. A multiscale scheme actually requires using a basis function that provides a local representation, hence locally supported or rapidly decaying. The contribution of this paper is to formulate a multiscale model, based on a spline parameterized by the scale value and on a pyramidal representation of images at different scales.

We consider the spline approximation problem with the 2nd order div-curl norm and either the luminance or the mass conservation equation, through the observation operators L i assessed on the n control points x i :

min J(w) = n i=1 (L i w -w i ) 2 + λ Ω α ∇div w 2 + β ∇curl w 2 (8)
Rather than exactly solving equation ( 8), which would lead to the thin-plate spline, the minimum is searched for among a set of spline functions suitable for the multiscale formalism and satisfying the two following properties. (1) The spline is defined from a unique bell-shaped radial basis function of unit support. The choice of this function is not critical as long as it is positive, decreasing and at least three times continuously differentiable in order to compute the 2nd order div-curl semi-norm. We make use of the basis function ψ proposed by [START_REF] Wendland | Piecewise polynomial, positive definite and compactly supported radial basis functions of minimal degree[END_REF] and defined as

ψ(r) = (1 -r) 6 (35r 2 + 18r + 3) for |r| ≤ 1. ( 2
)
The spline is a linear combination of translates of the basis function over a regular lattice of m grid points, whose sampling defines the scale parameter h. These translates are dilated by a factor γ proportional to h. The parameters defining the spline are the m weights q = (q j ) (each weight q j being homogeneous to a motion vector with u and v components) applied to the translates of the basis function. The parametric expression of the vector spline is thus:

w q,h (x) = vj ∈Z 2 ,hvj ∈Ω q j ψ( x -hv j γ ) (9)
where v j spans a regular lattice of unit spacing in the image domain Ω.

A new expression of the functional J is defined by substituting, in equation [START_REF] Bergen | Hierarchical model-based motion estimation[END_REF], w by its parametric form w q,h (9). Let us first consider the first term of J. If the observation operator is based on the luminance conservation equation, its new expression becomes:

I x Ψ q u + I y Ψ q v -I t 2 = A l q -I t 2 (10)
I t being the n-dimensional vector of the temporal derivatives at the control points; Ψ being the n × m matrix of general term ψ((x ik j )/γ, (y il j )/γ) with i indexing the n control points and j the m grid points (k, l); I x and I y are the n × n diagonal matrices of the image spatial derivatives at the control points. In the case of mass conservation, the first term of J becomes:

I x Ψ q u + I y Ψ q v + ID x Ψ + ID y Ψ -I t 2 = A m q -I t 2 ( 11
)
where I is the n × n diagonal matrix formed by the image values at control points, D x Ψ and D y Ψ are the matrices of the spatial derivatives of Ψ . Whatever the conservation equation, the first term of J is then rewritten as a quadratic function of q.

Let us now analyze the second term of J. By introducing the matrix of differential operators Q(D):

Q(D) = √ α∂ xx + √ β∂ yy ( √ α - √ β)∂ xy ( √ α - √ β)∂ xy √ α∂ yy + √ β∂ xx (12)
J then factorizes as:

α ∇div w 2 + β ∇curl w 2 = Q(D)w 2 (13)
The second term of J is finally rewritten as the quadratic expression Rq 2 , with:

R = √ α∂ xx Ψ + √ β∂ yy Ψ ( √ α - √ β)∂ xy Ψ ( √ α - √ β)∂ xy Ψ √ α∂ yy + √ β∂ xx Ψ . (14)
The substitution of w by the parametric expression w q,h allows J to be rewritten as a quadratic function of q:

J(q) = Aq -I t 2 + λ Rq 2 (15)
with A being either A l or A m depending on the conservation equation chosen.

Finding the minimum of J with respect to q is now a linear optimization problem.

The matrices A and R, in [START_REF] Corpetti | Dense estimation of fluid flows[END_REF], have a band structure since ψ has a compact support of size γ. The width of the band depends on the ratio of γ to the scale parameter h. If γ is smaller than h, the matrices A and R are diagonal and the vector spline is zero everywhere except in the vicinity of the grid points. If γ is large compared to h, the resulting vector spline can accurately approximate the thin-plate spline, but the A and R matrices are dense and require a heavy computational load. γ = 3h has been empirically chosen as a good compromise between the computational speed and the accuracy of the spline. The band structure allows an efficient numerical solving to be implemented.

Hierarchical motion estimation

A multiscale scheme is required for two main reasons. [START_REF] Korotaev | Retrieving Ocean Surface Current by 4D Variational Assimilation of Sea Surface Temperature Images[END_REF] The parametric spline allows the image motion to be assessed, given a spatial scale parameter h, and provided that the conservation equation can be computed. On satellite image sequences, a too strong motion and/or a too coarse time sampling cause large displacements between successive frames, preventing the linearization of the conservation equation. ( 2) Turbulent flows are associated with a large spectrum of spatial and temporal scales. We therefore make use of a pyramidal scheme, in which motion is hierachically computed from the coarsest to the finest scale.

Let I 0 and I 1 be two successive images of the sequence. Both are represented using a pyramid, from the full resolution I 0 (0) and I 1 (0) to the coarsest scale I 0 (p max ) and I 1 (p max ). To each index p corresponds a scale parameter h(p). The motion is initially computed at the coarsest scale with the parametric spline at scale h(p max ), yielding the motion field w(p max ). This initial coarse motion field is then progresively refined at each scale h(p) by first compensating the image I 0 (p) with w(p + 1) and computing the motion increment δw(p) between the compensated image and I 1 (p). The finest scale motion (p = 0) is thus expressed as the sum of the coarse scale motion w(p max ) and of the increments describing the finer resolutions:

w(0) = w(p max ) + 0 p=pmax-1 δw(p) (16)
The link between the scale parameter h(p) and the real spatial scale of the evolving image structures is not obvious: at one level of the pyramid, the motion is computed using a scale parameter h(p) corresponding to a basis function of support γ = 3h(p). The basis function is thus able to represent motion patterns with spatial size less than 3h(p); but there is no guarantee that all motion patterns of that size will be represented: this will occur only if enough control points have been selected in the existing patterns. The first result intends to demonstrate the efficiency of accounting for the conservation only at control points. For this purpose, the motion is computed using the multiscale vector spline and compared to the result of Corpetti's method [START_REF] Corpetti | Dense estimation of fluid flows[END_REF]. Both methods minimize the second order div-curl regularity constraint, make use of either luminance or mass conservation and are solved in a multiscale scheme. The two methods differ in the data confidence term of the minimized energy (computed on control points selected by double thresholding for the multiscale spline, on the whole image domain for Corpetti's method) and in the numerical minimization scheme (multiscale vector spline vs variational minimization). Two comparisons are displayed. First, the motion is computed using the luminance conservation equation on the synthetic 'OPA' sequence (on the left in figure 1), obtained by numerical simulation with the OPA ocean circulation model 3 . The OPA sequence consists of simulated images of sea surface temperature, used for computing motion. Additionaly the corresponding surface currents are available and used as the reference field for validation purposes. The results are displayed on figure 2. The mean angular error between the estimated and reference motion fields is 28 degrees for the multiscale spline and 42 degrees for Corpetti's method. The qualitative inspection of the motion field's streamlines and vorticity suggest that the motion of vortices is better assessed by the multiscale spline. A similar comparison on a Meteosat-5 sequence 4 acquired in the water vapor band is displayed on figure 3. The mass conservation equation is used as the 2D atmospheric flow can be considered as compressible to accomodate the effects of vertical motion. For this sequence, a sole qualitative assessment of results is possible. The multiscale spline is more accurate with respect to the location of the central vortex. It furthermore succeeds in capturing a rotating motion in the lower left part of the image, whereas Corpetti's method incorrectly computes a smooth laminar field.

Results

The second comparison is intended to demonstrate that the 2nd order div-curl regularity must be preferred to L 2 regularity for fluid motion assessment. The luminance conservation equation is considered and the motion is computed on the OPA sequence by the multiscale spline and the Horn and Schunck method [START_REF] Horn | Determining optical flow[END_REF]. The results are displayed on figure 4. Three different results are presented corresponding to different values of the λ coefficient assigned to the regularity component, so that both methods are tested with low, medium and high regularization. The angular errors for the multiscale spline are 30, 29 and 28 degrees (respectively for low, medium and high regularity), for the Horn and Schunk method 43, 47 and 49 degrees. The spline method is much more efficient as far as the detected location of eddies is concerned: only one vortex is detected by H&S method with low regularity, and none with medium and high regularity. This is a consequence of the L 2 regularization which favours laminar fields.

Figure 5 displays the motion fields estimated on the OPA sequence at three different scales. At the coarsest scale, the main vortices appear in the upper part of the image, and the large vortex in the bottom part is not detected at all. At the intermediate scale, more vortices appear. At finest resolution the location of vortices is improved and the large vortex in the bottom part of the image is even detected. This illustrates that the multiscale scheme actually links the size of the spatial structure with the spatial scale of the spline, although this link is not easy to interpret.

Conclusion and future work

This paper proposes an innovative approach for assessing motion on atmospheric and oceanographic satellite image sequences of fluid flows. Its characteristics are the following. (1) The data confidence is only taken into account at specific locations, the control points, where the conservation equation is workable, i.e. valid in 2D, or not degenerated to 0=0. The comparison made with a method using the same conservation equation and the same regularity constraint, but computing the conservation equation over the whole image domain, proves that using control points is essential for assessing the confidence in data. [START_REF] Horn | Determining optical flow[END_REF] The regularity of the motion field is controlled by the 2nd order div-curl constraint. This is especially adapted to fluid and turbulent flows, as it allows the variations of divergence and vorticity of the retrieved motion field to be controlled. This is highlighted by comparisons with methods based on L 2 regularization, which are unable to correctly assess the motion field in eddies. (3) A multiscale scheme is used to to allow estimation even in the case of large displacements and to provide a hierarchical representation of the motion field, distinguishing between large and fine scale structures such as vortices.

Being based on control points and on the 2nd order div curl constraint, the motion field should theoretically be expressed by a thin-plate spline. The theoretical contribution of this paper is to formulate a vector spline in a multiscale scheme, using a parametric spline model. This approach has two main advantages: (1) the minimum of the energy is obtained by solving a linear system; [START_REF] Horn | Determining optical flow[END_REF] the matrix to invert has a band structure, allowing an efficient implementation that makes the method affordable even for large sequences.

The approach requires further investigation on two issues. First, what is the link between the spatial scale of basis functions and the extent of image structures such as vortices? There is no obvious link according to the results and we believe that the control points must also be hierarchically organized for that purpose. Second, it is well known that spatial and temporal scales of turbulent flows are strongly linked, large vortices being stabler than smaller ones. It is thus necessary to consider motion estimation not only from a pair of successive images, but from a longer sequence in order to establish a multiscale representation both in space and time.

Fig. 1 .

 1 Fig. 1. Extract of the test sequences: left, OPA; right: Meteosat.

Fig. 2 .

 2 Fig. 2. Motion fields estimated on the OPA sequence using luminance conservation. Left to right: reference motion, multiscale spline, Corpetti and Mémin. Top to bottom: motion field, streamlines, vorticity.

Fig. 3 .

 3 Fig. 3. Motion fields estimated on the Meteosat sequence using mass conservation. Left: multiscale spline, right: Corpetti and Mémin. Top to bottom: motion field, streamlines.

Fig. 4 .

 4 Fig. 4. Comparison of motion estimates on the OPA sequence, luminance conservation. Top: multiscale spline with low, medium and high regularity; bottom: Horn and Schunk with low, medium and high regularity.

Fig. 5 .

 5 Fig. 5. Motion field estimated on the OPA sequence, from the coarsest to the finest (full) resolution.

Thanks to Marina Levy, LOCEAN, IPSL, France