F Gauthier

P.-E Roche

Shot noise of thermal plumes : Evidence of a boundary layer instability consistent with the onset of Kraichnan's Regime of convection

Introduction and Motivation

In 1997, an enhancement of heat transfer has been reported [START_REF] Chavanne | Observation of the ultimate regime in Rayleigh-Bénard convection[END_REF] in Rayleigh Bénard convection for Rayleigh numbers Ra 10 12 . This new regime -named "Ultimate Regime" in the original paper-was interpreted as the regime predicted by Kraichnan in 1962 [START_REF] Kraichnan | Turbulent thermal convection at arbitrary Prandtl numbers[END_REF] and which is characterized by the turbulence of the boundary layers of the convection cell. Although this interpretation has been indirectly supported by specific experimental tests, including the observation of the N u ∼ Ra 1/2 scaling predicted by Kraichnan for asymptotically high Ra [START_REF] Roche | Observation of the 1/2 power law in Rayleigh-Bénard convection[END_REF], a direct evidence of fluctuations in the boundary layer was still missing. We report such an observation based on the measurement of the shot-noise induced by thermal plumes in the heating plate [START_REF] Gauthier | Evidence of a boundary layer instability at very high Rayleigh number[END_REF].

Experimental set-up

Probing fluctuations directly within the boundary layer itself raises instrumentation challenges due to its thinness, estimated as a tenth of a millimeter in cells such as the one used in [START_REF] Chavanne | Observation of the ultimate regime in Rayleigh-Bénard convection[END_REF]. To overcome this difficulty, we measured the low-frequency temperature fluctuations of the heating plate of the convection cell, which was Joule heated at a constant flux. As explained in [START_REF] Gauthier | Evidence of a boundary layer instability at very high Rayleigh number[END_REF], operation at cryogenic temperatures (around 6 K) makes this measurement possible thanks to the thermal characteristic times of the fluid (He) compared to those of the bottom plate material (OFHC annealed copper).

The cell (insert of Fig. 1) is similar to the one used in [START_REF] Chavanne | Observation of the ultimate regime in Rayleigh-Bénard convection[END_REF]. The bottom plate temperature is monitored by 2 Ge thermistors (letter "T" on Fig. 1), distant from each other by a plate's radius . At the frequency of interest (of order 10-100 mHz) they give consistent signal. The thermistors' time response was measured with the 3ω method. As displayed on Fig 1 : their caracteristic frequency response is significantly larger than the frequencies of interest. 3 Results : a new signatures of the transition to the Ultimate Regime

!" !! !" " !" ! !" # !" $ !$" !#% !#" !!% !!" !% " !#"&
The insert of Fig. 2 shows temperature spectra for 5×10 10 < Ra < 6×10 13 and for Prandtl numbers of order 1 (1.0 ≤ P r ≤ 2.8). At low frequency, spectral density plateaus E plateau can be fitted on data. The main figure shows the spectra rescaled by E plateau and by arbitrary frequencies F c . The spectral densities E plateau can be made dimensionless with the temperature difference across the cell ∆, its height h and the kinematic viscosity of the fluid ν : P ⋆ = E plateau ν/(h∆) 2 . Fig. 3 (black symbols on the main plot) shows that P ⋆ roughly scales like Ra -2/3 below Ra ≃ 2 × 10 12 but increases faster with Ra above. This threshold Ra is the same as the one for which the heat transfer across the cell improves significantly, compared to the "hard turbulence" regime present at lower Ra. This is illustrated by the insert of Fig. 3 showing the Nusselt number N u (compensated by Ra 1/3 ) versus Ra.

As illustrated by the histograms and the skewness of the temperature fluctuations (see Fig. 4), the plate's fluctuations differ significantly from a gaussian distribution above Ra = 2 × 10 12 .

Interpretation and Conclusion : a boundary layer instability

Interestingly, for Ra < 2 × 10 12 , the low-frequency spectral density P ⋆ of the plate temperature spectra coincides with the corresponding quantity in the bulk of the flow (measured previously in a similar experiment) but not above

!" !! !" " !" ! !" !# !" !$ !" !! !" " F c p e nte -4 !" !# !" !! !" " !" ! !" !!# !" !!" !" !$ !" !% !" !&
f (e n Hz ) Ra ∼ 2 × 10 12 (see Fig. 3). A straightforward interpretation, consistent with the time scales of the system [START_REF] Gauthier | Evidence of a boundary layer instability at very high Rayleigh number[END_REF], consists in assuming that the temperature of the plate follows the slow temperature fluctuations of the well-mixed core region of the flow, with an offset corresponding to the temperature drop across the thermal boundary layer. On average, this temperature drop is ∆/2. The present results indicate that the boundary layer becomes significantly "noisy" -compared to the residual fluctuations in the core of the flow-above a threshold of order Ra ∼ 10 12 . This is consistent with the occurence of a instability in the thermal boundary layer, which would cause an increase of the fluctuations of the temperature drop across it.

Eplaq u e( f ) (e n K 2 . Hz -1 ) ' ' E(f) / E f / Fc f [Hz] E(f) [K^2/Hz] slope -4 plateau
!" !" !" !! !" !# !" !$ !" !% !! !"&' !"&( !"&% !"&# " "&# "&% ζ !! !" !# $ # " %$ !" %$ !# %$ $ %$ # %$ " %$ ! Ra
In 1962, R. Kraichnan predicted that the thermal boundary layers in Rayleigh-Bénard convection should become unstable at high enough Rayleigh number, leading to an improved heat transfer. The increase in heat transfer N u(Ra) reported for Ra ∼ 10 12 -10 13 in several experiments (since 1997) has been interpreted following Kraichnan's prediction. The present result confirms the occurrence of a boundary layer instability and is consistent with this interpretation.

Fig. 1 .

 1 Fig. 1. Frequency response of the various thermistances of the set-up, measured in-situ by the 3 ω technique. Insert : the cylindrical cryogenic convection cell. The letters "T" represents the thermometers used to measure temperature fluctuations.

Fig. 2 .Fig. 3 .

 23 Fig.2. Temperature power spectra of the bottom plate (continuous line) and of a local probe located in the core of the flow (circles)[START_REF] Chavanne | Observation of the ultimate regime in Rayleigh-Bénard convection[END_REF]. In the main plot, the spectra are rescaled by 2 ad-hoc parameters E plateau and Fc.

Fig. 4 .

 4 Fig.[START_REF] Gauthier | Evidence of a boundary layer instability at very high Rayleigh number[END_REF]. Skewness of the temperature fluctuations of the plate versus Ra with same symbols as Fig.3. Insert : Histograms of the (centered) temperature fluctuations normalized by the standard deviation. Ra from top to bottom : 5.6×10 10 , 5.5×10 11 , 1.7 × 10 12 , 4.5 × 10 12 , 1.6 × 10 13 and 5.7 × 10 13 . The continuous line corresponds to a gaussian.

Acknowledgements : We thank F. Chilla, Y. Gagne and B. Castaing for stimulating discussions.


Jane Tournois

Rahul Srinivasan

Pierre Alliez

Perturbing Slivers in 3D Delaunay Meshes

Isotropic tetrahedron meshes generated by Delaunay triangulations are known to contain a majority of well-shaped tetrahedra, as well as spurious sliver tetrahedra. As the slivers hamper stability of numerical simulations we aim at removing them while keeping the triangulation Delaunay for simplicity. The solution which explicitly perturbs the slivers through random vertex relocation and Delaunay connectivity update is very effective but slow. In this paper we present a perturbation algorithm which favors deterministic over random perturbation. The added value is an improved efficiency and effectiveness. Our experimental study applies the proposed algorithm to meshes obtained by Delaunay refinement as well as to carefully optimized meshes.

Introduction

Delaunay refinement algorithms [START_REF] Chew | Guaranteed-quality Delaunay meshing in 3D (short version)[END_REF][START_REF] Shewchuk | Tetrahedral mesh generation by Delaunay refinement[END_REF][START_REF] Li | Generating well-shaped Delaunay meshes in 3D[END_REF][START_REF] Rineau | Meshing 3D domains bounded by piecewise smooth surfaces[END_REF] have been extensively studied in the literature. They are amenable to analysis, and hence are reliable algorithms. In addition, the robust implementations of Delaunay triangulations which are now available greatly facilitate the implementation of Delaunay-based mesh refinement algorithms. However, most Delaunay refinement algorithms fail at removing all badly-shaped tetrahedra, and a special class of almost-flat tetrahedra (so-called slivers) may remain in the triangulation. These slivers, with dihedral angles close to 0 and to π, are problematic for many numerical simulations.

Slivers

Many finite element methods require discretizing a domain into a set of tetrahedra. These applications require more than just a triangulation of the domain for simulation and rendering. The accuracy and the convergence of these methods depend on the size and shape of the elements apart from the fact that the mesh should conform to the domain boundary [START_REF] Shewchuk | What is a good linear element? interpolation, conditioning, and quality measures[END_REF]. Both the bad quality and the large number of the mesh elements can negatively affect the execution of a simulation. It is required that all elements of the mesh are well-shaped as the accuracy of the simulations and computations can be compromised by the presence of even a single badly shaped element. In general it is desirable to bound the smallest dihedral angle in the mesh, from below. The Delaunay refinement technique guarantees a bound on the radiusedge ratio of all mesh elements, which is the ratio of the circumradius to the shortest edge length of a tetrahedron. Although in 2D this translates into a lower bound on the minimum angle in the mesh, in 3D it does not: a bound on the radius-edge ratio is not equivalent to a bound on the smallest dihedral angle.

The only bad elements that remain after Delaunay refinement are slivers. A sliver tetrahedron is formed by almost evenly placing its 4 vertices near the equator of its circumsphere (see Figure 1), and has a bounded radius-edge ratio. In such a sliver the smallest dihedral angle can be very close to 0 • , and a numerical simulation may be far from accurate in the presence of slivers.

Tetrahedron Quality

Several tetrahedron quality criteria have been defined and used in the literature depending on the application. The radius edge ratio ρ of a simplex is defined as the ratio of its circumradius to the length of the shortest edge. This measure, which is minimal for the regular tetrahedron, unfortunately cannot detect slivers, though it is used in Delaunay refinement algorithms to define bad simplices. The radius ratio, defined as the ratio of the inradius (insphere radius) to the circumradius (circumsphere radius), is another popular measure of tetrahedron quality. It is desired to ensure that radius ratio of all tetrahedra are bounded from below by a constant.

Another criterion for mesh generation is the minimum dihedral angle θmin. It can be shown that a lower bound on the radius ratio is equivalent to a lower bound on the minimum dihedral angle. In the sequel we choose this measure to evaluate the mesh quality as it is more intuitive and geometrically meaningful than, e.g., the radius ratio, which combines the six dihedral angles of a tetrahedron.

Consider an arbitrary tetrahedron τ with triangular faces T1, T2, T3, T4. Let the areas of these triangles be denoted by S1, S2, S3, S4 respectively, the dihedral angle between Ti and Tj by θij and the length of the edge between Ti and Tj by lij. The volume V of τ is given by

V = 2 3lij SiSjsinθij for i = j in {1, 2, 3, 4}. (1)
Let rC , rI be the circumradius and inradius of τ and ri be the circumradius of Ti for i in {1, 2, 3, 4}. We know that for any tetrahedron, ri ≤ rC . This gives Si ≤ πr 2 i ≤ πr 2 c , and we also have a bound on the volume V ≥ 4 3 πr 3 I . Using Equation 1, for i = j we have

4 3 πr 3 I ≤ 2 3lij SiSjsinθij ≤ 2 3lij π 2 r 4 c sinθij,
and we get

sinθij ≥ 2 π • r 3 I r 3 c lij rC ≥ 2 π • a 3 0 ρ0 ,
where a0 is the radius-radius ratio and ρ0 is the radius-edge ratio. Finally,

θij ≥ sin -1 2 π • a 3 0 ρ0 .
Li [START_REF] Li | Sliver-free 3-Dimensional Delaunay mesh generation[END_REF] uses a different parameter of tetrahedron quality to define a sliver. Denote the volume of tetrahedron pqrs by V and its shortest edge length by l. The volume per cube of shortest edge length (σ = V l 3 ) is used as a measure of the shape quality along with the radius-edge ratio ρ, or on its own [START_REF] Cheng | Sliver exudation[END_REF]. According to Li, a tetrahedron pqrs is called sliver if ρ(pqrs) ≤ ρ0 and σ(pqrs) ≤ σ0, where ρ0 and σ0 are constant.

Previous Work

The problem of removing slivers from a 3D Delaunay mesh has received some attention over the last decade. Delaunay refinement gets so close to providing a perfect output that removing the leftover slivers is generally performed as a post-processing step that is worth it. Previous work on removing and avoiding the creation of slivers can be classified into three parts: The Delaunay-based methods, the weighted Delaunay-based methods, and the non-Delaunay methods. For each part, postprocessing steps and complete mesh generation algorithms can be studied. This paper focuses on a post-processing step, devised to take as input a Delaunay mesh and to improve its quality in terms of dihedral angles.

Delaunay-based

Vertex Perturbation

Li [START_REF] Li | Sliver-free 3-Dimensional Delaunay mesh generation[END_REF][START_REF] Edelsbrunner | Smoothing and cleaning up slivers[END_REF] proposes to explicitly perturb the vertices incident to a sliver in an almost-good mesh, by locally relocating them so as to remove the incident slivers. The idea is based on the fact that, for any triangle qrs, the region of locations of the vertex p such that the tetrahedron pqrs is a sliver, is very small. Moving the point p out of this region ensures that the tetrahedron is not a sliver anymore, or has disappeared once the Delaunay connectivity is updated. This is achieved by moving the point p to a new location inside a small ball centered at p, whose radius is proportional to the distance from p to its the nearest neighbor. The author shows that for certain values of the involved parameters, there always exists some points in this ball which are outside all regions that form slivers with nearby triangles. Li uses the union graph concept to avoid circular dependencies on vertex perturbations. The following theorem [START_REF] Li | Sliver-free 3-Dimensional Delaunay mesh generation[END_REF] proves the existence of such a point that makes the mesh locally sliver-free.

Theorem 1 (Sliver theorem). If every simplex in a Delaunay triangulation has radius-edge ratio of at least ρ0, then there is a constant σ0 > 0 and a very mild perturbation S with σ(τ ) ≥ σ0 for each tetrahedron τ in the perturbed triangulation.

Based on this theorem, Li proposes an algorithm that applies mild random perturbations to the mesh until one which removes slivers is found. One drawback of the above result is the pessimistic theoretical estimate of the bounds on the involved parameters. These bounds are either too small or too large to have any significance. In practice, though this technique is very effective, when targeting a large bound on the minimum dihedral angle (e.g.15 • ), the average number of trials of random perturbations required is very large. In our experiments, it is not rare to apply hundreds of random perturbation trials on a single vertex before succeeding in removing a sliver. This number is not surprising when seeking a high minimum dihedral angle such as 15 • as the corresponding tetrahedron is not a real sliver anymore. However the fact that the perturbation succeeds even for a high minimum dihedral angle is at the core of our motivation. Finally, the fact that this method always maintains the mesh as a true Delaunay triangulation makes it both robust and practical.

Sliver-free mesh generation.

Some mesh generation algorithms are designed to avoid creating slivers. For example, Delaunay refinement can be modified by choosing a new type of Steiner point which does not create any sliver [START_REF] Li | Spacing control and sliver-free Delaunay mesh[END_REF][START_REF] Li | Generating well-shaped Delaunay meshes in 3D[END_REF][START_REF] Liu | New strategy for unstructured mesh generation[END_REF]. As an example, Chew's algorithm [START_REF] Chew | Guaranteed-quality Delaunay meshing in 3D (short version)[END_REF] inserts Steiner points in a randomized manner, to avoid the creation of slivers. This method has a theoretical lower bound of arcsin 1/4 ≈ 14.5 • on the angles of the triangular faces of the mesh.

Weighted Delaunay-based

Sliver exudation

First described by Cheng et al. [START_REF] Cheng | Sliver exudation[END_REF], sliver exudation is a technique based on turning a Delaunay triangulation into a weighted Delaunay triangulation [START_REF] Boissonnat | Curved Voronoi diagrams[END_REF], devised to trigger flips so as to increase the minimal angle. Edelsbrunner and Guoy [START_REF] Edelsbrunner | An experimental study of sliver exudation[END_REF] provide an experimental study of sliver exudation, and show that it works pretty well in practice as a post-treatment applied to a triangulation obtained by Delaunay refinement [START_REF] Rineau | Meshing 3D domains bounded by piecewise smooth surfaces[END_REF]. The main strategy of the algorithm consists of assigning a weight to each vertex so that the weighted Delaunay triangulation is free of any slivers after connectivity updates, without any changes over the vertex locations. This method successfully increases all dihedral angles above 5 • in the best configuration (see Section 3), but as admitted in [START_REF] Edelsbrunner | An experimental study of sliver exudation[END_REF], the theoretical bound on the dihedral angle is too small to be of any practical significance.

Beside being not strictly Delaunay anymore, the main disadvantage of sliver exudation is that the process often ends with leftover slivers near the boundary [START_REF] Edelsbrunner | An experimental study of sliver exudation[END_REF]. This is mainly due to the fact that sliver exudation is not allowed to modify the topology of the boundary of the mesh. Hence, weight assignments close to the boundary are constrained and do not always manage to remove the slivers.

Complete algorithm.

Cheng and Dey [START_REF] Cheng | Quality meshing with weighted Delaunay refinement[END_REF] propose a complete Delaunay refinement algorithm, combined with the sliver exudation technique. This type of weighted-Delaunay algorithm is also used to handle input domains containing sharp creases subtending small angles [START_REF] Cheng | Weighted Delaunay refinement for polyhedra with small angles[END_REF].

Non-Delaunay

Local combinatorial operations.

Though a Delaunay-refined triangulation is known to have nice properties on its angles in 2D [START_REF] Edelsbrunner | Algorithms in Comb. Geometry[END_REF], there is no theoretical guarantee on the dihedral angles in 3D. One valid choice consists of leaving the Delaunay framework by flipping some well-chosen simplices [START_REF] Shewchuk | Two discrete optimization algorithms for the topological improvement of tetrahedral meshes[END_REF][START_REF] Klingner | Aggressive Tetrahedral Mesh Improvement[END_REF], either as a post-processing step to the meshing process [START_REF] Krysl | Variational Delaunay approach to the generation of tetrahedral finite element meshes[END_REF], or during the whole process [START_REF] Cutler | Simplification and improvement of tetrahedral models for simulation[END_REF]. As long as the triangulation remains valid, flips can be performed on its edges and facets. Joe gives a description of all possible flips [START_REF] Joe | Construction of 3-dimensional improved-quality triangulations using local transformations[END_REF] that can be made in a triangulation, and a triangulation improvement algorithm through these flips. Although each improvement in this algorithm is local, the complete algorithm succeeds in improving the overall quality of the mesh.

Dealing with non-Delaunay meshes can also be combined with optimization steps, such as Laplacian smoothing [START_REF] Forsman | Tetrahedral mesh generation in convex primitives by maximizing solid angles[END_REF], which relocates each vertex to a new location computed as an average of the incident vertex positions. Laplacian smoothing can be applied to any valid triangulation.

Complete algorithm.

Some other types of triangulations, such as for example max-min solid angle triangulations [START_REF] Joe | Delaunay versus max-min solid angle triangulations for 3-dimensional mesh generation[END_REF] can be computed to improve the solid angles as compared to that in a Delaunay triangulation. This method generates a set of well-distributed points in the input polyhedral domain and first computes a Delaunay triangulation of these vertices. Then, local combinatorial transformations are applied to satisfy the local max-min angle criterion. These local transformations can in fact be applied to any triangulation as a post-processing step.

Instead of performing local improvements through flips in a Delaunay mesh, Labelle and Shewchuck [START_REF] Labelle | Isosurface stuffing: fast tetrahedral meshes with good dihedral angles[END_REF] propose a fast lattice refinement technique which constructs a triangulation based on two nested regular or adapted grids. In its graded version this algorithm provides a theoretical bound on the dihedral angles which is much more practical than provided by other algorithms.

Contribution

We present a sliver removal algorithm inspired by Li's random perturbation algorithm [START_REF] Li | Sliver-free 3-Dimensional Delaunay mesh generation[END_REF]. Our algorithm is made more deterministic by choosing a favored perturbation direction for each vertex incident to one or more slivers, before resorting to Li's random perturbation if the favored perturbation fails at removing the incident slivers. Our experiments show that the chosen deterministic directions are sufficient to remove more than 80% of the slivers of a mesh, leading to shorter computational times. In addition, our approach reaches higher minimum dihedral angles in practice.

Algorithm

We describe a sliver perturbation algorithm which improves in a hill-climbing manner the dihedral angles of an input isotropic Delaunay mesh. This algorithm can be used as a post-processing step after refinement or optimization.

To improve the dihedral angles of the mesh tetrahedra, the rationale behind our approach is as follows: each vertex v incident to at least one sliver is repeatedly relocated through a perturbation vector pv such that when v moves to v + pv, the incident sliver gets flipped. More specifically, the chosen direction for pv is not devised to improve the shape of the slivers, but rather to worsen them instead, so that they get flipped. Two directions are favored by the algorithm: the incident squared circumradius gradient ascent (see Section 2.1) and the sliver volume gradient descent (see Section 2.2). The length of the perturbation vector is heuristically chosen as a fraction (usually between 0.05 and 0.2) of the minimum incident edge length. If neither of these two perturbation vectors succeed in flipping a sliver we resort to random perturbations (see Section 2.3). If the whole sequence does not improve the local minimum dihedral angle then we restore the vertex to its original location before perturbation.

By construction, our combined perturbation algorithm is hill-climbing in the sense that the dihedral angles in the output mesh must be higher than the ones in the input mesh. The theoretical proofs of Li's method [START_REF] Li | Sliver-free 3-Dimensional Delaunay mesh generation[END_REF] concerning random perturbation apply to this combined perturbation method as we resort to it in case of failure of the deterministic perturbation.

When more than one sliver is incident to a vertex v, all perturbation vectors must be compatible (i.e., pushing in a similar direction) to be effective. In our current algorithm a set of perturbation vectors are said to be compatible if all their pairwise dot products are positive. The perturbation vector pv is then set to be the average of these vectors. When not compatible, v is perturbed only using random perturbations. The algorithm relies on a modifiable priority queue, built in a way such that vertices incident to fewer slivers are processed first. Hence, any "chain" of slivers (set of slivers sharing at least one vertex) is treated starting from its endpoints thereby minimizing the need to process vertices incident to more than one sliver. Insert all vertices affected by relocation into P, with their new perturbation vector. end while Note that each vertex move is conditional, as we want our algorithm to be hillclimbing in terms of dihedral angles. We need to check that the minimum dihedral angle of the triangulation does not decrease, and that the topology of the boundary is not affected. Otherwise, the relocation is canceled.

Algorithm 1 Sliver perturbation

Each time a vertex is effectively relocated, the priority queue is updated. Moving v to v in a Delaunay mesh makes combinatorial changes (and, hence, changes on incident dihedral angles) on the vertices incident to v before its removal, and the ones incident to v after its insertion. We first compute the perturbations associated with all these vertices, and insert them into the priority queue.

The order in which the vertices are processed in the priority queue is related to the vertex type. Interior vertices are processed first, since they are more likely to be easily perturbable than boundary vertices. The boundary vertices are constrained to be located on the boundary, and their move must not break the topology of the mesh. These constraints make them more difficult to perturb. The other ordering criteria is discussed in Section 3.

Circumsphere Radius

Fig. 3. Circumsphere of a sliver. Before perturbation (left), the sliver is close to the equatorial plane of its circumsphere. A very mild perturbation of one of the sliver vertices (right) makes its circumradius increase considerably.

In an almost-good isotropic tetrahedron mesh the distribution of the mesh vertices is locally uniform. Hence, perturbing the vertex locations so as to make the radius of the sliver's circumsphere explode triggers many flips as the empty circumsphere property must hold after Delaunay connectivity update.

Let τ be the sliver, and {pi}i=0,1,2,3 its vertices. Without loss of generality, and since the sequel remains true by translation, we can assume that p0 = 0 R 3 . We also assume that this vertex is fixed. Let c be τ 's circumcenter. We have ||c|| = R the radius of τ 's circumsphere. Then, ∇R 2 = ∇||c|| 2 . We aim at computing ∇R 2 . Let {pi}i=0,1,2,3 be the vertices of tetrahedron τ , with p0 = 0 R 3 and pi = (xi, yi, zi) for i = 1, 2, 3. Also, let p 2 i be ( Thus we have, ∇p

x 2 i + y 2 i + z 2 i ).
1 ||c|| 2 =     ∂||c|| 2 ∂x 1 ∂||c|| 2 ∂y 1 ∂||c|| 2 ∂z 1     with ∂||c|| 2 ∂x1 = ∂ ∂x1 D 2 x + D 2 y + D 2 z 4a 2 = 1 2a 3 • a • (Dx • ∂Dx ∂x1 + Dy • ∂Dy ∂x1 + Dz • ∂Dz ∂x1 ) - ∂a ∂x1 • (D 2 x + D 2 y + D 2 z ) , ∂||c|| 2 ∂y1 = 1 2a 3 • a • (Dx • ∂Dx ∂y1 + Dy • ∂Dy ∂y1 + Dz • ∂Dz ∂y1 ) - ∂a ∂y1 • (D 2 x + D 2 y + D 2 z ) , ∂||c|| 2 ∂z1 = 1 2a 3 • a • (Dx • ∂Dx ∂z1 + Dy • ∂Dy ∂z1 + Dz • ∂Dz ∂z1 ) - ∂a ∂z1 • (D 2 x + D 2 y + D 2 z ) .
and

∇p 1 a =    ∂a ∂x 1 ∂a ∂y 1 ∂a ∂z 1    =   y2z3 -y3z2 -(x2z3 -x3z2) x2y3 -x3y2   , ∇p 1 Dx =    ∂Dx ∂x 1 ∂Dx ∂y 1 ∂Dx ∂z 1    =    -2x1 ∂a ∂x 1 -2y1 ∂a ∂x 1 + p 2 2 z3 -p 2 3 z2 -2z1 ∂a ∂x 1 -p 2 2 y3 + p 2 3 y2   , ∇p 1 Dy =    ∂Dy ∂x 1 ∂Dy ∂y 1 ∂Dy ∂z 1    =    -2x1 ∂a ∂y 1 -p 2 2 z3 + p 2 3 z2 -2y1 ∂a ∂y 1 -2z1 ∂a ∂y 1 + p 2 2 x3 -p 2 3 x2   , ∇p 1 Dz =    ∂Dz ∂x 1 ∂Dz ∂y 1 ∂Dz ∂z 1    =    -2x1 ∂a ∂z 1 + p 2 2 y3 -p 2 3 y2 -2y1 ∂a ∂z 1 -p 2 2 x3 + p 2 3 x2 -2z1 ∂a ∂z 1   .
Following a gradient ascent scheme, the vertex position pi evolves this way:

p next i = pi + ∇piR 2 τ /||∇piR 2 τ ||
, where the step length is taken as a fraction of the minimum incident edge length to pi. A relocation is performed only if the new minimal dihedral angle in the tetrahedra impacted by the relocation is not smaller than it was before relocation. As shown by Figure 3, the squared radius of τ 's circumsphere increases very fast for a small perturbation of one of its vertices' positions. The circumsphere, now huge, most probably includes other mesh vertices, which triggers a flip to maintain the empty sphere Delaunay property.

Volume

One of the main characteristics of a sliver is that its volume is strictly positive albeit small with respect to its smallest edge length, and possibly arbitrarily small. This property can be exploited in order to apply a perturbation devised to generate a sliver with negative volume and hence to trigger a combinatorial change.

Let {pi}i=1,2,3 be the three fixed points of τ , and p0 the vertex to be perturbed. The volume of τ is

Vτ = 1 6 x0 y0 z0 1 x1 y1 z1 1 x2 y2 z2 1 x3 y3 z3 1
.

Then, we get the volume gradient:

∇p0Vτ = 1 6   y2z3 + y1(z2 -z3) -y3z2 -z1(y2 -y3) -x2z3 -x1(z2 -z3) + x3z2 + z1(x2 -x3) x2y3 + x1(y2 -y3) -x3y2 -y1(x2 -x3)   .
Following a gradient descent scheme, the vertex position pi evolves this way:

p next i = pi -∇piVτ /||∇piVτ ||,
where the step length is taken as a fraction of the minimum incident edge length to pi. A relocation is performed only if the new minimal dihedral angle of the tetrahedra impacted by the relocation is not smaller than it was before relocation. A negative tetrahedron volume triggers a flip to maintain a valid Delaunay triangulation.

Random perturbation

When both ∇V and ∇R 2 fail at flipping the considered slivers by vertex perturbation, we use a random perturbation based on Li's approach [START_REF] Li | Sliver-free 3-Dimensional Delaunay mesh generation[END_REF]. A perturbation satisfying three conditions (flip sliver, improve minimum dihedral angles, preserve restricted Delaunay triangulation) is searched for randomly inside a sphere centered at v. In accordance with Li's algorithm, the magnitude of the perturbation vector is set to fraction of the minimum incident edge length.

Experiments and Results

The algorithm presented has been implemented with the 3D Delaunay triangulation of the Computational Geometry Algorithms Library [1]. Our implementation of Li's random perturbation algorithm is based upon Algorithm 1, with one single perturbation type: the random one, described in Section 2.3. For each of the following experiments we set 100 trials of random perturbations (in our combined version as well as in the purely random algorithm).

The order in which the vertices are processed in the priority queue has been chosen empirically as a result of many experiments. Interior vertices are processed first, with priority over boundary vertices. Boundary vertices are constrained so as to remain on the domain boundary and their relocation is invalid if they modify the local restricted triangulation. This makes boundary vertices more difficult to perturb than interior vertices. The second order criterion is the number of incident slivers to the processed vertex. The idea behind this choice is that a chain of slivers (several incident slivers) is more difficult to perturb than an isolated sliver as the directions of gradients may not be compatible. However, if the endpoints of the chain are successfully perturbed, we ideally would not have to deal with vertices incident to more than one sliver. Thirdly, the vertex incident to a smaller dihedral angle is processed first, as our first goal is to remove the worst tetrahedra.

In our experiments, the ∇R 2 direction turns out to be more effective than ∇V at perturbing a sliver. On average, this perturbation is responsible for about 80% of all sliver flips. The ∇V perturbation accounts for about 15% of the flips while the random perturbation counts for the remaining 5%. The priority given to ∇R 2 over ∇V and random while picking the perturbation vector can be blamed for distorting these statistics, but we have chosen this order because it turns out to be the most effective. Giving priority to ∇V results in an overall slowdown. Random perturbation always remains the last resort in the combined perturbation algorithm as the deterministic directions are favored.

The following experiments show what our combined algorithm can achieve on meshes generated by Delaunay refinement alone and on some meshes which have been optimized after refinement. A mesh optimization algorithm is in general devised to improve the mesh quality [START_REF] Amenta | Optimal point placement for mesh smoothing[END_REF] while simpler algorithms aim at evenly distributing the vertices in accordance to a given mesh sizing function. Note that a mesh with well-spaced vertices does not mean an absence of slivers inside the mesh [START_REF] Talmor | Well-spaced points for numerical methods[END_REF], and hence sliver removal is still required. The mesh optimization schemes used in our experiments are the centroidal Voronoi tessellation [START_REF] Du | Centroidal Voronoi tessellations: applications and algorithms[END_REF] using the Lloyd iteration, and the Optimal Delaunay triangulation (ODT for short) [START_REF] Chen | Optimal Delaunay triangulation[END_REF]. Both of these optimization methods have been implemented in a way that respects the local density of the mesh. It is important to not modify the density of a graded mesh, and to not decrease its quality.

Figures 4 and5 provide the computation times and the best minimum dihedral angles obtained in our experiments. The same experiment has been carried out on many other models (not shown), giving similar results. Figures 4 and5 emphasize that, for the same definition of a sliver (in terms of smallest dihedral angle), the combined algorithm is faster in removing all slivers by explicit perturbation compared to using Li's random perturbation alone. Moreover the combined algorithm reaches higher minimum dihedral angles.

The algorithm obtains fairly high minimum dihedral angles when the input is a mesh obtained by Delaunay refinement. Figures 4 and5 illustrate that when the mesh is optimized prior to perturbation, the time taken for the algorithm to succeed in removing all slivers is shorter and that it can reach a higher minimum dihedral angle. As shown by histograms of Figure 4, the algorithm takes 611 seconds to perturb the mesh obtained after Delaunay refinement so that no dihedral angle is below 17 • . If the same mesh is optimized prior to perturbations the time taken goes down to 76 seconds for Lloyd and even further down to 11 seconds for ODT. Overall the same histograms show that a mesh optimized by ODT is easier to perturb and can reach a higher minimum angle (25 • ) than a mesh optimized by Lloyd (21 • ). However, optimization can be costly. The optimizations performed on Figure 4 meshes before applying perturbation took about 200 seconds. In spite of this additional cost, the combined perturbation algorithm remains more efficient than the random one. The same comments apply to Figure 5. The gradation of the mesh in Figure 5, along with the numerous high curvature regions, make it more difficult to perturb in a way that still preserves the gradation, even after optimization. Even in this case, ODT reaches a higher minimum angle.

For comparison we have also performed sliver exudation on meshes generated by Delaunay refinement and on meshes optimized after refinement. As expected sliver exudation performs better on the optimized meshes.

We performed two other experiments that were abandoned since they rarely succeeded in improving the mesh quality. While computing the perturbation of a vertex incident to more than one sliver, we tried combining ∇V vector of one of the slivers and ∇R 2 of the other by using their average as perturbation direction if they were compatible. In practice such a combination was almost never successful removing the slivers. The other aborted experiment consisted of removing from the mesh the vertices that every explicit perturbation failed to perturb. In practice this never resulted in improving the minimum dihedral angle.

Moreover, our experiments show that successively applying our combined algorithm to the mesh several times while progressively increasing the angle bound that defines a sliver provides higher minimal dihedral angles at the price of higher computation times. This amounts to giving priority to vertices incident to the worst slivers, cluster by cluster of minimum dihedral angles.

Table 1 summarizes the best angles obtained in this way using combined perturbation, random perturbation and sliver exudation. In this labor-intensive experiment we only measure how far we can go in terms of dihedral angles and do not consider timing. Finally, Figure 6 shows some Delaunay meshes obtained by Delaunay refinement followed by ODT optimization and perturbed with the combined algorithm along with their dihedral angle histograms.

Conclusion and discussion

We have presented a practical vertex perturbation algorithm for improving the dihedral angles of a 3D isotropic Delaunay triangulation. The key idea consists of performing a gradient ascent over the sliver circumsphere radius as well as a gradient descent over the sliver volume. All vertices incident to slivers are processed, in an order devised to improve effectiveness and computation times. We compare our approach with pure random perturbation and sliver exudation. 1. Angles. Minimum dihedral angles obtained by the different perturbation algorithms (combined perturbation, random perturbation, and sliver exudation). To achieve these maxima, combined perturbation takes about twice the exudation time, and random perturbation takes about six times the exudation time. Our experiments show that we are both faster and able to reach higher minimum dihedral angles. Our scheme is particularly well suited as a post-processing step after mesh optimization [START_REF] Tournois | Interleaving delaunay refinement and optimization for practical isotropic tetrahedron mesh generation[END_REF]. We also plan to use it in the context of mesh generation from multi-material voxel images [START_REF] Boltcheva | Mesh generation from 3d multi-material images[END_REF].

In the cases where all vertices of a sliver are on the domain boundary, the perturbation can fail in removing a sliver as the boundary vertices are too constrained. One way to extend our approach would be to also perturb the vertices incident to the sliver and whose relocation can impact the sliver. Future work will focus on obtaining a proof of termination of our combined perturbation algorithm, and some tighter lower bounds on output dihedral angles.

Fig. 1 .

 1 Fig. 1. Tetrahedron shapes. A sliver (left) has its four vertices close to a circle, four very small dihedral angles (close to 70.5 • ), and two very large (close to 180 • ). A regular tetrahedron (right) is well shaped and has its dihedral angles close to 70.5 • . Each of the other tetrahedra (middle) present a different type of degeneracy.

Fig. 2 .

 2 Fig. 2. Sphere. (Left) Graded mesh with 3195 vertices, output angles are in [23.5; 142.5]. (Right) Uniform mesh with 7041 vertices, output angles are in [30.02; 138.03].

  The center c of the circumsphere of τ is given by c =

Fig. 4 .

 4 Fig. 4. Dinosaur. Comparison of the timings for our perturbation and random perturbation (in seconds) w.r.t. the sliver angle bound α on the Dinosaur model meshes obtained by Delaunay refinement (left), followed by Lloyd optimization (middle) and ODT optimization (right).

Fig. 5 .

 5 Fig. 5. Bimba. Comparison of the timings for our perturbation and random perturbation (in seconds) w.r.t. the sliver angle bound α on the Bimba model meshes obtained by Delaunay refinement (left), followed by Lloyd optimization (middle) and ODT optimization (right).

Fig. 6 .

 6 Fig. 6. Delaunay meshes perturbed with combined perturbation algorithm after ODT optimization.

  Input: T : a Delaunay triangulation, α: the angle bound defining slivers, and Nmax: the maximum number of random trials.

	end if
	else
	revert every move of v and continue.
	end if
	end while
	if the number of slivers did not decrease, then
	Compute a new perturbation vector (another type, if possible),
	and re-insert v into P.
	end if

Let P be a priority queue of Delaunay vertices. Fill P with vertices incident to slivers, Compute perturbation vector pv for each vertex v in P, while P non-empty do Pop v from P, while no combinatorial change occurs, and #loops < Nmax do

Compute v = v + pv, if conditionally relocating v to v succeeds, then v ← v , if pv is random, then compute a new pv.

Acknowledgments

The authors wish to thank Mathieu Desbrun for helpful discussions and suggestions. We also thank Mariette Yvinec and Laurent Rineau for providing us with their implementation of sliver exudation.


Pierre-Frédéric Villard
email: p.villard@imperial.ac.uk

Wesley Bourne
email: wesley.bourne@morganstanley.com

Fernando Bello
email: f.bello@imperial.ac.uk

Interactive Simulation of Diaphragm Motion Through Muscle and Rib Kinematics

Modelling of diaphragm behaviour is of relevance to a number of clinical procedures such as lung cancer radiotherapy and liver access interventions. The heterogeneity in tissue composition of the diaphragm, as well as the various physiological phenomena influencing its behaviour, requires a complex model in order to accurately capture its motion. In this paper we present a novel methodology based on a heterogeneous model composed of mass-spring and tensegrity elements. The physiological boundary conditions have been carefully taken into account and applied to our model. Thus, it incorporates the influence of the rib kinematics, the muscle natural contraction/relaxation and the motion of the sternum. Initial validation results show that the behaviour of the model closely follows that of a real diaphragm.

Introduction

Diaphragm motion has a crucial influence on surrounding organs. Real time knowledge of its behaviour is of great significance in the context of various treatments for procedure rehearsal, planning and guidance. We are particularly interested in lung cancer radiotherapy and liver access procedures. In the first case, the lung is subject to the up and down movement of the diaphragm, which can in turn produce tumour movement of up to 5cm [START_REF] Shirato | Speed and amplitude of lung tumor motion precisely detected in fourdimensional setup and in real-time tumor-tracking radiotherapy[END_REF]. Liver motion is an important factor to consider in liver access procedures. It has been estimated to be up to 3cm [START_REF] Rohlfing | Modeling liver motion and deformation during the respiratory cycle using intensity-based free-form registration of gated MR images[END_REF]. Such motion is directly related to diaphragm movement, which is transmitted to the liver through its ligamentous attachments.

A possible solution to predict diaphragm motion is to apply displacement vectors obtained from non-rigid registration in combination with a statistical approach as in von Siebenthal et al [START_REF] Siebenthal | Intersubject modelling of liver deformation during radiation therapy[END_REF]. A disadvantage of this approach is that it tends to assume a reproducible breathing cycle, which has been demonstrated to be not always a valid assumption [START_REF] Shirato | Speed and amplitude of lung tumor motion precisely detected in fourdimensional setup and in real-time tumor-tracking radiotherapy[END_REF]. Our aim is to build a model that can be customized to the patient to comply with the accuracy requirements of radiotherapy. Physicallybased models of the diaphragm have already been proposed in the past (e.g. [START_REF] Promayon | Modelling and simulation of the respiration[END_REF][START_REF] Zordan | Breathe easy : Model and control of simulated respiration for animation[END_REF], etc.), but they are mainly based on mass-spring systems, with their motion dictated by a single parameter.

We proposed in [START_REF] Villard | Modelling Organ Deformation Using Mass-Springs and Tensional Integrity[END_REF] a model that takes into account the heterogeneity of the diaphragm, considering its main tissues types, as well the relevant physiological phenomena and various forces involved during normal diaphragm movement, thus introducing additional parameters to fine tune / control the motion. It consists of mass-spring elements for the muscular part and Tensegrity elements [START_REF] Ingber | Opposing views on tensegrity as a structural framework for understanding cell mechanics[END_REF] for the tendinous part. Here we focus on the respiration process. Deformations are produced by three different phenomena: the diaphragm contraction, the ribs kinematics and the sternum motion. This paper presents our methodology and initial validation results on a commercial virtual anatomy model, two radiotherapy patients and a normal volunteer.

Materials and methods

We start with a careful study of the relevant anatomy and physiology as it is important to understand the real behaviour of the diaphragm in order to be able to establish adequate boundary conditions and select suitable physical models. This section also indicates the integration technique used and the method to monitor the external forces.

Anatomy and physiology

Research on previous work about diaphragm physiology was focused on its anatomical, histological and mechanical components. The data obtained were analyzed and compared to determine viability and validity for its use. Based primarily on the work presented in [START_REF] Whitelaw | Shape and size of the human diaphragm in vivo[END_REF][START_REF] Gauthier | Threedimensional reconstruction of the in vivo human diaphragm shape at different lung volumes[END_REF][START_REF] Aladin | Inferences on passive diaphragm mechanics from gross anatomy[END_REF], we obtained the main physiological components involved in displacement of the diaphragm and studied how each one of them influences diaphragm motion. The thoracic diaphragm is a domed fibromuscular sheet that separates the thoracic and abdominal cavities. It is composed of two domes and the right dome is higher than the left dome (Fig. 1 left). The peripheral part of the diaphragm consists of muscular fibres that converge on a central tendon at the level of the xiphisternum. There are three different groups of muscular fibres:1) Sternal, attached to the back of the xiphoid process 2) Costal, attached on the inner surfaces of cartilages and adjacent portions of lower six ribs on either side and 3) Lumbar, attached to the lumbar vertebrae (Fig. 1

right).

Given the spatial dimension, topology and tissue composition of the diaphragm, it can be modelled as a heterogeneous surface composed of muscle and tendon. The origins of the muscle fibres define the attachment points of the diaphragm. They are used to define the boundary conditions for our model as follows: 1) The xiphoid process of the sternum, 2) The costal margin of thoracic wall (lower 6 ribs and their costal cartilages), 3) The end of ribs XI and XII, and 4) The vertebrae of lumbar region. We call them the attachment points and the details of the implementation will be described in the following section.

Deformation Modelling

In their survey of deformable models for surgery simulation, Meier et al. [START_REF] Meier | Real-time deformable models for surgery simulation: a survey[END_REF] claim that mass-spring systems and Boundary Element Models (BEM) are the best models for surgical simulators because they allow a realistic behaviour while supporting near real time performance. Thalmann and Nedel [START_REF] Nedel | Real time muscle deformations using mass-spring systems[END_REF] proposed a model of muscle based on a surface mass-spring system oriented according to the muscular fibre calling it the "action line". We borrow from this formulation the use of a predefined unique line to parameterise the action of a muscle, as well as the massspring formulation, in order to model the muscular component of the diaphragm. Concerning the tendinous part, which can be described as a tough band of fibrous connective tissue, we use a tensegrity model as proposed by Inger to model certain parts of the body [START_REF] Ingber | Opposing views on tensegrity as a structural framework for understanding cell mechanics[END_REF].

The word "tensegrity" comes from a contraction of "Tensional integrity". It is used to define a mechanical system with components that combine tension and compression in such a way as to enable the whole system to receive and apply forces, tensions and pressures. Tensegrity systems are composed of two kinds of elements: the elastic elements, which give the system tension, and the rigid elements, which have a constant length and that will exert compression forces.

The evolution of node positions is simulated by taking into account all of its links that could either be elastic or rigid. Each point is represented at a time t by the triplet (a(t), v(t), p(t)) representing the acceleration, velocity and position. As we define the forces to apply on the nodes with elastic and rigid links, we can calculate their acceleration with the fundamental law of dynamics:

∑ = a f . m (1)
It is then possible to compute the velocity and the position at time (t+Δt):

⎪ ⎪ ⎩ ⎪ ⎪ ⎨ ⎧ + = Δ + + = Δ + ∫ ∫ ∑ Δ + Δ + t t t t ) ( ) ( t) p(t ) ( / t) (t t t t dt t t mdt p v v f v (2)
The elastic links are simulated as in the mass-spring system with two parameters: their initial length l 0 and their elasticity k. A force F e exerting on an elastic link to elongate its length from l 0 to l is thus expressed as follows:

) - -k( e 0 l l = F (3)
We add a damping force F d to systematically and progressively decrease the velocity of each node.

v F d . -γ = (4)
The rigid links in the tensegrity part exert reaction forces to avoid changes in their length. We assume that there is no collision between rigid links. For each time step:

1. The elastic forces F e are computed using Eq. (3) as well as the new positions of the nodes as if there were no rigid links;

2. The rigid constrains are applied to ensure that the distance between two nodes remains constant; 3. Given two nodes A and B linked by a rigid link, A' and B' are their respective positions after applying the classical mass-spring algorithm (as if there was no rigid link). The real positions A 1 and B 1 considering the rigid links are given by the following:

A' A B' B BB AA1 1 1 1 / / =
(5) such that A 1 and B 1 remains on the line A'B'. [START_REF] Promayon | Modelling and simulation of the respiration[END_REF]. The new velocity is computed at each node.

The resulting set of differential equations is solved using a second order Runge-Kutta method:

⎪ ⎪ ⎩ ⎪ ⎪ ⎨ ⎧ Δ Δ - + + Δ + = Δ + Δ + + Δ Δ - = Δ + + Δ Δ - = Δ + = Δ + ∑ 2 ) ( 2 / 1 ) ( ) ( ) ( ) 2 / ( ) ( 2 / 1 ) ( ) ( ) ( 2 / 1 ) 2 / ( ) ( t t
p v p v a v v a v f a (6)
We will now focus on the movement of the diaphragm during inhalation. First, the anterior-posterior diameter increases because the sternum moves forwards as ribs are raised. Secondly, the transverse diameter increases due to a "pump handle" movement which elevates the ribs. Both movements can be simulated by the combined motion of the ribs, the sternum and the cartilage. We first assume that the ribs are rigid and their motion can be modelled by a kinematics law based on the finite helical axis method as in [START_REF] Didier | Breathing Thorax Simulation based on Pleura Behaviour and Rib Kinematics[END_REF]. The cartilage-sternum follows the motion of the ribs.

The attachment points are modelled by a mass-spring system with a very low elasticity in order to ensure that it remains attached. A fall-off function insures that the deformation field is continuous and smooth. The ones belonging to the diaphragm are automatically obtained such that the distance to the group {ribs, sternum, cartilage} is under a given threshold. Then, the vertical dimension increases due to diaphragm movement during inhalation. This motion can be simulated by the muscle action generation. We define an action line as a vertical line (inferiorsuperior) passing through given points and with a certain radius of action. The forces of contraction during inhalation or relaxation during exhalation are computed such that their intensity varies according to the distance from the action line.

Results

We now present the purposely written simulation software. We have applied our methodology to a virtual anatomy model and three patients. The focus for the former was on the realism of the motion and especially on the influence of the parameters described previously, whereas for the later we compared our simulated displacements with the real motion measured on a 4D CT scan. Our framework is coded using Java 3D API 4 . It allows the real-time fine tuning of simulation parameters in order to study various types of breathing, breath hold, ribs rotation, diaphragm muscle force intensity, etc. The deformation framework enables the visualization of the following elements: 1) External skin by transpar-ency (rigid model), 2) The spine (rigid model), 3) The ribs (modelled following a kinematics law), 4) The sternum and the ligaments attached to the ribs (modelled by a very stiff mass-spring system), and 5) The diaphragm (modelled as described above). Fig. 2 shows an illustration of the organ deformation display with examples of anatomy: a commercial virtual anatomy model and three real patients. The parameters of the simulation are: Damping=-0.02, Stiffness= 15, mass point=0.1 and time step=0.01. Results of the simulation and validation using both models are now presented.

Diaphragm Deformation Framework

Simplified model

We used a simplified model from a detailed commercial anatomy set (Anatomium 3D by CF Lietzau 3D Special Service). The anatomical parts are shown in Fig. 1. The diaphragm was remeshed to obtain a surface mesh composed of 11359 nodes and 22722 triangles. The tendinous tissue and the muscle tissue were separated by a plane. The rigid links of the tensegrity part were created by linking the upper surface to the lower surface with a heuristic algorithm. The reasoning behind it is as follows. Once a suitable distance has been chosen, a list of close masses is found. These masses could either be neighbouring masses on the same surface or ideally those on opposing surfaces. Testing the angle removes those on the same surface, leaving a sorted list of masses on the opposing surface.

Fig. 3 is an illustration of the simulation results. The diaphragm geometry at the beginning of the simulation (end of inhale) is represented by a wireframe. The final geometry (end of exhale) is represented by the coloured surface. It can be noticed that the tendon has kept its smooth geometry due to the tensegrity elements, while the muscle part has extended more due to their elastic elements. The influence of the three phenomena mentioned before can be seen here: the muscle relaxation that increases the height of the domes, the ribs rotation that presses the diaphragm on each side and the sternum going backward that pushes the diaphragm.

Now that the realism of the simulations has been evaluated, we focus on the computational cost of the simulation parameters. The tests below were performed on a dual core 2.4 ghz machine with 2 gigabytes of ram and a 256 megabyte graphics card. The operating system was Windows Vista Service Pack 1 with Java 1.6 installed. Table 1 provides a succinct summary of all the data. It gives the average iterations per second after five trials. The relationship between the parameters of the simulation and the mesh size are available. The computational time of the rib kinematics is shown to be quite high as when they are removed in the "Diaph.muscle" trial the average iterations per second hits its peak. Regarding the cost on the speed of the simulator by the addition of tensegrity, it shows that the largest change is visible in the high quality perfect model. This is due to the very high number of connections created in this case.

Patient specific models

Our method can be used on real patients. To test it, we first took 2 sets of 4D CT scan data acquired during normal tidal breathing, therefore the diaphragm motions were not as large as in the virtual model described above. We then acquired 2 MRI datasets from a volunteer subject at two different levels of breath holding.

The diaphragms were manually segmented using ITK-SNAP 5 . A mesh was then generated using the CGAL library6 with a mesh angle of 5, mesh distance of 1, and mesh radius of 1. The mesh was smoothed using vtkWindow-edSincPolyDataFilter() with a pass band of 0.2, followed by a 50\% decimation, followed by another 200 iterations of Windowed Sinc smoothing. An example of a segmented and meshed diaphragm is presented on Fig. 4. The mesh is composed of 6633 nodes and 13267 triangles. The rib cage and the lungs were segmented by the levelset method. Rib motion can be neglected in this case because of the normal tidal breathing and the position of the patient lying down on the scanning table. This assumption has been validated by computing the displacement of the ribs. Therefore, only the diaphragm relaxation needs to be taken into account from the initial exhale to the final exhale. The diaphragm will be attached (as in the virtual model) to the locations described in §2. The simulation is applied with small forces on the previously defined action line and the lungs are linked to the diaphragm with springs.

The simulation was compared to the imaging data at the same point in the breathing cycle by exporting the resulting lung geometry from the fine-tuned simulation and creating a mesh from the segmented 4D CT diaphragm. The differences were then measured and studied with the MESHDEV 7 software that com-putes the distance between two triangular meshes. The result can be seen on Fig. 5. The image on the left represents the distance between the surface of the lungs from the CT scan at beginning and end of inhale. The image on the right represents the distance between the simulated end of inhale and the end of inhale extracted from medical images. Finally, Table 2 provides a succinct summary of the computing time. It can be seen that it is slower than the virtual simulator, attributed to the high count of vertices involved in the simulation.

Discussion

The physiological behaviour of this generic model was discussed with our clinical collaborator, who validated that the antero-posterior, transverse and vertical dimensions of the thorax are increased as illustrated by the model. It was also pointed out that only the domes descend when breathing lightly, whereas with a deep breath, the domes descend further and the central tendon can move from the level of vertebrae T8 to T9. In full exhale, the right dome reaches the 4th intercostal space and the left dome reaches the 5th rib. This complex behaviour further justifies the need for a heterogenous model such as that presented here.

Concerning the patient specific model, it can be seen that the main output of the simulation is the relaxation of the two domes. The contribution is only slightly significant on Subject 1. This is mainly due to the diaphragm displacement (beginning to end of inhale) which is not large enough to correctly appreciate the advantage of the simulation. This displacement is much larger with the 4D CT scan of Subject 2 and with the MRI of Subject 3. The average diaphragm displacements are, respectively, 5.8mm with a standard deviation of 5mm and 16mm with a standard deviation of 9.4mm. The aim of our simulator is to obtain a geometry resulting from the simulation that is as close as possible to the final position of the diaphragm. Fig. 5 shows a clear correlation between both meshes, with average displacements of 4.29mm with a standard deviation of 3.6mm for Subject 2, and 12mm with a standard deviation of 6.7mm for Subject 3. The diaphragm contraction/relaxation modelling is good enough to simulate a realistic behaviour for liver access procedure training, but it would need to be further improved to be used in radiotherapy where precision is crucial.

The performance evaluations of Tables 1 and2 show that the addition of the tensegrity constraint makes a slightly slower simulation, but it is still more than adequately fast. The recorded times were very similar. This demonstrates that, even in large models with hundreds of rigid links, the computational time is not massively increased.

Rib kinematics makes a considerable difference to the speed. This is a result of some large equations being solved to update the locations at each iteration. If a patient is observed to only be using the diaphragm and no rib action during respiration (as subject 1, 2 and 3), a large computational cost can be avoided. Another influencing factor on the speed of a simulation is the size of the mesh. Both virtual and patient models were tested with different resolutions. The evaluation shows an intuitive linear relationship between the number of vertices in the scene and the time taken in computation.

Conclusions

We have presented here a novel methodology to simulate the behaviour of the diaphragm. It is based on a heterogeneous model composed of mass-spring and tensegrity elements. These latter elements are used to rigidify the tendinous tissue. The physiological boundary conditions have been carefully studied and applied to our model: the diaphragm motion follows more than one parameter, thus our model also takes into account the influence of the rib kinematics, the muscle natural contraction/relaxation and the motion of the sternum.

Our implementation runs in real-time and can be fine-tuned by adjusting the simulation parameters, e.g. if the patient is holding her breath or suddenly having a normal tidal breathing. The anatomy and the physiology can also be customized by generating patient-specific segmentations, mesh models and individual diaphragm actions. In the future, we aim to further validate our method and apply it to liver biopsy planning and rehearsal, as well as to monitoring lung motion in a radiotherapy context.

Fig. 1 .

 1 Fig. 1. Left: Thorax anatomy and right: Attachment points of the diaphragm.

Fig. 2 .

 2 Fig. 2. Interface display of ribcage, lungs and diaphragm a) model b,c and d): patients.

Fig. 3 .

 3 Fig. 3. Illustration of diaphragm motion during inhale.

Fig. 4 .

 4 Fig.4. Left: Segmented diaphragm inside CT scan and right: Combination of segmentation and the mesh.

Fig. 5 .

 5 Fig. 5. Distance error measurement on the diaphragm between left: beginning and end of real inhale and right: simulated end of inhale and real end of inhale.

Table 1 .

 1 Iterations per second of all the simulations on the virtual anatomy model. The trials have been averaged into a single value.

	Resolution Vertices	Heterogenous (All-Forces)	Homogenous (All-Forces)	Ribsmotion (Heterogenous)	Diaph.muscle (Heterogenous)
	High	20909 8.13	8.56	8.24	9.68
	Medium 16023 10.41	10.69	10.93	11.76
	Low	12047 12.79	13.07	13.27	14.70

Table 2 .

 2 Iterations per second of all the simulations on a real patient model (subject 1). The trials have been averaged into a single value.

	Resolution Vertices	Heterogenous (All-Forces)	Homogenous (All-Forces)	Ribsmotion (Heterogenous)	Diaph.muscle (Heterogenous)
	High	20909 8.13	8.56	8.24	9.68
	Medium 16023 10.41	10.69	10.93	11.76
	Low	12047 12.79	13.07	13.27	14.70

Pierre-Frédéric Villard , Wesley Bourne and Fernando Bello

http://java.sun.com/products/java-media/3D/

http://www.itksnap.org/

http://www.cgal.org/

http://www.meshdev.sourceforge.net/


Alexis Ballier

Bruno Durand

Emmanuel Jeandel

Tilings robust to errors

Introduction

Tilings are a basic and intuitive way to express geometrical constraints. They have been used as static geometrical models of computation since Berger proved the undecidability of the so-called domino problem [START_REF] Berger | The undecidability of the domino problem[END_REF] by capturing geometric aspects of computation [START_REF] Robinson | Undecidability and nonperiodicity for tilings of the plane[END_REF][START_REF] Berger | The undecidability of the domino problem[END_REF][START_REF] William | Nonrecursive tilings of the plane[END_REF][START_REF] Myers | Nonrecursive tilings of the plane[END_REF][START_REF] Durand | Complex Tilings[END_REF][START_REF] Gurevich | A remark on Berger's paper on the domino problem[END_REF][START_REF] Börger | The Classical Decision Problem[END_REF]. The model assumes the reliability of local color-constraints of tilings, hence several research tracks were aimed at constructing tilesets that are reliable to errors, be it as a computing model [START_REF] Durand | Fixed point and aperiodic tilings[END_REF][START_REF] Durand | High complexity tilings with sparse errors[END_REF] or as a model for DNA self-assembly [START_REF] Winfree | Proofreading tile sets: Error correction for algorithmic self-assembly[END_REF][START_REF] Reif | Compact error-resilient computational DNA tiling assemblies[END_REF].

In this paper we are interested in error robustness of tilesets. We can see an erroneous tiling as an usual tiling by Wang tiles [START_REF] Hao | Dominoes and the ∀∃∀ case of the decision problem[END_REF] where we allow a small proportion of colors on adjacent edges of tiles not to match. We give a more formal definition of this notion in Section 1. Some constructions using fixed point methods construct tilesets that are robust to a small proportion of errors [START_REF] Durand | Fixed point and aperiodic tilings[END_REF][START_REF] Durand | High complexity tilings with sparse errors[END_REF]. Our goal here is to focus on general properties that imply error-robustness or non-error-robustness of some classes of tilesets.

In Section 2 we give a construction of "error-cleaning functions" for tilesets that allow only periodic tilings. We prove that for this class of tilesets it is possible to apply well known hierarchical constructions [START_REF] Gacs | Reliable cellular automata with self-organization[END_REF][START_REF] Durand | On stability of computations by cellular automata[END_REF] so that we can repair every erroneous tiling with probability one granted the probability of errors is sufficiently small.

On the other hand, we prove that the family of tilesets that produce only a countable number of tilings are not locally robust: the correction of a finite number of errors in those tilings may always require a modification of an infinite number of tiles. This is incompatible with local error correction and thus robustness. Locally robust tilesets can be expressed as properties of their ground state configurations as used to model crystals [START_REF] Radin | Global order from local sources[END_REF][START_REF] Radin | Z n versus Z actions for systems of finite type[END_REF].

Finally, by using classical constructions [START_REF] Börger | The Classical Decision Problem[END_REF][START_REF] Robinson | Undecidability and nonperiodicity for tilings of the plane[END_REF][START_REF] Ollinger | Two-by-two substitution systems and the undecidability of the domino problem[END_REF] for encoding Turing machines into tilesets, we prove in Section 4 that these two classes of tilesets are recursively inseparable, which shows that it is not possible to obtain a simple (recursive) characterization of tilesets robust to errors.

Definitions

We present notations and definitions for tilings; we focus our study on tilings of the plane but most of the results naturally extend to higher dimensions. We define tilings by local constraints as they give the most straightforward definitions for tiling errors. Different models are used in literature, such as Wang tiles [START_REF] Hao | Proving theorems by pattern recognition II[END_REF] or subshifts of finite type [START_REF] Lind | Multidimensional Symbolic Dynamics[END_REF], but one can easily transform one formalism into another (see [START_REF] Cervelle | Tilings: recursivity and regularity[END_REF] for more details and proofs).

In our definition of tilings, we first associate a color to each cell of the plane. Then we impose a local constraint on them, that is we describe which colorings are allowed and which are not. More formally, Q is a finite set, called the set of colors. A configuration c consists of cells of the plane with colors, thus c is an element of Q Z 2 . We denote by c(i) the color of c at the cell i ∈ Z 2 . For an element x of Z 2 we denote by c + x the configuration i → c(i + x).

Definition 1 (Patterns).

A pattern P is a finite restriction of a configuration i.e., an element of Q V for some finite domain V of Z 2 . A pattern appears in a configuration c if it can be found somewhere in c; i.e., if there exists a vector t ∈ Z 2 such that c(x + t) = P (x) on the domain of P .

Definition 2 (Local constraints).

A local constraint is a pair τ = (V, δ). V is a finite domain of Z 2 and is called the neighborhood. δ is the constraint function from Q V to {0, 1}.

The idea behind this formalism of local constraints is to define which patterns are allowed and which are not. A pattern is allowed if and only if it maps to 0 by the constraint function. The local constraint τ = (V, δ) naturally extends to a global constraint function ∆ τ : Q Z 2 → {0, 1} Z 2 by applying it uniformly on every cell of the plane: ∆ τ (c)(x) = δ((c + x) |V ). This corresponds to a sliding block code [13, Chapter 1, § 5].

Definition 3 (Tilings).

A configuration c ∈ Q Z 2 is said to be valid for τ = (V, δ) (or a tiling by τ , or allowed by τ ) when it satisfies the local constraint everywhere, that is for every cell x, ∆ τ (c)(x) = 0.

The set of tilings by τ is denoted by T τ . In this paper we only consider tilesets that can tile the plane, thus T τ = ∅ is an implicit condition of all the results.

As we want to study tilings with some errors, this definition of classical tilings naturally extends to tilings with errors by considering an error repartition where some cells are not correctly tiled: Definition 4 (Tiling with errors). Let e be an element of {0, 1} Z 2 and c a configuration, we say that c is a tiling by τ with error repartition e if ∆ τ (c) = e.

One may remark that with this definition, repairing an error is different from replacing the erroneous tile with a correct one as errors may have consequences that require replacing other tiles, even if those tiles are locally correct.

In this paper we first prove that the consequences of such a correction are not problematic in the case of tilesets that allow only periodic tilings (Section 2) but may have consequences on infinitely many cells in the case of tilesets that allow countably many tilings (Section 3).

Before entering the core of the problem of tilings with errors we need to recall a couple of structural results on regular tilings of the plane. We embed Q with the discrete topology, Q Z 2 with the induced product topology. A classical result on the set of configurations is its compactness, as a direct application of Tychonoff's theorem:

Proposition 1. Q Z 2 is a compact perfect metric space (a Cantor space).
The metric we consider, as induced by the product topology, is defined by d(x, y) = 2 -|i| where |(a, b)| denotes |a| + |b| and i is the point closest to (0, 0) (i.e., with minimal norm) such that x(i) = y(i). The ball of center x and radius k is the set of points y such that |x -y| ≤ k. In this paper we reformulate the compactness of sets of tilings in a way that suits our needs: Lemma 1. For every finite subset D of Z 2 there exists C, a finite extension of D such that if P is a pattern defined on C that respects the local constraints then P |D appears in a valid tiling of Z 2 .

This means that if we can tile sufficiently large but finite patterns we are sure that a small part of it will appear in a valid tiling of the whole plane. The function that is given a tileset τ and outputs the set C, even with a fixed D = {(0, 0)}, is uncomputable since this would allow us to decide the domino problem which is a well known non-recursive problem [START_REF] Berger | The undecidability of the domino problem[END_REF][START_REF] Ollinger | Two-by-two substitution systems and the undecidability of the domino problem[END_REF][START_REF] Robinson | Undecidability and nonperiodicity for tilings of the plane[END_REF].

Proof. We prove that for every pattern P there exists a finite domain C P of Z 2 such that there exists a valid pattern defined on C P that contains P if and only if P appears in a tiling of the plane. This proves the lemma by taking C to be the finite union of all C P where P is a pattern of domain D.

If P appears in a tiling of the plane take C P to be the domain of P . Suppose that there is a pattern P that does not appear in a tiling of the plane but such that there exists arbitrary large valid extensions of it. By compactness, this sequence of patterns can be extracted to converge towards a configuration. This configuration contains P and is a valid tiling of the plane.

⊓ ⊔

Robustness

In this section we show how, in the case of tilesets that allow only periodic tilings, it is possible to reconstruct a valid tiling from one with a small proportion of errors by a local modification. The method is the same that was used to prove the error robustness of strongly aperiodic tilesets in [START_REF] Durand | Fixed point and aperiodic tilings[END_REF][START_REF] Durand | High complexity tilings with sparse errors[END_REF]; since sets of periodic tilings have a very simple structure we are able to easily apply the same methods.

We consider error distributions e ε such that each cell of Z 2 has probability ε to be equal to 1 independently of other points, i.e., a Bernoulli distribution.

We first describe a generic process [START_REF] Durand | On stability of computations by cellular automata[END_REF][START_REF] Durand | Fixed point and aperiodic tilings[END_REF] of sorting errors into "islands of errors" that given a repartition taken with a Bernoulli distribution transforms it into an error free repartition with probability one (Section 2.1). Then we present some structural results on tilesets that will allow us to apply this generic process to these tilesets (Section 2.2).

Iterative cleaning process

When ε is small the intuition is that cells with value 1 will be sparse; we will not have big clusters of 1's; however big clusters have a non-null probability to appear thus they almost surely appear on the infinite plane Z 2 . Even if there is such an unavoidable problem, we will see that we can decompose error repartitions into different layers then repair each layer incrementally so that eventually everything gets repaired. First of all, let us define what a layer is: Definition 5. For an error repartition e ∈ {0, 1} Z 2 and a point x ∈ Z 2 such that e(x) = 1, x is said to be in an (i, j)-island of e if there exists no point at distance between i and j in e that has value 1, that is:

∀y ∈ Z 2 , e(y) = 1 ⇒ |x -y| ∈ [i; j]
We denote by I i,j (e) the set of points of Z 2 that are in an (i, j)-island of e. Figure 1 depicts some islands. These islands can be seen as isolated clusters of errors. The idea now is to remove the islands, hence our definition of an erasing function:

Definition 6 (Erasing function). The function E i,j from {0, 1} Z 2 into itself erases the (i, j)-islands of an error repartition:

{0, 1} Z 2 → {0, 1} Z 2 e → c : c(x) = 0 if x ∈ I i,j (e) e(

x) otherwise

Now if we consider integer sequences (α n ) n∈N and (β n ) n∈N , we can think about applying successively the functions E αn,βn to an error repartition so that we obtain an iterative cleaning process [START_REF] Durand | On stability of computations by cellular automata[END_REF][START_REF] Durand | Fixed point and aperiodic tilings[END_REF] that erases small islands of errors, then bigger ones, then even bigger ones, etc. Definition 7 (Iterative cleaning process). Let e be an error repartition and consider two sequences α n and β n . We define the iterative cleaning of e by (α, β) denoted (e α,β n ) n∈N by:

e α,β 0 = e e α,β n+1 = E αn,βn (e α,β n ).

α n β n Figure 1. Islands of rank n
We call "islands of rank n" the points that we corrected at the n th iteration. This operation is pictured on Figure 1 where the black points are the points with value 1, those that will be "cleaned" at this step are marked by a cross, and we can see that we have a kind of security belt between α n and β n . The important part of this process is that cleaning the islands of rank n at the n th iteration creates more islands that we will catch at the n+1 th step. The following theorem catches this phenomenon:

Theorem 1 ( [START_REF] Durand | On stability of computations by cellular automata[END_REF][START_REF] Durand | Fixed point and aperiodic tilings[END_REF]). If the sequences α n and β n match the conditions:

-∀i, 8(β 1 + ... + β i ) < α i+1 ≤ β i+1 - i log(βi) 2 i < ∞
Then there exists ε > 0 (sufficiently small) such that, almost surely, for any element of {0, 1} Z 2 taken with the Bernoulli distribution of probability ε, this iterative cleaning process removes all the ones.

We do not give a proof of this result and refer the reader to the original [START_REF] Durand | On stability of computations by cellular automata[END_REF][START_REF] Durand | Fixed point and aperiodic tilings[END_REF]. It is easy to check that the following sequences match the conditions of Theorem 1 for n sufficiently large:

α n = 2 n 2 β n = nα n
Moreover, what interests us is the fact that lim n→∞ β nα n = ∞, which means that we can have arbitrary large security belts.

Error-robustness of periodic tilesets

In this section we consider tilesets such that their only valid tilings are periodic. We describe how to use Theorem 1 for tilings.

With the process described in Section 2.1, when n is sufficiently large, α n -β n is also large. In our application to tilings with errors, this means that the islands of rank n are surrounded by a belt of width α nβ n that is correctly tiled. However, this zone is only locally tiled correctly. Lemma 1 ensures that when we have a sufficiently large correctly tiled pattern then a smaller part of it appears in a tiling of the whole plane. Hence, if n is sufficiently large, the islands of rank n are surrounded by a belt of large width, say k, where each pattern of size k × k appears in a tiling of the whole plane. One may remark that k may be much smaller than α nβ n but k can still be as big as we would like it to be if we take n sufficiently big. This is depicted on Figure 2. Now we want to replace the erroneous zone by other tiles so that the whole zone is tiled correctly but we are facing a problem: even if we can ensure that the islands of rank n are surrounded by a belt of patterns that appear in a valid tiling of the whole plane, like on Figure 2, how can we be certain that the zone surrounded by the belt can be properly filled by tiles so that no error remains?

A pattern P defined on D is said to be k-extensible if there exists a decomposition of D in (D i ) 1≤i≤n such that:

-∪ 1≤i≤n D i = D -For every i, P |Di appears in a valid tiling of the plane -For every i < n, D i ∩ D i+1 contains a ball of radius k.

With the help of the previous remarks, it is clear that for any k there exists a sufficiently large n such that the islands of rank n are surrounded by a belt that is k-extensible.

It is already known that a tileset allows only a finite number of tiling if and only if it allows only periodic tilings [1, Theorem 3.8], hence we can focus on this simpler case: Lemma 2. If a tileset τ admits only a finite number of tilings then there exists k such that any k-extensible pattern P can be found in a tiling of the whole plane.

Proof. There exists ε > 0 such that the distance between two different tilings is greater than ε. Let k be an integer such that any two configurations that are equal at their center on the pattern defined on [-k; k] × [-k; k] are at distance strictly smaller than ε. We now prove that any k-extensible pattern can be found in a tiling of the whole plane.

Let n and (D i ) 1≤i≤n be such that P is k-extensible over (D i ) 1≤i≤n . We may assume that n = 2, the result for any n being obtained by an easy induction of this simpler case.

Let a be a tiling containing P |D1 and b one containing P |D2 . There exists a ′ and b ′ that are shifted forms of, respectively, a and b such that a ′ and b ′ are equal on P |D1 ∩ P |D2 at their center. a ′ and b ′ are both tilings by τ . Since P |D1 ∩ P |D2 contains a ball of radius k, by our choice of k we obtain that a ′ = b ′ because the distance between them is strictly lower than ε. Therefore P appears in a ′ which is a tiling of the whole plane.

⊓ ⊔ Now, with n sufficiently large, Lemma 2 tells us a bit more on the belt surrounding islands of rank n: there exists such a surrounding belt that appears in a valid tiling of the plane. On Figure 2 this means that the outer (blue) belt is in fact part of a tiling of the plane. This allows us to state our main theorem of robustness: Theorem 2. If a tileset allows only periodic tilings then it is robust to a small probability of errors.

Proof. By Theorem 1 from [START_REF] Ballier | Structural aspects of tilings[END_REF], tilesets that allow only periodic tilings are exactly the tilesets that allow only a finite number of tilings. Hence Lemma 2 applies to these tilesets. By Lemmas 1 and 2, there exists N such that any belt of width N contains a smaller belt that appears in a tiling of the plane.

If we take α n = 2 (n+N ) 2 and β n = (n + N )α n , for every n, the belt has width at least N , hence the (finite) set of points surrounded by this belt can be filled by tiles from the tiling where the belt appears. Moreover these sequences match the conditions of Theorem 1 thus we can repair every tiling taken with a small probability of errors.

⊓ ⊔

We remark that this method of obtaining a valid tiling by surrounding the errors with a security belt also works obviously for a tileset that allows everything (T τ = Q Z 2 ). We can complicate it a little bit by considering Q = {0, 1, 2} and the tilings may contain only 1's and 2's or only 0's. Moreover it is possible to use this method to obtain aperiodic tilings robust to errors [START_REF] Durand | Fixed point and aperiodic tilings[END_REF][START_REF] Durand | High complexity tilings with sparse errors[END_REF].

Non robustness

In the previous section we proved that there exist ways of correcting errors for some tilesets. It would not make much sense to define exactly what we call a "tileset robust to errors" since there may exist other methods for correcting errors but it seems natural that every way of correcting errors has to be local, hence our definition of local robustness: Definition 8 (Locally robust). We say that a tileset is locally robust if for any finite repartition of errors (that contains only finitely many 1's), any tiling with this error repartition can be transformed in a tiling without error by modifying only finitely many cells.

This definition can be related to the ground state configurations used by C. Radin to model crystalline order [START_REF] Radin | Global order from local sources[END_REF][START_REF] Radin | Z n versus Z actions for systems of finite type[END_REF]. Recall that a configuration is said to be ground state if whenever we modify a finite part of it we do not decrease the number of tiling errors in it. In that sense, locally robust tilesets are the tilesets for which ground state configurations are either tilings without any error or have an infinite number of tiling errors.

The structural results from Section 2 ensure that tilesets which allow only periodic tilings are locally robust: find n such that all the errors are in the same island of rank n and replace the finite area surrounded by the belt by a valid one.

Theorem 3. Tilesets that allow countably many tilings are not locally robust.

First recall a structural result about such tilesets: Theorem 4 ( [START_REF] Ballier | Structural aspects of tilings[END_REF]). If a tileset allows a countable number of tilings then it allows one with exactly one direction of periodicity.

Corollary 1. If a tileset τ allows a countable number of tilings then it allows one with exactly one direction of periodicity that can be seen as a bi-infinite word ω xyz ω with |x| = |y| = |z| and y = z.

Proof. Take a configuration c with exactly one direction of periodicity v from Theorem 4. We can represent c as a bi-infinite word w on a finite alphabet Σ. By representation we mean that we can decode letters of Σ in blocks of tiles such that when the blocks from w are repeated along v we obtain c. It is easy to prove that we can obtain from τ a tileset τ ′ over Z such that any tiling by τ ′ can be decoded in such a way in a tiling by τ . Without loss of generality we may also assume that the neighborhood of τ ′ is {-1, 0, 1}: it suffices to code the tileset τ ′ into Wang tiles [START_REF] Cervelle | Tilings: recursivity and regularity[END_REF]. Even after all these codings, Σ is still finite.

Since Σ is finite, there exists j > i > 0 such that w(i) = w(j). Therefore, the word w is of the form xayaz where a is a letter of Σ and x and z are infinite words.

If x is equal to ω (ay), w is equal to ω (ay)z ′ . If x is not equal to ω (ay), x(ay) ω is also allowed by τ . In both cases, we obtain a bi-infinite word w = x ′ (z ′ ) ω that is not periodic and x ′ is an infinite word (the case where x is equal to ω (ay) being symmetric). Since w is not periodic it can be written as w = xy(z) ω where |y| = |z| and y = z. We repeat the same argument on the infinite x part of w to write it in the desired form.

⊓ ⊔

This corollary gives us a configuration that we prove to be incompatible with local robustness. Such a configuration is depicted on Figure 3. Imagine now that we shift half a plane horizontally by one cell. We obtain a configuration like on Figure 4(a). This configuration has only a finite number of tiling errors: around the cell where we broke the vertical line of y's. Now, what can we do in order to correct it? The only solution seems to shift back half a plane in order to restore the vertical line of y's. We now prove it is indeed the only solution by obtaining a contradiction with our hypothesis that the set of all possible tilings is countable.

x x x x x y z z z z x x x x x y z z z z x x x x x y z z z z x x x x x y z z z z x x x x x y z z z z x x x x x y z z z z Proof (of Theorem 3). Let c be a tiling as in Corollary 1: c is of the form ω xyz ω with |x| = |y| = |z| = p and y = z. We modify c by shifting half a plane by p to obtain a configuration like depicted on Figure 4(a). Since this transformation breaks only finitely many tiling rules, suppose that we can correct this by modifying only a finite number of cells. We obtain a tiling like on Figure 4(b) where we have a semi-infinite line of y's in one direction and a shifted semi-infinite line of y's in the other direction, separated by a pattern that repaired the error. We obtained a transformation of the tiling and can repeat it on every sufficiently long vertical line of y's. This transformation gives us a different tiling each time we apply it since y = z. Such vertical lines of y's appear infinitely many times, therefore we obtain 2 ℵ0 valid tilings for τ , a contradiction.

⊓ ⊔

Recursive inseparability of robust and non-locally robust tilesets

Tilesets that allow countably many tilings are not locally robust (Theorem 3) while tilesets that allow only periodic tilings are robust to errors (Theorem 2). In this section we prove that those classes of tilesets are not recursively separable, hence neither are robust and non locally robust tilesets.

Theorem 5. Robust tilesets and non-locally robust tilesets are recursively inseparable.

Proof. We will assume the reader familiar with the encoding of Turing machines into tilesets à la Berger [START_REF] Berger | The undecidability of the domino problem[END_REF], if not please refer to the detailed constructions in [START_REF] Börger | The Classical Decision Problem[END_REF][START_REF] Ollinger | Two-by-two substitution systems and the undecidability of the domino problem[END_REF]. Every tiling by such a tileset contains arbitrary large squares on which we can force the bottom left corner. This corner is where we put the start of a Turing machine computation. We can see the rows of these squares as the Turing machine's tape, time is going bottom-up. Then we do the following:

-If the machine halts and outputs 1 then we force a periodic tiling: when this halting state reaches the border of the square we force a new border such that the only way to tile the plane is to repeat periodically this square. This is exactly what is done in [3, Appendix A]. -If the machine does not halt, the tileset tiles aperiodically with an infinite computation of the Turing machine inside. -If the machine halts and outputs 2 then force the periodicity vertically but allow only a new color, blue, that forces a monochromatic half plane at its left and another color at its right, green, that also forces a monochromatic half plane.

This new tileset always tiles the plane with the uniform blue and green tilings. If the machine halts and outputs 1 then the tileset allows only periodic tilings. If it halts and outputs 2 then it allows countably many tilings because the vertical computation line can appear at a countable number of positions. If it does not halt then it allows an uncountable number of tilings.

The class M 1 of Turing machines that halt and output 1 is recursively inseparable of the class M 2 of TM that halt and output 2: consider the sets C i (i ∈ {1, 2}) of Turing machines that halt on i with their code as input, then a Turing machine M that outputs 2 if its input is not in C 2 and 1 if it is not in C 1 ; there exists such a machine that always halts if we suppose C 1 and C 2 to be recursively separable and M ( M ) gives the contradiction.

With the previous tileset construction, the class M 1 is recursively encoded into tilesets that are robust to errors and the class M 2 into tilesets that are not locally robust.

⊓ ⊔

Conclusions and open problems

In this paper we have shown how it is possible to repair tiling errors (as defined in Section 2) on very simple tilesets (the periodic ones). This correction process relies on the fact that we can surround errors by correct zones which we are sure appear in a valid tiling of the plane; while this is true for periodic tilesets we remark that there exists some other tilesets for which it is also true, it would be interesting to obtain a characterization of such tilesets. We also proved that tilesets that allow countably many tilings are not locally robust; being locally robust is a necessary condition for being able to apply our iterative correction process. While this iterative correction process allows to repair periodic tilesets and even some aperiodic ones [START_REF] Durand | Fixed point and aperiodic tilings[END_REF][START_REF] Durand | High complexity tilings with sparse errors[END_REF] we believe there exists some tilesets that would be locally robust but on which this iterative correction process will not work.

The recursive inseparability between tilesets that we are able to repair and the ones that are not locally robust shows that there is no simple characterization of error-robustness for tilesets.

:

  Locally correct. : Possibly incorrect. : Locally correct and appears in a tiling of the plane.

Figure 2 .

 2 Figure 2. Example of what happens with a sufficiently large security belt.

Figure 3 .

 3 Figure 3. Example of a tiling from Corollary 1.

  x x x y z z z z z z x x x y z z z z z z x x x y z z z z z z x x x x y z z z z z x x x x y z z z z z x x x x y z z z z z (a) Small modification of c x x x x x y z z z z z z z z x x x x x y z z z z z z z z x x x x x y z z z z z z z z x x x x x y z z z z z z z z x x x x x y z z z z z z z z x x x x x x y z z z z z z z x x x x x x y z z z z z z z x x x x x x y z z z z z z z x x x x x x y z z z z z z z x x x x x x y z z z z z z z x x x x x x y z z z z z z z (b) Corrected version of a small modification of c

Figure 4 .

 4 Figure 4. Transformations of c