Cédric Montero
email: cedric.montero@univ-montp2.fr

Bruno Clair

Tancrède Alméras

Arie Van Der Lee

Joseph Gril

Relationship between wood elastic strain under

Keywords: wood, cellulose, X-ray diffraction

HAL is

Introduction

Wood as a composite material

Wood is one of the oldest and the most abundant natural composite material on earth. During hundred of millions years, it has been optimized by biological evolution to provide trees with their needs for mechanical support, leading to a very sophisticated structure with high mechanical performance [START_REF] Salmén | Cell wall features with regard to mechanical performance . A review COST Action E35 2004 -2008 : Wood machining -micromechanics and fracture[END_REF]. Its use as a raw or transformed material is still increasing today in the context of growing concerns for environmental issues. Studying the origin of its mechanical properties basically aims at promoting and improving the use of this cheap and abundant resource. Since many principles of composite material design (cellular, multi-layered, fiber-reinforced, prestressed), exist in wood at different levels and are applied in an infinity of varieties depending on the tree species and the wood type, understanding the relationships between its ultra-structure at different length scales and its mechanical behaviour is also of fundamental interest for "bioinspiration" [START_REF] Martone | Mechanics without Muscle: Biomechanical Inspiration from the Plant World[END_REF]. Wood structure and properties Wood is a cellular material made of multi-layered elongated cells (typically 1 mm long and 20 µm diameter). In between the cells, a middle lamella essentially made of lignins and pectins acts as a glue joint. The cell wall is divided into a primary wall and a secondary wall sub-divided into 3 layers (S1, S2, S3). The material of which wall layers are made can be described as fibre-reinforced nanocomposite. Its fibre phase, made of partly crystalline cellulose microfibrils having nanometric width, is embedded in a matrix of lignin and hemicellulose. As the microfibrils represent an important fraction of the wall (~50%) and have very high axial stiffness (>130 GPa), their organisation is determinant for the mechanical performance of the cell-wall. In each layer, cellulose microfibrils are generally not parallel to the cell axis, but organized as a helicoid around the cell, oriented at a specific angle. The S2 layer is the more cellulosic and by far the thickest of the secondary layers, thus its microfibril angle (MFA) mainly governs the mechanical properties in the fibre direction [START_REF] Cave | The anisotropic elasticity of the plant cell wall[END_REF]. A detailed description of wood multi-scale structure including graphical representation can be found in literature [START_REF] Eitelberger | Prediction of transport properties of wood below the fiber saturation point -A multiscale homogenization approach and its experimental validationPart I: Thermal conductivity[END_REF][START_REF] Salmén | Micromechanical understanding of the cell-wall structure?[END_REF]. Wood micro-mechanics Several models have been proposed to predict the mechanical properties of wood from the organization and properties of cell wall components [START_REF] Salmén | Micromechanical understanding of the cell-wall structure?[END_REF][START_REF] Salmén | A model for the prediction of fiber elasticity[END_REF][START_REF] Harrington | Modelling the elastic properties of softwood -Part I: The cell-wall lamellae[END_REF][START_REF] Hofstetter | Micromechanical modeling of solid-type and plate-type deformation patterns within softwood materials. A review and an improved approach[END_REF]. These models generally account for the cellular nature of wood, the composite nature of the cell wall and the effect of microfibril orientation. As wood cells have a very high aspect ratio, the behaviour of wood in its longitudinal direction is generally approximated as a parallel association of cell wall layers. Each cell wall layer is considered as a composite medium where microfibrils are also associated in parallel with the matrix phase. These assumptions implicitly consider infinitely long and straight fibres and/or microfibrils, and perfect cohesion between the components. In order to evaluate these assumptions, it is interesting to quantify the contribution of the individual components on the macroscopic behaviour. The behaviour of lignins and hemicelluloses was for a long time inaccessible from strain measurement until recent progress in spectroscopy tools made it possible to monitor strain of specific components [START_REF] Åkerholm | Interactions between wood polymers studied by dynamic FT-IR spectroscopy[END_REF]. Cellulose behaviour on the other hand is easier to measure thanks to its crystalline structure. Indeed, X-ray diffraction (XRD) patterns provide information about the orientation of cellulose in the cell wall (MFA) [START_REF] Yamamoto | Method of Determining the Mean Microfibril Angle of Wood over a Wide-Range by the Improved Cave's Method[END_REF][START_REF] Cave | Theory of X-ray measurement of microfibril angle in wood[END_REF][START_REF] Ogurreck | Analytical description of the scattering of cellulose nanocrystals in tracheid wood cells[END_REF], but also about the state of strain of the cellulose crystal, through the analysis of changes in lattice spacing in response to mechanical loads [START_REF] Suzuki | Mechanical deformation of crystal lattice of cellulose in Hinoki wood[END_REF][START_REF] Sobue | X-Ray measurement of lattice strain of cellulose Crystals during the shrinkage of wood in the longitudinal direction[END_REF]. XRD, crystal strain and micromechanics XRD has been used in composite materials to quantify the respective contributions of the matrix and the reinforcing fibres to the apparent strain (e.g. [START_REF] Young | Analysis of stress transfer in two-phase polymer systems using synchrotron microfocus X-ray diffraction[END_REF]). A number of studies reported the strain of cellulose measured by XRD during mechanical tests performed on wood. Because cell-walls are fibre-reinforced composites, the axial strain of the microfibril is of particular interest. It can be observed through the change in position of the (004) reflection of the cellulose crystal. This method has been used to localize elastic strain during water sorption or desorption [START_REF] Sobue | X-Ray measurement of lattice strain of cellulose Crystals during the shrinkage of wood in the longitudinal direction[END_REF][START_REF] Clair | Mechanical behavior of cellulose microfibrils in tension wood, in relation with maturation stress generation[END_REF][START_REF] Zabler | Moisture changes in the plant cell wall force cellulose crystallites to deform[END_REF], to measure the crystal strain during internal stress release [START_REF] Clair | Mechanical behavior of cellulose microfibrils in tension wood, in relation with maturation stress generation[END_REF], and to monitor the change in mechanical stress during wood formation [START_REF] Clair | Maturation stress generation in poplar tension wood studied by synchrotron radiation microdiffraction[END_REF]. Experiments based on simple tensile tests [START_REF] Suzuki | Mechanical deformation of crystal lattice of cellulose in Hinoki wood[END_REF][START_REF] Suzuki | Relation between Young's modulus and the cell wall structures of Sugi (Cryptomeria japonica D. Don)[END_REF][START_REF] Nakai | The relationship between macroscopic strain and crystal lattice strain in wood under uniaxial stress in the fiber direction[END_REF][START_REF] Nakai | Mechanical behavior of the crystal lattice of natural cellulose in wood under repeated uniaxial tension stress in the fiber direction[END_REF][START_REF] Peura | The effect of axial strain on crystalline cellulose in Norway spruce[END_REF] provide estimations of cellulose contribution to wood deformation. These studies revealed that the ratio between crystal strain and wood strain is always lower than 1, with a value generally ranging from 0.5 to 0.95. Only one study reports the behaviour of cellulose during compression tests [START_REF] Nakai | The relationship between macroscopic strain and crystal lattice strain in wood under uniaxial stress in the fiber direction[END_REF], where the strain ratio ranged from 0.25 to 1. The analysis of this ratio can bring useful information about the deformation process of a material at the ultra-structural scale. In the case of relatively simple composites, interpretation of these data can be rather straightforward in terms of strain localization, because only one change in scale is involved between the crystal scale and the macroscopic scale (e.g. [START_REF] Young | Analysis of stress transfer in two-phase polymer systems using synchrotron microfocus X-ray diffraction[END_REF]). In the case of wood, however, many level of organisation are involved between the crystal (nanometre scale) and the macroscopic material (centimetre scale). Interpretation of the relation between crystal strains and wood strains can be given only under explicit assumptions about the mechanisms involved at each scale. For example, under the classical assumptions made in wood micro-mechanical models (infinitely long fibres and microfibrils, perfect cohesion between components and between layers) the strain ratio is expected to be related to the MFA by simple geometric considerations. Departure from this relation would indicate that some of the above-mentioned assumptions are not correct. Similarly, the strain ratio is expected to be the same under tensile and compressive loading, provided wood is loaded in its elastic domain.

Objectives

In order to check the real contribution of crystalline cellulose to wood elastic behaviour, a precise estimate of this strain ratio must be obtained on the same specimen both in tension and compression. This has been for the first time achieved here through the simultaneous observation of the two sides of a wood specimen submitted to a bending load. Here we present a four-point bending device designed for in-situ XRD, allowing to measure with sufficient accuracy during a single test the behaviour of crystalline cellulose both under tension and compression. It aims at quantifying the contribution of cellulose to wood elastic strain and checking its linearity.

Material and method

Tested material

Experiments were performed on poplar (Populus spp.) wood (air-dry specific gravity approx. 0.55). Specimens were taken far enough from the pith to neglect ring curvature and machined at a size of 2×3.2×50 mm 3 along radial, tangential and longitudinal direction (R, T, L) respectively. Experiments were performed in ambient conditions (22 ± 0.4 °C and 49 ± 7%RH). Specimens were taken in mature wood far enough from the pith to neglect ring curvature. Their orientation was chosen to get the direction of deflection parallel to the T direction and the X-ray beam perpendicular to the T direction, in order to minimize heterogeneity along strain profiles. The thickness along R was optimized to maximize the diffraction signal (trade-off between the amount of diffracting material and signal attenuation). The specimen height along T was set to maximize the measurement zone while keeping acceptable R/T ratio for bending tests.

Bending tests

Specimens were tested in four-point bending in joist orientation with 20 mm inner span and 40 mm outer span in a device specifically designed for the X-ray goniometer (fig. 1). The screw-induced displacement of the inner holders applies pure bending on the specimen, inducing a linear strain gradient with compression on the upper side and tension on the lower side. Strain was recorded using two 5 mm-long strain gages (Kyowa KFG-5-120-C1-11-L1M3R) glued with cyanoacrylate adhesive at the mid-span of the specimen, on the tension and compression surfaces. Preliminary tests made on matched specimens showed that the linear elastic limits of the material were ε - = -0.22% in compression and ε + = 0.3% in tension. Different strain levels (0.2%, 0.3%, 0.4%) were applied to the specimens in this experiment. Experimental set-up for X-Ray diffraction Diffraction patterns were obtained with an X-ray diffractometer (Agilent Technologies Gemini S) using the CuKα source (wavelength λ = 0.154 nm, operated at 50 kV and 30 mA) with a 0.8 mm collimator. The exposure time was set to 40 seconds to optimize the diffraction signal and the total test duration (approx. 30 min each profile). The detector was positioned at a distance D = 135 mm to the specimen, as to focus the signal acquisition on the first (004) reflection and maximize the resolution of the obtained data. The specimen was asymmetrically positioned according to the Bragg condition at 17° to maximize the intensity of the (004) reflection. The detector resolution was 2048 x 2048 pixels (31 µm/pixel) (fig. 2). The specimen was scanned from top to bottom both in the initial and deformed state with a step Δy = 0.2 mm. A diffraction pattern focused on the (004) plane was recorded at each step. The beam being fixed, the scan is performed by vertical translation of the bending device. The relative position y of the device is measured by a displacement sensor (Mitutoyo Absolute Digimatic Inticator) with 1 µm accuracy. Because the specimen is also translated during the bending, the position of the X-ray beam relative to the border of the specimen has to be precisely determined for each scan to compute a coherent strain field. This was done using a specific method based on the interpretation of the signal near the specimen border (see Suppl. Mat.). As a validation, the specimen height h is determined with this interpretation and compared to measurements made with the calliper. The specimen deflection after loading was also measured and found consistent with the measured strain under the assumptions of beam theory.

X -r a y b e a m

Determination of lattice spacing d 004

The cellulose crystal strain is measured by computing the displacement of the first reflection of the (004) crystal plane (which is normal to the cellulose microfibril axis). The theoretical (004) lattice spacing of cellulose is d 004 th = 0.258 nm, corresponding to a scattering angle of 2θ 004 th = 34.73° for the considered wavelength. Slight changes in d 004 due to mechanical strain induce a slight change in θ 004 and therefore a displacement of the diffraction peak. The lattice distance was calculated using Bragg's equation linking the wavelength of the X-ray source λ, to the lattice distance d 004 and the scattering angle θ 004 : λ = 2d 004 sinθ 004 (Eq. 1)

The exact scattering angle θ 004 is obtained from its geometric relation with H 004 (distance between the diffraction peak and beam centre measured on the detector), and D (distance between the specimen and the detector): tanθ 004 =H 004 /D (Eq. 2)

In this analysis, the precision of the lattice spacing estimation critically depends on the precision obtained on these geometrical parameters. The relevant distance D to be taken into account is the distance between the detector and the centre of mass of the area of the specimen crossed by the beam. This distance must be calibrated each time a diffractogram is recorded, because the surface of the wood specimen is never perfectly plane and vertical and also because the wood specimen may slightly move out of plane during the experiment. For this purpose, a reference crystal material (spherical gold powder, particle size diameter: 0.5-0.8 µm, d Au 111 = 0.235 nm) was laid on the external surface of the specimen. The precise distance D g between the specimen surface and the detector was determined from the position of the gold reflection. The distance D was defined as D = D g + R/2 (with R the specimen width). The azimuth distribution of the (200) cellulose reflection was used for the determination of the microfibril angle of the specimen using Cave's method [START_REF] Cave | Theory of X-ray measurement of microfibril angle in wood[END_REF]. The radial position of the (004) reflection was used to determine parameter H 004 . In these experimental conditions, the displacement of peak position H 004 was only a few pixels. The following procedure was applied to get sub-pixel determination of this parameter: the diffractograms were integrated along the azimuth to obtain 1-D radial profiles between 2θ=33° and 2θ=40° (fig. 3) and filtered to remove outlying data (Zinger spots due to cosmic radiation). The abscissa was converted into lattice spacing using equations ( 1) and ( 2), and a second-order polynomial was fit to the data around the peak position. The precise position of the maximum, directly giving the lattice spacing d 004 , was then computed from the polynomial coefficients. Computation of strain profiles The macroscopic strain field along specimen height h is assumed to be plane according to beam theory. Therefore the strain at any distance y from the upper surface of the specimen can be linearly interpolated from the gage measurements ε + and ε -:

ε wood (y) = ε --(ε + -ε -)y/h (Eq. 3)
The profile of crystal strain ε cryst (y) is obtained from the estimations of the lattice spacing in initial and deformed states at position y:

ε cryst (y) = d 004 def (y)/d 004 ini (y) -1 (Eq. 4)
The occurrence of a small amount of stress relaxation during the test was minimized by the short duration of the test and neglected in further analyses.

Results

Profiles of lattice spacing and crystal strain

The recorded profiles of lattice spacing in initial and deformed states are shown in figure 4. Some variations of lattice spacing can be noticed along the initial profile. Variations near the border of the specimen are due to a weaker signal obtained at this level and other possible artefacts due to the state of wood near the surface and the presence of the strain gage. Data at these positions are removed in further analysis. Variations within the specimen are clearly observed but these sources of disturbance are independent of the applied stress so that they are also present in the deformed state, and will be eliminated when subtracting the lattice distances to obtain the strain. The lattice profile in the deformed state clearly differs from that in initial state: the lattice distance is reduced on the upper part of the specimen corresponding to the compressed side, and increased on the lower part corresponding to the part under tension. Cellulose crystal strain profiles computed from these data are shown in figure 5, together with the profile of macroscopic strain. The proximity between the middle of the specimen height and its neutral line (for which the strain is zero) confirms symmetrical specimen loading in tension and compression. The crystal strain is clearly negative in the upper part, positive on the lower part, and zero in the middle of the specimen. Although local variations can be noticed, the profile appears approximately linear. A comparison with the macroscopic strain shows that, at any position, the cellulose crystal strain is slightly lower in magnitude than the wood strain. Cellulose lattice strain versus macroscopic strain The relationship between wood macroscopic strain and cellulose crystal strain is shown in figure 6. This relationship is close to linear, showing that when wood is submitted to a strain, the cellulose deforms proportionally. In order to detect a possible difference in response to compressive and tensile stresses, the average ratio Γ of crystal strain to wood strain was determined by linear regression independently for the compressed side (Γ -=0.81) and the tension side (Γ + =0.65). A numerical difference was found for this specimen, and this experiment was repeated on different poplar specimens with similar MFA to check the significance of this result. Another example of result is shown in Figure 7. The results of five tests (table 1) show that the difference between the ratios obtained from the two sides is not significant. Additional specimens were tested with larger strain up to the elastic limit of the material (see Table 1). Results for a same specimen at two strain levels are shown in figure 8. At higher strain levels, the results are similar to those shown in figures 4-5-6: the profile of crystal strain is linear along the specimen height, and its magnitude is proportional to that of wood strain, both on the compression and on the tension sides. The obtained Γ ratios were close to those obtained at a strain of 0.2% (table 1).

Fig. 8: Relationship between cellulose crystal strain and wood macroscopic strain (specimen 1) at two levels of applied strain; firstly at Δε = 0.1903% (✕) and secondly increased at Δε =0.3957% (◊). The grey lines (⎯) indicates linear regressions on data relative to the larger strain level.

Discussion

Measuring lattice strain Our results show that the strain of crystalline cellulose can be measured with sufficient precision to evaluate its contribution to macroscopic wood strain. It is noted that the profile of cellulose apparent lattice spacing is often found non-uniform inside the unloaded wood specimen. This heterogeneity of lattice spacing at the initial state is due to at least two factors: (1) A possible artefact due to heterogeneity of specimen density. Indeed, heterogeneity (e.g. due to the presence of the large cavities of vessels) induces slight changes in the position of the centre of mass of the specimen. Error propagation analysis shows that an error of 100 µm in the estimation of the position of centre of mass induces an error of 0.00014 nm on the estimation of lattice spacing (equivalent to a strain of 0.05%). ( 2) The presence of eigenstresses in the material (e.g. due to heterogeneous drying shrinkage of the specimen), leading to a heterogeneous state of cellulose strain prior to loading. In order to cope with this heterogeneity, it is necessary to carefully monitor the position of each exposition, in order to compute crystal strain by subtracting the lattice spacing measured on the same material points before and after loading. To do so, we designed a specific procedure based on the interpretation of the signal at the specimen border (see Suppl. Mat.) and calibration with gold powder.

Linearity of cellulose strain

We tested wood specimens up to their proportional stress/strain limit, and checked that, within this range, cellulose strain is proportional to macroscopic strain. This proportionality was checked in three ways: comparing different positions corresponding to different strain levels within the bent specimen; comparing the tensile and compressed sides of the specimen; comparing the strain ratio obtained for different applied strain levels.

Our conclusions are in accordance with most in situ tensile tests performed at low strain levels [START_REF] Suzuki | Mechanical deformation of crystal lattice of cellulose in Hinoki wood[END_REF][START_REF] Suzuki | Relation between Young's modulus and the cell wall structures of Sugi (Cryptomeria japonica D. Don)[END_REF][START_REF] Nakai | The relationship between macroscopic strain and crystal lattice strain in wood under uniaxial stress in the fiber direction[END_REF][START_REF] Nakai | Mechanical behavior of the crystal lattice of natural cellulose in wood under repeated uniaxial tension stress in the fiber direction[END_REF]. However, in contrast with earlier reports [START_REF] Nakai | The relationship between macroscopic strain and crystal lattice strain in wood under uniaxial stress in the fiber direction[END_REF], we show that the strain ratio does not differ significantly in tension and compression. This was made possible by the use of a bending test, where the specimen is symmetrically loaded in tension and compression at the same time.

Coincidence of the neutral line between wood strain and cellulose crystal strain confirms that they are proportional within the specimen, suggesting that any possible localised effects due to e.g. the compression at the loading points are of second-order.

Micromechanical interpretation of strain ratio

The strain ratio we obtained is nearly 0.75, which is in the range of values reported in tensile tests: Suzuki [START_REF] Suzuki | Mechanical deformation of crystal lattice of cellulose in Hinoki wood[END_REF] estimated a ratio of 0.65, and Nakai [START_REF] Nakai | Mechanical behavior of the crystal lattice of natural cellulose in wood under repeated uniaxial tension stress in the fiber direction[END_REF] reported a ratio between 0.55 and 0.85. Only Peura [START_REF] Peura | The effect of axial strain on crystalline cellulose in Norway spruce[END_REF] reports a much lower ratio (approx. 0.1), but this was obtained on very thin specimens for which the estimated specific modulus of elasticity was also unusually low.

Interpretation of this strain ratio involves at first order the geometrical effect of the MFA. Indeed, cellulose microfibrils are not aligned with the fibre axis, so that the strain in the direction of cellulose ε C is related to the strain in the principal directions of the fibre wall through the microfibril angle µ. The rigorous formula for this relation is:

ε C = c 2 ε W L + s 2 ε W T -csε W S
, where s=sinµ, c=cosµ, ε W L and ε W T are the strains in the longitudinal and transverse directions of the cell wall, respectively, and ε W S is the shear strain in the wall. Because the wood is loaded only in its longitudinal direction parallel to the cell wall axis, the transverse stress in the cell wall is neglected, so that the strain of the cell wall is related to the longitudinal strain through the Poisson's ratio (ν W LT ) of the wall:

ε W T = -ε W L ν W LT .
Following a number of authors [START_REF] Hofstetter | Micromechanical modeling of solid-type and plate-type deformation patterns within softwood materials. A review and an improved approach[END_REF][START_REF] Schniewind | Wood as a linear orthotropic viscoelastic material[END_REF][START_REF] Alméras | Modelling anisotropic maturation strains in wood in relation to fibre boundary conditions, microstructure and maturation kinetics[END_REF], we assume that shear strains are fully restrained at the cell-wall level (ε W S =0), because they are prevented by the cohesion and anti-parallel configuration of adjacent cell walls. Under this assumption, the relationship between cellulose strain and wall longitudinal strain reduces to:

ε C = ε W L (c 2 -s 2 ν W LT ) (Eq. 5)
Note that, for the strain level considered in this experiment, the change in the MFA due to the deformation of the specimen is negligible (it can be easily shown from geometrical consideration and was experimentally demonstrated by Peura [START_REF] Peura | The effect of axial strain on crystalline cellulose in Norway spruce[END_REF]), so that a unique value of µ can be used. The Poisson's ratio of the cell wall can be estimated from micro-mechanical models [START_REF] Alméras | Modelling anisotropic maturation strains in wood in relation to fibre boundary conditions, microstructure and maturation kinetics[END_REF]. For the MFA used in our experiment (µ=17°), the model estimates ν W LT ≈0.8, so that the estimation of the ratio ε C /ε W L is 0.85. This quantitative estimation of the effect of MFA shows that it is not sufficient to explain the observed values of strain ratio. Nakai [START_REF] Nakai | Mechanical behavior of the crystal lattice of natural cellulose in wood under repeated uniaxial tension stress in the fiber direction[END_REF] reports strain ratios estimated from tensile tests on specimens with different MFA, ranging from 0.3 (MFA=21°) to 0.95 (MFA=9°). Our calculations show that, for that set of observation too, the effect of the MFA was not sufficient to explain the observed values. The difference between experimental estimates of strain ratio and those expected from the geometrical effect of the MFA can be ascribed to other deformation mechanisms occurring in wood. The effect of amorphous cellulose associated in series with the crystalline material has been recognized as a possible factor by previous authors [START_REF] Nakai | Mechanical behavior of the crystal lattice of natural cellulose in wood under repeated uniaxial tension stress in the fiber direction[END_REF][START_REF] Peura | The effect of axial strain on crystalline cellulose in Norway spruce[END_REF]. Indeed, this mechanism is able to generate macroscopic strain without generating crystal strain, and could therefore partly explain why the strain ratio is lower than expected from purely geometrical considerations. This effect can be represented by an additive contribution ε A to the total strain, independent of the crystal strain ε C , as illustrated on figure 9.

Figure 9: Schematic representation of the sources of longitudinal deformation in wood. Total strain ε can be divided into a contribution from the cell-wall ε W L , that depend on the crystal strain ε C through the MFA µ, and an additive contribution ε A that does not involve the deformation of crystalline material.

The relation between total strain, crystal strain, MFA and this additive contribution is expressed by:

ε = ε C /(c 2 -s 2 ν W LT ) + ε Α (Eq. 6)
Actually, due to the complex hierarchical structure of wood, many alternative mechanisms may contribute to macroscopic strain without contributing to crystal strain. (i) At the macroscopic level, an effect of grain angle, a slight difference between the direction of wood fibres and the longitudinal direction of the specimen, would induce a projection effect similar to that of the MFA.

(ii) At a sub-millimetric scale, wood is not only made of fibres but also of a small fraction of different types of cells called ray parenchyma, oriented in the radial direction of wood, i.e. perpendicular to fibres. The loading of wood along the fibre direction probably generates complex strain fields near the areas where rays are present, possibly involving bending and shear effects that would contribute to ε A . (iii) The wood fibre itself is actually hollow and tapered, so that some bending of the wall is also likely to happen and have a similar effect. (iv) Wood fibres are glued together by an amorphous component called the middle lamella. At that level, wood can be considered as a fibre-reinforced composite where the matrix is the middle lamella and the reinforcing elements are wood fibres. In such composites, stress redistribution induces a so-called "shear-lag" effect, resulting in strain concentrations near the fibre ends and therefore smaller strain in the fibres fraction than in the bulk material, as has been shown on artificial composites (e.g. [START_REF] Young | Analysis of stress transfer in two-phase polymer systems using synchrotron microfocus X-ray diffraction[END_REF]). (v) At the level of the cell-wall, we here assumed a full shear restraint, following previous authors [START_REF] Schniewind | Wood as a linear orthotropic viscoelastic material[END_REF], but the occurrence of a small amount of shear has been considered by other authors and cannot be ruled out without further verification. Such an effect would indeed also contribute to decreasing the microfibril strain / wall strain ratio. (vi) Moreover, the fibre wall itself has a multilayer structure, and shear strain could possibly occur at the interface between layers, so that the strain at the level of a given layer could also be lower than the strain of the bulk wall. (vii) Each layer is constituted of cellulose microfibrils reinforcing a matrix of lignin and hemicelluloses. The cellulose microfibrils are neither perfectly straight nor infinitely long. As a consequence, a slight initial curvature of the microfibrils could result in bending effects that would partly allow the extension of the wall material without extending the microfibril itself. (viii) Finally, the finite length of the microfibrils induces a shear lag effect, which may also contribute to explain the difference between crystal strain and wall strain.

In order to discriminate between these mechanisms or evaluate their relative contribution, complementary observations at the relevant length scales would be necessary. Our experiment however provides an estimation of their cumulated effect through the estimation of ε A . The total contribution of these "additive" effects due to localized shear, bending deformations or extension of amorphous material is evaluated (using equation 6) near 12% of the total strain for the poplar wood we studied. This can be considered as a second-order effect regarding elastic strain, justifying why most models of elastic properties neglect these factors and concentrate on the effect of cell-wall layered structure and microfibril orientation. However, the behaviour of these amorphous compartments may be of primary importance regarding wood hygroscopic behaviour (see e.g. [START_REF] Zabler | Moisture changes in the plant cell wall force cellulose crystallites to deform[END_REF]), and possibly other non-elastic wood behaviour, such as visco-elastic creep, mechanosorptive creep or post-elastic response.

Conclusion

A methodology was developed to precisely evaluate the crystal strain of cellulose during a bending test. Results show that during elastic deformations, cellulose crystal strain remains proportional to wood macroscopic strain. The mean strain ratio measured on four specimens was close to 0.75 and did not differ significantly between tension and compression. The effect of cellulose orientation was taken into account through classical projection formula. This analysis showed that the values of strain ratio can be only partly explained by this geometrical effect, and evidenced a contribution of strains located in non-crystalline components of wood. This contribution was evaluated to be approximately 12% of macroscopic strains, and hypotheses for its origin are proposed.

Supplementary material

Measurements of lattice strain are performed by subtracting the lattice spacing recorded in deformed and initial state. Because the state of cellulose is substantially heterogeneous within a piece of wood, it is critical to record lattice spacing before and after loading at the same position within the specimen. For this, the displacement induced by the bending must be precisely accounted for. We designed a specific procedure to do so. The principle of this procedure is to use variations in signal intensity occurring near the specimen border (Fig. S1-a) to locate the relative position y C of the beam centre. Let R b be the beam radius and Δy the distance between successive shoots (with Δy<2R b , so that successive shoots overlap, see Fig. S1-b). We assume that the signal intensity is, at first order, proportional to the amount of crystalline cellulose crossed by the beam, i.e. to the area of wood S illuminated by the beam (Fig. S1-c). When a shoot is made outside of the specimen (y C >R b ) the diffracting wood area is zero (S=0). When the shoot is entirely within the specimen (y C <-R b ), the diffracting wood area constant (S=S 0 =πR b 2 ). In between, the diffracting area is directly related to the distance y C between the beam centre and the specimen border: (Eq. S1) The intensity for full diffraction I 0 corresponding to an illuminated wood area S 0 can be estimated from the data in the core of the specimen. Then, the illuminated area S for any shoot can be estimated from the corresponding intensity I: S = I*S 0 /I 0 . By solving eq. S1 numerically, the value of y C can be deduced from the value of S.

Finally, the position to be associated to a given shoot is the centre of mass of the diffracting area, y G . It can be deduced from y C using eq. S2: This procedure is based on the assumption that the amount of crystalline material (and therefore diffraction intensity) is uniform within the sample. Fig. S1-a shows that this is only approximately true within the core of the sample. The procedure happens nevertheless to be quite robust to this assumption, since the height of the specimens determined with this method was very close to that determined with a calliper (relative difference ≈ 1%). Moreover, variations in signal intensity are highly reproducible (as attested by comparison of the two curves of figure S1-a), so that if this heterogeneity induces a slight error on the determined position, the same error will be made on the two profiles, so that the calculation of strain will refer to the same material points and remains valid.

y G = f (y C ) = 2 S(y C ) y R b 2 -(y -y C ) 2

Fig. 1 :

 1 (a) Schematic view of the four-point bending device and specimen dimensions (all dimensions in mm). ε -and ε + is the strains measured with a strain gage on the upper (compressive) and lower (tensile) side, respectively. (b) Photography of four-point bending device specially designed for the X-ray goniometer head.

Fig. 2 :

 2 Fig.2: Schematics of the experimental set-up for XRD in-situ bending tests. The specimen is oriented at the Bragg angle (17°) from the normal to the x-ray beam position. D g : distance between the outer surface of the specimen and the detector. H: distance between the center of the transmitted beam and peak position of the diffraction peak.

Fig. 3 :

 3 Fig. 3: Example of radial profile of cellulose (004) and gold (111) obtained from diffractogram integration

Fig. 4 :

 4 Fig. 4: Profiles of cellulose lattice spacing along the height of specimen 1, at the initial state (⎯✕⎯) and deformed state (⎯♢⎯). The macroscopic strain applied is ε -= -0.194% and ε + = 0.184%.

Fig. 5 :

 5 Fig. 5: Profile of cellulose crystal strain (✕) and wood macroscopic strain (---) along the height of specimen 1. The gray line (⎯) indicates linear regression of the cellulose crystal strain.

Fig 6 :

 6 Fig 6: Relationship between cellulose crystal strain and wood macroscopic strain (specimen 1). ✕ experimental results, ⎯ linear regression of experimental results, ⎯ 1:1 line

Fig 7 :

 7 Fig 7: Relationship between cellulose crystal strain and wood macroscopic strain (specimen 2) ✕ experimental results, ⎯ regression line of experimental results, ⎯ 1:1 line

0

  Fig. S1: (a) Profiles of signal intensity along specimen height y, at initial and deformed states. (b) beam overlap for successive shoots (c) relationship between beam centre position and diffracting wood area.

Table 1 :

 1

	Specimen MFA (°)	Δε (%)	Г -	Г +
	1	17.7	0.1903	0.81	0.65
	1*	17.7	0.3957	0.73	0.73
	2	17.3	0.2175	0.76	0.78
	3	18.0	0.2048	0.68	0.79

Results of bending tests on different specimens. MFA: mean microfibril angle, Δε=(ε + -ε -)/2: macroscopic strain level applied, Г -(resp. Г + ): mean strain ratio in the compressive zone (resp. tensile zone).

*Same specimen at two strain levels applied.

Acknowledgements

The authors would like to thanks Mr. Gilles Camp (LMGC) for the realization of the bending apparatus. This work has been performed in the framework of the "MechWood" project supported by the Austrian Research Promotion Agency (FFG, 815234/12791) and regularly discussed during meetings of European COST Action FP0802 "Experimental and Computational Micro-Characterization Techniques in Wood Mechanics".


Yves Martinez-Maure

Rigidity and Bellows-type Theorem for hedgehogs

published or not. The documents may come

Rigidity and Bellows-type Theorem for hedgehogs

Introduction

In 1813, A. L. Cauchy proved (almost rigorously) his famous rigidity theorem: Any convex polyhedron of R 3 is rigid (that is, no convex polyhedron of R 3 can be continuously deformed so that its faces remain rigid) [START_REF] Cauchy | Sur les polygones et polyèdres, second mémoire[END_REF]. First examples of flexible polyhedra were discovered by R. Bricard in 1897 [START_REF] Bricard | Mémoire sur la théorie de l'octaèdre articulé[END_REF], but these « Bricard's flexible octahedra » are self-intersecting. The question of rigidity of embedded nonconvex polyhedra remained open until 1977 when R. Connelly discovered a first example of flexible sphere-homeomorphic polyhedron [START_REF] Connelly | A counterexample to the rigidity conjecture for polyhedra[END_REF]. In the late seventies, R. Connelly and D. Sullivan formulated the so-called « Bellows conjecture » stating that whenever we perform a rigid deformation of a flexible polyhedron P (that is, a continuous deformation of P that changes only its dihedral angles), the volume of P remains constant. This conjecture was proved by I. Sabitov for sphere-homeomorphic polyhedra [START_REF] Kh | The volume of a polyhedron as a function of length of its edges. (Russian)[END_REF] and by R. Connelly, I. Sabitov, and A. Walz for general orientable 2-dimensional polyhedral surfaces [START_REF] Connelly | The Bellows conjecture[END_REF].

In 1927, E. Cohn-Vossen proved that smooth closed surfaces of R 3 with everywhere positive Gaussian curvature are rigid [START_REF] Cohn-Vossen | Zwei Satze uber die Starrheit der Eilflchen[END_REF]. Smooth closed surfaces of R 3 with everywhere positive Gaussian curvature are also Gauss infinitesimally rigid [START_REF] Stoker | Differential geometry[END_REF], that is rigid with respect to the Gaussian curvature regarded as a function of the outer unit normal [START_REF] Izmetiev | Infinitesimal rigidity of convex surfaces through the second derivative of the Hilbert-Einstein functional II: Smooth case[END_REF]Section 2]. In this paper, we consider rigidity and Gauss infinitesimal rigidity for hedgehogs of R 3 (regarded as Minkowski differences of closed convex surfaces of R 3 with positive Gaussian curvature). As noticed by I. Izmestiev, Gauss infinitesimal rigidity can be interpreted as « infinitesimal » uniqueness in the Minkowski problem (prescribing Gauss curvature as a function of the unit normal) of which the author already studied the uniqueness part [START_REF] Martinez-Maure | Uniqueness results for the Minkowski problem extended to hedgehogs[END_REF]. In this previous work [START_REF] Martinez-Maure | Uniqueness results for the Minkowski problem extended to hedgehogs[END_REF], we presented different ways of constructing pairs of non-congruent hedgehogs that share the same curvature function (i.e., inverse of the Gaussian curvature). This will allow us to give examples of nontrivial (i.e., distinct from a point) hedgehogs that are not Gauss infinitesimally rigid. Assume we have a one parameter family of hedgehogs (H ht ) t∈[0,1] , all with the same curvature function. We do not know whether they must be congruent in R 3 . However, we shall give a bellows-type theorem for hedgehogs: under an appropriate differentiability condition of the family with respect to the parameter, we shall prove that all the hedgehogs of the family considered have the same algebraic volume.

Basic definitions C 2 -hedgehogs in R n+1
As is well-known, every convex body K ⊂ R n+1 is determined by its support function h

K : S n -→ R, where h K (u) is defined by h K (u) = sup { x, u |x ∈ K }, (u ∈ S n ),
that is, as the signed distance from the origin to the support hyperplane with normal vector u. In particular, every closed convex hypersurface of class C 2 + (i.e., C 2 -hypersurface with positive Gaussian curvature) is determined by its support function h (which must be of class C 2 on S n [15, p. 111]) as the envelope H h of the family of hyperplanes with equation x, u = h(u). This envelope H h is described analytically by the following system of equations x, u = h(u) x, . = dh u (.) .

The second equation is obtained from the first by performing a partial differentiation with respect to u. From the first equation, the orthogonal projection of x onto the line spanned by u is h (u) u and from the second one, the orthogonal projection of x onto u ⊥ is the gradient of h at u (cf. Figure 1). Therefore,

for each u ∈ S n , x h (u) = h(u)u + (∇h) (u)
is the unique solution of this system.

Figure 1. Hedgehogs as envelopes parametrized by their Gauss map

Now, for any C 2 -function h on S n , the envelope H h is in fact well-defined (even if h is not the support function of a convex hypersurface). Its natural parametrization x h : S n → H h , u → h(u)u + (∇h) (u) can be interpreted as the inverse of its Gauss map, in the sense that: at each regular point x h (u) of H h , u is a normal vector to H h . We say that H h is the hedgehog with support function h (cf. Figure 2). Note that x h depends linearly on h. In fact, we can introduce a more general notion of hedgehogs by regarding hedgehogs of R n+1 as Minkowski differences of arbitrary convex bodies of R n+1 [START_REF] Martinez-Maure | Geometric study of Minkowski differences of plane convex bodies[END_REF]. But in the present paper, we shall only consider hedgehogs with a C 2support function and we will refer to them as 'C 2 -hedgehogs'.

Gaussian curvature and algebraic volume of C 2 -hedgehogs

Let H n+1 denote the R-linear space of C 2 -hedgehogs defined up to a translation in the Euclidean linear space R n+1 and identified with their support functions. Analytically speaking, saying that a hedgehog H h ⊂ R n+1 is defined up to a translation simply means that the first spherical harmonics of its support function is not specified.

As we saw before, elements of H n+1 may be singular hypersurfaces. Since the parametrization x h can be regarded as the inverse of the Gauss map, the Gaussian curvature

K h of H h at x h (u) is given by K h (u) = 1/ det [T u x h ],
where T u x h is the tangent map of x h at u. Therefore, singularities are the very points at which the Gaussian curvature is infinite. For every u ∈ S n , the tangent map of x h at the point u is

T u x h = h(u) Id T u S n + H h (u)
, where H h (u) is the symmetric endomorphism associated with the Hessian of h at u. Consequently, if λ is an eigenvalue of the Hessian of h at u then λ + h (u) is (up to the sign) one of the principal radii of curvature of H h at x h (u) and the so-called 'curvature function' R h := 1/K h can be given by

R h (u) = det [H ij (u) + h (u) δ ij ] , (1)
where δ ij are the Kronecker symbols and (H ij (u)) the Hessian of h at u with respect to an orthonormal frame on S n . The case n = 2. From (1), the curvature function

R h := 1/K h of H h ⊂ R 3 is given by R h = (λ 1 + h) (λ 2 + h) = h 2 + h∆ 2 h + ∆ 22 h,
where ∆ 2 denotes the spherical Laplacian and ∆ 22 the Monge-Ampère operator (respectively the sum and the product of the eigenvalues λ 1 , λ 2 of the Hessian of h). So, the equation we shall be dealing with will be the following

h 2 + h∆ 2 h + ∆ 22 h = 1/K.
Note that the so-called 'mixed curvature function' of hedgehogs of R 3 , that is,

R : H 2 3 → C (S 2 ; R) (f, g) → R (f,g) := 1 2 (R f +g -R f -R g ) is bilinear and symmetric: (i) ∀ (f, g, h) ∈ H 3 3 , ∀λ ∈ R, R (f +λg,h) = R (f,h) + λR (g,h) ; (ii) ∀ (f, g) ∈ H 2 3 , R (g,f ) = R (f,g) . For any h ∈ H 3 , we have in particular R -h = R h . Note that R (1,f ) = 1 2 (∆ 2 h + 2h
) is (up to the sign) half the sum of the principal radii of curvature of H h ⊂ R 3 . The (algebraic) area and volume of a hedgehog H h of R 3 is defined by

s (h) = S 2 R h dσ and v (h) = S 2 hR h dσ.
where σ is the spherical Lebesgue measure on S 2 and R h = det (Hess (h) + h.Id 2 ) the curvature function (that is, the inverse of the Gaussian curvature K h of H h ) [START_REF] Martinez-Maure | De nouvelles inégalités géométriques pour les hérissons[END_REF]. The (algebraic) area s (h) of H h can be interpreted as the difference s + (h)s -(h), where s + (h) (resp. s -(h)) denotes the total area of the smooth regions of H h on which the Gaussian curvature is positive (resp. negative). The (algebraic) volume v (h) of H h can be regarded as the integral over R 3 -H h of the index i h (x) defined as algebraic intersection number of an oriented halfline with origin x with the surface H h equipped with its transverse orientation (number independent of the oriented half-line for an open dense set of directions).

Gauss infinitesimal rigidity in the context of hedgehogs

In this work, we shall use the Banach spaces C m that were introduced by L. Nirenberg in his study of the Minkowski problem in R 3 , (m ∈ N).

Definition 1 Let H h be a C 2 -hedgehog of R 3 . A smooth deformation of H h is the data of a differentiable map h : [0, 1] → C 2 , t → h t := h (t, .) such that h 0 = h. Definition 2 Let H f be a C 2 -hedgehog of R 3 . An infinitesimal isogauss defor- mation of H f is the data of a family (H f +tg ) t∈R of hedgehogs of R 3 , x f +tg : S 2 → H f +tg ⊂ R 3 u → x f (u) + tx g (u)
where H g is a hedgehog of R 3 such that the mixed curvature function R (f,g) :=

1 2 (R f +g -R f -R g ) be identically zero on S 2 . Definition 3 Let H f be a C 2 -hedgehog of R 3 . If every infinitesimal isogauss deformation (H f +tg ) t∈R of H f is trivial
, that is such that H g is reduced to a single point, then the hedgehog H f will be said to be Gauss infinitesimally rigid.

Remark 1. The hedgehog H g is reduced to a single point if, and only if, its support function g is the restriction to S 2 of a linear form on R 3 , which amounts to saying that its curvature function R g is identically zero on S 2 [7, Theorem 1]. Therefore, the hedgehog H f is Gauss infinitesimally rigid if, and only if,we have:

∀g ∈ C 2 S 2 ; R , R (f,g) = 0 =⇒ (R g = 0) . Remark 2. If the hedgehog H f ⊂ R 3 is trivial (that is, reduced to a point), then H f is not Gauss infinitesimally rigid. Indeed, for every regular C 2 -hedgehog H g ⊂ R 3 , we have R (f,g) = 0 although R g is not identically zero on S 2 .

Gauss infinitesimal rigidity of regular

C 2 -hedgehogs of R 3
Let us recall the proof of the Gauss infinitesimal rigidity (with respect to the curvature function) of regular C 2 -hedgehogs of R 3 (that are closed convex surfaces of class C 2 + in R 3 ). It is essentially a rewriting of the proof by J. Stoker [START_REF] Stoker | Differential geometry[END_REF]: Let H f be a regular C 2 -hedgehog of R 3 . Clearly, the regularity of H f is equivalent to the strict positivity of its curvature function R f := 1/K f . If (H f +tg ) t∈R defines an isogauss deformation of H f , then we have [START_REF] Martinez-Maure | Uniqueness results for the Minkowski problem extended to hedgehogs[END_REF]Lemma 5]:

R 2 (f,g) ≥ R f .R g
and hence R g ≤ 0 on S 2 . By taking the origin to be an interior point of the convex body bounded by H f in R 3 , we may assume without loss of generality that f > 0 so that f Rg ≤ 0 on S 2 . Now, by symmetry of the mixed volume of hedgehogs of R 3 [START_REF] Martinez-Maure | Hedgehogs and zonoids[END_REF], we get:

0 = S gR (f,g) dσ = S f R (g,g) dσ = S fR g dσ,
where σ is the spherical Lebesgue measure on S 2 . Therefore, R g is identically zero on S 2 which implies that H g is reduced to a single point by Remark 1.

Relation to Minkowski problem

There is a close connection between Gauss infinitesimal rigidity and the uniqueness question in the Minkowski problem extended to hedgehogs. This is due to the following equivalence:

∀ (f, g) ∈ C 2 S 2 ; R 2 , (R f = R g ) ⇐⇒ R (f +g,f -g) = 0 .
In [START_REF] Martinez-Maure | New notion of index for hedgehogs of R 3 and applications[END_REF][START_REF] Martinez-Maure | Uniqueness results for the Minkowski problem extended to hedgehogs[END_REF], the author gave examples of pairs of non-congruent hedgehogs of R 3 having the same curvature function. From each of these examples, we can deduce examples of nontrivial hedgehogs that are not Gauss infinitesimally rigid. It is for instance the case of the pair of hedgehogs of R 3 given by:

f (u) :=    0 if z ≤ 0 exp(-1/z 2 ) if z > 0 and g (u) :=    exp(-1/z 2 ) if z < 0 0 if z ≥ 0,
where u = (x, y, z) ∈ S 2 ⊂ R 3 . Indeed, we have clearly R (f,g) = 0. Therefore, these two nontrivial hedgehogs H f and H g are not Gauss infinitesimally rigid.

Only nonanalytic examples are known. The question of knowing whether there exists a pair of noncongruent analytic hedgehogs of R 3 with the same curvature function remains open (by 'analytic hedgehogs', we mean 'hedgehogs with an analytic support function'). As a consequence, the question of knowing whether there exist examples of nontrivial analytic hedgehogs that are not Gauss infinitesimally rigid is also open.

A bellows-type theorem for hedgehogs Lemma 4 The curvature function R :

C 2 → C 0 , h → R h is differentiable on C 2 , and: ∀ (f, g) ∈ C 2 × C 2 , dR f (g) = lim t→0 t =0 R f +tg -R f t = 2R (f,g) .
Proof of the lemma. Indeed, we have:

∀t ∈ R * + , R f +tg -R f = R f + 2tR (f,g) + t 2 R g -R f = t 2R (f,g) + tR g ,
and hence

lim t→0 t =0 R f +tg -R f t = lim t→0 t =0 2R (f,g) + tR g = 2R (f,g) . Theorem 5 Let H h be a C 2 -hedgehog of R 3 . If a smooth deformation of H h , say h : [0, 1] → C 2 , t → h t := h (t, .) ,
preserves the curvature function (that is, is such that R ht = R h for all t ∈ [0, 1]), then it also preserves the algebraic volume:

∀t ∈ [0, 1] , v (h t ) = v (h) .
Proof of Theorem 5. By assumption, the map R

• h : [0, 1] → C 0 is constant.
Since h is differentiable by assumption and R by Lemma 4, R • h is differentiable and the chain rule gives:

∀t ∈ [0, 1] , (R • h) ′ (t) = 2R ( h(t), ( ∂h ∂t ) (t)
) . Therefore, a differentiation yields:

∀t ∈ [0, 1] ,
R (h(t),( ∂h ∂t )(t)) = 0.

Now, for every t 0 ∈ [0, 1], we have : Besides, by symmetry of the mixed volume of hedgehogs, we have:

∀t ∈ [0, 1] -{t 0 } , v (
1 3 S 2 ∂h ∂t (t 0 ) R h(t 0 ) dσ = 1 3 S 2
h (t 0 ) R (h(t0),( ∂h ∂t )(t0)) dσ.

From (2), we then deduce that:

∀t 0 ∈ [0, 1] , (v • h) ′ (t 0 ) = lim t→t 0 t =t 0 v (h (t)) -v (h (t 0 )) t -t 0 = 0,
and thus all the hedgehogs of the family (H h t ) have the same (algebraic) volume.

Figure 2 .

 2 Figure 2. A hedgehog with a C 2 -support function Hedgehogs with a C 2 -support function can be regarded as the Minkowski differences of convex hypersurfaces of class C 2 + . Indeed, given any h ∈ C 2 (S n ; R), for all large enough real constant r, the functions h+r and r are support functions of convex hypersurfaces of class C 2 + such that h = (h + r)r.

Figure 3 .

 3 Figure 3. Hedgehogs as differences of convex bodies of class C 2 +

Remark 3 .

 3 Noncongruent hedgehogs that share the same curvature function may of course have different (algebraic) volume. It is for instance the case of the hedgehogs shown on Figure 4 whose support functions f , g are defined on S 2 by f (u) := ) f (u) if z = 0 0 if z = 0, where u = (x, y, z) ∈ S 2 ⊂ R 3 .

Figure 4 .

 4 Figure 4. Same curvature function and different algebraic volume


Xiaohong Jia
email: xhjia@cs.hku.hk

Yi-King Choi
email: ykchoi@cs.hku.hk

Bernard Mourrain
email: mourrain@sophia.inria.fr

Wenping Wang
email: wenping@cs.hku.hk

An Algebraic Approach to Continuous Collision Detection for Ellipsoids

Keywords: moving ellipsoids, characteristic equation, continuous collision detection, algebraic conditions, subresultants

We present algebraic conditions for characterizing three configurations formed by two ellipsoids in R 3 that are relevant to collision detection: separation, external touching and overlapping. These conditions are given in terms of explicit formulae expressed by the subresultant sequence of the characteristic polynomial of the two ellipsoids and its derivative. For any two ellipsoids, the signs of these formulae can easily be evaluated to classify their configuration. Furthermore, based on these algebraic conditions, an efficient method is developed for continuous collision detection for two moving ellipsoids under arbitrary motion.

Introduction

Collision detection finds many applications in computer graphics, computer animation, CAD/CAM as well as computational physics (see [START_REF] Eberly | 3D Game Engine Design[END_REF]). Since collision detection for general freeform moving objects is computationally very expensive, bounding volumes are often used to approximate the freeform objects in order to reduce computational cost. Ellipsoids are good candidates for such bounding volumes since they have low algebraic degree and are tight bounding volumes for a wide class of objects [START_REF] Bouville | Bounding ellipsoids for ray-fractal intersection[END_REF]; [START_REF] Lu | Variational 3d shape segmentation for bounding volume computation[END_REF]). Minimum bounding (or enclosing) ellipsoids and their computations have long been studied as a classical mathematical problem, see e.g., [START_REF] Welzl | Smallest enclosing disks (balls and ellipsoids)[END_REF]; [START_REF] Kumar | Minimum-volume enclosing ellipsoids and core sets[END_REF]; [START_REF] Todd | On khachiyan's algorithm for the computation of minimum-volume enclosing ellipsoids[END_REF]; [START_REF] Schröker | Uniqueness results for minimal enclosing ellipsoids[END_REF], and have found important applications to not only CAGD and computer graphics but also other areas such as data uncertainty analysis.

Collision detection for ellipsoids has been an active research topic over the past years, see [START_REF] Ju | A novel collision detection method based on enclosed ellipsoid[END_REF]; [START_REF] Rimon | Obstacle collision detection using best ellipsoid fit[END_REF]; [START_REF] Shiang | Estimate of minimum distance between convex polyhedra based on enclosed ellipsoids[END_REF]; [START_REF] Wang | Efficient collision detection for moving ellipsoids using separating planes[END_REF]; [START_REF] Choi | Continuous collision detection for two moving elliptic disks[END_REF][START_REF] Choi | Continuous collision detection for ellipsoids[END_REF]. [START_REF] Wang | An algebraic condition for the separation of two ellipsoids[END_REF] provide algebraic conditions for the three important configurations, i.e., separation, external touching and overlapping, of two static ellipsoids in R 3 , the three dimensional affine space. It is shown that the relative position of two ellipsoids is related to the root pattern of their characteristic equation. Specifically, two ellipsoids are separate if and only if their characteristic equation has two distinct positive roots; and they touch each other externally if and only if there is one positive double root. This result allows us to compute the configuration of two ellipsoids by simply determining the number of the positive roots of a polynomial, which can be derived by the classical Decart's rule and the modified sign variation number of the signed subresultant sequence (e.g., see [START_REF] Basu | Algorithms in Real Algebraic Geometry[END_REF]).The algebraic conditions given by [START_REF] Wang | An algebraic condition for the separation of two ellipsoids[END_REF] lays the theoretical foundation for the follow-up practical applications in collision detection for ellipsoids.

Collision detection mostly deals with moving objects. When the positions of the two ellipsoids are given in a sequence of discrete time frames, which is often the situation in computer animation, temporal and geometric coherence are exploited in [START_REF] Wang | Efficient collision detection for moving ellipsoids using separating planes[END_REF] for speeding up the computations. A separating plane for two non-intersecting ellipsoids at a frame is calculated which can then be used to quickly identify if the ellipsoids are still separate or not in the next frame. The more computationally intensive ellipsoidellipsoid intersection test needs to be carried out only when the separating plane fails to guarantee separation of the ellipsoids.

Continuous collision detection (CCD) in which object motions are given by continuous funcions of a time parameter t has been gaining increasing interests in the past decade (see e.g., [START_REF] Redon | Fast continuous collision detection between rigid bodies[END_REF]; [START_REF] Teschner | Collision detection for deformable objects[END_REF]). It focuses on determining the collision status of objects over a specific time span and is exact in the sense that no discretization of the time domain is needed. For CCD of two moving ellipoids with continuous motions, the method in [START_REF] Wang | Efficient collision detection for moving ellipsoids using separating planes[END_REF] mentioned above cannot easily be extended to deal with the problem since the computation of a separating plane requires solving for the roots of the characteristic equation. Based on the algebraic condition by [START_REF] Wang | An algebraic condition for the separation of two ellipsoids[END_REF], [START_REF] Choi | Exact collision detection of two moving ellipsoids under rational motions[END_REF] reduces the problem to an analysis of the zero set of the characteristic equation which is a bivariate polynomial equation; however, their algorithm involves the brute-force computation of the zero set and is therefore slow. Later [START_REF] Choi | Continuous collision detection for ellipsoids[END_REF] use the same CCD algebraic formulation, and further develop an efficient search scheme to determine the collision time instants from the characteristic equation in real time. The basis idea is to find the contact time instants by alternative searches in the two parameter domains of the characteristic equation. Unfortunately, the use of Bézier clipping in the search means that this method is only applicable to ellipsoids whose motions are expressed as rational functions.

In this paper, we aim at establishing algebraic conditions, based on those by [START_REF] Wang | An algebraic condition for the separation of two ellipsoids[END_REF], which are suitable for direct application to CCD of moving ellipsoids under arbitrary motions (such as the commonly used helical motions as in Example 5.2). Our algebraic conditions are based on five explicit formulae derived from the subresultant sequences of the characteristic polynomial of two ellipsoids and its derivative.

Essentially, these conditions distinguish the different root patterns of the characteristic equation, therefore distinguishing the configurations formed by two ellipsoids, without having to solve for the roots. We also develop an algorithm for CCD of two moving ellipsoids under arbitrary continuous motions which are not necessarily rational. In this case, the formulae will directly lead to five functions in the time parameter t, and CCD can simply be done by solving these five functions.

We note that [START_REF] Gonzalez-Vega | Solving the separation problem for two ellipsoids involving only the evaluation of six polynomials[END_REF] translate the algebraic conditions in [START_REF] Wang | An algebraic condition for the separation of two ellipsoids[END_REF] to a set of closed form formulae to characterize the separation of two ellipsoids, and their results are also based on the subresultant sequence of the characteristic polynomial and its first derivative. However, their conditions do not distinguish the two conditions of external touching and overlapping. Our present work is inspired by [START_REF] Gonzalez-Vega | Solving the separation problem for two ellipsoids involving only the evaluation of six polynomials[END_REF] and improves it by further distinguishing external touching and overlapping. The distinction of these two configurations is not only useful in collision detection but also of theoretical interest in its own right.

The remainder of the paper is organized as follows. In Section 2 we review the algebraic conditions given by [START_REF] Wang | An algebraic condition for the separation of two ellipsoids[END_REF] and present the algebraic tool of subresultant sequences which will later be used to derive our explicit formulae. We then derive in Section 3 the five explicit formulae for distinguishing the root patterns of the characteristic equation, thus characterizing the configurations of two ellipsoids. In Section 4, we examine the computational cost by optimizing the evaluation of the five explicit formulae. In Section 5, we present examples on applying our method to continuous collision detection for ellipsoids and compare the efficiency of our approach with [START_REF] Choi | Continuous collision detection for ellipsoids[END_REF]. We conclude the paper in Section 6.

Preliminaries

Given two ellipsoids A : X T AX = 0 and B : X T BX = 0 in R 3 , where X = (x, y, z, w) T are the homogeneous coordinates of points ( x w , y w , z w ) ∈ R 3 and A, B are 4 × 4 coefficient matrices with elements in R, the characteristic polynomial of the ellipsoids A and B is defined by

f (λ) = det(λA + B).
The characteristic polynomial f (λ) has degree 4 in R[λ], where R[λ] is the polynomial ring with real coefficients. Since A and B represent ellipsoids, we have det(A) < 0 and det(B) < 0; hence dividing each coefficient of f by the leading coefficient det(A), we get the monic form of f :

f (λ) = λ 4 + aλ 3 + bλ 2 + cλ + d,
where d = det(B)/ det(A) > 0. It also follows that zero cannot be a root of f (λ) = 0. We further assume that the interiors of the ellipsoids A and B are defined by X T AX < 0 and X T BX < 0. Now we define the three configurations of two ellipsoids: separate, external touching and overlapping.

An ellipsoid is regarded as a solid bounded by the boundary surface X T AX = 0. Two ellipsoids are separate if their boundary surfaces and interiors share no common points; otherwise, they are said to intersect. Furthermore, two intersecting ellipsoids are said to overlap if their interiors have a common point; otherwise they touch externally. That is, two intersecting ellipsoids may overlap or touch externally.

The following algebraic conditions are given by [START_REF] Wang | An algebraic condition for the separation of two ellipsoids[END_REF] on the configurations of two ellipsoids A and B.

Theorem 2.1. [START_REF] Wang | An algebraic condition for the separation of two ellipsoids[END_REF] 1. The characteristic equation f (λ) = 0 always has at least two negative roots; We shall derive explicit expressions whose signs indicate the root pattern as described in Theorem 2.1. Such a symbolic approach avoids computing the roots of the characteristic equation and makes it possible to develop an efficient method for continuous collision of two moving ellipsoids when there is a time parameter involved. Note that a symbolic treatment has been proposed for two ellipses in [START_REF] Choi | Continuous collision detection for two moving elliptic disks[END_REF], which provides a basis for continuous collision detection for two moving ellipses therein. However, it is a more difficult problem for two ellipsoids because, unlike the case of ellipses, the appearance of a double root of the characteristic equation does not necessarily mean any configuration change for two ellipsoid. See [START_REF] Choi | Continuous collision detection for two moving elliptic disks[END_REF] for a brief discussion about this difficulty.

We now introduce the concept of the subresultant sequence, an algebraic tool to be used in our derivation.

For more details, see for example, [START_REF] Geddes | Algorithms for Computer Algebra[END_REF], [START_REF] Basu | Algorithms in Real Algebraic Geometry[END_REF], von zur Gathen and Gerhard (1999), [START_REF] Kerber | Division-free computation of subresultants using Bezout matrices[END_REF].

Definition 2.2. [START_REF] Basu | Algorithms in Real Algebraic Geometry[END_REF] Let

f (λ) = n k=0 p k λ k , g(λ) = m k=0 q k λ k be two polynomials in R[λ] with degrees n = deg(f ) ≥ deg(g) = m.
The i-th Sylvester-Habicht matrix of f and g, denoted by SyHa i (f, g), is defined by

SyHa i (f, g) :=                    p n • • • • • • • • • • • • p 0 0 0 0 . . . . . . 0 . . . . . . p n • • • • • • • • • • • • p 0 . . . 0 q m • • • • • • • • • q 0 . . . . . . . . . . . . 0 0 . . . . . . . . . . . . q m • • • • • • • • • q 0 0 • • • 0                             m -i          n -i .
The i-th signed subresultant sr i (λ) of f and g is the determinant of the

(n + m -2i) × (n + m -2i) matrix whose first n + m -2i -1 columns are taken from the first n + m -2i -1 columns of SyHa i (f, g),
and
the last column is the polynomial sequence

λ n-i-1 f, λ n-i-2 f, • • • , f, g, λg, • • • , λ m-i-1 g.
Note that the first signed subresultant sr 0 (λ) is equal to the resultant Res(f, g) of f and g.

In the sequel, we consider the subresultant sequence of a polynomial f and its first derivative f . Denote the coefficient of the degree j term of the polynomial sr i by

sr ij , j = 0, • • • , deg(sr i ).
The following property of subresultant sequence will be important to our later analysis.

Lemma 2.3. Let f (λ) ∈ R[x]. Then deg(gcd(f (λ), f (λ))) = k if and only if sr kk = 0 and sr ii ≡ 0, i = 0, • • • , k -1. Furthermore, we have sr k (λ) = gcd(f (λ), f (λ)).

Explicit Formulae for Configurations of Two Ellipsoids

In this section we shall derive explicit formulae for characterizing the configurations of two static ellipsoids. Throughout we shall repeatedly apply the classical Descartes' rule of signs and the modified sign variations of the signed subresultants. See the details in, for example, [START_REF] Basu | Algorithms in Real Algebraic Geometry[END_REF].

Proposition 3.1 (Descartes' rule of signs). Let f (λ) = a n λ n + • • • + a 0 ∈ R[λ]. Then the number of positive roots of f (λ) = 0 is equal to Var(a n , • • • , a 0 ) -2k
for some non-negative integer k, where Var(s) is the number of sign variations in a sequence s.

Corollary 3.2. The number of negative roots of Table 1 summarizes Lemma 3.3 to Lemma 3.6.

f (λ) = 0 is equal to Var((-1) n a n , (-1) n-1 a n-1 , • • • , a 0 )) - 2k for some non-negative integer k. Let f (λ) = λ 4 + aλ 3 + bλ 2 + cλ + d
Next we are going to show how to determine whether the characteristic equation f (λ) = 0 has four real roots. Consider the characteristic polynomial f (λ) and its derivative f (λ), together with their first three subresultants sr 0 (λ), sr 1 (λ), sr 2 (λ). Denote the sequence

P = f (λ), f (λ), sr 2 (λ), sr 1 (λ), sr 0 , (1)
f (λ) = λ 4 + aλ 3 + bλ 2 + cλ + d f (λ) = 4λ 3 + 3aλ 2 + 2bλ + c sr 2 (λ) = (-8b + 3a 2 )λ 2 + (2ab -12c)λ + ac -16d sr 1 (λ) = (-6a 3 c + 2a 2 b 2 -12a 2 d + 28abc -8b 3 -36c 2 + 32bd)λ -9a 3 d + a 2 bc + 3ac 2 + 32abd -4b 2 c -48cd sr 0 = -192cd 2 a + 256d 3 + 144c 2 db + b 2 a 2 c 2 -6c 2 da 2 + 18c 3 ba + 144ba 2 d 2 -4b 3 a 2 d + 16b 4 d -4c 3 a 3 -128d 2 b 2 -4b 3 c 2 -27a 4 d 2 -80cb 2 ad + 18cba 3 d -27c 4 .
(2)

Using the same notation as in [START_REF] Basu | Algorithms in Real Algebraic Geometry[END_REF] Proposition 3.7. [START_REF] Basu | Algorithms in Real Algebraic Geometry[END_REF] Let f be a polynomial of degree n in R[λ], and let a < b be elements in R∪{-∞, +∞} that are not roots of f . Denote by P the signed subresultant sequence {f, f ,

sr n-2 (f, f ), • • • , sr 0 (f, f )}.
Then the number of roots of f counting without multiplicities in R is equal to MVar(P; a, b).

Lemma 3.8. The characteristic equation f (λ) = 0 has four real roots if and only if one of the following four cases holds.

1. sr 22 > 0, sr 11 > 0, sr 0 > 0; this occurs if and only if f (λ) = 0 has four distinct simple roots.

2. sr 22 > 0, sr 11 > 0, sr 0 = 0; this occurs if and only if f (λ) = 0 has one double root and two distinct simple roots.

3. sr 22 > 0, sr 11 = 0, sr 0 = 0; this occurs if and only if f (λ) = 0 has two double roots or a simple and a triple root.

4. sr 22 = 0, sr 11 = 0, sr 0 = 0; this occurs if and only if f (λ) = 0 has one quadruple root.

Proof. The above enumeration covers all possible cases under which the characteristic equation f (λ) = 0 has four real roots. We now establish the corresponding algebraic conditions for these cases. Hence the only possible choice is that sr 22 > 0, sr 11 > 0.

For case 3, since f (λ) = 0 has two double roots or a simple and a triple root, the degree of gcd(f (λ), f (λ))

is 2. By Lemma 2.3, we have sr 0 = 0 and sr 11 = 0. By Proposition 3.7, MVar(P; -∞, +∞) = 2. Note that MVar(P; -∞) = MVar +, -, sr 22 , 0, 0 MVar(P; +∞) = MVar +, +, sr 22 , 0, 0 .

Hence the only possible choice here is that sr 22 > 0.

For case 4, by Lemma 2.3, f (λ) = 0 has one quadruple root if and only if sr 22 = sr 11 = sr 0 = 0.

Next we shall provide algebraic conditions for characterizing the configurations of two ellipsoids. (a) since sr 0 = 0, f (λ) = 0 has no multiple root. Hence by Lemma 2.1 the two ellipsoids are separate.

(b) since sr 0 = 0 and sr 11 = 0, f (λ) = 0 has one double root and two simple roots. The double root λ 0 is the root of sr 1 = sr 11 λ + sr 10 , and hence negative because sr 11 > 0, sr 10 > 0. Therefore by Lemma 2.1 the two ellipsoids cannot be externally touching, and are therefore separate.

" =⇒ ":

Since A and B are externally touching, by Theorem 3.6 f (λ) = 0 has four real roots, two of which are a positive double root and the other two are negative. This also leads to two subcases: ii. f (λ) = 0 has one simple root and one triple root λ 0 . Then sr 20 /sr 22 = λ 2 0 > 0, which contradicts the fact that sr 22 > 0 and sr 20 < 0. Hence this subcase never happens.

Remark 3.1. According to Theorem 3.10, the subresultant coefficient sr 10 is crucial for distinguishing the two configurations of separation and external touching. Note that Gonzalez-Vega and Mainar (2008) uses principal subresultant sequences (i.e., sr 22 , sr 11 , sr 0 ) to derive the explicit formulae which do not distinguish separation from external touching of two static ellipsoids. We achieves this distinction by considering the sign of sr 10 .

Table 2 summarizes Theorem 3.10 on characterizing the configuration, i.e., separation, external touching or overlapping, of two static ellipsoids.

Computation Costs

Our method for determining the configuration of two ellipsoids involves only the evaluation of the explicit formulae sr 22 , sr 20 , sr 11 , sr 10 , sr 0 and Var( 1,a,b,c,d). Here we adopt the optimized evaluation of these five polynomials provided by [START_REF] Emiris | Real algebraic numbers and polynomial systems of small degree[END_REF]

> 0 > 0 > 0 separate 2 > 0 > 0 > 0 = 0 separate (2) > 0 > 0 < 0 = 0 externally touch (3) > 0 < 0 = 0 = 0 externally touch
The remaining cases correspond to overlapping ellipsoids. and let

∆ 2 = b2 -c W 1 = d -bc T = -9W 2 1 + 27∆ 2 ∆ 3 -3W 3 ∆ 2 ∆ 3 = c2 -b d W 2 = bē -c d A = W 3 + 3∆ 3 W 3 = ē -b d B = -dW 1 -ē∆ 2 -c∆ 3 T 2 = AW 1 -3 bB ∆ 1 = A 3 -27B 2 .
The two explicit formulae are then given by sr 22 := ∆ 2 , sr 20 := -W 3 , sr 11 := T, sr 10 := T 2 , sr 0 := ∆ 1 up to a positive constant multiple. The above expressions take 28 multiplications and 12 additions.

Application: Continuous Collision Detection for Two Moving Ellipsoids

Algorithm

Let M (t) be a 4 × 4 matrix, whose entries are arbitrary smooth functions in t, that represents an arbitrary continuous deformation and motion in R 3 . By applying such deformation and motion to an ellipsoid A : X T AX = 0, we obtain a moving ellipsoid A(t) :

X T M -1 (t) T AM -1 (t)X = X T A(t)X = 0.
f (λ; t) = λ 4 + a(t)λ 3 + b(t)λ 2 + c(t)λ + d(t), ( 3
)
where a(t), b(t), c(t), d(t) are all functions in t. We shall study the subresultants sr 0 (t) and sr i (λ; t), i = 1, 2 of f (λ; t) and f λ (λ; t), and denote the coefficient of the degree j term in λ of the polynomial sr i (λ; t) by

sr ij (t), i = 1, 2, j = 0, . . . , deg(sr i ). Theorem 5.1. Let 0 ≤ t 1 < t 2 < • • • < t m ≤ 1, m ≥ 0, be all the distinct contact time instants in [0, 1]
at which two given moving ellipsoids are in external touch. Let t 0 = 0, t m+1 = 1 and let δ i be an arbitrary

number in interval (t i , t i+1 ), i = 0, • • • , m.
Then the configuration of the two moving ellipsoids (i.e., whether they overlap or are separate) does not change during time interval (t i , t i+1 ), and therefore can be decided by their configuration at the time instant t = δ i , i = 0, . . . , m.

Proof. Since the center positions and the shapes of the two moving ellipsoids vary continuously, if the configuration of two ellipsoids changes from separation to overlapping, or vice versa, there must be a contact time instant t * at which the ellipsoids are in external touch. Therefore, within a time interval (t i , t i+1 ) that does not contain any other contact time instants, the configuration of the ellipsoids remains the same.

Hence we need only check the status of the two moving ellipsoids at a time instant δ i ∈ (t i , t i+1 ) to decide their configuration during the entire time interval (t i , t i+1 ).

We can see from the above theorem that the primary task for continuous collision detection for two moving ellipsoids is to determine the so-called contact time instants at which the ellipsoids are touching externally. The configuration of the ellipsoids in between two contact time instants can then be easily identified using the algebraic conditions as established in Theorem 3.10.

Next, we explain how these contact time instants can be found.

Theorem 5.2. Suppose that sr 0 (t) ≡ 0. If two ellipsoids A(t) and B(t) touch externally at t 0 , we have

sr 0 (t 0 ) = 0.
Proof. This follows immediately from Theorem 3.10.

Note 1. The special situation sr 0 (t) ≡ 0 happens when the characteristic equation f (λ; t) = 0 always has a double root λ(t) for any time instant t ∈ [0, 1], i.e., f (λ; t) = (λ -λ * (t)) 2 f (λ; t). Geometrically this occurs when two moving ellipsoids, at every time instant t, touch internally or externally, or have a reducible intersection in C 3 (see [START_REF] Tu | Using signature sequence to classify intersection curves of two quadrics[END_REF] and Example 5.1).

Lemma 5.3. Given two moving ellipsoids A(t) and B(t), t ∈ [0, 1], we have either one of the following cases: [START_REF] Bhatia | Matrix Analysis[END_REF]), the roots of f (λ; t) = 0 can be expressed as complex-valued polynomials α i (t), where i = 1, 2, 3, 4. Without loss of generality, let α 1 (t) ≡ α 2 (t) > 0 be the double root for t ∈ [t 0 , t 1 ]. Since α 1 (t) and α 2 (t) are polynomials, we must have α 1 (t) ≡ α 2 (t) > 0 for t ∈ [0, 1], which means that the moving ellipsoids are always in external touch. Hence, the moving ellipsoids can either be in external touch for all t ∈ [0, 1] or be in external touch only at some (if there is any) discrete contact time instants.

1. A(t
Theorem 5.4. Suppose that sr 0 (t) ≡ 0. If the moving ellipsoids A(t) and B(t) are in external touch only at some discrete contact time instants t i ∈ [0, 1], where i = 1, . . . , n and n ≥ 0, then at each t i , we have

sr 11 (t i ) = 0.
Proof. Since sr 0 (t) ≡ 0, we have f (λ; t) = (λ -λ * (t)) 2 f (λ; t) (Note 1). Since A(t) and B(t) are not in external touch for all t, λ * (t) cannot be a positive double root. Also, since λ = 0 cannot be a root of the characteristic equation f (λ; t) = 0 for any t, by continuity of the roots of a polynomial [START_REF] Bhatia | Matrix Analysis[END_REF],

λ * (t) < 0 for all t. Now, consider at a contact time instant t i , f (λ; t) = 0 has an additional positive double root and therefore f (λ; t i ) = 0 has two double roots. By Lemma 3.8, we therefore have sr 11 (t i ) = 0.

Note 2. Consider two moving spheres that are in external touch at only some discrete contact time instants.

Since they always have a reducible intersection in C 3 , no matter whether they are separate or not, their characteristic equation always contains a double root (which is negative). Hence, sr 0 (t) ≡ 0. Furthermore, sr 11 (t i ) = 0 if the two spheres are in external touch at t i .

Theorem 5.5. Suppose that sr 0 (t) ≡ 0 and sr 11 (t) ≡ 0. Note 3. Theorem 5.6 implies that when sr 0 (t) ≡ 0 and sr 11 (t) ≡ 0, the configuration of the ellipsoids throughout the time span can be determined by their configuration at any particular time instant, e.g., at t = 0.

We now summarize the above analysis in Algorithm 1 for continuous collision detection for two moving ellipsoids. Using Theorems 5.2, 5.4, 5.5 and 5.6, we obtain a set Z of time instants which captures all the contact time instants of two moving ellipsoids by solving for the roots of some functions under different conditions. The set Z may also contain other time instants not corresponding to any contact, which can be eliminated easily by checking against the algebraic conditions given by Theorem 3.10. Again by Theorem 3.10, the configuration of the ellipsoids at each interval defined by two consecutive contact time instants can then be determined. confirmed to be the contact time instants of the ellipsoids. We check the collision states at

δ 0 = 0.02 ∈ (0, t 1 ), δ 1 = 0.4 ∈ (t 1 , t 2 ), δ 2 = 0.95 ∈ (t 2 , 1),
and conclude that the two ellipsoids are separate in [0, t 1 ), overlap in (t 1 , t 2 ), and are separate again in (t 2 , 1] (Figure 2).

Comparison

In this section, we compare our method with [START_REF] Choi | Continuous collision detection for ellipsoids[END_REF] 3. Each test is run for 1,000 times and the average running time is taken.

The degree 2 rational rigid motion includes the simple yet commonly used motion in which an object assumes a degree 2 rotation plus a linear translation. Under this kind of rigid motion, the two methods have comparable performances and can both complete CCD in about 0.1 ms. The next pair of ellipsoids we study are under degree 2 rational affine motion, that is, one with deformation. The degree of f (λ; t) in t is 36 in this particular example. The method in [START_REF] Choi | Continuous collision detection for ellipsoids[END_REF] needs much longer (about 2.5 ms ) to compute CCD, because their method deals with a bivariate function f (λ; t) and basically needs to find a pathway λ(t) such that f (λ; t) > 0 for all t to declare that the ellipsoids are always separate. The time taken therefore depends not only on the degree of the motion but also on the topology of the zero set of f (λ; t). In this example, the two moving ellipsoids are in close proximity from time to time but remain separate within the entire time span. The approach in [START_REF] Choi | Continuous collision detection for ellipsoids[END_REF] therefore takes longer to find the pathway λ(t).

On the other hand, our method does not depend on the complexity of f (λ; t) and can solve the CCD in 0.7 ms.

In the last example, the moving ellipsoids are under rational motions of degree 4 with large deformations (Example 2 of [START_REF] Choi | Continuous collision detection for ellipsoids[END_REF]). The degree of f (λ; t) in t is 48. The [START_REF] Choi | Continuous collision detection for ellipsoids[END_REF] method takes about 1 ms while ours takes about 13 ms to complete CCD. The slower performance of our method is due to the high degree in the subresultant expressions. The degree of sr 0 (t) is 288; its composition and root finding are therefore time consuming. When only the first contact time instant of two moving ellipsoids is required (which is a common output for CCD), our method does not need to solve for all roots of the subresultants and can complete CCD in 0.4 ms.

We remark here that the above examples serve to demonstrate the efficiency of our method when time performance is of major concern. Both [START_REF] Choi | Continuous collision detection for ellipsoids[END_REF] and ours are exact continuous collision detection methods in the sense that no discretization of the time domain is needed. However, as we mentioned in Section 1, [START_REF] Choi | Continuous collision detection for ellipsoids[END_REF] solve a bivariate characteristic equation using numerical computations. We therefore use a float-point implementation of our method for comparison with [START_REF] Choi | Continuous collision detection for ellipsoids[END_REF]. Note that, when high numerical accuracy is desired, our method has the advantage that exact arithmetic can be used to achieve any required accuracy. We also listed in Table 4 the corresponding time costs for the same collision detection examples under degree 2 rigid motion and degree 2 affine motion using symbolic computation.

Conclusions

We use five explicit formulae to decide the geometric configurations, that are, separate, externally touching or overlapping of two ellipsoids. Our derivation is based on the algebraic conditions provided in [START_REF] Wang | An algebraic condition for the separation of two ellipsoids[END_REF], which shows the correspondence between the root patterns of the characteristic polynomial and the configurations of two ellipsoids. The explicit formulae are composed of the coefficients of the subresultant sequence of the characteristic polynomial and its first derivative. These algebraic formulae can be applied naturally to continuous collision detection for two moving ellipsoids whose relative positions as well as shapes vary along time governed by arbitrary continuous functions. In future work, we expect to apply simple techniques and develop explicit formulae to determine the relative geometric configuration of two general quadrics.

2.

  The two ellipsoids A and B are separate if and only if f (λ) = 0 has two distinct positive roots; 3. The two ellipsoids A and B touch each other externally if and only if f (λ) = 0 has a positive double root; 4. The two ellipsoids A and B overlap if and only if f (λ) = 0 has no positive root.

  be the characteristic polynomial of two ellipsoids A and B. Then we have the following lemma. Lemma 3.3. Var(1, a, b, c, d) = 4. Proof. Suppose on the contrary that Var(1, a, b, c, d) = 4. Then we shall have Var(1, -a, b, -c, d) = 0, which by Corollary 3.2 means that the characteristic equation f (λ) = 0 has no negative root. This contradicts Theorem 2.1. Lemma 3.4. If Var(1, a, b, c, d) = 0, the two ellipsoids A and B overlap. Proof. By Descartes' rule of signs, Var(1, a, b, c, d) = 0 implies that the characteristic equation f (λ) = 0 has no positive root. By Theorem 2.1, the two ellipsoids overlap. Hence given the characteristic polynomial f (λ) of two static ellipsoids, we can immediately tell that the two ellipsoids overlap if Var(1, a, b, c, d) = 0. However, the converse of Lemma 3.4 is not true. Lemma 3.5. If two ellipsoids A and B overlap, then Var(1, a, b, c, d) = 0 or 2; further if their characteristic equation f (λ) = 0 has four real roots, then Var(1, a, b, c, d) = 0. Proof. Since the ellipsoids A and B overlap, f (λ) = 0 has no positive root. By Descartes' rule of signs and Lemma 3.3, Var(1, a, b, c, d) = 0 or 2. Now if f (λ) = 0 has four real roots, since none of these roots can be positive, we have Var(1, -a, b, -c, d) = 4 and therefore Var(1, a, b, c, d) = 0. Using the above lemmas, we have Lemma 3.6. Two ellipsoids A and B are separate or externally touch if and only if their characteristic equation f (λ) = 0 has four real roots and Var(1, a, b, c, d) = 2. Proof. " =⇒ ": By Theorem 2.1, the configuration of separation or external touch suggests the characteristic equation has four real roots. On the other hand, By Proposition 3.1, Var(1, a, b, c, d) is an even number, hence can be 0, 2 or 4. But by Lemma 3.3 Var(1, a, b, c, d) = 4, and by Lemma 3.4 Var(1, a, b, c, d) = 0. So Var(1, a, b, c, d) = 2, and therefore f (λ) = 0 must have four real roots. " ⇐= ": Since Var(1, a, b, c, d) = 2, by Proposition 3.1 f (λ) = 0 has zero or two positive roots. If f (λ) = 0 has no positive root, the two ellipsoids overlap, which by Lemma 3.5 yields Var(1, a, b, c, d) = 0 and hence a contradiction. Therefore f (λ) = 0 has two positive roots, that is, the two ellipsoids are either separate or touch externally.

  , let MVar(P; a) denote the modified 1 number of sign variations in a sequence of polynomials P = P 0 , . . . , P n evaluated at a ∈ R ∪ {-∞, +∞}, that is, MVar(P; a) = MVar(P 0 (a), . . . , P n (a)). Also, denote MVar(P; a, b) = MVar(P; a) -MVar(P; b), where a, b ∈ R ∪ {-∞, +∞}.

For case 1

 1 , by Proposition 3.7, MVar(P; -∞, +∞) = 4. Note that MVar(P; -∞) = MVar +, -, sr 22 , -sr 11 , sr 0 MVar(P; +∞) = MVar +, +, sr 22 , sr 11 , sr 0 .Hence the only possible choice is that sr 22 > 0, sr 11 > 0 and sr 0 > 0.For case 2, since f (λ) = 0 has a double root, by Lemma 2.3 we have sr 0 = 0. Again by Proposition 3.7, MVar(P; -∞, +∞) = 3. Note that MVar(P; -∞) = MVar +, -, sr 22 , -sr 11 , 0 MVar(P; +∞) = MVar +, +, sr 22 , sr 11 , 0 .

Lemma 3. 9 .

 9 Let f (λ) = 0 be the characteristic equation of two ellipsoids A : X T AX = 0 and B : X T BX = 0. If sr 22 ≤ 0 then the two ellipsoids overlap.Proof. If sr 22 ≤ 0, then MVar(P; -∞, +∞) ≤ 2. By Proposition 3.7, f (λ) = 0 has at most 2 distinct real roots. As f (λ) = 0 has at least two negative roots (counted with multiplicity), considering f (-∞) > 0, f (+∞) > 0 and f (0) = d > 0, f has either no positive roots or two double roots of opposite sign. For the latter case, we must have both MVar(P; -∞, +∞) = 2 and sr 11 = sr 0 = 0. But sr 22 ≤ 0 yields MVar(P; -∞, +∞) ≤ 1, which is a contradiction. Thus f (λ) = 0 has no positive root and by Lemma 2.1 the two ellipsoids overlap.Theorem 3.10. Let f (λ) = 0 be the characteristic equation of two ellipsoids A : X T AX = 0 and B :X T BX = 0.1. The two ellipsoids A and B are separate if and only if Var(1, a, b, c, d) = 2 and (a) sr 22 > 0, sr 11 > 0, sr 0 > 0; or (b) sr 22 > 0, sr 11 > 0, sr 10 > 0, sr 0 = 0. 2. The two ellipsoids A and B touch each other externally if and only if (a) sr 22 > 0, sr 11 > 0, sr 10 < 0, sr 0 = 0; or (b) sr 22 > 0, sr 20 < 0, sr 11 = 0, sr 0 = 0. In the other cases, the two ellipsoids overlap. Proof. 1. " =⇒ ": Since A and B are separate, by Theorem 3.6, Var(1, a, b, c, d) = 2 and f (λ) = 0 has four real roots, two of which are distinct positive reals and the other two are negative. This leads to two subcases: (a) f (λ) = 0 has two distinct negative roots. By Lemma 3.8, we have sr 22 > 0, sr 11 > 0, sr 0 > 0; (b) f (λ) = 0 has one negative double root. By Lemma 3.8, sr 22 > 0, sr 11 > 0, sr 0 = 0. By Lemma 2.3, gcd(f (λ), f (λ)) = sr 1 = sr 11 λ + sr 10 = sr 11 (λ -λ 0 ), where λ 0 = -sr 10 /sr 11 is the negative double root of f (λ) = 0. Hence sr 11 > 0 yields sr 10 > 0. " ⇐= ": By Theorem 3.8, both a) and b) indicate that f (λ) = 0 has four real roots. Since Var(1, a, b, c, d) = 2, by Theorem 3.6 the two ellipsoids are separate or externally touch.

  (a) f (λ) = 0 has two distinct negative roots, which by Lemma 3.8 yields sr 22 > 0, sr 11 > 0, sr 0 = 0. By Lemma 2.3, sr 1 = sr 11 λ + sr 10 = gcd(f (λ), f (λ)) and therefore λ 0 = -sr 10 /sr 11 is the positive double root of f (λ) = 0. Now since sr 11 > 0, we have sr 10 < 0.(b) f (λ) = 0 has one negative double root. By Lemma 3.8 sr 22 > 0, sr 11 = 0, sr 0 = 0. By Lemma 2.3sr 2 = sr 22 λ 2 + sr 21 λ + sr 20 = gcd(f (λ), f (λ)) = sr 22 (λ -λ 0 )(λ -λ 1 ), where λ 0 is the positive double root and λ 1 is the negative double root. Hence sr 20 /sr 22 = λ 0 λ 1 < 0. Since sr 22 > 0, we have sr 20 < 0. " ⇐= ": (a) by Lemma 3.8, f (λ) = 0 has one double root and two simple roots. By Lemma 2.3, sr 1 = sr 11 λ + sr 10 = gcd(f (λ), f (λ)) and therefore λ 0 = -sr 10 /sr 11 > 0 is the double root of f (λ) = 0. By lemma 2.1 the two ellipsoids are externally touching. (b) by Lemma 3.8 there are two subcases: i. f (λ) = 0 has two double roots. By Lemma 2.3 sr 2 = sr 22 λ 2 +sr 21 λ+sr 20 = gcd(f (λ), f (λ)) = sr 22 (λ -λ 0 )(λ -λ 1 ), where λ 0 and λ 1 are the two double roots. Since λ 0 λ 1 = sr 20 /sr 22 < 0, one of the double root should be positive. By Lemma 2.1 the two ellipsoids are in external touch.

Figure 1 :

 1 Figure 1: Configurations of two moving ellipsoids under rational motion in Example 5.1. The ellipsoids touch externally at t = t 1 and t = t 2 .

  on continuous collision detection for two ellipsoids. Since the method in Choi et al. (2009) only deals with ellipsoids under rational motions, the examples we use here are also confined to rational motions. Both algorithms are implemented in C++ and the tests are run on a workstation with an Intel Xeon 3.33-GHz CPU. Double precision floating-point arithmetic is used for all computations in the comparison. Polynomials are represented in the Bernstein form in order to improve the robustness and accuracy of the computations. Root solving of polynomials are then done by subdivision using the de Casteljau algorithm. The algorithms are applied to three pairs of moving ellipsoids under different motion types to detect their collision states over a specific time span and the performance of the algorithms are listed in Table

Figure 2 :

 2 Figure 2: Configurations of two ellipsoids, one static and the other moves under non-rational helical motion, in Example 5.2. The ellipsoids touch externally at t = t 1 and t = t 2 .

Table 1 :

 1 Configurations of two ellipsoids and the corresponding values for Var(1, a, b, c, d) and the number of real roots for f (λ) = 0.

	Var(1, a, b, c, d) # of real roots for f (λ) = 0 Configuration of two ellipsoids
	0	2 or 4	overlap
	2	4 2	separate or externally touch overlap
	4		impossible
	where by Definition 2.2,

  obtained from the Bezoutian matrix. Let

				sr 2 (λ)		sr 1 (λ)
	Cases	Var(1, a, b, c, d)	sr 22		sr 20		sr 11	sr 10	sr 0	Ellipsoids configurations
	(1)	2
		b = -	a 4	, c =	b 6	, d = -	c 4	, ē = d,

Table 2 :

 2 Algebraic conditions for characterizing the configuration of two ellipsoids. The case numbers correspond to that of Lemma 3.8.

  Next we extend our method to solving continuous collision detection for two moving ellipsoids A(t) : X T A(t)X = 0

	and B(t) : X T B(t)X = 0, whose center positions vary and shapes deform with respect to a continuously
	varied parameter t ∈ [0, 1]. The two moving ellipsoids A(t) and B(t) are said to be collision-free if A(t) and
	B(t) are separate for all t ∈ [0, 1]; otherwise A(t) and B(t) collide.
	In this continuous setting, the characteristic polynomial associated with A(t) and B(t) becomes a bivari-
	ate polynomial both in parameter λ and parameter t which is given by f (λ; t) = det(λA(t) + B(t)). Clearly,
	f (λ; t) is of degree 4 in λ with coefficients as functions of t. We divide f (λ; t) by its leading coefficient in λ
	and get

  Let t 1 , t 2 , . . . , t n , n ≥ 0, be all the distinct real roots of sr 11 (t) = 0 in [0, 1]. Let δ be an arbitrary number in [0, 1] \ {t 1 , t 2 , . . . , t Proof. This follows from Lemma 5.3 and Theorem 5.4.Theorem 5.6. Suppose that sr 0 (t) ≡ 0 and sr 11 (t) ≡ 0. Then either the two ellipsoids A(t) and B(t) touch externally throughout [0, 1], or they do not have any external touch at all in [0, 1].Proof. Since sr 0 (t) ≡ 0 and sr 11 (t) ≡ 0, by Lemma 3.8, f (λ; t) = 0 always have two double roots or one quadruple root in λ. Due to the continuity of the roots of a polynomial (see e.g.,[START_REF] Bhatia | Matrix Analysis[END_REF]) and the fact that λ = 0, λ = ±∞ are not roots of f (λ; t) = 0, the signs of the roots remain the same in [0, 1]. Now, if f (λ; t i ) = 0 has a positive double root λ 0 (t i ) for some t i ∈ [0, 1], then λ 0 (t) remains positive throughout [0, 1],

	which means that A(t) and B(t) are in external touch throughout [0, 1]. On the other hand, if f (λ; t) = 0
	does not have a positive double root for any t ∈ [0, 1], then A(t) and B(t) do not touch externally at all in
	the interval.

n }. If the ellipsoids A(δ) and B(δ) touch externally, then A(t) and B(t) are in external touch throughout [0, 1].

  Algorithm 1: Collision Detection of Two Moving Ellipsoids Input : The characteristic polynomial f (λ; t) = λ 4 + a(t)λ 3 + b(t)λ 2 + c(t)λ + d(t) of two moving ellipsoids A(t) and B(t).

	Output: Three sets S, I and T containing time instants or intervals in which A(t) and B(t) are
	separate, overlap and touching externally, respectively.
	begin
	Z ← ∅, Z ← ∅, S ← ∅, I ← ∅, T ← ∅;
	if sr 0 (t) ≡ 0 then
	if sr 11 (t) ≡ 0 then	// Theorem 5.6
	if A(0) and B(0) are separate then S ← [0, 1];
	else if A(0) and B(0) overlap then I ← [0, 1];
	else if A(0) and B(0) are touching externally then T ← [0, 1];
	return;
	else
	Z ← {t | sr 11 (t) = 0};	// Theorem 5.5
	Select δ ∈ [0, 1] \ Z;
	if A(δ) and B(δ) touch externally then	// Conditions of Theorem 3.10(2)
	T ← [0, 1];
	return;
	else
	Z ← {t | sr 0 (t) = 0}
	foreach t i ∈ Z do
	if A(t i ) and B(t i ) touch externally then T ← T ∪ {t i }; // Conditions of Theorem 3.10(2)
	Z ← T ∪ {0, 1};
	Let Z = {t 1 , . . . , t n } where 0 = t 1 < • • • < t n = 1;
	foreach i = 1, . . . , n -1, select δ i ∈ (t i , t i+1 ) do	// Theorem 5.1
	if A(δ i ) and B(δ i ) are separate then	// Conditions of Theorem 3.10(1)
	S ← (t i , t i+1 );
	else
	I ← (t i , t i+1 )

Table 3 :

 3 Motion Type deg t (f (λ; t)) deg(sr 0 (t)) Run-time performance of our method against[START_REF] Choi | Continuous collision detection for ellipsoids[END_REF] to solve CCD of two moving ellipsoids under rational motions. The timing under the column "Solving f (λ; t)" is the time taken for determining the collision status over a given time span. The total time for CCD is the sum of the time taken for setting up and solving f (λ; t). All timings are averaged over 1,000 runs.

		Time (ms)
	Set up	Solve f (λ; t)	Total for CCD
	f (λ; t)	Choi 2009	Our	Choi 2009	Our

Table 4 :

 4 Run-time performance of our method to symbolically solve CCD of two moving ellipsoids under rational motions. All timings are averaged over 1,000 runs.

We first delete those polynomials identical to zero in the sequence P. Then the modified number of sign variations MVar(P, a) is similarly defined as the commonly used number of sign variations Var(P, a) except that the case {+, 0, 0, +} or {-, 0, 0, -} (exactly two zeros between the two nonzero number) counts the variation for two but not zero. See page 330 of Basu et al. (2006) for details.

Examples

Example 5.1. Let A(t) and B(t) be two moving ellipsoids defined by (x + 12t -11) 2 4 + y 2 + z 2 = 1 and (x -3) 2 4 + (y -4t + 2) 2 + (z -4t + 4) 2 = 1, respectively, where t ∈ [0, 1]. The characteristic polynomial associated with A(t) and B(t) is

up to a constant multiple. Since sr 0 (t) ≡ 0, we shall next compute

Solving the real roots for sr 11 (t) in [0, 1], we obtain t 1 = 0.5395042868, t 2 = 0.8722604191 which are both confirmed to be the contact time instants of the ellipsoids by checking against the conditions in Theorem 3.10.

Selecting

and by checking the collision states at t = δ 0 , δ 1 , δ 2 , we conclude that the two ellipsoids are separate during time interval [0, t 1 ), overlap during (t 1 , t 2 ), and are separate again during (t 2 , 1] (Figure 1).

The following example shows that our approach not only works for rational motions but also allows arbitrary functional motions, e.g., helical motions, of two moving ellipsoids.

Example 5.2. Let A and B be two static ellipsoids defined by

respectively. Let A(t) be a moving ellipsoid defined by applying to A a rotation about the axis (1, 0, 0) T by an angle 10t and then a translation along the helical curve P (t) = (cos 10t, sin 10t, 10t), where t ∈ [0, 1]. The coefficients of the characteristic polynomial associated with the moving ellipsoid A(t) and the static ellipsoid B are functions in t that contains trigonometric terms; the expressions are long and hence are omitted here.

We solve the transcedental function sr 0 (t) by the simple bracketing and bisection method. Other methods, such as the secant method, can also be used for root finding (see [START_REF] Press | Numerial Recipes 3rd Edition: The Art of Scientific Computing[END_REF] for more details).

The real roots of sr 0 (t) in [0, 1] are found to be t 1 = 0.0749830692 and t 2 = 0.8913371204, which are both