Ricardo Uribe Lobello
email: ricardo.uribe-lobello@liris.cnrs.fr

Florence Denis
email: florence.denis@liris.cnrs.fr

Florent Dupont
email: florent.dupont@liris.cnrs.fr

Adaptive surface extraction from anisotropic volumetric data: contouring on generalized octrees

Keywords: Surface extraction, Multi-resolution models, Volumetric representation, Image processing

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

Introduction

Volumetric data are already being used as a valuable source of information in multiple domains. It can be obtained from the segmentation of a set of images obtained from a MRI or from 3D rasterization in order to produce a voxel representation of the object. Simi-larly, the recent developments of time-of-flight camera technology makes it possible to create 3D models from objects by using depth information in a scene. Volumetric datasets are often characterized by a strong anisotropy caused by the increasing resolution of images (millions of pixels) produced compared with the small quantity of provided images (usually some hundreds). This will produce highly stretched voxels affecting the ability of conventional octree based algorithms to capture high frequency details of the surface. Voxel data can be represented as a scalar function F so that F (x, y, z) ∈ R | (x, y, z) ∈ Z 3 defined at the nodes of a regular grid Ω. Let v x,y,z be the voxel centred at the index (x, y, z). A labelled volumetric object is defined as the set of voxels V = {∀v x,y,z | F (x, y, z) = l} where l is a value usually established in a segmentation process. However, for rendering and geometrical modelling, it is more suitable to handle an explicit surface representation. This has increased the development of algorithms to extract a piecewise surface approximation from the discrete surface of V , noted ∂V . Generating accurate surfaces from voxels is challenging because many datasets can include degenerated configurations that will induce non-manifold surfaces. Moreover, ∂V can contain sharp features that will force to increase arbitrarily the sampling. In these cases, uniform sampling is not efficient because most of the samples are far from the surface. Therefore, we need a method that adaptively samples ∂V and generates an accurate surface. In order to be used by post processing algorithms as re-meshing or rendering, the extracted surface must be a topological 2-manifold, which means that every edge on the surface has to be shared by at most two polygons while each vertex in the surface must be surrounded by a 1-ring of edges. Additionally, for numerical applications, the quality of individual polygons is an important factor that can strongly affect the performance and convergence of numerical algorithms. In this article, we present a surface extraction method well adapted to anisotropic volumetric datasets. It produces visually appealing meshes while preserving the quality of most of the triangular elements.

Related Work

There is a vast bibliography concerning surface generation methods. Most methods used to extract surfaces from labelled volumetric data are based on the well known Marching Cubes algorithm [START_REF] Lorensen | Marching cubes: A high resolution 3D surface construction algorithm[END_REF]. Dual Contouring methods have been mainly used with implicit functions but some of them have been partially adapted to volumetric data. We briefly summarize relevant methods hereafter. The Marching Cubes (MC) algorithm is based on a regular division of the volumetric data in cubical cells that can be processed separately. All possible intersection patterns have been reduced to 14 pre-calculated cases resumed in a look-up table. However, MC extracts very dense surfaces with many bad-shaped triangles. Several algorithms have made improvements to MC [START_REF] Wilhelms | Octrees for faster isosurface generation[END_REF][START_REF] Shu | Adaptive marching cubes[END_REF][START_REF] Nielson | The asymptotic decider: resolving the ambiguity in marching cubes[END_REF]. Kobbelt et al. [START_REF] Kobbelt | Feature sensitive surface extraction from volume data[END_REF] replace the regular grid by an octree and only process the cells which intersect the surface. Varadhan et al. [START_REF] Varadhan | Topology preserving surface extraction using adaptive subdivision[END_REF] control the octree subdivision using a topology preserving criterion in order to make every cell MC-compatible and obtain consistent adaptive meshes. Kazhdan et al. [START_REF] Kazhdan | Unconstrained isosurface extraction on arbitrary octrees[END_REF] propose an algorithm to extract closed manifold surfaces with MC from unrestricted octrees. However, all these algorithms apply MC triangularion inside cells and generate a lot of bad quality triangles. In addition, as surface nodes are located on the edges of octree cells, sharp features are lost during reconstruction. Dual contouring (DC) methods (see Ju et al. [START_REF] Ju | Dual contouring of hermite data[END_REF][START_REF] Schaefer | Dual contouring: The secret sauce[END_REF]) are able to build adaptive meshes and to capture sharp features in the presence of Hermite Data (surface normals). In these methods, surface nodes are not created on cell's edges but one node is generated inside every cell. However, this approach can produce topological artefacts and non-manifold configurations as illustrated in figure 2. By limiting DC octree adaptivity, Varadhan et al. [START_REF] Varadhan | Feature-sensitive subdivision and isosurface reconstruction[END_REF] have proposed an improved DC that avoids non-manifold configurations. Then, in Dual Marching Cubes (DMC), Nielson [START_REF] Nielson | Dual marching cubes[END_REF] used a look-up table to produce surfaces that are dual to those created by Marching cubes by allowing cells to contain more than one dual node in order to solve the non-manifold vertex configurations. However, we observe that by simply applying the DMC algorithm, non-manifold edge configurations can still appear between cells sharing an ambiguous face (see figure 2). Concerning non-manifold configurations, several approaches propose to create dual grids extracted from an octree data structure. These grids are aligned to features of the surface by using Hermite data on the octree cells. Schaefer et al. [START_REF] Schaefer | Dual marching cubes: Primal contouring of dual grids[END_REF] have proposed to use a dual hexahedral grid combined with marching cubes (Dual Marching Cubes DMC2) to extract the surface. However, as dual cells are not necessarily convex, this algorithm can produce self-intersected surfaces (see figure 3). In addition, aligning the grid to shape features can make many degenerated triangles appear. In order to avoid self intersections, Manson and Schaefer [START_REF] Manson | Isosurfaces over simplicial partitions of multiresolution grids[END_REF] built a dual tetrahedral grid that do not contain inverted tetrahedra. Then, they applied a marching tetrahedra (MT) to extract the surface. This solution produces intersection-free surfaces, unfortunately, using MT, combined with a feature aligned grid produces meshes with many bad-shaped triangles. Finally, Manifold Dual Contouring (MDC) [START_REF] Schaefer | Manifold dual contouring[END_REF] produces adaptive manifold meshes by using DMC over a regular octree and a topology preserving criterion that maintains the topology of the original object. Nevertheless, this approach can still generate self-intersected surfaces and its bottom-up strategy for mesh simplification is not suitable to handle strongly anisotropic datasets.

In this article, we are mainly interested in adaptive surface extraction from the voxel representation of an object. We want to obtain manifold meshes at several resolution levels without creating too many degenerated triangles. Therefore, we have considered a dual approach as the most suitable option. We present a dual algorithm to produce multi-resolution and adaptive surfaces that is appropriate to volumetric data. It produces meshes with good-shaped triangles and is able to approximate sharp features. The dual approach is inspired by Nielson's Dual Marching Cubes (DMC) [START_REF] Nielson | Dual marching cubes[END_REF]. This algorithm allows us to place several surface nodes inside each octree cell and gives us more freedom in nodes localization in order to better approximate the underlying surface ∂V . Our solution extends DMC over a generalized octree data structure to extract surfaces adapted to the geometry of ∂V . We start by defining some preliminary concepts and the notation that is going to be used. Let be V the volumetric object as defined in the introduction and ∂V its boundary. In this paper, a cell is defined as an axis-aligned hexahedral block containing voxels and it is equivalent to the cell of an octree. Hereafter, our octree cells are going to be noted as cells. Cell volume will indicate the voxels contained inside a cell and we intend to check only the values of the eight voxels in the corners of every cell. The vertices of every cell are labelled to indicate whether they are inside the object.

Definition 1 An ambiguous face exists when the face contains two vertices inside ∂V interleaved with two vertices outside ∂V .

Our algorithm is able to consider at most four surface components inside a single cell. Every component must have a single connected boundary and be well approximated by a plane.

Octree data structure

In order to adapt our space subdivision approach to the anisotropy of volumetric data, our octree data structure is not the usual one: the octree cells do not have to be always subdivided in the tree main axial directions. In consequence, we can choose to subdivide a cell in only one or two axial directions as illustrated in the figure 4 for a quadtree. For the multi-resolution data structure, we have implemented a hashed octree that uses a Morton space filling code [START_REF] Lewiner | Fast generation of pointerless octree duals[END_REF] as index. A Morton code can compactly encode hierarchical and spatial information about cells. However, this code have been conceived to work in a regular grid, so, to be usable in a generalized and adaptive data structure, it must be completed with some spatial coordinates information. Our approach is based on the idea that, even if an octree cell c is not divided in all directions, the sub-cells S contained inside it must fulfil, or tile, all the space inside c. This will also mean that all Morton codes indexing sub-cells inside c must be pointing to cells in S. This distribution has to be done with respect to the space occupied for every subcell as illustrated in figure 4 for cell indexed by 101. This characteristic of our octree construction algorithm combined with the spatial coordinates information contained inside cells allows us to access any neighbouring cell even if the octree has not been divided regularly. Our algorithm is shown in 1.

Algorithm 1: Get Neighbour algorithm input : Cell c and a face f . output: A face adjacent cell a through the face f .

level ← GetLevel (c); coincidentFace ← GetCoincidentFace (f ); offset ← 1; while offset ≤ level do a ← GetAdjacentCellByFace (c,f ,offset); if Exists (a) then cellLimitC ← GetCellLimitByFace (c,f ); cellLimitA ← GetCellLimitByFace (a,coincidentF ace); if cellLimitC == cellLimitA then return adjacentCell; else offset = offset +1; end end end
Our function GetAdjacentCellByFace returns a cell face adjacent to the provided cell. This is done by obtaining the spatial coordinates of the current cell, adding an offset in the corresponding direction and converting the spatial coordinates into a Morton code. As anisotropic subdivision is usually done at the deepest levels of the octree, our offset will not be larger than the current depth of the octree. If no neighbouring cell exists, our method looks her parent cell by shifting the three least significant bits of its Morton code. The method GetCel-lLimitByFace returns the geometrical limit of a cell in a direction deduced from the face provided. This ensures that two cells are really face adjacent because the irregular subdivision of the octree can produce cells that are close in terms of Morton encoding but far geometrically. Finally, as our algorithm is strongly dependent on the access to neighbouring cells, we use a hash table to index the octree cells and make almost constant access times.

Octree construction

The octree construction uses an iterative top-bottom algorithm that builds an octree based on ∂V topological and geometrical characteristics. First, to capture ∂V general features, we regularly divide the octree until some user defined minimal depth (usually 2 or 3) is reached. Then, we apply several criteria in order to decide if a cell must be subdivided: we check if the octree cell is inhomogeneous, then, we apply a complex cell criteria as explained below.

Definition 2 Let C be a cell. C is complex if at least one of its faces is complex. A cell face is complex if:

-It intersects ∂V and all its vertices are either inside or outside V (see figure 5 left). -At least one of its edges intersects ∂V more than once (see figure 5 right). A Complex face criterion allows us to detect cells that contain a piece of the surface that is smaller than the current level of subdivision. In the octree construction context, this criterion is implemented by extracting the connected components of the voxels contained in every homogeneous cell face (all vertices are inside or outside ∂V ). In addition, a complex edge criterion detects configurations that can induce non-manifold surfaces when a given cell edge intersects ∂V more than once. This strategy allows us to only check voxels that belongs to cell faces to ensure that all cells intersected by ∂V are going to be processed.

Rule 1 Subdivision Rule: If a cell C is complex, it has to be subdivided. Otherwise, it is no longer subdivided and it is marked as an octree leaf.

The application of the subdivision rule 1 will build an octree where no leaf cell that intersects ∂V is complex. From this octree, we can extract closed manifold surfaces at any level of resolution.

Adaptive surfaces and manifold verification: in order to adapt our surfaces to the curvature of ∂V , we propose to use a second rule based on the normals at the intersection points of cell edges with ∂V . To calculate the normals, we have used a surface approach based on the measure of the average of unitary normals in the neighbourhood of an edge intersection with ∂V [START_REF] Flin | Adaptative estimation of normals and surface area for discrete 3d objects[END_REF]. The previous rule gives valuable information about the shape of ∂V inside the cell. As consequence, δ close to 1 will produce an almost regular mesh. On the contrary, δ close to 0 will generate a highly simplified surface. The curvature rule 2 is applied once a cell has succeed the complex cell test. In this way, a normal can be calculated at the unique intersection point of every edge that traverses ∂V and the curvature rule 2 can be checked. Cells that have not been subdivided because the curvature of ∂V exceeds δ inside them are marked as compact.

Our subdivision rule 1 ensures that our octree will be subdivided until there are no more complex cells. The limit case will be a regular divided octree, without any non-manifold vertex configuration because of the dual equivalence between the Marching Cubes and the Dual Marching Cubes look-up tables. However, the addition of the rule 2 for the curvature can lead to nonmanifold configurations. This is due to the local application of subdivision rules 1 and 2. A cell c is not subdivided if it is not complex and its curvature is less than a user defined threshold. That means that the topology of ∂V inside c can be well approximated by using one of the cases of the DMC look-up table. However, c neighbouring cells can keep subdividing until revealing a more complex topology that was not visible from c. This can include ∂V multiple traverses through c edges (see figure 6

top).

In order to solve this, we propose to extract information about ∂V shape in face adjacent cells over the faces of compact cells. Even if Westermann et al. [START_REF] Westermann | Real-time Exploration of Regular Volume Data by Adaptive Reconstruction of Iso-Surfaces[END_REF] have proposed a similar approach to fix surface cracks by using an extension of the marching cubes triangulation, their approach is limited to neighbouring cells with at most one level of difference in the octree. Our approach does not have this limitation and it is just guided by the topological complexity of ∂V on the faces of compact cells.

First, for every compact cell c, we get the adjacent cells at the same octree level. Then, for every adjacent cell through a face f , we extract all its leaf cells S that are face adjacent with c. By using the cells in S we build a quadtree with all the faces of cells in S that are adjacent to c. This quadtree contains the values of F (x, y, z) over its vertices and can represent the intersection between V and f produced by the octree subdivision (see figure 6 middle). So, we can apply a complex face and complex edge test and determine if the face is complex. If so, we divide the adjacent compact cell and we add its sub-cells to the compact cells list (see figure 6 middle). Our algorithm is entirely based on the information provided by the evaluation of F (x, y, z) in the octree vertices. We have implemented this procedure recursively and its stop guarantee is provided by the fact that Fig. 6 2D illustration of our local refinement and manifold verification algorithm. First row shows an object in black contained inside two cells (left). Then, we apply our refinement algorithm over the right cell (center). This octree cells configuration will generate a non-manifold configuration (right) with dual nodes (in red) contained in not refined (compact) cells (in light blue). Middle row, we extract the value of octree vertices that are contained inside grey cells edges (left), if those edges are complex, we refine the compact cells (center) and we update the list of compact cells (in grey) to check (right). Bottom row, we repeat this test recursively (left) until verified all compact cells (center). Then, we can extract the surface (right).

neighbouring cells that have the same size will be covered by the DMC look-up-table. In practice, our algorithm stops when all compact cells having complex faces or edges have been eliminated by subdivision (see figure 6 bottom). Proof of correctness: Let state that in Dual Marching Cubes, all the patches have to pass through intersected edges. In the case of a complex edge, let be C a compact cell at depth n and S a set of cells at level {n + x | x > 0}. Let be E an edge of S and e i the edges of cells in S so that {∀e i | e i ⊂ E} and {∀j = i | e j ∩ e i = ∅}. Without lost of generality, we can say that if we force every edge in C to traverse ∂V only once, we also assure that only one cell in S is traversed by ∂V . This comes from the fact that all edges e i ∈ S are disjoint. Then, ∂V can traverse just one cell edge in S. In the case of a complex face, let F be a face of a compact cell. Then, if we assure that no isolated connected component exits inside F , there cannot be intersected edges to connect with the dual nodes created for the simple connected components around the edge of the compact cell. Our local refinement algorithm can be used to produce surfaces that are not only adapted to the curvature but guided by external criteria. For some applications, a model has to be refined locally in some areas corresponding to the interest of the user. An example of this feature is illustrated in figure 7, where our algorithm is applied to the "Lucy" surface model that is progressively refined to obtain a more detailed surface on the top of the model. Our algorithm guarantees that the refined model is still a closed valid manifold. Finally, in order to improve the general quality of triangles in our surfaces, we can keep the level differ-ence between octree leaf cells in 1, even though, our algorithm is not limited by this.

Dual nodes generation

Once we have applied the manifold verification algorithm, we use the Dual Marching Cubes look-up-table (LUT) (see figure 8) to create the surface nodes inside our octree cells. However, the application of the DMC LUT does not guarantee that all edges connected by two dual nodes are shared by at most two polygons. If we see a cell as a graph with eight vertices and twelve edges, DMC inserts dual nodes based on the number of connected components that are inside ∂V and generates non-manifold edge configurations in cells sharing an ambiguous face. In figure 9 left, there is only one connected component inside ∂V (in blue) and, in consequence, only one dual node is created.

We have solved this problem by considering the connectivity of the graph of vertices that conforms a cell (eight vertices and twelve edges). As showed in the figure 9, blue vertices are inside the surface ∂V and white ones are outside. Then, a blue connected component is composed by one or more blue vertices that are connected by edges which two endpoints are blue. Respec- tively, white connected components are defined in the same way for white vertices. For each cell c with a single blue component and one or multiple ambiguous faces, we enumerate the adjacent cells A c through the ambiguous faces that also have a single blue component. Then, for each cell in A c , we extract the white connected components. In these cells, the existence of an ambiguous face implies the existence of at least two such white connected components. Therefore, we use the white components to generate the dual nodes in c and in all adjacent cells in A c . This strategy produces two vertices instead of one and solves the problem recovering the tunnel in ∂V (see figure 9 In the cases (a) and (b), there is only one ambiguous face per cell. The use of blue connected components will produce a non-manifold edge. Therefore, we use the white connected components to extract the surface patches. In (c), there is a cell with three ambiguous faces (right) and more that one blue connected component, then, it is possible to use the blue components to produce the patches. In (d), two cells with three ambiguous faces containing more that one blue component that can be used to produce the dual nodes and extract the surface patches.

Locating dual nodes inside cells instead of cell edges produces surfaces having less waving and reconstruc- Fig. 11 Comparison of our algorithm with DMC2 [START_REF] Schaefer | Dual marching cubes: Primal contouring of dual grids[END_REF] and MC [START_REF] Lorensen | Marching cubes: A high resolution 3D surface construction algorithm[END_REF] on a volumetric "Block" dataset of 512 × 512 × 128 voxels with a curvature parameter λ of 0.9 and a maximal octree depth of 7. DMC2 has problems to capture some of the sharp edges and produces more degenerated triangles. tion artefacts. In addition, dual nodes inside cells provide more freedom in order to capture sharp edges. To illustrate this, in figure 11, we compare the result of our algorithm with DMC2 [14] on a mechanical "Block" volumetric object of 512 × 512 × 128 voxels. As it can be seen, our algorithm produces smoother surfaces while avoiding reconstruction artefacts especially visible in sharp edges. This is mainly explained because, even if DMC2 aligns the dual grid to the features of the volume, as it uses MC to extract the surface, it localizes the surface nodes on the edges of topological cubes reducing its ability to capture sharp edges. On the contrary, even if our algorithm uses an octree, as it localizes nodes inside cells, it has more freedom to place surface nodes close to the features of the V .

Connectivity generation

A minimal edge is one that does not contain any other edge in a neighbouring octree cell (see figure 12). In dual contouring algorithms, there is a one to one relationship between the surface patches and the minimal edges that are intersected by the surface. Then, in order to extract the connectivity of the surface, we used an algorithm proposed in Dual Contouring [START_REF] Ju | Dual contouring of hermite data[END_REF] that traverses the octree to enumerate the minimal edges intersected by the surface. This method is based on three recursive methods: CellProc, FaceProc and EdgeProc. CellProc receives a cell c as parameter, FaceProc receives two face adjacent cells and EdgeProc receives four edge adjacent cells (in the regular case, three in the adaptive case). The algorithm works as follows, CellProc receives a cell c and calls itself recursively in every sub-cell of c. Then, it makes twelve calls to Face-Proc with every pair of face adjacent sub-cells. Finally, it makes six calls to EdgeProc with all four sub-cells sharing an edge (see figure 13 top). FaceProc receives two cells sharing a common face f and calls itself four times with every pair of subcells sharing a face contained in f . Then, it makes four calls to EdgeProc with every four sub-cells that share an edge contained in f (see figure 13 middle). Finally, EdgeProc receives four edge adjacent cells and makes two recursive calls to EdgeProc with all four sub-cells sharing a half-edge contained in the edge (see figure 13 bottom). Only EdgeProc generates surface patches traversing minimal edges. When EdgeProc receives four cells that share a minimal edge, it connects the dual nodes in those cells that are related to the minimal edge to form a quad in the uniform case (see figure 14 left) and a triangle in the adaptive one as illustrated in figure 14 right.

Dual nodes localization

Dual nodes localization is important to improve the geometrical approximation of the surface. A good localization algorithm will strongly avoid, or at least reduce, aliasing on the surface. In addition, our method must be robust enough to be applied over the usually noisy volumetric data surfaces. We propose a dual node localization algorithm based on the centroid b s of every connected component of V inside a cell. Let be C a leaf cell in the octree where V ∩ C = ∅. As we know the dimensions of C and its mass M , we can obtain its centroid b c . Then, we calculate the centroid of V ∩ C named b s and its mass m s . Finally, the mass and barycentre of the complementary b e of (V ∩ C) can be calculated using equation M b c = m s b s + m e b e . Then, we use b s , b e and their respective masses to estimate an initial position v p for the dual node that will lie on the segment between b s and b e . However, in some configurations, this initial position v p can lie far from the surface (either inside or outside ∂V ). To approach the dual node to the surface, we use the segment defined by b s and b e and noted as S. First, we move the dual node from its initial position v p to the intermediate point of the segment connecting v p with the endpoint of S that traverses ∂V , noted v e . Then, we test if this intermediate point is over the surface, if it is the case, we stop, if not, we apply recursively the same operation to the half-segments of the segment v p to v e . In order to validate if the node v is over the surface, we use the voxels spacing as an interval and we check if the current node position falls into an interval that makes part of the discrete surface of ∂V . Dual node localization process is illustrated in figure 15. Our dual node localization method has two main characteristics: it sticks dual nodes to the volumetric surface of the component and, as each surface component inside a cell must have a simple connected boundary, it tends to localize the node inside the kernel or center area of the surface which points can be connected to any other point of the surface by using a line [START_REF] Varadhan | Topology preserving surface extraction using adaptive subdivision[END_REF]. These properties are useful to avoid self intersections and to improve the general quality of triangles.

In figure 16, we compare our approach with the Dual Marching Tetrahedra (DMT) algorithm [START_REF] Manson | Isosurfaces over simplicial partitions of multiresolution grids[END_REF]. In figure 16a, DMT generates a good approximation but the mesh and the histograms of angles show a lot of triangles degenerated having very small (under 5 • ) minimal angles or huge maximal angles (over 100 • ). This is because DMT uses a dual tetrahedral grid aligned to the features of ∂V combined with Marching Tetrahedra triangulation. On the other side, as it can be seen in figure 16b, our approach produces a more homogeneous min- imal and maximal angle distributions around 30 • and 85 • respectively, producing far less degenerated triangles.

We considered other localization methods such as Quadric Error Functions [START_REF] Garland | Surface simplification using quadric error metrics[END_REF], however, these methods can easily localize surface nodes outside the original volumetric object (generating self-intersections) and be strongly affected by noisy segmentation processes that will not produce a smooth discrete surface. Our experiments have shown that our method works well on smooth or noisy segmented surfaces.

Results

One of the main characteristics of our solution is its ability to generate closed manifold surfaces at any resolution level. It can be used to produce multi-resolution models based on a user defined resolution. In figure 17, several "Armadillo" models are built based on a maximal octree level from 4 to 8. These models do not necessary have the same topology (see red circles around fingers) because thin features cannot be captured at lower resolutions and some sections of surface can be separated in the final mesh. However, this feature can be useful in visualization where the speed of transmission is more important than the topology of the model. On the contrary, when the topology of the model is an important factor, the curvature threshold can be used to obtain simplified models that will keep the same topology of the original volumetric object. This is shown in figure 18 where several topological equivalent surfaces are extracted from an "Armadillo" model of 512× 512×256 voxels. Observe that all surfaces keep the same topology and no thin components are separated.

RMS distance (measured with Metro [1]

) in the different models are not strongly affected by the curvature simplification and a good approximation is still reached with a high level of simplification. Times in the table include the octree construction. The most expensive step is the normals calculation and our dual node localization algorithm because they are calculated directly on volumetric data.

In order to compare our solution with other relevant methods in the literature, we used the Adaptive Marching Cubes (AMC) of Kazhdan et al. [START_REF] Kazhdan | Unconstrained isosurface extraction on arbitrary octrees[END_REF] and the Dual Marching Cubes (DMC2) of Schaefer et al. [START_REF] Schaefer | Dual marching cubes: Primal contouring of dual grids[END_REF]. We applied these algorithms over a set of volumetric datasets extracted from classic polygonal models (see our recon-Fig. 18 Simplified surfaces with respect to the curvature threshold. Surfaces are extracted from an "Armadillo" volume of 512 × 512 × 256 voxels using an octree of maximal level 7. Observe how thin features as the fingers are connected in all the surfaces (see red circles [START_REF] Schaefer | Dual marching cubes: Primal contouring of dual grids[END_REF]. All surfaces are extracted with an octree of maximal depth of 8 in order to capture the finer components of the surface and a curvature λ parameter of 0.9 to reduce the number of triangles in the flat areas of the volume. * A degenerated Triangle is defined as having a minimal angle of at most 2 degrees.

structions in figure 19). We set the curvature parameter λ in 0.9 to force a cell subdivision if there is a slight curvature in the surface ∂V , and the maximal depth of the octree in 8 so to guarantee that we capture the main features of the models. All these results are summarized in table 2.

As it can be seen, even if all three algorithms are almost equivalent in mesh size and minimal angles, our algorithm generates consistently far less bad-shaped triangles defined as triangles containing minimal angles smaller that 2 degrees. In the case of AMC, MC triangulation is used inside cells and its technique to close the surface where the resolution changes usually produces a lot of bad-shaped triangles. On the other side, DMC2 builds a dual hexahedral grid from an octree using a dual node inside every octree cell (see [START_REF] Schaefer | Dual marching cubes: Primal contouring of dual grids[END_REF]) but it still uses the MC triangulation. In addition, dual cells are not necessarily convex and their shapes can be very twisted. On the contrary, our method relies on the octree cells, uses multiple dual nodes inside every cell and provides more freedom in surface nodes localization.

Our experiments over multiple curvature and octree depth parameters showed that our method keeps generating surfaces with far less degenerated triangles. This is understandable because the main reasons for creating bad shaped triangles in AMC and DMC are their restrictions in the nodes localization which is not fundamentally affected by a higher resolution (higher octree depth), neither a more regular space subdivision (bigger curvature threshold).

Conclusions and perspectives

We have presented a robust surface extraction algorithm that can be used to obtain adaptive surfaces from highly anisotropic volumetric data. Our method operates directly on volumetric data and does not need any pre-calculated information. We have seen that almost all methods tackle one or two particular issues but they do not solve all the iso-surface extraction problems at once. It is clear that is not easy to find a way to address all limitations but our method provides a good trade-off between surface quality and geometrical accuracy. As our solution can obtain manifold meshes at any level of resolution, it can be used for data exploration, visualization and transmission applications. In future work, we would like to improve the surface approximation by using more precise curvature estimators inside octree cells. Another interesting possibility is to use curvature information to move the division planes to the more curved zones of the surface instead of dividing octree cells at middle points to better capture the surface features. In addition, we are implementing an out-of-core version of our solution combined with parallel strategies in order to process huge datasets.

Fig. 1 (

 1 Fig. 1 (Left) Non-manifold vertex configuration and (right) non-manifold edge configuration in original Dual Contouring algorithm.

Fig. 2 (

 2 Fig. 2 (Left) Two DMC look-up table cases 17 connected by an ambiguous face. (Right) Non-manifold edge configuration (line in red).

Fig. 3

 3 Fig. 3 Schaefer et al. DMC2 self intersection problem. Vertices in blue are inside ∂V (white are outside), the surface is the blue dotted line, red circles are dual nodes and grey lines are dual cell edges. Yellow lines are (self-intersected) surface edges that are generated connecting surface nodes located on dual cell edges.

Fig. 4

 4 Fig. 4 Generalized quadtree data structure (left) and its representation as Morton codes (right).

Fig. 5 (

 5 Fig. 5 (Left) a complex cell with a tunnel and (right) a complex cell with a complex edge.

Rule 2

 2 Let C be a cell and I i=0...n its intersection points with ∂V . Let be n i=0...n the normals on I i=0...n . If Min(n i • n j |i = j} ≤ δ, C must be subdivided. δ is a user provided parameter with values in the closed interval [0, 1].

Fig. 7

 7 Fig.7Progressive "Lucy" (mesh and surface) model extracted from a 512 × 256 × 128 voxels volume. A more refined mesh is produced adaptively on the top of the model.

Fig. 8

 8 Fig. 8 DMC Look-up table of surface patches generated by using connected components on the graph defined by the cell vertices. Blue vertices are inside ∂V , external faces are orange, internal faces are grey. Dual nodes are noted in red.

Fig. 9 (

 9 Fig. 9 (Left) Non-manifold edge generated through an ambiguous face since two portions of the surface share a common edge (in red). (Right) Connectivity inversion, two dual nodes are created in each cell and a tunnel in ∂V is recovered. Dual nodes are noted in red.

  ). Some examples of patches generated by our algorithm are shown in figure 10.

Fig. 10

 10 Fig.10Four cases where two cells share an ambiguous face. In the cases (a) and (b), there is only one ambiguous face per cell. The use of blue connected components will produce a non-manifold edge. Therefore, we use the white connected components to extract the surface patches. In (c), there is a cell with three ambiguous faces (right) and more that one blue connected component, then, it is possible to use the blue components to produce the patches. In (d), two cells with three ambiguous faces containing more that one blue component that can be used to produce the dual nodes and extract the surface patches.

  (a) Dual marching cubes (DMC2) with 54.888 triangles. (b) Our algorithm with 46.888 triangles.

Fig. 12

 12 Fig. 12 Examples of minimal edges. The minimal edges does not contain any other smaller edge. Some examples are showed in blue. Non minimal edges contain more that one edge in edge adjacent cells of smaller size (in orange). Not all minimal (or not minimal) edges are highlighted in this figure.

Fig. 13

 13 Fig. 13 CellProc, FaceProc and EdgeProc methods as they are used to extract the topology.

Fig. 14

 14 Fig. 14 EdgeProc patch creation. In the regular case, Edge-Proc receives four cells (left) and creates a patch with the four minimizers (in red) of edge adjacent cells. (Right) In the adaptive case, EdgeProc receives one of the cells repeated, so, only three cells and it creates a triangular patch. Blue vertex is inside ∂V , all others are outside.

Fig. 15

 15 Fig. 15 2D illustration of the algorithm to estimate a good dual node localization by using connected components. Cell centroid (red), V ∩ C centroid (green), (V ∩ C) c centroid (blue), dual node (yellow).

  (a) Dual marching Tetrahedra with 133.512 triangles and average min/max angles of 16.5/103.4. (b) Ours with 127.188 triangles and average min/max angles of 33.1/89.6.

Fig. 16

 16 Fig.[START_REF] Uribe Lobello | Génération de surfaces adaptatives à partir de données volumiques binaires[END_REF] Comparison between the DMT algorithm of Manson and Schaefer and our algorithm on a volumetric "Block" dataset of 512 × 512 × 256 voxels with a curvature parameter λ of 0.9 and a maximal octree depth of 8. Histograms show that our algorithm obtains much better min/max angle distributions.

Fig. 17

 17 Fig. 17Multi-resolution "Armadillo" model (512 × 512 × 256 voxels) with a curvature threshold of 0.9 and five different resolutions based on a maximal octree depth from 4 (1.458 triangles) to 8 (347.922 triangles). The topology of the surface changes through different resolutions (see red circles).

  Fig. 17Multi-resolution "Armadillo" model (512 × 512 × 256 voxels) with a curvature threshold of 0.9 and five different resolutions based on a maximal octree depth from 4 (1.458 triangles) to 8 (347.922 triangles). The topology of the surface changes through different resolutions (see red circles).

Fig. 19

 19 Fig. 19 From left to right: Dragon, Horse, Buddha, Armadillo, Block and FanDisk models. All models are extracted from a 512 × 512 × 128 voxels representation.

Algorithm 2 :

 2 Octree constructioninput : Root cell c. Minimal minLevel and maximal maxLevel level subdivision. output: A regular octree until level minLevel and an adaptive octree subdivision between minLevel and maxLevel.

	Add c to the list toP rocess;
	while toP rocess is not empty do
	currentCell ← First(toProcess)
	n ← Level (currentCell)
	if (n≥ minLevel) and (n ≤ maxLevel) then
	divide ← no
	if Complex (currentCell) then
	divide ← yes
	else
	if ( Inhomogeneous (currentCell) and
	( HighCurvature (currentCell)) then
	divide ← yes
	else
	Mark currentCell as a compact leaf.
	end
	end
	if divide then

Comment: currentCell has to be divided; S = {q1, q2, q3, ..., q8} ← Divide(currentCell); foreach cell qi in S do Add qi in toP rocess end end end end

Table 1

 1 resumes sizes, approximation and timing statistics for the surfaces of figure18. Hausdorff and

	Curvature #Triangles Hausdorff RMS Time(secs)
	0.9	89094	0.0043	0.15	9.1
	0.7	42822	0.0091	0.21	8.3
	0.5	27274	0.012	0.33	7.3
	0.3	20740	0.014	0.41	6.5
	0.1	16922	0.014	0.45	6.3

Table 1

 1 Statistics obtained for the surface extraction of a volume "Armadillo" (512 3 voxels) with different curvature values. Surfaces are always a closed manifold.

Table 2

 2 Comparison of our method with Adaptive Marching Cubes (AMC) [Kazhdan et al. 2007] and Dual Marching Cubes (DMC2)

	Criterion		# Triangles		# Degenerated* Triangles	Minimal Average Angles
	Models/Methods	AMC	DMC2	Ours	AMC	DMC2	Ours	AMC	DMC2	Ours
	Dragon	122074	107342	133586	2360	11942	140	31.4	28.2	33.1
	Horse	36408	42316	61262	850	5020	162	31.23	27.6	32.4
	Buddha	119124	134342	165921	1411	12459	271	31.9	28.8	32.84
	Armadillo	137098	156820	90773	3160	15741	173	31	27.88	36.8
	Block	55250	71144	57796	1136	8659	58	33.94	29.86	31.4
	FanDisk	43094	52344	46896	1252	7040	63	32.53	27.6	33.4

). Statistics of the extraction process for these surfaces are provided in table 1.

Acknowledgements

The authors would like to thank the Stanford 3D Scanning Repository and Georgia Tech Large Geometric Models Archive for the models used in this article. We would also like to thank the reviewers for their

valuable commentaries. This article is an extended version of an article [START_REF] Uribe Lobello | Génération de surfaces adaptatives à partir de données volumiques binaires[END_REF] presented in the french conference CORESA 2012.


Georges-Henri Cottet
email: georges-henri.cottet@univ-grenoble-alpes.fr

High Order Semi-Lagrangian Particle Methods

published or not. The documents may come

masses give exact weak solutions to advection equations written in conservation form. For simple linear advection equations, the approximations of exact solutions by particles, measured on distribution spaces, therefore only relies on quadrature estimates for initial conditions and right hand sides. A typical error estimate for the solution U of an advection equation reads

(U -U h )(•,t) W -m,p ≤ C 1 exp (C 2 T ) h m
for t ≤ T , where h is the initial inter-particle spacing, m is the order of the quadrature rule using particle initial locations as quadrature poins, C 1 and C 2 are positive constants depending on the flow regularity, and W -m,p is the dual of W m,q , with 1/p + 1/q = 1). When one wishes to recover smooth quantities U ε h from the Dirac masses carried by the particles, one needs to pay for the regularization involved in the process and a typical error estimate becomes :

(U -U ε h )(•,t) L p ≤ C( ε r + h m /ε m ),
where r is the approximation order of the regularization used to mollify the particles (the reference [START_REF] Cottet | A new approach for the analysis of vortex methods in 2 and 3 dimensions[END_REF] provide detailed proofs of the above estimates in the context of vortex methods). This estimate immediately shows the dilemma of particle methods : the regularization size must be small, and controls the overall accuracy of the method, but it must contain enough particles so that the term h m /ε m does not compromise the convergence of the method. This constraint is even more stringent if particles are involved in additional terms of the model, in particular pressure gradient terms (in compressible flows) or diffusion terms. In the later case error terms of the form h m /ε m+2 arise. In any case, proper convergence would require, on top of ε → 0, that h/ε → 0, the so-called overlapping condition. This condition is in practice difficult to satisfy, in particular for 3D flows. Instead, a constant ratio h/ε is most often chosen in refinement studies.

On the other hand, even if the overlapping condition is not satisfied, particle methods still enjoy conservation properties and some kind of adaptivity which goes with the belief that "particles go where they are needed". This belief is however often more a hand-waving argument than a reasonable assumption based on solid grounds. Figure 1 shows very simple examples which illustrate the shortcomings of particle methods in the simulation of 2D inviscid vortex flows. In this case the 2D Euler incompressible Euler equations in vorticity form

∂ ω ∂t + div (uω) = 0
are discretized by particles of vorticity. In the left picture, a typical particle distribution is shown for an initial vorticity field with support in an ellipse. This figure shows that particles tend to align along directions related to the flow strain, creating gaps in the vorticity support. The right picture corresponds to an axisymmetric initial vorticity field, leading to a stationary solution. In this example ω 0 (x) = (1 -|x| 2 ) 3 . Error curves, in the energy norm for the particle velocities, corresponding to h = h 0 = 0.1 and h = h 0 /2 are plotted, with a constant ratio h/ε. This figure shows that, despite the smoothness of the solution, the expected initial gain in accuracy is almost completely lost after a short time due to the distortion of the particle distribution. Whenever point-wise values are required with some accuracy (for instance to recover satisfactory spectra in turbulent flows or local pressure or vorticity values on an obstacle) the overlapping condition cannot be ignored.

Remeshing and semi-lagrangian particles

Although several methods have been considered to overcome the lack of overlapping of particles while keeping their grid-free nature, particle remeshing is to our knowledge the only tool which so far allowed to deliver in a clear-cut way accurate results for complex two and three-dimensional dynamics, in particular in incompressible flows. Remeshing was already used in early simulations for some specific flow topology, like vortex sheets [START_REF] Krasny | Desingularization of periodic vortex sheet roll-up[END_REF] or filaments [START_REF] Martin | Numerical investigation of three-dimensional evolving jets subject to axisymmetric and azimuthal perturbation[END_REF], but its first systematic use goes probably back to [START_REF] Koumoutsakos | High Resolution simulations of the flow around an impulsively started cylinder using vortex methods[END_REF][START_REF] Koumoutsakos | Inviscid Axisymmetrization of an Elliptical Vortex[END_REF], where pioneering results where obtained for flow past a cylinder at challenging Reynolds numbers and for the problem of axisymmetrization of elliptical vortices. The first numerical comparisons of these methods with spectral methods in 3D turbulent flows were performed in [START_REF] Cottet | A comparison of spectral and vortex methods in three-dimensional incompressible flows[END_REF].

Remeshing consists of redistributing particle masses on nearby grid points, in a way that conserves as many moments of the particles as possible. The number of moments, and hence the accuracy, dictates the size of the remeshing kernel. Con-serving the 3 first moments (including mass) in some sense guarantees that remeshing has not a diffusive net effect. This has been considered as a minimal requirements in the references already cited and in all following works.

The remeshing frequency can be a point of debate. However there are two aspects to consider. On the one hand remeshing a particle distribution which has already been highly distorted is likely to produce numerical noise. On the other hand, the time scale on which particles are distorted is the same as the one on which the particle advection should be discretized, namely 1/|∇a| ∞ , where a is the advecting velocity field. It is therefore natural to remesh particles every few time-steps and numerical truncation errors coming from remeshing must be accounted for in the numerical analysis. Recently [START_REF] Cottet | High order Semi-Lagrangian particles for transport equations: numerical analysis and implementation issues[END_REF] particle methods have been analyzed from this point of view -in other words as semi-lagrangian methods.

To describe the method and discuss its accuracy, let us consider the 1D model linear advection problem -which is somehow the engine of particle methods in all applications :

θ t + (aθ ) x = 0, x ∈ R,t > 0, ( 1
)
where a is a given smooth velocity field. A particle method where particles are remeshed at each time step can be recast as

θ n+1 i = ∑ θ n j Γ x n+1 j -x i ∆ x , i ∈ Z d , n ≥ 0. (2)
In the above equation ∆ x is the grid size on which particles are remeshed (assuming a regular grid), x j are the grid points and Γ is the remeshing interpolating kernel.

x n+1 j is the result of the advection at time t n+1 of the particle located at x j at time t n . Note that to generalize the method to several dimensions one may use similar formulas with remeshing kernels obtained by tensor products of 1D kernels (this is the traditional way) or, following [START_REF] Magni | Accurate, non-oscillatory remeshing schemes for particle methods[END_REF], one can alternate the advection steps in successive directions, with classical recipes to increase the accuracy of the splitting involved in this process. This later method is economical when one uses high order kernels (with large supports) as, in 3D, its computational cost scales like O(3M) instead of O(M 3 ) for a kernel involving M points in each direction.

The moment conservation properties mentioned earlier to be satisfied by the remeshing kernel Γ can be expressed as

∑ k∈Z (x -k) α Γ (x -k) = 1 if α = 0 0 if 1 ≤ α ≤ p , x ∈ R, (3)
for a given value of p ≥ 1. An additional requirement is that Γ is globally in W r+1,∞ and of class C ∞ in each integer interval (in practice Γ is a polynomial in these intervals), and satisfies the interpolation property : Γ (i -j) = δ i j . In the simple case of an Euler explicit scheme to advect particles, x n+1 j = x j + a(x j ,t n ) and when the time step sastisifes the condition

∆t < |a | -1 L ∞ , (4)
on can prove [START_REF] Cottet | High order Semi-Lagrangian particles for transport equations: numerical analysis and implementation issues[END_REF] that the consistency error of the semi-lagrangian method is bounded by O(∆t + ∆ x β ) where β = min (p, r). Using higher order Runge-Kutta schemes increase the time accuracy, as expected. Moreover, at least for kernels of order up to 4, under appropriate decay properties for the kernel Γ one can prove for a large class of kernels the stability of the method under the sole assumption (4).

Let us give a sketch of the consistency proof for the case r = p = 1 if a is only a function of x and the Euler scheme is used to advance particles. We start from [START_REF] Berger | Adaptive mesh refinement for hyperbolic partial differential equations[END_REF] and assume that θ n j = θ (x j ,t n ) where θ is the exact smooth solution to the advection equation and we want to prove that

θ n+1 i = θ (x i ,t n+1 ) + O(∆t 2 + ∆ x 2 ).
We write j = i + k and, in (2) ,

θ n j = θ n i + k∆ x θ x (x i ,t n ) + O(∆ x 2 ).
Particle advection with the Euler scheme gives

x n+1 j = x j + a j ∆t = x i + k∆ x + a i ∆t + [a(x i + k∆ x) -a(x i )]
and thus

Γ x n+1 j -x i ∆ x = Γ (k + λ i + ν[a(x i + k∆ x) -a(x i )]) = Γ (k + λ i ) + k∆t a (x i )Γ (k + λ i ) + O(∆ x 2 ),
where we have used the notations

ν = ∆t/∆ x , λ i = a i ν.
We thus obtain

θ n+1 i = ∑ k θ (x i ,t n ) + k∆ x θ x (x i ,t n ) + O(∆ x 2 ) Γ (k + λ i ) + k∆t a (x i )Γ (k + λ i ) + O(∆ x 2 ) .
The moment properties of order 0 and 1 yield

∑ k k Γ (k + λ i ) = -λ i , ∑ k k Γ (k + λ i ) = -1.
We therefore have

θ n+1 i = θ (x i ,t n ) -∆t a(x i ) θ x (x i ,t n ) -∆t a (x i ) θ (x i ,t n ) + O(∆t 2 + ∆ x 2 ) = θ (x i ,t n ) -∆t(a θ ) x (x i ,t n ) + O(∆t 2 + ∆ x 2 )
and, by (1),

θ n+1 i = θ (x i ,t n+1 ) + O(∆t 2 + ∆ x 2 )
which proves our claim.

Fig. 2 Refinement study for the flow [START_REF] Cohen | Biorthogonal bases of compactly supported wavelets[END_REF]. CFL value is equal to 12. Black-circle curve : kernel Λ 2,1 ; red-square : kernel Λ 4,2 ; blue-triangle : kernel Λ 6,4 ; dashed lines indicate slopes corresponding to second and fourth order convergence.

In the general consistency result mentioned above, one can also check that the order of spatial accuracy is p whenever one can ensure that after the advection step each grid cell contains exactly one particle. [START_REF] Cottet | High order Semi-Lagrangian particles for transport equations: numerical analysis and implementation issues[END_REF] contains explicit expressions of kernels of order up to 6, denoted by Λ p,r where p and r measure the moment and regularity properties of the kernel. It also contains a number of refinement studies which suggests that in practice for a kernel Λ p,r the observed order of accuracy is between min (p, r) and p. Figure 2 shows the results of a typical refinement study on a 2D level set benchmark. This case consists of a level set function corresponding to a disk of radius 0.15 centered at (0.5, 0.15) in the periodic box [0, 1] 2 . The velocity field is given by

a(x 1 , x 2 ,t) = f (t) -sin 2 (πx 1 ) sin(2πx 2 ), sin(2πx 1 ) sin 2 (πx 2 ) , (5)
with f (t) = cos(πt/12). This field produces a strong filamentation of the solution culminating at t = 6 then drives the solution back to its initial state at t = 12, where numerical errors can be recorded. The actual convergence rate observed in this example for the kernels Λ 2,1 , Λ 4,2 and Λ 6,4 are respectively 1.87, 3.17 and 5.92.

It is interesting to note that several authors have recently advocated the use of particle methods to correct dissipative effects of finite-difference of finite-volume level set methods. Roughly speaking the idea is to seed particles at sub-grid levels near the interface and use these particles to rectify the location of the interface (see [START_REF] Lossaso | Two-way coupled SPH and particle level set fluid dynamics[END_REF] for instance). However [START_REF] Cottet | High order Semi-Lagrangian particles for transport equations: numerical analysis and implementation issues[END_REF] shows that, both in 2D and 3D, plain semilagrangian particle methods, with appropriate remeshing kernels (second order is actually enough) deliver better results with fewer points and larger time-steps.

The possibility of combining high accuracy with stability non constrained by the grid size makes semi-lagrangian particle methods appealing tools for turbulent transport. In [START_REF] Lagaert | Hybrid spectral particle method for the turbulent transport of a passive scalar[END_REF] this was exploited to investigate universal scaling laws for passive scalars advected in turbulent flows. In this study the accuracy of particle methods was first compared to classical spectral methods. Figure 3 shows a typical comparison of scalar spectra and of the pdf of the scalar dissipation for a turbulent flow. In this experiment a second order kernel was used for the particle method. These results, and several other diagnostics, indicate that except for the very tail of the spectra most of the scales are well captured by the particle method with the same resolution as for the spectral method. A factor 1.2 between the grids was found sufficient to resolve satisfactorily the scalar also in the dissipative range. In the case of high Schmidt numbers (ratio between flow viscosity and scalar diffusivity,) even with this requirement for slightly increased resolution, the gain in CPU time over the spectral method resulting from the use of large time steps in the particle method reached a factor 80.

Adaptive semi-lagrangian particles

Remeshing particles somehow detracts particle methods from self adaptivity (however illusive that notion might be, as we have seen). To restore some kind of adaptivity in the method it is natural to rely on the grid on which particles are remeshed. Like for grid-based methods, one may envision several ways of doing so. One way is to assume a priori that grid refinement is desirable in some parts of the flow, typically zones which are close to fixed boundaries. Another way is to adapt the grid to the smoothness of the solution itself.

Semi-lagrangian particle methods on non-uniform grids

In that case we assume that a non-uniform grid (the physical space) is obtained by a predefined mapping from a cartesian uniform grid (the reference space). The method works as follows [START_REF] Cottet | Vortex methods with spatially varying cores[END_REF]. At each time step, particles are pushed in the physical space, then mapped in the reference space. In this space regular remeshing formulas are used on cartesian grids, and particle locations are meshed back to the physical space. Figure 4 is an illustration of this method in the context of 2D vortex particle methods for the simulation of the rebound of a vortex dipole impinging on a wall. In this example the grid is stretched in both directions in an exponential manner with respect to the distance to the wall. Applications in 3D flows of similar technics can for instance be found in [START_REF] Ploumhans | Vortex Methods for Direct Numerical Simulation of Three-Dimensional Bluff Body Flows: Application to the Sphere at Re = 300, 500, and 1000[END_REF]. Fig. 3 Comparison of spectral and particle methods in a turbulent flow with a Schmidt number equal to 50 [START_REF] Lagaert | Hybrid spectral particle method for the turbulent transport of a passive scalar[END_REF]. Top picture : spectra of the scalar variance. Bottom picture: pdf of the scalar dissipation. Red curves : spectral method with 1024 3 points; green curves : particle method with 1024 3 points ; blue curves : particle method with 1280 3 points.

Particle methods with adaptive mesh refinement

Particle remeshing also enables to incorporate Adaptive Mesh Refinement (AMR) finite difference techniques [START_REF] Berger | Adaptive mesh refinement for hyperbolic partial differential equations[END_REF] to adapt the particle discretization to the solution itself. [START_REF] Bergdorf | Multilevel Adaptive Particle Methods for Convection-Diffusion Equations[END_REF] describes how to define, move and remesh patches of particles at different resolution in a consistent way. Figure 5 illustrates how the method allows to capture filaments ejected by an elliptical vortex in a 2D inviscid flow. In this example the refinement was based on the vorticity gradients. We do not enter more into the details of this method as we believe that it is outperformed by the more recent wavelet-based method described below.

Wavelet-based multi-level particle methods

The concept of multi-resolution semi-lagrangian particle methods was recently pushed further in [START_REF] Bergdorf | A Lagrangian Particle-Wavelet Method[END_REF], using wavelet tools. The idea of combining particle and wavelet methods can actually already been found in [START_REF] Basdevant | Methode des ondelettes mobiles[END_REF]. In this reference wavelet served as particle shapes in a grid-free method. The method however did not find practical ways to address the issue of interacting and recombining wavelets. In [START_REF] Bergdorf | A Lagrangian Particle-Wavelet Method[END_REF] instead, because particles are remeshed at every time-step on regular grids, particle methods can inherit concepts and techniques used in the context of finite-difference methods (see for instance [START_REF] Vasilyev | Solving Multi-dimensional Evolution Problems with Localized Structures using Second Generation Wavelets[END_REF]). Nonetheless the semi-lagrangian character of the method introduces original and interesting features. To describe the method and understand how it works, let us consider again the 1D advection equation on the real line, first with constant velocity then in the general case. The multi-dimensional case and the application to incompressible flows will be next outlined. We consider equation ( 1). The method is based on the following classical wavelet decomposition, in the framework of interpolating bi-orthogonal wavelets [START_REF] Cohen | Biorthogonal bases of compactly supported wavelets[END_REF] :

θ (x) = ∑ j c 0 k φ 0 j (x) + ∑ j L-1 ∑ l=l 0 d l j ψ l k (x) (6)
where l 0 (resp. l = L) corresponds to the coarsest (resp. finest) level. In the above Fig. 6 Advection of top-hat profile in a uniform field by a 2-level particle method. Left picture : black curve : result for uniform coarse resolution; blue curve : result at the coarsest level for a 2level method; red squares active particle at the fine level. Right picture : number of active particles as a function of time.

equation φ l j and ψ l j are respectively the scaling and wavelet (or detail) basis functions, centered around the grid point x l j = j∆ x l , where ∆ x l = 2 L 0 -l ∆ x 0 is the grid size at level l. This decomposition allows to define the solution at the different scales :

θ l (x) = ∑ k c l j φ l j (x). (7)
Scale and detail coefficients are related between successive levels through classical filtering operations :

c l i = ∑ j h(2i-j)c l+1 j , d l i = ∑ j g(2i-j)c l+1 j , c l+1 i = ∑ j [h(2 j -i)c l j +g(2 j -i)d l j ] (8)
where the filter functions g, h, g , h depend on the particular wavelet system chosen.

Let us now describe one time-step of the algorithm, first for a constant velocity value. The methods advances the solution scale by scale in the following manner. Assume at time t n = n∆t the solution is known on grid points belonging to the nested grids (x l j ) j,l . For each level l ∈ [L 0 , L]

• a wavelet analysis selects grid points which correspond to detail coefficients d l j above a given threshold • particles are initialized on these grid points with, for grid values, the corresponding scale parameters c l j

• particles are pushed and remeshed on the grid (x l k ) k . For the remeshing step above to be consistent one needs to ensure that active particles are always surrounded by "enough" particles carrying consistent values of θ l . This is done in a way similar to what would be done in a finite-difference method by creating ghost particles in the neighborhood of the active particles. The values of the solution assigned at time t n to these grid points is obtained by interpolation from θ l-1 . In order to make sure that only relevant values at level l are retained at the end of the iteration, a tag value equal to 0 is assigned to the ghost particles while the value 1 is assigned to the active particles. These tag values are pushed and remeshed along with the particles. After remeshing, only particles with tag values different from 0 are retained.

Finally, when all active particles have been moved and remeshed at all scales, values of the function at a given scale are updated by values at the next scale on even points when available.

Figure 6 illustrates the method in the case of a translating top-hat profile. The remeshing kernel is Λ 2,1 which for a constant velocity corresponds to a second order method. The scaling function is a piecewise linear function and the detail coefficients correspond to central finite-differences of the second derivatives. In the left picture we compare, after the time needed to travel 15 times the width of the top-hat, the results obtained with uniform one-level resolution at ∆ x 0 = 0.01 and ∆ x 1 = 0.005 and the wavelet-based method using these two levels. For the twolevel case, only the coarse resolution is showed together with active particles at the higher resolution. One can see that all solutions exhibit overshoots, as expected from a second order method, but the two-level solution limits the overshoot even at the coarse level and is very close to the higher resolution solution. The number of active particles, shown on the right picture of the figure, increases to respond to the oscillations created by the remeshing, then stabilizes.

Let us now consider the case of a non uniform velocity. In this case, scales are not advected in an independent fashion but interact as a result of compression and dilatation. To allow fine scales to appear from coarse scales, there are two options. The first one is similar to what would be done in a finite difference method [START_REF] Vasilyev | Solving Multi-dimensional Evolution Problems with Localized Structures using Second Generation Wavelets[END_REF]. It consists of considering for each scale at the beginning of each step additional ghost particles at the level immediately above on which values of solutions are interpolated from the coarse scale. Another option, simpler and more specific to particle methods, is to remesh each scale on a scale twice smaller. At the end of the timestep, the value of the solution on this smaller scale is chosen to be either the result of the remeshing either of the coarse grid or the fine grid, depending on the value of the flag described earlier.

It is important to note that, in both options, the time scale on which the scale l + 1 appears from scale l is governed by 1/|a | ∞ which is consistent with the maximum time step allowed for the particle method. In other words for time steps corresponding to large CFL numbers, only scale l + 1 can appear form scale l is the condition (4) is satisfied. Figure 7 illustrates the method for two scales for an initial condition consisting on a sine wave advected in a velocity given by a(x) = 1 + cos(2πx). As the wave travels to the right it is subject to successive compressions/amplification and dilatation/dampening. As a result, the number of active particles oscillates (right picture). The left picture shows that the adaptive method allows to recover the high resolution results. To go beyond the 1d toy problem just considered, there are again two options. One is to use the same ideas but to rely on multidimensional tensor product wavelets, with the added complexity that one grid point is associated to several wavelets. This is the option followed in [START_REF] Bergdorf | A Lagrangian Particle-Wavelet Method[END_REF]. The other one, presumably simpler and which only uses 1D wavelets, would be to split advection into three successive directional advections, along the lines of [START_REF] Magni | Accurate, non-oscillatory remeshing schemes for particle methods[END_REF].

In vortex flows, where these methods have been primarily applied, an additional work is to compute velocities created by multi-level vortices. The solution chosen in [START_REF] Bergdorf | A Lagrangian Particle-Wavelet Method[END_REF] is to use a fast multipole method, where each particle associated to the value of the solution and the associated grid size is accounted for in the Biot Savart law. The left picture of Figure 8 illustrates the method for the flow around swimming fishes [START_REF] Gazzola | Reinforcement learning and wavelet adapted vortex methods for simulations of self-propelled swimmers[END_REF]. The power of the method is here fully apparent. The complex geometries and the associated boundary conditions are dealt with by a penalization method an the nested multi-level cartesian grids allow to capture at a minimal cost the fine vortices created by the swimmers, despite the low accuracy of the penalization method. This method has been implemented for both 2D and 3D geometries [START_REF] Rossinelli | MRAG-I2D: Multi-resolution adapted grids for remeshed vortex methods on multicore architectures[END_REF]. It has been in particular applied in combination with optimization technics to determine efficient swimming strategies [START_REF] Gazzola | Reinforcement learning and wavelet adapted vortex methods for simulations of self-propelled swimmers[END_REF].

Conclusion

Particle methods with particle remeshing at each time-step can be analyzed as semilagrangian conservative methods. The accuracy of the method is governed by moment and regularity properties of the remeshing kernel and high order kernels can be derived in a systematic fashion. Adaptivity in the method can be reinforced in particular by using wavelet-based multi-resolution analysis. In any case, the semilagrangian nature of the method allows to use time-steps which are not constrained by the grid size. In many applications this can lead to substantial computational savings.

Fig. 1

 1 Fig.1Simulation with a grid-free vortex particle methods[START_REF] Cottet | Vortex methods[END_REF]. Left picture: particle distribution in the simulation of an elliptical vortex patch. Right picture: error curves for a stationary axisymmetric vortex. Solid line h = 0.005, dashed line h = 0.01. ε/h = 1.5.

Fig. 4

 4 Fig.4Vortex dipole impinging on a wall[START_REF] Cottet | Vortex methods with spatially varying cores[END_REF]. Top pictures : results at two successive times with remeshing on an exponentially stretched grid. Bottom pictures : results with uniform grids at the later time for the coarsest (left) and finest (right) grids. Only a small percentage of the active particles are shown by dots.

Fig. 5

 5 Fig.5Blocks of refined particles in an AMR implementation of particle methods for the simulation of an elliptical vortex[START_REF] Bergdorf | Multilevel Adaptive Particle Methods for Convection-Diffusion Equations[END_REF].

Fig. 7

 7 Fig. 7 Multilevel particle method in a compression/dilatation flow. See figure 6 for captions.

Fig. 8

 8 Fig.8Multi-level vorticity particle method with penalization around complex geometries . Top figure : 2D calculation around multiple fishes. Bottom picture : 3D vorticity passed a wind turbine (curtesy of P.Koumoutsakos)


IPDM doctoral workshop 11-12/06/2016 Laure-Anne Parpaleix, 1st year PhD candidate MINES ParisTech, PSL Research University, Center for Management Science (CGS) Innovation and growth potential: managing investment in middle market companies

Growing challenges and research gap

The search for levers to sustain economic growth, which is a key and topical challenge for public policy and private actors, has turned into an intensive research topic in both economy and management. Sustaining growth through investments remains at the heart of public policies insofar as growth is assumed to be correlated to economic drive notably through increased firms' profit and employment rate. Therefore, search for a better understanding of growth determinants has been extensively discussed in numerous theoretical and empirical studies [START_REF] Coad | Empirical investigations into the characteristics and determinants of the growth of firms[END_REF].

Throughout the continuous refinement in the economic models investigating the dynamics of investment for sustained growth, technological progress and, further, innovation have taken an increasing role. In Adam Smith first theory on macro-economic growth [START_REF] Smith | An Inquiry into the Nature and Causes of the Wealth of Nations[END_REF] technical progress emerged as a key factor driving macro-economic growth. Over time, several theories further developed causal links between technological progress and both macro and micro growth rate. First defined as an exogenous and undefined variable it has gradually been endogenised since neoclassical models, and refined as, among others, a productivity factor [START_REF] Solow | A Contribution to the Theory of Economic Growth[END_REF] or a stock of knowledge [START_REF] Arrow | The economic implications of learning by doing[END_REF]. Overall, innovation stands out as one of the main growth driver [START_REF] Ahlstrom | Innovation and Growth: How Business Contributes to Society[END_REF], a result that is thought to help steering the investment choices. Yet, the underlining processes linking investment to innovation on the one side, and innovation to growth on the other remain ill-understood.

On the one hand, despite voluminous and diverse literature [START_REF] Cameron | Innovation and Growth: a survey of the empirical evidence[END_REF] produced by the study of innovation impact on economic growth, a precise relationship has yet to be unequivocally established [START_REF] Demirel | Survey of the literature on innovation and economic performance[END_REF]. As for now, correlations based on static and in retrospect input-output economic models remains inconclusive (Gupta, Guha et al. 2013). In these studies, the several growth indicators (e.g. turnover, employment rate, productivity) follow independent paths and are not correlated to any innovation factor (e.g. R&D expenditure, patents, propensity to innovate) [START_REF] Del Monte | R&D and the growth of firms: empirical analysis of a panel of Italian firms[END_REF]. Thus characterizing the relationship between firm innovativeness and growth raises the stake of identifying accurate and appropriate measurements for both growth and innovative activities.

On the other hand, the classical thought that R&D spending is statistically linked to innovation, thus far roughly correlating the issue with the right amount of resources invested, has been repeatedly proven false over large sets of data. The absence of a mechanical impact of investment on innovative output, known as "the R&D paradox" (Le Masson, Weil et al. 2010), shows that beyond the initial financial inputs, appropriate design management and governance models seem necessary to support sustainable innovative activities.

A wide range of researches in the field of innovation management has built on this perspective to offer refined patterns of innovative activities process (e.g. spillovers theories, resource based view, dynamic capabilities (O'Connor 2008). However, their impact on firms' growth dynamics has not been pinpointed. Because of this failure in understanding growth drivers, investors lack of managing strategies allowing them to guide their investing policies. Thus, my researches aim at describing new growth dynamics in order to shape investment guiding tools.

My researches will be conducted in partnership with Bpifrance venture capital fund.

Bpifrance is a state-owned organism backing national public policies by assisting companies through different types of investments. Bpifrance doctrine states that it should "behave like a prudent investor operating under market conditions to serve public interest". Therefore, the venture capital fund has both to meet financial performances and to carry out a public service mission. This dual goal is assumed to be reached by ensuring firms forthcoming growth. Bpifrance share the view that innovation capabilities are an accurate indicator of firm's growth, as a breeding ground for both economic value and public benefit. Thus, Bpifrance investors challenge lays analyzing firms innovative capabilities to identify highgrowth potential.

Research goals

Based on recent management researches finding, my thesis aims at answering to the following questions raised by the above assessment:

 How to characterize innovation capabilities in order to assess firms' growth potential (i.e. indicators, measurement methodology, and dynamics)?

 How can investment reach greater impact on firm sustainable growth? How to shape investments tools (i.e. innovative firms' selection mechanisms, financial vehicles) and steering process accordingly?

As regards the empirical basis, my research will focus on middle market companies as there are assumed to be one of the main growth engines. This firm category is of particular interests as it is a quite new one so the relationship between innovation and growth has been little examined and stays even more mysterious. In 2008 a French law created this new business category to reach more precise analysis on these firms considered as valuable assets for French economy. It gathers companies employing between 250 and 5 000 workers with an annual turnover below 1.5 billion euro. In-between SMEs and large firms, French middle market companies' category, as German Mittelstand, gathers long lasting firms, half of them being family firms, with entrenched innovation capabilities that should drive value generation and employment opportunities. However, this single category hides a large diversity of innovation and growth patterns. These various potential growth dynamics could explain why middle market companies growth lack of steadiness and persistence. Growth process remains tricky to characterize and monitor with current statistic tools used for SMEs and large firms, hence the need to adjust current tools and investigate new indicators.

Theoretical framework

Despite abundant literature, the relationship between innovation and growth remains illunderstood.

Historically, economists were the first to pinpoint and explore technological progress as a driver of economic growth through macro-economic frameworks [START_REF] Smith | An Inquiry into the Nature and Causes of the Wealth of Nations[END_REF]. Over time, various theories refined causal links between technological progress and growth rate at different level of analysis ( [START_REF] Solow | A Contribution to the Theory of Economic Growth[END_REF], [START_REF] Arrow | The economic implications of learning by doing[END_REF], (Romer 1986)). Gradually, a shift happened from theories developed at an aggregate scale to a firm and an individual scale [START_REF] Nelson | An evolutionary theory of economic change[END_REF]. Various disciplines contribute to model economic growth including economy (macro, micro and political economy) [START_REF] Cameron | Innovation and Growth: a survey of the empirical evidence[END_REF], management [START_REF] Teece | Explicating dynamic capabilities: the nature and microfoundations of (sustainable) entreprise performance[END_REF], finance [START_REF] Timmons | Venture capital's role in financing innovation for economic growth[END_REF], law [START_REF] Levine | Law, Finance, and Economic Growth[END_REF]) and psychology (entrepreneurship literature). Nevertheless, the relationship between innovation and growth at the firm level remains an intellectual puzzle.

I conducted an analysis of 29 articles selected among the most relevant papers (i.e. taken among the most cited and published in highly ranked papers) mainly extracted from three literature review ( [START_REF] Ernst | Patent applications and subsequent changes of performance: evidence from timeseries cross-section analyses on the firm level[END_REF]  Despite some studies highlighting a positive correlation between innovative activities and growth, its conclusions often cannot be duplicated. Besides, others (e.g. (Brouwer, Kleinknecht et al. 1993)) demonstrate that innovation and growth remains unrelated (Bottazzi, Dosi et al. 2001).

 Those economic studies mainly analyze innovative outcomes in retrospect (counting patents or new products and processes) more than innovative capabilities hence a need for further investigation.

It reflects both a technical issue (i.e. time period, firm size, type of industry) and a debate about the nature of targeted growth (i.e. sales, employment, market share, goodwill). Establishing a correlation between innovation capabilities and growth potential raises stakes in identifying accurate and appropriate measurements for both growth potential and innovative capabilities beforehand. Management literature should contribute to overcome these issues.

Management literature on design functions and ecosystems enables to refine the concept of growth

Theories on public action often assume a correlation between R&D investments, innovation outcomes and growth. However, current economic studies, aiming at defining in retrospect a correlation between inputs and outputs, are not sufficient to demonstrate such a relationship. The difficulty to demonstrate the impact of R&D investment on companies' innovation and growth is known as "the R&D paradox" (Le Masson, Weil et al. 2010). An analysis of top 700 international companies' investments in industry demonstrated no correlation between R&D expenditures and firms growth rates of turnover. Other studies confirm those results whatever timescales, performance indicators or sectors used [START_REF] Jaruzelski | Money Isn't Everything[END_REF]. For example, among the most innovative companies, some, like Apple in 2010, have reached high-growth rates spending less in R&D than the average firm in their field. Therefore, taking into account the intensity of R&D investment only seems too restrictive. Beyond the initial financial inputs, appropriate design management and governance models seem necessary to support firms' innovation capabilities and assess innovation potential beforehand. Recent researches on innovative design regimes [START_REF] Hatchuel | A new approach of innovative design: an introduction to C-K theory[END_REF] should help reviewing traditional economic indicators, develop new ones tailored for measuring beforehand firms innovation dynamics and thus guide future growth and investments strategies.

Based on these works I build two research hypotheses:

 Micro-economic growth can be fruitfully described as the increase of innovation capabilities. Literature on design functions [START_REF] Hatchuel | Growth of the firm by repeated innovation: towards a new micro-economics based on design functions[END_REF] shows that production function [START_REF] Ricardo | On the Principles of Political Economy and Taxation[END_REF][START_REF] Ricardo | On the Principles of Political Economy and Taxation[END_REF]) is insufficient to account for the development of firms, especially innovative firms. They show that considering a single financial output (such as turnover) does not take into the account the capacity to create new products, concepts and knowledge/skills (Hatchuel and Masson 2001). Instead, most of the recent developments in innovation management show that innovative capabilities building, absorptive capacity, etc. are pivotal elements to account for firms' growth potential [START_REF] Teece | Explicating dynamic capabilities: the nature and microfoundations of (sustainable) entreprise performance[END_REF]. We hypothesize that growth regimes can be differentiated based on whether companies are able to expand their knowledge, goods and ideas or not  Growth can be better observed at an ecosystem level: as literature on knowledge spillover [START_REF] Audretsch | R&D Spillovers and the Geography of Innovation and Production[END_REF] has shown effects beyond a single company, the challenge lies in identifying appropriate patterns to describe partnerships and all kinds of social returns [START_REF] Griliches | The search for R&D spillovers[END_REF]) (Griffith).

To summarize, my researches aim at refining the concept of growth:

 By investigating growth process as the expansion of design capabilities instead of the usual economic growth patterns. The empirical case studies should enable identifying several growth path linked to several investment patterns, contribution to economic value and public interest.

 By providing criteria that would allow taking into account firms ecosystems dynamics. For instance through the analysis of the diversity of partnerships, platform-based designs diffusion, spillover impacts …

New approaches for research on relationships between investment and innovation

Literature discusses the coupling between investment and innovation by questioning how a firm capital structure influences its innovative projects performance. Several studies advocate that venture capital funded R&D generates substantial increase in innovative outputs compared to corporate R&D without debating on causes of this productivity differential [START_REF] Kortum | Assessing the contribution of venture captial to innovation[END_REF]. Compared to corporate R&D or relationship lending, venture capital main asset relies in its ability to combine an initial financial input and support throughout the project lifetime. Capital venture, by interfering in firms at critical stages of development affects their design management and governance models. Then, different hypothesis could be investigated to explain R&D investment performance differential, for example whether it is due to strict project selection mechanisms ex-ante, to closer ex-post monitoring (Gerasymenko, Clercq et al. 2015) or to specific innovation management methods [START_REF] Engel | Accelerating Corporate Innovation: Lessons from the Venture Capital Model[END_REF]. Effects of growth support policies through investments can now be reviewed in the light of recent literature on design theories and innovation capabilities. The issue lies in identifying the right investment models according to the aimed growth dynamic.

Research approach

My research is following multiple paths:

 Contextualization through historical cases of public policies strategies that have supported firm growth. The aim is to size a public operator's specific assets for sustaining firms' growth compared to external private partners or internal corporate capabilities.

 Identify different growth dynamics according to the several innovation capabilities paths.

o Characterization of innovation capabilities through multiple bodies of literatures. Both economic literature and management literature on dynamic capabilities theories (O'Connor 2008) argue on innovative capabilities dynamics.

o Thanks to Bpifrance records, cases study of financed firms following several growth models will contribute, in the light of design theories (e.g. Concept-Knowledge theory [START_REF] Hatchuel | Innovation répétée et croissance de la firme: mico-économie et fonctions de conception[END_REF])), to model innovation capabilities and suggest news indicators.

Qualitative studies enable the focus on surprising cases that are essential to single parameters out before proceeding to quantitative tests on databases.

 Understanding of the impact of investments on innovation and growth to adjust investments tools and steering processes through finance, innovation management and public action policies literatures. At the same time, I will conduct empirical analysis of investors' reasoning (notably when assessing firms' innovation capabilities, forecasting potential growth or shaping strategic guidance).

 In parallel of the above mentioned qualitative analysis, quantitative tests based on multiple databases provided by Bpifrance, completed if necessary by external sets of data (INSEE, CASD Community innovation survey and social data statements) will contribute to refine hypothesis on the coupling between investment, innovation and growth.

Expected outcomes and implications for future research and practitioners

On one hand, I expect the following theoretical implications:

 An in-depth study of design capabilities expansion determinants leading to a better understanding of the determinants of growth and the role of innovation capabilities in firms' growth trajectories.

 Original model with parameters linked to the description of innovation capabilities.

 Clarification of diverging firms growth objectives which are not always in line with public interest and its potential impact on innovation and investments strategy.

On the other hand, I aim at managerial implications regarding tools to steer the investment strategies (i.e. innovative firms' selection mechanisms, financial vehicles). Besides, I also expect findings regarding public policies as this research will question the hypothesis that the best strategy to serve public interest, as a state-owned venture capital fund, is maximizing a few individual firms' growth.

Central issues to be advised on during the doctoral workshop

The theoretical background of my research project builds on several bodies of literature such as: innovation management, economic growth theories, finance. Therefore, one challenge is to identify which research communities to discuss my work with. The doctoral workshop would be a great opportunity to discuss this issue.

I am also interested in feedbacks on the methodological approach. What are the potential drawbacks of analysing a few empirical cases before going on with quantitative tests?

Finally, I would be pleased to receive any comments on the following topics at the heart of my thesis:

-Insights on literature bodies dealing with the indicators measuring the impact of externalities special features of middle market companies' (including family firms) innovation capabilities management and economic growth dynamics

  ,[START_REF] Del Monte | R&D and the growth of firms: empirical analysis of a panel of Italian firms[END_REF],[START_REF] Demirel | Survey of the literature on innovation and economic performance[END_REF]). Most of these quantitative studies are based on correlation tests either to investigate if growth follows a random walk (literature on Gibrat's Law[START_REF] Mansfield | Entry, Gibrat's Law, innovation and the growth of firm[END_REF])) or to evaluate the impact of endogenous (i.e. patent granted, R&D expenditure, new product) and exogenous indicators (i.e. spillovers, clusters) on past firms' growth. While theoretically, economic approach assumes that innovation fosters firm performance, quantitative researches has failed to generate cumulative results:  Heterogeneity of growth measure is a first obstacle[START_REF] Delmar | Measuring Growth: Methodological Considerations and Empirical Results[END_REF], Delmar, Davidsson et al. 2003) as the several growth indicators (e.g. turnover, employment rate, firm survival, productivity, market value) are not correlated (i.e. an increased in turnover can go along either with a fall or a rise in employment rate). As growth is heterogeneous in nature[START_REF] Delmar | A taxonomy of high-growth firms[END_REF], a definition of growth including several dimensions (i.e. instead of one single indicator) is needed to distinguish several kinds of growth. Similarly, measuring innovation is one of central methodological issues. The question is: how to find accurate measures of innovation and examine its impact on firm's performance? Innovation is often measured in retrospect through R&D spending or numbers of patents. Both are quantitative indicators failing in differentiating how disruptive innovations are.