G Cherchi

P Alliez

R Scateni

M Lyon

D Bommes

Selective Padding for Polycube-based Hexahedral Meshing

Keywords: Descriptors (according to ACM CCS): I.3.5 [Computer Graphics]: Computational Geometry and Object Modeling-Curve, surface, solid, and object representations

Figure 1: Overview. Our pipeline takes as input a model and its polycube mapping (a); we compute the relative hex-mesh and locate the surface areas in need of padding analyzing the mapping quality (b); we set and solve a binary problem to find a set of facets to extrude in order to create a selective padding layer (c); we compute and analyze the mapping with the new hex-mesh structure (d).

Introduction

Volumetric meshes are an ideal tool for physically-based simulations, especially in the fields of mechanical simulation and computational fluid dynamics. While tetrahedral meshes are very simple to obtain and flexible in adapting to the computational domain, hexahedral meshes have better numerical properties and lower complexity. One straightforward way to generate a hexahedral mesh on a given domain passes through polycube mapping.

Polycubes are three-dimensional domains formed by connecting axis-aligned cuboids face to face. The simplicity of their structure ex-plains their popularity for several computer graphics problems such as quadrilateral and hexahedral meshing (e.g., [HXH10, GSZ11]), efficient texture mapping (e.g., [THCM04, CL ⇤ 10]), morphing (e.g., [FJFS05]), spline fitting (e.g., [WHL ⇤ 07]), etc.

When the volumetric domains to partition in hexahedra are regular enough (e.g., with limited range of levels of details like mechanical and CAD models), polycubes are an ideal tool for generating hexahedral meshes, by computing volumetric mappings with their corresponding domains. However, the quality of the hexahedra directly depends on the mapping distortion. While the distortion is typically negligible inside the domain (because of the regularity of the lattice we use to grid the polycube), it can be excessive near the boundary. For instance, when a convex boundary edge of the polycube with a folding angle of 90 degrees, maps to a flat part of the domain boundary, with 180 degrees.

A common solution consists in "padding" the entire domain by adding extra layers of hexahedra, in a way that pushes inside all 90-degree edges and replaces them with edges incident to two hexahedra. Padding offers a means to trade deformation for number of elements, and increases the overall quality of the domain.

In this paper, we present an automated method able to insert padding elements only where they are needed to increase the overall quality. Our algorithm works selectively on the volumetric domain, keeping untouched the parts that are directly derived from the polycube mapping having already a good quality. It reaches this goal operating selective insertion of sheets of hexahedra inside the domain. We review next the previous work relating to padding strategies, and we position our solution in terms of complexity and quality of the final hexahedral mesh.

State of the Art

Polycubes and Parameterization-based hex-meshing. Introduced in 2004 [THCM04], polycubes attracted a great attention from the Computer Graphics community, for several applications. An important field in which polycubes have immediately spread is the one of volumetric applications and hexahedral meshing, in particular for domains with limited range of levels of details [LVS ⇤ 13, HJS ⇤ 14, YZWL14, GSZ11, FXBH16]. Shortly after their formalization, it was clear that the use of polycubes for purposes different from texture mapping needed an algorithmic strategy for their generation [START_REF] Fan | Automatic PolyCube-Maps[END_REF]HWFQ09,[START_REF] Gregson | All-Hex Mesh Generation via Volumetric PolyCube Deformation[END_REF]. Recent approaches improve the polycube quality by reducing the mapping distortion, lowering the number of corner singularities [LVS ⇤ 13, HJS ⇤ 14] and preserving axis-rotation invariants [FBL16]. However, most of these approaches do not address the quality of polycube-generated hexmeshes.

A wide range of approaches aim at generating well-structured volumes, based on skeletons (e.g., [START_REF] Livesu | Skeletondriven Adaptive Hexahedral Meshing of Tubular Shapes[END_REF]), gridding (e.g., [LJLJ15]) or expanding methods (e.g., [TBM96]). Parameterization-based methods are another important class of hex-meshing methods. They map the input volume to another parametric space where the final mesh connectivity is generated. With the polycubes spread, new approaches were proposed for meshing [GSZ11, LVS ⇤ 13] and optimization [CLS16, GDC15]. The added value of these approaches is clear: reasoning with a simple structure (together with a mapping function with the original shape) yields several advantages in terms of efficiency and simplicity of implementation, as a polycube can be trivially hex-meshed with a regular grid. Nevertheless, most algorithms devised for polycube generation overlook the quality of the elements of the final mesh.

The structure of the final mesh is defined by the shape of the polycube and the lattice in which it has been gridded (the corners of the polycube become singularities in the mesh). Such a structure is regular with high quality elements inside the domain, but it may exhibit very low quality elements near the domain boundary. For this reason, polycube-based hex-meshes often require a post-processing step referred to as padding (or pillowing [START_REF] Mitchell | Pillowing doublets: refining a mesh to ensure that faces share at most one edge[END_REF]) on the entire domain boundary. As detailed by Shepherd [She07,[START_REF] Shepherd | Hexahedral mesh generation constraints[END_REF], the padding operation starts with an initial mesh from which a subset of hexahedra is defined to create a shrink set. The shrink set is separated from the original mesh and shrunk. The void left by such a shrinking process is filled by adding a new layer of hexahedra.

The motivation for our approach stems from the observation that, for specific shapes, the padding operation is not necessary for the whole boundary. In addition, it may in some cases worsen the mesh quality. We explore the possibility to perform such a padding operation via local sheet insertion.

Hex-mesh refinement A wide range of local refinement algorithms have been proposed both for quadrilateral and hexahedral meshes [SDW ⇤ 10]. A common objective is to change the mesh resolution, to decrease the valence of inner vertices or to adapt the mesh density in specific areas as required by FEM simulations.

Zhu et al. [ZCWG14]

proposed a method to improve the quality of CAD-based hex-meshes. While the user deforms the CAD model, the associated hex-mesh is automatically improved by adding or removing hexahedral sheets by using dual operations, in order to keep the resolution and the quality of the mesh constant. In our setting the regularity of the mesh is guaranteed by the polycube properties. We can hence limit our topology changes to the mesh boundary.

Chen et al. [START_REF] Gao | An approach to achieving optimized complex sheet inflation under constraints[END_REF] introduced an approach to achieve complex sheet inflation under various constraints, in order to improve the mesh quality. The method takes as input a set of user-defined boundary mesh edges and a set of hexahedra. The edges specify the boundary position where the new sheets should be inserted, and then the algorithm computes, through an iterative solving of a Max-Flow and optimization steps, the whole layer position. In our meshes, even if we allow the sheet insertion in all the mesh, we restrict the quality analysis to the mesh surface, because polycube-based hex-meshes always have regular and good-quality inner elements. For this reason, it is easier, in our case, to automatically detect the position where the insertion of a new hexahedral sheet can improve the mesh structure. 17] proposed an automated block decomposition method based on sheet operations, which generates a block decomposition from which a high-quality hex-mesh stems. They start from a B-rep solid model, compute the relative tet-mesh and finally extract a hex-mesh. They insert and collapse whole sheets of hexahedra to improve the obtained hex-mesh by solving an integral linear problem. We instead rely on the good polycube-based inner structure and insert sheets of hexahedra to improve the quality of the near-surface elements in the final hex-mesh. It is interesting enough that they agree on the importance to develop an algorithm for a robust insertion of boundary hexahedral layers, which is what we propose in this paper.

Wang et al. [WSC

Owen et al.

[OSE17] contributed a template-based approach for generating locally refined all-hex meshes. Using a restitched set of split configurations, a local refinement of the hex-mesh structure is performed, yielding elements with minimum Scaled Jacobian of 0.3. Our approach hinges upon a similar set of templates. Since we apply the sheet insertions in polycube based hex-meshes, we can afford to use a restricted set of dual operators to achieve the desired result. Indeed, since we only have edges belonging to 1 to 4 facets in polycube-based hex-meshes, we can obtain our goal just by using the set of operators that increase by one the valence of the hex-mesh edges, where it is needed.

Wang et al. [START_REF] Wang R | Hex mesh topological improvement based on frame field and sheet adjustment[END_REF] presented a method to improve the topology of hex-meshes via frame field optimization and sheet operations. Starting from a hex-mesh in which they build an initial frame field, they optimize the field in order to obtain a high-quality one, that can be used to identify the most problematic areas in the mesh. Then, they adjust the structure of the mesh via a set of sheet operations. Even if we work with different topologies and perform other types of analysis, Section 6 compares our approach with this work as we are adopting different methods to achieve similar goals. Indeed, working with polycube-based hex-meshes, we can afford to focus our analysis on the mesh boundary, while they need to analyze the whole structure to perform the optimization.

Positioning and Contributions

In this work we present an algorithm to perform selective and localized padding into polycube-based hex-meshes. The introduction of distortion when mapping an object with curved surfaces or nonright angles into an axis-aligned shape like a polycube (e.g., a sphere mapped to a cube) is inevitable. We, thus, aim here to improve the quality of the mapping by adding hexahedral elements only in selected and limited areas of the hex-mesh topology.

State-of-the-art padding operation -global padding in the whole hex-mesh surface -can sometimes have the opposite effect of locally worsening the quality of the mesh instead of improving it. Our idea is to analyze the quality of the mapping between the polycube space and the object space, and then to identify where a selective insertion of sheets of hexahedra can improve the quality of the final mesh. We change the topology of the mesh only where it is needed, and we leave it untouched in areas where the quality is already acceptable.

Rationale and Overview

The starting point of our approach is the mapping quality of the hexahedral mesh elements. We choose as error metric the Scaled Jacobian (SJ) which computes, for each hexahedron, a value between 1 and 1. When SJ = 1 the hexahedron is a perfect cube with highest possible quality. When SJ  0 the hexahedron is flipped and thus unfit to further processing. The quality of a hexahedron is judged low when SJ is smaller than a user-specified threshold.

We now select a simple example to convey the intuition behind our approach. Assume a box-cylinder object formed by a right box and a cylinder on top, as depicted by Figure 2. When computing the hex-mesh of such an object by gridding its polycube, the distortion is localized in the cylinder portion of the object, where the polycube forms right angles.

Performing global padding on this model decreases the quality of elements which are of high quality in the input model, as shown by second row of Figure 3. By analyzing the local mapping distortion and performing selective padding instead, we insert a sheet of It is forthwith evident the local negative effect of global padding where it is not necessary. Performing global padding increases the quality around the cylinder surface but also decreases the quality of the right box, while adding 850 elements. Our selective padding instead increases the quality only where it is needed (improving both the minimal SJ and the average SJ) and adds only 336 elements.

Figure 4 depicts a sub-set of three template configurations used to extrude the facets of a hexahedron element. To minimize the number of added hexahedra, we allow the new layers of hexahedra to turn, around an edge (edge turn) or a vertex (vertex turn). Note that an edge turn introduces 2 singular edges, while a vertex turn introduces 7 singular edges. For this reason, in the following sections, we propose a way to balance between introduced elements and introduced singularities. Overview. Our approach takes as input a volumetric mesh and its polycubes generated by the Polycut algorithm [LVS ⇤ 13]. It proceeds in two main steps: mapping analysis (Section 5) followed by selective padding (Section 4). The selective padding step refines the mesh structure to provide additional elements, and hence degrees of freedom, for existing mesh optimization approaches.

More specifically, the role of the mapping analysis step is to identify a set of facets delineating hexahedra with high mapping distortion. According to an analysis of the dihedral angles between facets of bad quality hexahedra, we determine the set of facets used as padding constraints for the global solver.

Starting from these facets, our goal is to selectively pad the mesh with just-enough hexahedra to reduce distortion where needed. Our solution to preserve the structure of the hexahedral mesh is to proceed by sheet insertion, the sheets being decomposed into a series of consistent facet extrusion operators.

We formulate a model with a constrained objective function with binary variables: one variable per facet; one variable per edge (to count edge turns) and one variable per vertex (to count vertex turns). The goal of the optimization is to find a satisfactory balance between quality, number of elements and number of singularities.

Selective Padding

We pose the local padding problem as a binary all-linear problem with a set of constraints that preserve the consistency of the topological hex-mesh structure. Given a set of Hard-constrained Facets (HF) which should be padded, the solution of the binary problem yields a set of facets which allow a consistent padding that includes at least the facets from HF. Our formulation enables the user to trade the number of additionally padded facets for the number of singularities introduced via padding.

Simple Binary Problem

Let M = (V, E, F, H) be the polycube-based input hex-mesh composed of vertices, edges, facets and hexahedra. We create a binary variable x f i for each facet specifying whether this facet should be padded (through extrusion) or not. By definition, every inner edge of a polycube-based hex-mesh has four incident facets. If we need to change the mesh by extruding facets (as explained in Section 4.3), we can extrude, for an inner edge, only two or four facets without creating topological inconsistencies. The outer edges, on the other hand, can be incident to 1, 2 or 3 hexahedra. We denote by E1H, E2H and E3H these sub-sets of edges.

Our first objective is to pad as few facets as possible:

E padding = |H| 2 3 Â fi2F\HF x f i (1)
where the term |H| 2 3 is added to achieve resolution independence.

We enforce two constraints during optimization. Firstly, all facets in HF must be padded:

x f i = 1 8 f i 2 HF (2)
Secondly, in order to achieve a valid, hex-topology preserving padding, we require the number of padded facets around each edge, but the ones belonging to E1H and E2H, to be even:

 fi2F(e j ) x f i = 2k j 8e j 2 E \ (E1H [ E2H) (3)
where k j is an integer variable defined for each edge involved in this constraint, and F(e j ) is the set of the facets incident to edge e j . We need to treat them in different ways with respect to the constraint.

1. Edges in E1H can have 0, 1 or 2 selected incident facets. All these configurations are suitable. 2. Edges in E2H can only have three incident facets, two on the surface and one inside. Almost all the possible paddings are legal: (i) the padding of the inner facet, (ii) the padding of both surface facets, (iii) the padding of one of the surface facets along with the inner one, causing an edge turn, and (iv) the padding of all three facets. We insert in the model a custom constraint to avoid the case of padding only one of the surface facets that would cause a topological inconsistency. 3. Edges in E3H are covered in the general constraints of Eq 3.

We are aware that in this way we are excluding legal cases but, considering how we define the insertion of new layers via facet extrusion, they would not generate valid solutions. We are, moreover, not limiting the possibility to reach an optimal solution. Figure 5 illustrates an example solution of this simple binary problem. We extend next our formulation in order to reduce the number of inserted singularities.

Binary Problem Extension

Minimizing E padding alone under constraints 2 and 3 yields a valid solution with the lowest number of extra elements. In practice, however, this may not always be the desired solution. As depicted by c 2018 The Author(s) . We thus extend our problem formulation with a binary variable te i for every edge e i 2 E, recording a turn configuration at the location of e i (cf. Figure 4 middle) and thus the introduction of a pair of valence 3 and 5 edges. Additionally, we add a binary variable tv l for every vertex v l 2 V , recording a vertex turn configuration (cf. Figure 4 bottom). Using these edge and vertex variables we can formulate the additional objective of keeping the number of introduced singularities low (e.g., minimizing the number of layer turns) via:

E complexity = |H| 1 3 Â e j 2E ⇤ \E1H te j + Â vl 2V ⇤ \V 1H tv l (4)
where E ⇤ and V ⇤ are the sub-sets of edges and vertices that are incident to two orthogonal facets of HF since these layer's turns are unavoidable and will be in the final solution. As for edges, V 1H denotes the sub-sets of vertices incident to only one hexahedron.

As before we add the term |H| 1 3 to the left sum to render the formulation independent from the hex-mesh resolution.

In order to ensure that the indicator variables te i are 1 if and only if an edge turn configuration is present, we add the following constraint:

te j = |x f i x f k | 8e j 2 E ⇤ \ E1H, ~fi = fk and f i , f k 2 F(e j ) (5)
where F(e j ) denotes the set of facets incident to edge e j . According to Equation 5, to detect a possible edge turn in the edge e j we consider a pair of facets f i , f k 2 F(e j ) having the same orientation ( ~fi = fk ). As can be observed in Figure 6, the value of the subtraction |x f i x f k | determines whether an edge turn is present in e j .

Similarly, to ensure that the indicator variables tv l are 1 exactly when a vertex turn configuration is present we add the following constraint:

tv l = |te i te k | 8v l 2 V ⇤ \V 1H, ẽi = ẽk and e i , e k 2 E(v l ) (6)
where E(v l ) denotes the set of edges incident to vertex v l . According to Equation 6, to find a vertex turn we consider a pair of edges e i , e k in the set of edges incident to v l in E(v l ) having the same orientation (ẽ i = ẽk ). We subtract, in absolute value, the te variables of the selected edges and can thus detect whether a vertex turn is present in v l .

Finally, under constraints 2, 3, 5 and 6, we optimize a linear combination of E padding and E complexity :

min E = E padding + l • E complexity (7)
Adjusting the coefficient l provides the user with a means to trade the number of extra elements for the number of extra singularities. Figure 12 illustrates the output solution for three different values of l.

Note also that, satisfying Equation 2 for all surface facets is a feasible solution, so we can say that a solution satisfying the aforementioned constraints always exists.

Sheet Insertion

The output of the solver is a set of facets PF representing the areas where one or several hexahedral sheets must be inserted. We generate the padding layer by extruding each facet f i 2 PF and transforming it into a hexahedron. The set of extruded facets forms the new layer. We extrude each facet in both directions, considering a fraction of the edge lengths of incident hexahedra as a reference, except for the surface facets which are extruded only towards the inside. The final structure of the mesh, including the new elements, is deduced from analyzing the global configuration of facets in PF. As shown by Figure 7, the extra singularities are decided in accordance to the adjacent facets.

Mapping Analysis

In the previous section we described a method which, given a set of facets HF that need to be padded, finds a complete set of facets which can be extruded while preserving the hex-mesh structure. We now detail how we find a suitable set HF such that the completed padding improves the quality of the hex-mesh.

Distortion per Facet

As polycube-based hex-meshes are commonly well shaped inside the domain, our main idea is to analyze the boundary hexahedra in order to identify the ones which would benefit from being padded.

The quality of a hex-mesh is commonly measured via the Scaled Jacobian (SJ). For each hexahedron in the hex-mesh we compute the minimum determinant of the Jacobian matrix, evaluated at each of its 8 corners, and the center of the element, divided by the corresponding edge lengths. According to the Verdict manual [SEK ⇤ 07], a "good" quality hex-mesh should have only hexahedra h i such that SJ(h i ) 0.5. As we observed in our experiments, improving the quality of the mesh beyond this value is sometimes possible. Therefore, we decided to limit our analysis to the boundary hexahedra whose Scaled Jacobian is lower than T = 0.6, considering the remaining ones as already good. This value may be, of course, adjusted to fit application-specific requirements.

The quality of a boundary hexahedron can often be improved by padding. However, padding all its facets is not suitable (cf. Figure 3). We must therefore transition from a distortion measure per hexahedron to a distortion measure per facet, where a high facet distortion indicates that padding is required. The set HF of facets which must be padded can then be simply defined as those facets whose distortion measure exceeds a user-specified threshold T derived from the quality requirement of a specific application.

We do the transition in two steps. First, we define a distortion measure per boundary edge based on the dihedral angle between incident facets. The distortion of a facet is then defined as the maximum of the distortion of the four facet edges. In Section 4.1, we introduce the sub-sets of edges E1H, E2H and E3H as the edges respectively incident to 1, 2 and 3 hexahedra. Since padding the facets incident on the edges of E2H introduces two new hexahedra which again share a E2H edge, distortion is unlikely to improve. We therefore only consider boundary edges in E1H and E3H, both shown in Figure 8. For each boundary singular edge, belonging to either E1H or E3H, we compute the dihedral angles between its incident boundary facets as q = ! n i • ! n j (see Figure 8 to identify n i and n j ), and we measure how much they deviate from their ideal values (90 for the E1H edges or 270 for the E3H ones). For E1H surface edges, we define D(e) = q if q 0.5 and 0 otherwise. For E3H surface edges, we compute q based on the two incident inner facets f i and f j , because they determine the angle we need to split in case of high distortion. We define D(e) = |q| if q  0.5 and 0 otherwise. We finally assign a value D( f ), between 0 and 1, to the surface facets, to record how much it is necessary to extrude them in order to improve the quality of the mesh. We define D( f ) as the maximum distortion of the four incident edges: D( f ) = max e2 f D(e). In this step we ignore areas of the mesh composed only by hexahedra with SJ T .

Padded Facets

Each facet of the polycube hex-mesh surface is now assigned a distortion value. In our experiments, we observed that the direct use of these values to determine the set of constrained padding facets HF is not ideal as they may form a fragmented set with isolated low distortion facets in the middle of high distortion patches, and viceversa (see Figure 10). To obtain a more consistent set of uniform patches of facets HF we first smoothly propagate the distortion values, then define HF via the solution of a Max Flow -Min Cut problem. To propagate the distortion values associated with the surface facets, we apply a simple iterative flooding algorithm that starts from facets with D( f ) 6 = 0 and fills the adjacent empty ones. The process is iterated for a maximum number of steps denoted by k in the following formula:

D( f i ) n+1 = D( f j ) n • e n 2 8 f i 2 N b ( f j ), 0  n  k (8) c 2018 The Author(s) .
where n is the current iteration, e n 2 is a term to favor a soft propagation of the deformation values (see Figure 9), and N b ( f i ) is the set of neighboring boundary facets of f i . The described formula is applied only to facets with D( f ) = 0, while the facets with a value different from 0 are not changed.

The Max Flow -Min Cut graph is defined by two nodes for the two used labels (L 1 for facets f 2 HF and L 0 for the other facets), a node for each surface facet, an arc between adjacent surface facets and an arc between facet nodes and label nodes. We then formulate the Max Flow -Min Cut problem as follows:

E(L) = Â f 2SH P f (L f ) + Â h fp, fqi Ppq(Lp, Lq) ( 9
)
where SH is the set of hexahedra with at least one face on the boundary. P f (L f ) represents a penalty for cutting an arc between a facet f i and the label L fi . We define P fi as:

P fi (L i ) = 1 |l i D( f i )| ( 10
)
where l i = 0 for P fi (L 0 ) and l i = 1 for P fi (L 1 ). Assume D( f i ) is close to 1: P fi (L 0 ) ' 0 and P fi (L 1 ) ' 1. Therefore, it is convenient to cut the arc between f i and L 0 and to assign to f i the label L 1 . Note that assigning the facet f i to the label L 1 corresponds to inserting f i in the HF set of hard constraints.

The right sum represents the penalty for cutting an arc between two adjacent facets fp and fq. In other words, it assigns a price to assigning two labels Lp and Lq to two adjacent facets fp and fq and Lp 6 = Lq. We define Ppq as: The solution of the Max Flow -Min Cut returns a collection of facets ready to be used as hard constraints (HF) in the binary problem formulation described in Section 4. As shown by Figure 10 they are organized in consistent patches with neither holes nor isolated facets.

Ppq(Lp, Lq) = [1 |D( fp) D( fq)|] • [D( fp) + D( fq)] (11
While the above steps may appear ad-hoc, they offer a robust way to produce homogeneous patches of hard constraints. We also considered the alternative solution to compute one distortion value per facet then use such values as soft constraints by modifying the objective function. However, considering the negligible time required by the graph-cut solve step (fractions of a second), we decided to keep this method in order to provide the solver with a set of hard constraints. By using hard constraints obtained as described above, we achieve better results (concerning both quality and singularity count) in considerably lower time. Moreover, by visual inspection of the hard constraints, users are provided with the guarantee that they will be part of the final solution.

Results

All our experiments are conducted on a computer equipped with an Intel Core i7 4GHz processor, 16GB of RAM and Linux OS. We used the data structures of the Cinolib [Liv17] library for developing our C++ suite of code, the Max-Flow algorithm proposed in [START_REF] Boykov | An experimental comparison of min-cut/max-flow algorithms for energy minimization in vision[END_REF] to compute the hard constraints set, and Gurobi [GO16] as numerical solver. We used Polycut [LVS ⇤ 13] to produce polycube-maps, implemented the meshing pipeline described in [START_REF] Gregson | All-Hex Mesh Generation via Volumetric PolyCube Deformation[END_REF] to produce all the polycube-based hex-meshes, and used the hex-mesh optimizer described in [LSVT15] to optimize the showed results. For a fair comparison, we apply the optimizer in the no-padded version and in the ones with global and selective padding, to show the differences in using it to the different meshes. During optimization we set the surface attraction as high as possible to preserve the original shape of the model. In this way, we always obtain, for each one of the models in this paper, a maximum Hausdorff distance lower than 0.009 w.r.t. the bounded box diagonal.

In all the figures depicting results we use the following colorcode: the mesh separatrices are depicted in red, and the quality of the mesh elements ranges from green (good) to red (bad). For each model we compare the hex-mesh directly derived from the polycube with the hex-meshes with global and selective padding applied. Table 1 reports the most relevant data for our models: hexahedra count; quality of hexahedra in terms of minimum and average Scaled Jacobian; number of singularities (vertices and edges); time required by the solver to identify the set of facets to extrude. The sheet insertion operation takes less than a few seconds and thus we do not report it, since it is negligible. 1: Statistics. We compare three meshes: input, with global padding and with our selective padding method. For each model we report the number of hexahedra (#H), the number of singular vertices and edges (#Sv and #Se), and the minimum and average Scaled Jacobian (mSJ and aSJ). For the Chamfer and Lego models we record the results obtained with different l values, to measure the impact of l over the final mesh. For the "Double hinge" model we report the version without padding the holes (NH) and the one with padded holes (WH). It shows that padding concave angles, sometimes, does not significantly improve the final quality. We apply our algorithm on the models Test 1, 2, 4 and 5 taken from [START_REF] Wang R | Hex mesh topological improvement based on frame field and sheet adjustment[END_REF].

Model

As we commented in the introduction, the padding operation usually improves the global quality of a volumetric mesh. We are able to show that our selective padding allows to obtain a substantial quality improvement over the global padding. In Table 1 we also show that, on the top of obtaining a better quality, our method adds fewer extra elements than global padding. The only trade-off to pay is, sometimes, the increased number of singularities.

Table 2 records a brief comparison between the results of our algorithm and the one described in [START_REF] Wang R | Hex mesh topological improvement based on frame field and sheet adjustment[END_REF], applied to the same domains. As we have not been granted access to the original software, we performed comparisons on a series of results produced by their algorithm. Therefore, starting from the same shape, we extract and optimize a polycube-based structured hex-mesh, while they analyze and optimize an unstructured one. This allows us to start from a regular structure of good quality inside the shape, and apply our algorithm just on the surface. We obtain comparable or better results for both the minimum and average Scaled Jacobian.

Extra elements vs extra singularities. The improvement of the hex-mesh structure requires to find the right trade-off between the number of extra elements and the number of extra singularities. The user-specified l parameter controls the number of turns of the inserted sheets and therefore the addition of singular vertices to the mesh structure. Figure 12 shows the differences in the final result for the Lego model, depending on the used l value. With a low l value the solver can insert turns everywhere in the mesh structure, with the goal of padding as few facets as possible. With a high l value the number of turns -and new singularities -is limited at the price of adding more extra elements. Finding the well-balanced l value allows to reach a reasonable compromise between extra elements, extra singularities and final quality, as it is shown in Figure 12d. Timing. The solver's timing depends on the number of elements in the mesh and its shape. The solution space varies depending on the input structure of the mesh. We usually use "coarse" meshes to solve the padding problem since a refinement step is always possible in post-processing (e.g., splitting each hexahedron in eight sub-hexahedra). When a model is complex enough not to allow to start from coarse hex-meshes, the solver can require up to several minutes of computation to produce the set of facets to extrude.

Model

Mechanical parts vs. organic shapes. As mentioned in Section 1, polycubes are an ideal tool for generating hexahedral meshes of quite regular shapes with a limited set of details. In this class of objects, mechanical pieces and CAD models are relevant subsets.

As it is clear from our results, we focus our attention on this class of objects. Indeed, the application of our selective padding on organic and free-form polycube-based shapes produces the same result of global padding. In Figure 14, we show the use of our algorithm on the Bunny model. It is evident that, to improve the mesh quality, we need to push inside all the singular edges of the mesh.

Limitations

While our experiments show excellent results, improving the quality of both the original and the global padding models, we cannot prove that our approach reaches the maximum possible quality for the chosen application field. Searching for the best trade-off between complexity and distortion would require either trying all possible l values and then selecting the one yielding the best results, or proceeding by dichotomy. According to our experiments, the value of the l parameter which leads to the best results depends on the shape of the model. It is thus not possible to suggest the silver bullet value of l that could work for any model. The choice of l can be a fine-tuning task and we let the user set it interactively.

Our method can detect distortion in the proximity of all concave tunnels in models with genus greater than 0. For simplicity, we refer to these particular shapes as "Holes". However, padding the holes is not always relevant for improving the overall quality. Experiments carried out on the Double hinge model (Figure 13) show that the quality can even decrease.

In the end, we produce a new hex-mesh topology suitable for obtaining a good quality final mesh, but the final quality is strictly dependent on the hex-mesh optimizer. We use the approach from [LSVT15] but, considering that the untangling is still an unsolved task in 3D, our results may not have the desired quality despite their good structure. Timing. The solver's timing depends on the number of elements in the mesh and its shape. The solution space varies depending on the input structure of the mesh. We usually use "coarse" meshes to solve the padding problem since a refinement step is always possible in post-processing (e.g., splitting each hexahedron in eight sub-hexahedra). When a model is complex enough not to allow to start from coarse hex-meshes, the solver can require up to several minutes of computation to produce the set of facets to extrude.

Mechanical parts vs. organic shapes. As mentioned in Section 1, polycubes are an ideal tool for generating hexahedral meshes of quite regular shapes with a limited set of details. In this class of objects, mechanical pieces and CAD models are relevant subsets.

As it is clear from our results, we focus our attention on this class of objects. Indeed, the application of our selective padding on organic and free-form polycube-based shapes produces the same result of global padding. In Figure 14, we show the use of our algorithm on the Bunny model. It is evident that, to improve the mesh quality, we need to push inside all the singular edges of the mesh.

Limitations

While our experiments show excellent results, improving the quality of both the original and the global padding models, we cannot prove that our approach reaches the maximum possible quality for the chosen application field. Searching for the best trade-off between complexity and distortion would require either trying all possible l values and then selecting the one yielding the best results, or proceeding by dichotomy. According to our experiments, the value of the l parameter which leads to the best results depends on the shape of the model. It is thus not possible to suggest the silver bullet value of l that could work for any model. The choice of l can be a fine-tuning task and we let the user set it interactively.

Our method can detect distortion in the proximity of all concave tunnels in models with genus greater than 0. For simplicity, we refer to these particular shapes as "Holes". However, padding the holes is not always relevant for improving the overall quality. Experiments carried out on the Double hinge model (Figure 13) show that the quality can even decrease.

In the end, we produce a new hex-mesh topology suitable for obtaining a good quality final mesh, but the final quality is strictly dependent on the hex-mesh optimizer. We use the approach from [LSVT15] but, considering that the untangling is still an unsolved task in 3D, our results may not have the desired quality despite their good structure.

Conclusion and Future Work

We introduced in this paper a novel pipeline for the generation of quality hexahedral meshes, where high quality of a hexahedron refers to low deviation from the perfect cube. Our pipeline utilizes a polycube mapping for decomposing the input 3D domain into portions which are simple to discretize into a hexahedral mesh. In this pipeline our main contribution is a selective padding step which automatically adds sheets of hexahedra only where they are needed, in order to increase the global quality of the output hexahedral meshes. These sheets can form turns inside the domain, which induce extra edge and vertex singularities. Compared to global padding or greedy straight sheet insertion, our approach improves the global quality while generating fewer hexahedra. It works at its best when applied to input domains bounded by a piecewise planar surface, which is typical of mechanical parts. On concave holes the gain is only moderate.

In the future we plan to explore an automatic parameter selection approach in order to find the best balance between quality and complexity. We also wish to explore a global optimization approach to select the optimal set of constrained facets that yields the maximum quality. We then intend to analyze whether a combination of padding and inverse padding (remove hexahedral layers from the mesh) can provide better results. Furthermore, we plan to extend our approach to more general polycube-based hex-mesh structures, like those used by Fang et al. [START_REF] Fang | All-hex meshing using closed-form induced polycube[END_REF].

Figure 2 :

 2 Figure 2: Introducing distortion. From left to right: the tet-mesh of the input model, its polycube (tet-mesh) and the final hex-mesh with distortion highlighted. Here and in the other figures we use the color ramp on the right as a quality indicator where red stands for highly distorted elements, yellow for medium quality, and green stands for good quality elements.

Figure 3 :

 3 Figure 3: Global vs selective padding applied to the box-cylinder model. Top: global (left) and selective (right) padding. Bottom: we show the quality of the elements induced by the two padding strategies.

  c 2018 The Author(s)

Figure 4 :

 4 Figure 4: Padding via facet extrusion. Top to bottom: padding a single facet (one added hexahedron), padding two facets (two added hexahedra and one edge turn) and padding three facets (three added hexahedra and a vertex turn).

Figure 5 :

 5 Figure 5: An example of our padding strategy. In (a) the red set of facets has to be extruded; in (b) the solution obtained with the basic formulation; the mesh in (c) is similar but, in this case, padding straight to the bottom would insert many hexahedra. Our simple formulation finds the solution in (d) with fewer added hexahedra by introducing edge and vertex turns inside the mesh.

Figure

  Figure 5d, the padding may introduce singularities inside the mesh. While adding few more hexahedra is often better than introducing extra singularities, turns of the padding layer provides us with a means to avoid generating many extra elements. A trade-off between extra elements and extra singularities is required.

Figure 6 :

 6 Figure 6: Detection of an edge turn. Left: two selected facets with similar orientation. Middle: |x f i x f k | = 0 hence no edge turn is detected. Right: |x f i x f k | = 1 hence an edge turn is detected.

Figure 7 :

 7 Figure 7: Example of padding configurations and associated sheet insertions. When four facets are selected (a) the resulting new layer of hexahedra introduces new singularities, depicted with red circles (b). If the same four facets are selected, together with four more ones (c) the final layer is different (d), and there are no extra singularities added.

Figure 8 :

 8 Figure 8: Example of E1H edge (left) and of E3H edge (right).

Figure 9 :

 9 Figure 9: Example of propagation of the D( f ) values on the domain boundary, with k = 4.

  ) where [1 |D( fp) D( fq)|] measures the difference between D( fp) and D( fq), and [D( fp) + D( fq)] favors cuts between arcs in low distortion areas.

Figure 10 :

 10 Figure 10: Computing hard constraints. Left: D( f ) values associated to each facet after the first step; notice, on the top model the isolated high-distortion elements, on the bottom model the isolated low-distortion elements (pointed by the arrows). Middle: D( f ) values after the propagation. Right: final constraints resulting from solving the Max Flow -Min Cut formulation, which have filled the gaps.

Figure 11 :

 11 Figure11: A gallery of results obtained with our method compared to the input mesh and global padding. GP stands for global padding and SP for selective padding (ours). In red are singular edges. Colors indicates quality as described in Figure2. In the front half of the mesh we show only the hexahedra with a Scaled Jacobian below 0.75, the others are transparent.

Figure 12 :

 12 Figure 12: Padding the "Lego" mesh. Global padding (b), and selective padding (c, d, e). The padding layers are shown in blue in (f, g, h).

Figure 13 :

 13 Figure 13: Padding the "Double hinge" mesh. Global padding (b), and selective padding, including holes (c) and without holes (d). The padding layers are shown in blue in (e , f).

Figure 12 :

 12 Figure 12: Padding the "Lego" mesh. Global padding (b), and selective padding (c, d, e). The padding layers are shown in blue in (f, g, h).

Figure 13 :

 13 Figure 13: Padding the "Double hinge" mesh. Global padding (b), and selective padding, including holes (c) and without holes (d). The padding layers are shown in blue in (e , f).

Figure 14 :

 14 Figure 14: Selective padding in organic shapes. On the left the model without padding, on the right the selective padding (equivalent to the global one) applied on the same model.

Table 2 :

 2 A comparison between results obtained with our approach and those obtained with[START_REF] Wang R | Hex mesh topological improvement based on frame field and sheet adjustment[END_REF].

	[WGZC18]		Ours
	mSJ aSJ mSJ aSJ
	Test 1 .35	.94	.66	.98
	Test 2 .34	.88	.73	.96
	Test 4 .64	.96	.68	.99
	Test 5 .39	.89	.41	.98

c 2018 The Author(s)

Acknowledgments

We thank the authors of [WGZC18] for sharing with us the models used in their article. Gianmarco Cherchi gratefully acknowledges the European Erasus+ program for funding his period at INRIA. This work is partly financed by the DSURF PRIN 2015 (2015B8TRFM) project.


Davide Scaramuzza
email: davide.scaramuzza@ieee.org

Nicolas Criblez
email: nicolas.criblez@gmail.com

Agostino Martinelli
email: agostino.martinelli@ieee.org

Roland Siegwart
email: r.siegwart@ieee.org

Robust Feature Extraction and Matching for Omnidirectional Images

This paper presents a new and robust method for extracting and matching visual vertical features between images taken by an omnidirectional camera. Matching robustness is achieved by creating a descriptor which is unique and distinctive for each feature. Furthermore, the proposed descriptor is invariant to rotation. The robustness of the approach is validated through real experiments with a wheeled robot equipped with an omnidirectional camera. We show that vertical lines are very well extracted and tracked during the robot motion. At the end, we also present an application of our algorithm to the robot simultaneous localization and mapping in an unknown environment.

Introduction

One of the most important problems in vision based robot navigation systems is the search for correspondences in images taken from different viewpoints. In the last decades, the feature correspondence problem has been largely investigated for standard perspective cameras. Furthermore, some works have provided robust solutions for wide-baseline stereo matching, structure from motion, ego-motion estimation, and robot navigation (see [START_REF] Matas | Robust wide baseline stereo from maximally stable extremal regions[END_REF], [START_REF] Kadir | An affine invariant salient region detector[END_REF], [START_REF] Mikolajczyk | Indexing based on scale invariant interest points[END_REF], [START_REF] Lowe | Distinctive image features from scale-invariant keypoints[END_REF], [START_REF] Lindeberg | Feature detection with automatic scale selection[END_REF], [START_REF] Mikolajczyk | An affine invariant interest point detector[END_REF], [START_REF] Tuytelaars | Matching widely separated views based on affine invariant regions[END_REF], [START_REF] Baumberg | Reliable feature matching across widely separated views[END_REF], and [START_REF] Tuytelaars | Matching widely separated views based on affine invariant regions[END_REF]). Some of these works normalize the region around each detected feature using a local affine transformation, which attempts to compensate for the distortion introduced by the perspective projection. However, such methods cannot be directly applied to images taken by omnidirectional imaging devices because of the non-linear distortions introduced by their large field of view. In order to apply those methods, one needs first to generate a perspective view out of the omnidirectional image, provided that the imaging model is known and that the omnidirectional camera possesses a single effective viewpoint [START_REF] Shree | Catadioptric omnidirectional camera[END_REF]. An application of this approach can be found in [START_REF] Mauthner | Region matching for omnidirectional images using virtual camera planes[END_REF]. There, the authors generate perspective views from each region of interest of the omnidirectional image. This image unwrapping removes the distortions of the omnidirectional imaging device and enables the use of state-of-theart wide-baseline algorithms designed for perspective cameras. Nevertheless, other researchers have attempted to apply to omnidirectional images standard feature detectors and matching techniques, which have been traditionally employed for perspective images. In [START_REF] Micusik | Structure from Motion with Wide Circular Field of View Cameras[END_REF], for instance, the authors check the candidate correspondences between two views using RANSAC algorithm. Finally, other works have been developed, which extract one-dimensional features from new images called Epipolar plane images, under the assumption that the camera is moving on a flat surface [START_REF] Briggs | Robot Navigation Using 1D Panoramic Images[END_REF]. These images are generated by converting each omnidirectional picture into a 1D circular image, which is obtained by averaging the scan lines of a cylindrical panorama. Then, 1D features are extracted directly from such kinds of images. In this paper, the features we want to track from omnidirectional images are real world vertical features, which are predominant in structured environments. In our experiments, we used a wheeled robot equipped with an omnidirectional camera, which had the camera axis perpendicular to the direction of motion of the robot. Because of this settings and assuming the environment to be flat, all world vertical lines project into radial lines on the image plane. The novelty of this paper consists of a robust method to match vertical lines between images taken by an omnidirectional camera during the motion of the robot. Matching robustness is achieved by creating a descriptor which is unique and distinctive for each feature. Furthermore, the proposed descriptor is invariant to rotation. This descriptor is based on the image gradients. The robustness of the approach was validated through real experiments by using a robot equipped with an omnidirectional camera. In this paper, we show that vertical lines are very well extracted and tracked during the robot motion. At the end, we also present an application of our algorithm to the robot Simultaneous Localization And Mapping (SLAM) in an unknown environment. This paper is organized as follows. In section 2, we will discuss the vertical line extraction. In section 3, we will illustrate how to mark a feature by means of a descriptor, and, in section 4, we will describe the feature matching process. Finally, in section 5, we will present our experimental results and the application of our algorithm to SLAM.

Vertical Line Extraction

Our platform consists of a wheeled robot equipped with an omnidirectional camera looking upwards. The main advantage of such kind of camera is that it provides a 360 • field of view of the scene, which gives a very rich and sparse information. Figure 1 shows a sample picture taken by our omnidirectional camera. As the circular external boundary of the mirror is visible in the image, we use a circle detector to determine the location of the center. In our arrangement, we set the camera-mirror system perpendicular to the floor where the robot moves. This guarantees that all vertical lines of the environment converge towards the image center. To extract the vertical lines, we first compute the image gradients (e.g. we used a Sobel filter), and then we keep only those gradients whose orientation looks towards the image center up to ±5 • . This 10 • tolerance allows us to deal with the effects of the floor irregularities in the projections of the 3D vertical lines. After this filtering, we apply non-maxima suppression and we end up with a binary edge map (Fig. 2). The next step will be identifying the most reliable vertical lines. To this end, we divide the omnidirectional image into 720 predefined uniform sectors, which give us an angular resolution of 0.5 • . By summing up all binary pixels that vote for the same sector, we obtain the histogram shown in Fig. 3. As observed in Fig. 2, there are many potential vertical lines in a structured environment. To keep the most reliable and robust features, we choose only those lines whose length covers at least half of the angle of view (Fig. 3). Finally, we use non-maxima suppression to avoid the features to be too close (Fig. 4).

Building the Descriptor

In section 4, we will describe our method for matching vertical lines between consecutive frames while the robot is moving in an unknown environment. To make the feature correspondence robust to false positives, each vertical line is given a descriptor, which is unique and distinctive for each feature. Furthermore, this descriptor is invariant to rotation. In this way, finding the correspondent of a vertical line can be done by looking for the line with the closest descriptor. In the next subsections, we will describe how to build this descriptor.

Rotation Invariance

To build the descriptor, we extract a predefined number of circular areas (namely 3 areas) in fixed position along a given radial line (Fig. 5). The centers of these circular areas are equally spaced and the radius is chosen such that the circles touch without overlapping. Then, each area is smoothed by a Gaussian window of variance σ and the image gradients (magnitude and phase) are computed within each of these areas. The rotation invariance is achieved redefining the gradient phase relatively to the radial line's angle. (Fig. 5).

Orientation Histograms

To make the descriptor robust to false matches, we split each circular area into two parts and consider each one individually (Fig. 6). In this way, we preserve the information about what we have on the left and right sides of the feature. For each side of each circular area, we compute an orientation histogram (Fig. 7) of all gradient vectors. The whole orientation space (from -π to π) is divided into N b equally spaced bins. In order to decide how much of a certain gradient magnitude m belongs to the adjacent inferior bin b and how much to the adjacent superior bin, each magnitude m is weighted by the factor (1 -w), where

w = N b • o -b 2π (1)
with o being the observed orientation in radians. Thus, m(1 -w) will vote for the adjacent inferior bin, while mw will vote for the adjacent superior bin.

According to what has been mentioned so far, each bin contains the sum of the weighted gradient magnitudes which belong to the correspondent orientation interval. We observed that this weighted sum made the orientation histogram more robust to image noise. The reader observe that the orientation histogram is already rotation invariant because the gradient angles have been referred to the radial line's angle.

Building the Feature descriptor

The computed orientation histograms help to build the feature descriptor. Indeed, the descriptor is an N -element vector containing all histogram values of the circular areas. For instance, by extracting 3 circular areas for each vertical feature and choosing 30 bins for each histogram, the length of the feature descriptor will be

N = 3areas • 2parts • 30bins = 180 (2)
Finally, it is important to know that all feature descriptors have the same length N . To achieve slight illumination invariance, every descriptor is normalized to 1. This choice relies on the hypothesis that the image intensity changes linearly with illumination. Although this is not true in nature, this approximation proved to work properly.

Feature Matching

As every vertical feature has its own descriptor, the correspondent of a vertical line in the consecutive images can be searched among the features with the closest descriptor. As a distance measure between two vector descriptors A and B, we use the Euclidean distance:

d(A, B) = N k=1 (A(k) -B(k)) 2 (3)
As a consequence, the correspondent of a feature, in the current image, is expected to be the one, in the consecutive image, with the minimum distance. However, if a feature is no longer present in the next image, there will be a closest feature anyway. For this reason, we define three tests to decide whether a feature correspondent exists and which the correspondent is. Before describing the three criterions, let us introduce some definitions. Say {A 1 , A 2 , ..., A N A } and {B 1 , B 2 , ..., B N B } two sets of feature descriptors extracted at time t and t -1 respectively (where N A and N B are the number of features in the first and second image). Then, say

D i = {d(A i , B j ), j = 1, 2, ..., N B )} ( 4
)
the set of all distances between a given A i and all B j (j = 1, 2, , N B ). Finally, say minD i = min(D i ) the minimum of the distances between A i and all B j .

First Test

The first test checks that the distance from the closest descriptor is smaller than a given threshold. As the threshold depends on the length of the descriptor, we set

minD i = F 1 • N ( 5
)
where N is the descriptor length. By this criterion, we actually set a bound for the maximum acceptable distance to the closest descriptor.

Second Test

The second test checks that the distance from the closest descriptor is smaller enough than the mean of the distances from all other descriptors, that is:

minD i = F 2 • < D i > (6)
where < D i > is the mean value of D i and F 2 clearly ranges from 0 to 1. This criterion comes out of experimental results. In Table 1, we show a real comparison among the distances between the descriptor A 1 at time t and all the descriptors at time t -1. There, the descriptor B 1 is the correct correspondent of A 1 . The reader might also observe that its distance is smaller than the mean of all other distances.

Third Test

Finally, the third test checks that the distance from the closest descriptor is smaller than the distance from the second closest descriptor:

minD i = F 3 • SecondSmallestDistance, (7)
where F 3 clearly ranges from 0 to 1. As in the previous test, the third test raises from the observation that, if the correct correspondence exists, then there must be a big jump between the closest and the second closest descriptor. Factors F 1 , F 2 , and F 3 are to be determined experimentally. The empirical values we used for these parameters are shown in Table 2.

Table 2. The parameters used by our algorithm with their empirical values F1 = 0.0075 F2 = 0.55 F3 = 0.85

Experimental Results

In this section, we present some experimental results obtained by moving our robot in a real indoor environment. In the first subsection, we show the performance of our feature tracker during the motion of the robot, while in the second one, we present the results of our feature tracker applied to SLAM.

In these experiments, the robot was moving at about 0.15 m/s. The image size was 640x480 pixels and the frame rate was 3 Hz.

Performance of the Feature Tracking

In this experiment, we guided our robot through an office-like environment for about 70 meters. The results of the feature tracking are shown in Fig. 8. The plot refers to a short path of the whole trajectory while the robot was coping with an L-shaped trajectory. As the reader may observe, many features are detected and tracked over the time. Indeed, all lines of the plot appear to be smooth and homogeneous. Furthermore, the lines do not intersect, meaning that there was no false matching. We can also notice that the algorithm was able to match vertical elements even when the correspondent features were not observed in the previous image (e.g. observe the large gap between the dots pointed to by the arrow in Fig. 8). Indeed, when the correspondence is not found in the last frame, our algorithm starts looking into all previous frames (actually up to the twentieth frame), and stops when the correspondence is found. By zooming into the plot of Fig. 8, we found that some lines are given different numbers. For instance, feature number 24 is labeled as 31 after some frames. And the same happens with features 42 and 49. When this happens, it means that the algorithm found no correspondence for the current feature, and thus, the feature is labeled as a new entry, but in fact this is a false new entry. After having visually checked every single frame of the video sequence, we found 6 false matches and 22 false new entries. Comparing these errors to the 2631 corresponding pairs detected by the algorithm over the whole video sequence, we had 1.06% of mismatches. Furthermore, we found that the false matches occurred every time the camera was facing objects with repetitive texture. Thus, the ambiguity was caused by the presence of vertical elements which were almost identical. On the other hand, a few of false new entries occurred whenever the displacement of the robot between two successive images was big. However, the reader should observe that when a feature matches with no other previous feature of the last frames, it is better to believe this feature to be new rather than commit a false matching.

Application to SLAM

We applied our feature tracker to two important problems in autonomous navigation, that is, sensor self-calibration and SLAM. Regarding the former, the results can be found in [START_REF] Scaramuzza | A Flexible Technique for Accurate Omnidirectional Camera Calibration and Structure from Motion[END_REF], [START_REF] Scaramuzza | A Toolbox for Easily Calibrating Omnidirectional Cameras[END_REF], and [START_REF] Martinelli | Automatic Self-Calibration of a Vision System during Robot Motion[END_REF]. In this section, we show only the results we obtained for SLAM. We implemented the standard EKF based SLAM. In particular, the EKF estimates the vector:

X = [x r , y r , θ r , X 1 , Y 1 , ..., X N O , Y N O ] T (8)
where [x r , y r , θ r ] is the robot configuration, X i , Y i are the Cartesian coordinates of the i-feature in the map, and N O is the number of observed features. Our mobile robot is equipped with wheel encoders and with the same omnidirectional camera adopted in the experiments described in section 5.1. The bearing observations provided by the omnidirectional camera consist of the vector z whose components are:

z 1 = arctan( Y 1 -y r X 1 -x r ) . . . (9)
z NO = arctan( Y NO -y r X NO -x r )

To initialize a new feature in the map, consecutive bearing observations as in [START_REF] Martinelli | Automatic Self-Calibration of a Vision System during Robot Motion[END_REF], which refer to the same feature, are integrated with the odometry. Then, the estimation is improved by integrating the information coming from all the bearing observations through the EKF. The result is shown in Fig. 9 where both the robot trajectory and the position of the features are shown.

Conclusion

In this paper we introduced and discussed a new and robust method to extract and match vertical lines between images taken by an omnidirectional camera. The basic idea to achieve robust feature matching consists of creating a descriptor which is unique and distinctive for each feature. Furthermore, this descriptor is invariant to rotation. To evaluate the performance of our approach, we performed real experiments where we evaluated the quality of the matching. We conclude that the proposed approach is very robust and precise. Finally, we adopted the proposed method to implement an EKF based SLAM.

Fig. 1 .

 1 Fig. 1. An image taken by our omnidirectional camera.

Fig. 2 .

 2 Fig. 2. The final binary edge map after non-maxima suppression.

Fig. 3 .

 3 Fig.3. The number of binary pixels voting for a given orientation angle. The orientation axis ranges from 0 to 720 half degrees.

Fig. 4 .

 4 Fig. 4. Extraction of the most reliable vertical features from an omnidirectional image.

Fig. 5 .

 5 Fig. 5. Extraction of the circular areas.To have rotation invariance, the gradient phase is referred to the orientation of the vertical line.

Fig. 6 .

 6 Fig. 6. Two parts of a circular area.

Fig. 7 .

 7 Fig. 7. An example of gradient orientation histograms for the left and right sides of a circular area.

Fig. 8 .

 8 Fig. 8. Feature tracking during the robot motion. In y-axis, we have the angle of sight of each feature, and in the x-axis, the frame number. Each dot represents a feature detected in the current frame. The lines represent the tracked features. The number reported on some dot appears only when a new feature is detected.

Fig. 9 .

 9 Fig. 9. The results obtained by implementing a simple EKF-SLAM, which uses the proposed feature tracker. The black line is the trajectory estimated by using the odometry alone. The red line is the trajectory estimated by the EKF using both the odometry and vertical lines. The blue points represent the map ground truth provided by a laser range finder. The red circles are the detected verticals features.

Table 1 .

 1 The distances between the descriptor A1 at time t and all descriptors B j , j = 1, 2, .., NB at time t -1

	B1 B2 B3 B4 B5 B6 B7
	2.38 5.42 4.55 5.79 5.66 6.17 5.43

Acknowledgment

This work was conducted within the EU Integrated Projects COGNIRON (The Cognitive Robot Companion) and BACS (Bayesian Approach to Cognitive Systems). It was funded by the European Commission Division FP6-IST Future and Emerging Technologies under the contracts FP6-IST-002020 and FP6-IST-027140 respectively.


Keywords: nonlocal, damage, negative hardening Résumé nonlocal, endommagement, écrouissage négatif

come

Introduction

The effect and the formulation of boundary conditions -such as free edges, notches and initial cracks -remain an open question for nonlocal models. The main drawback of the classical nonlocal integral theory [1] consists in the nonphysical interaction, through the nonlocal averaging process, of points across a crack or a hole. The definition of natural boundary conditions of vanishing strain normal derivative at a free edge is still under discussion for gradient formulations [2,3]. The continuous nucleation of a crack of zero thickness is not so simple as the thickness of a localization band is more or less proportional to the internal length introduced. Local behavior along free edges -i.e. with a vanishing internal length -has been obtained by some authors [START_REF] Pijaudier-Cabot | chapter in Nonlocal Modeling of Materials Failure[END_REF]. The consideration of an internal length evolving with damage [START_REF] Geers | [END_REF]6] seems a way to properly bridge Damage Mechanics and Fracture Mechanics as the internal length may then vanish for large values of damage.

In the present note, one proposes a solution -bringing also questions -to these main difficulties. The idea is to keep the nonlocal averaging process but to quantify the distance between points as an effective distance, i.e. as a distance function for instance of the geometry and the matter encountered beween interacting points. One proposes to define such an effective distance with respect to a dynamic process: with respect to how information or wave propagates between interacting points. This can be made through the introduction of an internal time, constant, instead of a internal length, measured as evolving. Dynamics is important to define a link between a characteristic time and a characteristic length, either when viscosity is introduced [7,9] or when the physical defects obscuration phenomenon encountered in high speed dynamics and multi-fragmentation is taken into account [8]. Some authors even introduce the classical nonlocal theory by comparing the characteristic wavelength of the deformation field to an intrinsic length of the material [10], still a dynamics vocabulary. Remind also that the physics is often dynamic at small (discrete) scales as in quantuum mechanics, as in dislocations theories or as at the scale of the macromolecular chains in elastomers, even if the continuum macroscopic behavior is quasi-static. Let us then accept then the importance of the dynamics effects and of wave propagations. But as just mentioned in order to address quasi-static non-viscous structural failure.

Chapter 1

Wave propagation in plain, notched and damaged media

This part illustrates the wave propagation in different media. It contributes to the understanding of the differences or similarities between a notched and a damaged zone.

Finite Element modelling

In this chapter one considers the problem of an impact on a square plate of the Figure 1.1 (1.5 m×1.5 m×0.01 m). An impulse force F (t) according to the out of plane z-axis (close to a Dirac) is applied to the plate (Figure 1.2) in order to generate a wave propagation in this plate.

In order to numerically simulate this problem, one uses the Finite Element code Dyna3D to obtain the transient solution of a wave propagation. The plate is modeled with 3D under-integrated brick elements with an elastic material (E=35 GPa, ν=0,2, ρ=2400 kg/m 3 ). No boundaries conditions are applied in the simulations. The discrete equations of motion are solved with an explicit scheme so that the time step is directly linked to the element size (Courant's condition satisfied). In order to gain a good description of the wave propagation, one choses to have small elements (0.015 m×0.015 m×0.01 m). One obtains then approximately 360000 elements with a time step equal to 5.10 -4 ms. This time step will be the accuracy of the different computed times in the problem. The length of the notch is 0.5 m, its width 0.1 m. Figure 1.6 shows the map of the z-acceleration in the plate at t=100 ms. Notice again that one can see the wave reflection on the free faces of the plate and in particular the reflection due to the notch (upper lip). This means here Ẽ = E/100. Figure 1.8 shows the map of the z-acceleration in the plate at t=100 ms. One sees the wave reflection on the free edges of the plate as in the previous cases. It is important to see that in this case, the damaged zone acts almost as a free surface for the wave propagation. This is due to the large differences in impedance between the initial undamaged material and the damaged one.

Propagation in plain media

Propagation in notched media

Propagation in damaged media

Discussion

These simple examples illustrate the wave propagation without and with notches (nothing new). First, wave propagation is affected by the free edges (classical reflection). Second, the comparison between the real notch and the damaged zone allows to conclude that highly damaged material (low Young's modulus) behaves as a notch.

CHAPTER 1. WAVE PROPAGATION

These two properties are the basis of the nonlocal integral formulation proposed next.

Chapter 2

Nonlocal theory with internal time

Softening constitutive equations classically lead to spurious dissipation modes and to mesh dependency. The need of the definition -and the introduction -of an internal length in the models is now established. But in which form ? In a gradient form ? In an inegral form ? From an internal viscosity (delay-damage) combined then with dynamics ? The main idea of such regularizations is to average the "engine" -i.e. the thermodynamics force, denoted next V in the general case or Y or ˆ for damage -responsible for the strain localization. The procedure to define a nonlocal variable V nl from its local expression V introduces a characteristic length l c considered as a material parameter.

One do not intend here to solve the problem in the general case and focus only on the nonlocal integral theories.

Nonlocal integral theories -Boundary effect

The classical nonlocal theory [1] uses the integral

V nl (x x x) = 1 V r Ω ψ x x x -s s s l c V(s s s)ds V r (x x x) = Ω ψ x x x -s s s l c ds (2.1)
over the whole domain Ω in order to define a nonlocal quantity V nl -built from the local variable V -to be used in thermodynamics state or evolution laws. In Eq. (2.1) ψ is the nonlocal weight function, positive and decreasing with respect to the distance x x x -s s s between the considered point x x x and all the structure points s s s. The normalizing factor is denoted V r (x x x) and l c is an internal or characteristic length. Two classical ψ-functions are

ψ(ξ) = e -1 2 ξ 2 or ψ(ξ) = 1 -ξ 2 2
(2.

2)

The expressions (2.1) do then define the same contribution to V nl of points x x x and s s s

13

• across a crack than any points x x x and s s s separated by the same distance x x x -s s s (for instance path x x x 1 -x x x 6 versus path x x x 2 -x x x 5 of figure 2.1),

• across a damaged zone: for example the interaction between points x 3 -x 4 is not affected by the zone at a damage level D 0 , i.e. at a Young's E(1 -D 0 ) much lower than the Young's modulus of the virgin material.

This illustrates the need to deal with such a boudary condition effect but also to make equivalent a real crack and a highly damaged zone.

Specimen (a) Specimen (b)

x 1 x 2 x 3 x 6 x 5 x 4 x 1 x 2 x 3 x 6 x 5 x 4
Specimen (c)

x 1 x 2 x 3 x 6 x 5 x 4 D=D 0

Time instead of distance

In order to solve the problem, one proposes in this work to keep the nonlocal integral framework and to consider the information time propagation τ xs between points x x x and s s s instead of the classical distance x x x -s s s . A nonlocal theory with internal time is then simply defined by replacing Eq. (2.1) by:

V nl (x x x) = 1 V r Ω ψ τ xs τ c V(s s s)ds V r (x x x) = Ω ψ τ xs τ c ds (2.3)
with τ xs the information propagation time taken next as the time for a wave to propagate from point x x x to point s s s and τ c a material parameter. As the wave time propagation τ sx from point s s s to point x x x is identical to τ xs the nonlocal weight function thus built is symmetric. Such a nonlocal averaging process applies to physical laws of different nature.

Concerning continuum mechanics and elasticity, plasticity and/or damage, the variables V are often equivalent strains, accumulated plastic strains or strain energy densities.

Example: nonlocal damage models

For example, local isotropic damage theories for quasi-brittle materials define

• the damage as a loss of stiffness D = 1 -Ẽ/E with Ẽ (resp. E) the damaged (resp. initial) Young modulus, equation also rewritten in 3D as the elasticity law coupled with isotropic damage

σ σ σ = E E E(1 -D) : (2.4)
with σ σ σ, and E E E respectively the stress, the strain and the Hooke tensors.

• the damage evolution as a function of a local variable V either equal to the thermodynamics force Y = 1 2 : E E E : associated with damage (Marigo model) or to an equivalent strain as Mazars strain ˆ = + : + ,

D = g(V) (2.5)
with g a nonlinear function.

The non local damage law is simply written

D = g(V nl ) (2.6)
instead of Eq. (2.5) with the nonlocal averaging process (2.3).

For concrete, the microcracks due to tension are mainly orthogonal to the loading direction, when the microcracks due to compression are mainly parallel to the loading direction. The damage state has then to be represented by a tensorial variable D D D either a fourth rank tensor or a second rank tensor. The use of a second order damage tensor is more convenient for practical applications (as well as for the material parameters identification) and this is the choice which has been made. The damage anisotropy induced by either tension or compression is simply modeled by the consideration of damage evolution laws ensuring a damage rate proportional to the positive part of the strain tensor, i.e. a damage governed by the principal extensions [START_REF] Mazars | The unilateral behavior of damage concrete[END_REF].

The full set of constitutive equations for the local anisotropic damage model reads [START_REF] Desmorat | Modèle d'endommagement anisotrope avec forte dissymétrie traction/compression, 5è journéees du Regroupement Francophone pour la Recherche et la Formation sur le Béton (RF2B)[END_REF][START_REF] Desmorat | Nonlocal anisotropic damage model and related computational aspects for quasi-brittle materials[END_REF] • Elasticity,

= 1 + ν E σ σ σ - ν E tr σ σ σ 1 1 1 (2.7)
with E the Young modulus and ν the Poisson ratio.

• Effective stress,

σ σ σ = (1 1 1 -D D D) -1/2 σ σ σ D (1 1 1 -D D D) -1/2 D + 1 3 tr σ σ σ 1 -tr D D D --tr σ σ σ 1 1 1 (2.8)
where(.) D stands for the deviatoric part and . for the positive part of a scalar.

• Damage criterion,

f = ˆ -κ(tr D D D) (2.9)
so that the condition f < 0 → elastic loading or unloading, f = 0, ḟ = 0 → damage growth, where ˆ = + : + is Mazars equivalent strain built from the positive part + of the strain tensor and where

κ(tr D D D) = a • tan tr D D D aA + arctan κ 0 a (2.10)
• Induced damage anisotropy governed by the positive extensions,

D = λ 2 + (2.11)
In the rate independent formulation, the damage multiplier λ is determined from the consistency condition f = 0, ḟ = 0.

There are 5 material parameters introduced: E, ν for elasticity, κ 0 as damage threshold and A and a as damage parameters.

The model is simply made nonlocal, either from the classical integral theory or from the new integral nonlocal with internal time theory, by replacing Mazars equivalent strain ˆ by its nonlocal form in the damage criterion function, becoming

f = ˆ nl -κ(tr D D D)
(2.12) instead of Eq. (2.9) with the nonlocal averaging process (2.3).

Effective or "dynamic" distance -Vanishing internal length

In a plain and uncracked medium the internal length and internal time concepts are equivalent as x x x -s s s = c τ xs (2.13) with c the information celerity taken as the wave speed. Then, if the internal time is related to the internal length as l c = c τ c the weight functions are equal,

ψ x x x -s s s l c = ψ τ xs τ c (2.14)
In a homogeneously damaged medium at D = D 0 , the wave speed is proportional to the square root of the damaged Young's modulus and depends on the damage level as

c = c 1 -D 0 (2.15)
One has in this case

τ xs = x x x -s s s c (2.16) and τ xs τ c = x x x -s s s c τ c = x x x -s s s l c √ 1 -D 0 > x x x -s s s l c (2.17)
which shows that the effective or "dynamic" distance

x x x -s s s eff = x x x -s s s √ 1 -D 0 (2.18)
between two points increases, as expected, with damage. Eq. (2.17) defines in an equivalent manner an effective internal length

lc = l c 1 -D 0 (2.19)
which tends to zero when D 0 tends to unity in accordance with Pijaudier-Cabot et al. results [START_REF] Pijaudier-Cabot | chapter in Nonlocal Modeling of Materials Failure[END_REF] of a material behavior becoming local on free edges (at least in the direction normal to the edge). In 1D and in the non uniform case, the effective distance is defined as the integral over the path [s, x] as

x s (1 -D(x )) -1/2 dx , in 3D as c τ xs Chapter 3

Equivalence between a crack and a damaged zone

In order to illustrate the formulation ability, consider the specimens of figure 3.1: (a) a square plate, (b) a notched specimen and (c) a specimen with a damaged zone at D = D 0 . Vertical z-acceleration fields are also plotted at the same instant for the 3 specimens. It is clearly noticed that the information initiated as an impulse at x x x = x x x 1 takes longer to reach point x x x 6 for both notched and damaged specimens (the wave generated in the example has to turn around the notch). Note the waves reflexion on upper notch and on upper damaged zone sides. When the damage D 0 becomes large (0.99 in the figure), notched and damaged specimens are found equivalent.

Connectivity matrices

To quantify the approach, consider six points x x x p in these specimens. The distances between these points, used in the classical nonlocal theory, are calculated in table 1 for the plate (a) (with l = 0.5 m).

x x xp-x x xq lc

x x x 1 x x x 2 x x x 3 x x x 4 x x x 5 x x x 6 x x x 1 0 l/l c 2l/l c l √ 5/l c l √ 2/l c l/l c x x x 2 l/l c 0 l/l c l √ 2/l c l/l c l √ 2/l c x x x 3 2l/l c l/l c 0 l/l c l √ 2/l c l √ 5/l c x x x 4 l √ 5/l c l √ 2/l c l/l c 0 l/l c 2l/l c x x x 5 l √ 2/l c l/l c l √ 2/l c l/l c 0 l/l c x x x 6 l/l c l √ 2/l c l √ 5/l c 2l/l c l/l c 0 Table 3.1: Connectivity table built from the classical distances x x x -s s s = x x x p -x x x q
The connectivity table is rewritten in a matrix form as a connectivity matrix ( √ 2 and √ 5 replaced by their numerical values to simplify further comparisons),

C C C = l l c •         0 1 2 2.
        (3.1)
The classical connectivity matrix C C C is the same (and is symmetric) for the 3 structures (a), (b) and (c). It does not take into account the presence of notches nor the occurence of damage as

C 16 = C 25 = C 34 = l/l c .
In order to determine the connectivity matrix with the new internal time formulation, one proceeds as follows:

• an impulse force according to the out of plane z-axis (close to a Dirac) is applied at each point x x x = x x x p in a dynamic finite element analysis of a 3D model of a plate (thickness of 0.1 m, 720000 nodes and 360000 elements, free boundary conditions),

• the times τ xs at which the first pic of the z-acceleration information arrives at point s s s = x x x q are recorded. Divided by τ c they are put in the form of connectivity matrices C C C τ ,

• the expressions for τ xs /τ c are synthesized by use of the relation l c = c τ c in the form C C C τ = l lc • A A A and compared to table 1 results for the different specimens. For the concrete square plate and notched specimen one ends up to (E = 35000 MPa, ν = 0.2, ρ = 2400 kg/m 3 ):

C C C (a) τ = l l c •         0 1 2.008 2.233 1.420 1 1 0 1 1.420 1 1.420 2.008 1 0 1 1.420 2.233 2.233 1.420 1 0 1 2.008 1.420 1 1.420 1 0 1 1 1.420 2.233 2.008 1 0         (3.2) C C C (b) τ = l l c •         0 1 2.
        (3.3)
where the notch presence is taken into account (boxed terms), leading for the notched specimen to C (b)

τ 16 = τ (b) 16 /τ c = 1.494 l/l c intead of C (a) τ 16 = τ (a)
16 /τ c = l/l c for the plate with no notch.

For the specimen with the damaged zone at D = D 0 = 0.99, the connectivity matrix corresponding to the nonlocal internal time analysis reads

C C C (c) τ = l l c •         0 1 2.
        (3.4)
Again ony the terms 16 and 61 of the connectivity matrices are changed and found close to the value 1.494 l/l c obtained with reals notches: the damaged zone behaves as a notch, damage and notch being both taken into account by the proposed nonlocal with internal time analysis.

Straight 1D wave propagation

For a better understanding, focus on times τ 16 , τ 25 , τ 34 wave propagation in specimen (c) from x x x 1 to x x x 6 , from x x x 2 to x x x 5 , from x x x 3 to x x x 4 . The corresponding distances are equal in the 3 cases (equal to l). Due to the presence of the notch the time τ 16 is of course larger than the time τ 25 correspnding to a straight path with no notch (and τ 25 /τ c = l/l c ).

A wave propagating along the path x x x 1 -x x x 6 in specimen (c) slows down to the celerity c when meeting the damaged zone. If one only consider the straight path x x x 1 -x x x 6 for the wave propagation, one has

τ 16 = l c + e c • 1 - √ 1 -D 0 √ 1 -D 0 or τ 16 τ c = l l c 1 + e l • 1 - √ 1 -D 0 √ 1 -D 0 (3.5)
where e is the thickness of the damaged zone at D = D 0 . The figure 3.2 shows the normalized increase of time τ 34 = τ (c) 34 with respect to the damage value D 0 and for different ratios e/l. The ratio (τ 34 /τ c )/(l/l c ) in case (c) is equal (resp. close) to unity for a zero (resp. small) damage, the nonlocal theory with internal time recovering then the classical nonlocal theory. The very large increase obtained for large values of the damage, enhanced by a large damaged zone thickness, proves that both a real crack and a highly damaged zone are equivalent in the proposed nonlocal framework. This property is emphasized when the ratio of the weight functions ψ(τ 34 /τ c )/ψ(l/l c ) is drawn (figures 3.3 to fig: ratioexp02) for l/lc = 5, 1 and 0.2 and for the Gaussian weight function ψ(ξ) = e -1 2 ξ 2 : the nonlocal spatial interaction between points across a damaged zone strongly diminishes with damage increase, the larger the points distance l with respect to the internal length l c , the stronger the effect.

Chapter 4

Bar with a damaged zone The Gaussian weight function ψ(ξ) = e -1 2 ξ 2 is considered. The normalizing factors V r (x) are compared for both analyses in figure 4.1a where D 0 = 0.99, l = 2.5 cm, L = 50 cm and where the characteristic length is taken as l c = 10 cm (twice the size 2l of the damaged zone). The classical normalization does not "see" the damaged zone and averages across it when the new nonlocal with internal time approach behaves for the undamaged domains almost as for two independent bars, as expected. The nonlocal strains obtained with both approaches are compared in figure 4.1b. In this particular piecewise constant strain field, the formulation with internal time gives as expected a nonlocal strain field closer to the local strain field than the classical nonlocal approach for which too much importance is given to points over the damaged zone when performing the integral (4.1). Note that in case of structural failure, strain localization leads to non homogeneous fields. The proposed nonlocal averaging then fully acts and makes the solution regular.

= 1 V r L -L ψ • (s)ds V r (x) = L -L ψ ds (x)

Conclusion

A new nonlocal integral formulation is proposed. An internal time is introduced leading to the equivalent definition of an effective or "dynamic" distance and of an evolving internal length. The important feature is that the distance bewteen points is not the argument of the averaging weight function anymore. It is replaced by the information or wave time propagation between these points. The nonlocal weight function build is then symmetric, even in non homogeneous bodies. Pre-computations of wave propagation in the considered structure allow to build the corresponding nonlocal connectivity matrix, with of course the open question of the wave type.

The cracks and notches presence -and if necessary their closure -are naturally taken into account within the wave propagation study and Pijaudier-Cabot and Bazant nonlocal theory is recovered far from the boundaries. The proposed approach makes equivalent a crack and a highly damaged zone, as points across a notch have a small contribution in the nonlocal averaging.

Important point, no assumption on the medium isotropy is made. The proposed nonlocal framework includes anisotropy, either initial or induced.

Last, when no damage is considered as in plasticity with negative hardening models, the internal time concept may still be used, for example by making nonlocal -through Eq. (2.3) -the accumulated plastic strain. The presence of existing notches and cracks will be naturally taken into account if propagation of elastic waves defines the time τ xs but the equivalence between a localized zone and a crack will not be obtained. To gain this last feature, two possibilities are:

• to consider an elasticity coupled with damage even if the plasticity part of the model remain not affected by damage,

• to define τ xs from plastic waves propagation (no need of damage then). 29

FFigure 1 . 1 :Figure 1 . 2 :

 1112 Figure 1.1: Impact on a square plate.

Figure 1 .

 1 Figure 1.3 shows the specimen used to model the wave propagation in the plain plate. Figure1.4 shows the map of the z-acceleration in the plate at a time t=100 ms. Notice that one can see the wave reflection on the free faces of the plate (up and left).

Figure 1 . 3 :

 13 Figure 1.3: Plain specimen.

Figure 1 . 4 :

 14 Figure 1.4: Wave propagation in the plain specimen.

Figure 1 .

 1 Figure 1.5 shows the specimen used to model the wave propagation in the notched plate.The length of the notch is 0.5 m, its width 0.1 m. Figure1.6 shows the map of the z-acceleration in the plate at t=100 ms. Notice again that one can see the wave reflection on the free faces of the plate and in particular the reflection due to the notch (upper lip).

Figure 1 . 5 :

 15 Figure 1.5: Notched specimen.

Figure 1 . 6 :

 16 Figure 1.6: Wave propagation in the notched specimen.

Figure 1 .

 1 Figure1.7 shows the specimen used to model the wave propagation in a plate where previous notch is replaced by a damaged material at a damage level D=0.99. In the computation one uses for the damaged zone a constant effective Young's modulus Ẽ = E(1-D) instead of the initial modulus E. This means here Ẽ = E/100. Figure1.8 shows the map of the z-acceleration in the plate at t=100 ms. One sees the wave reflection on the free edges of the plate as in the previous cases. It is important to see that in this case, the damaged zone acts almost as a free surface for the wave propagation. This is due to the large differences in impedance between the initial undamaged material and the

Figure 1 . 7 :

 17 Figure 1.7: Damaged specimen.

Figure 1 . 8 :

 18 Figure 1.8: Wave propagation in the "damaged" specimen.

Figure 2 . 1 :

 21 Figure 2.1: Notched and damaged specimens.

Figure 3 . 1 :

 31 Figure 3.1: Square plate, notched and damaged specimens. Geometry and pictures of wave propagation at the same instant.

Figure 3 . 2 :Figure 3 2 cFigure 3 2 cFigure 3 2 c

 32323232 Figure 3.2: Ratio τ 34 τc / l lc versus damage D 0 .

  ) with as strain (x) = σ/E if |x| > l, (x) = σ/E(1 -D 0 ) else and where ψ = ψ |x-s| lc for the classical nonlocal analysis, ψ = ψ τxs τc for the nonlocal analysis with internal time.

Figure 4 . 1 :

 41 Figure 4.1: a) Comparison of normalizing factors V r (x) ; b) Comparison of normalized nonlocal strains (lengths in cm, black: nonlocal with internal time, grey: classical nonlocal, dot: local strain).

  Introduction of an internal timein nonlocal integral theories .

	Internal report LMT n˚268
	ISBN : 978-2-11-097305-4
	R. Desmorat and F. Gatuingt
	LMT-Cachan
	ENS Cachan/CNRS/Université Paris 6/PRES Universud Paris
	61, av. du Président Wilson
	F-94230 Cachan, France
	2007

  As an example emphasizing how a damaged zone is taken into account, consider a bar in tension at a stress level σ. The length of the bar is 2L, a damaged zone at large D = D 0 corresponds to x ∈ [-l, l]. The non local analyses perform the integrals (2.1) over the whole bar,

	nl

X 1 X 2 X 3 X 6 X 5 X 4 Specimen (a) X 1 X 2 X 3 X 6 X 5 X 4 Specimen (b) X 1 X 2 X 3