Mireille El Gheche
email: el.gheche@univ-mlv.fr

Jean-Christophe Pesquet
email: pesquet@univ-mlv.fr

Joumana Farah
email: joumanafarah@usek.edu.lb

Caroline Chaux
email: chaux@univ-mlv.fr

Béatrice Pesquet-Popescu

DISPARITY MAP ESTIMATION UNDER CONVEX CONSTRAINTS USING PROXIMAL ALGORITHMS

Keywords: Proximity operator, total variation, tight frame, proximal algorithms, convex optimization, ℓ p -norm, Kullback-Leibler divergence

In this paper, we propose a new approach for estimating depth maps of stereo images which are prone to various types of noise. This method, based on a parallel proximal algorithm, gives a great flexibility in the choice of the constrained criterion to be minimized, thus allowing us to take into account different types of noise distributions. Our main objective is to present an iterative estimation method based on recent convex optimization algorithms and proximal tools. Results for several error measures demonstrate the effectiveness and robustness of the proposed method for disparity map estimation even in the presence of perturbations.

INTRODUCTION

Given any two images of the same scene acquired by stereoscopic cameras, it can be seen that a degree of similarity exists between the two views. Indeed, the pixel at position x in the left image I L corresponds to a pixel at position xu in the right image I R : the disparity of those pixels is denoted by u. The matching problem amounts to searching for the disparity field u which minimizes an error measure. Consequently, we can express the stereo matching problem as :

minimize u J(u), (1)
where,

J(u) = (x,y)∈D φ(I L (x, y) -I R (x -u(x, y), y)) (2)
and φ is assumed to belong to Γ 0 (R) which is the class of a proper lower-semi continuous convex function from R to ] -∞, +∞]. D ⊂ Z 2 is the considered finite image domain.

The recovery of the depth estimation has been increasingly used in a variety of fields such as movies, web networks, 3D reconstruction and computer games. During the last decades, various methods have been introduced in disparity estimation, especially in computer vision [START_REF] Scharstein | A taxonomy and evaluation of dense two-frame stereo correspondence algorithms[END_REF]. Stereo matching algorithms are generally classified into two categories: feature based and area based ones. Feature based methods match the feature elements between two images [START_REF] Medioni | Segment-based stereo matching[END_REF], while area based methods perform matching between pixels. Currently, the associated problems are often solved using discrete optimization techniques based on a global approach. Algorithms such as dynamic programming [START_REF] Veksler | Stereo correspondence by dynamic programming on a tree[END_REF], graph cuts [START_REF] Kolmogorov | Computing visual correspondence with occlusions using graph cuts[END_REF] and variational methods [START_REF] Miled | A convex optimisation approach for depth estimation under illumination variation[END_REF] were proposed. Note that variational-based disparity estimation methods demonstrated excellent performance compared with the state-of-the-art.

In [START_REF] Miled | A convex optimisation approach for depth estimation under illumination variation[END_REF], a convex energy function approximation of J is derived and minimized subject to various convex constraints (modeled by convex sets (C i ) 1≤i≤m ) arising from prior knowledges and observed data. The optimization problem (1) to be solved then becomes minimize Liu∈Ci,i∈{1,...,m}

J(u), (3)
where (L i ) 1≤i≤m are linear operators. However, the employed algorithm requires not only the strict convexity of the criterion to be considered but a quadratic form of it. This paper describes a new DDE (Dense Disparity Estimation) approach that generalizes the convex optimization approximation presented in [START_REF] Miled | A convex optimisation approach for depth estimation under illumination variation[END_REF]. The proposed method is based on some proximity operator (prox) properties and then, it is no longer limited to strictly convex quadratic criteria. Therefore, explicit expressions of the prox of some norms such as ℓ 1 -norm and more generally ℓ p -norm with p > 1 or the Kullback-Leibler divergence allow us to apply such proximal algorithms to our disparity map estimation problem. By allowing a rich choice of distance functions, the proposed technique is well-suited for dealing with different types of noise corrupting the observed data, such as Poisson noise, salt and pepper noise or Gaussian noise.

The remainder of the paper is organized as follows: Section 2 describes the basis tools employed in our approach. Section 3 lays out the problem statement. Results on the Middlebury dataset 1 are provided in Section 4. Finally, some conclusions are drawn in Section 5.

BACKGROUND

We will now present some tools which are useful in the solution of our minimization problem. In what follows, we will focus on the definition of the proximity operators, total variation and tight frames.

Proximity operator

The proximity operator of a convex function is a natural extension of the projection operator P C onto a nonempty closed convex set C ⊂ R N [START_REF] Combettes | The foundations of set theoretic estimation[END_REF][START_REF] Combettes | Proximal splitting methods in signal processing[END_REF]8]. The projection P C y of a point y ∈ R N onto C is the solution to the problem:

minimize x∈R N ι C (x) + 1 2 x -y 2 (4)
where

ι C ∈ Γ 0 (R N ) is the indicator function of C.
In a seminal paper [START_REF] Moreau | Fonctions convexes duales et points proximaux dans un espace hilbertien[END_REF], Moreau proposed the following extension of the notion of projection by replacing the function ι C by an arbitrary function f ∈ Γ 0 (R N ). Then, the problem can be rewritten as

minimize x∈R N f (x) + 1 2 x -y 2 . ( 5
)
This problem admits a unique solution which is the proximity operator prox f y of f at y. Proximity operators have very attractive properties [START_REF] Combettes | Proximal splitting methods in signal processing[END_REF] that make them particularly well-suited to design iterative minimization algorithms. In Table 1, we summarize some of them that will be needed to obtain closed form expressions of the required proximity operators.

Total variation

Let H and V be two spatial convolution operators corresponding respectively to the computation of horizontal and vertical discrete gradients. Then, a discrete version of the total variation (TV) for which proximity tools can be applied is the following [START_REF] Combettes | A proximal decomposition method for solving convex variational inverse problems[END_REF][START_REF] Pustelnik | Parallel ProXimal algorithm for image restoration using hybrid regularization[END_REF]:

TV(u) = (x,y)∈D |Hu(x, y)| 2 + |V u(x, y)| 2 . ( 6
)
It should be noticed that, in our simulations, we apply periodic convolutions, which means that the operators H and V are diagonalized by a Fast Fourier Transform.

1 http://cat.middlebury.edu/stereo/newdata.html

f (x) prox f x φ(x -z), z ∈ R N z + prox φ (x -z) φ(x/ρ), ρ ∈ R \ {0} ρ prox φ/ρ 2 (x/ρ) φ(Lx), semi-orthogonal, x + ν -1 L ⊤ (prox νφ (Lx) -Lx) L ∈ R M ×N , LL ⊤ = νI, ν > 0 ι C (x) = 0 if x ∈ C +∞ if x / ∈ C P C x χ > 0, α > 0 -χ ln(x) + αx, if x > 0 +∞ if x ≤ 0 x -α + |x -α| 2 + 4χ 2
Table 1. Some proximity operator properties (I denotes the identity matrix) [START_REF] Combettes | Proximal splitting methods in signal processing[END_REF].

Tight frame

Frame representations [START_REF] Han | Frames, bases, and group representations[END_REF] and more precisely tight frame representations have been widely used during the last decade. Such transforms can be described by an analysis frame operator F and a synthesis frame operator F ⊤ . F is said to be a tight frame when F ⊤ F = νI, where ν > 0.

A simple example of a tight frame is the union of m orthogonal wavelet bases [START_REF] Mallat | A wavelet tour of signal processing[END_REF], thus leading to a tight frame representation with ν = m.

In this paper, we will restrict our attention to the union of 4 shifted orthonormal Haar bases over 1 resolution level.

PPXA+ Algorithm

Many algorithms have been formulated in the literature to solve optimization problems like (1). In our case, we use the Parallel Proximal Algorithm (PPXA+) [START_REF] Pesquet | A parallel inertial proximal optimization method[END_REF] which is a flexible tool. PPXA+ (algorithm 1) allows us to minimize a convex criterion J on some closed convex constraint sets (C i ) 1≤i≤m . It consists of computing, in parallel, the projections onto the different convex sets (C i ) 1≤i≤m and the proximity operator of the criterion J. Suppose that the following assumptions hold.

• m i=1 L ⊤ i L i + I is an invertible matrix. • ∃ ǔ ∈ R K such that (∀i ∈ {1, . . . , m}) L i ǔ ∈ ri(C i ) and ǔ ∈ ri(dom J). 2 • (λ n ) n∈N is a sequence of relaxation parameters such that (∀n ∈ N) λ ≤ λ n+1 ≤ λ n < 2, where λ ∈]0, 2[.
Then, the sequence (u n ) n∈N generated by Algorithm 1 converges to a solution to Problem (3), provided that such a solution exists.

[Initialization] (w 1 , ..., wm

) ∈]0, +∞[ m ,γ > 0 (z i,0 ) 1≤i≤m+1 ∈ R N 1 × . . . × R Nm × R K Q = ( m i=1 w i L ⊤ i L i + γI) -1 u 0 = Q( m i=1 w i L ⊤ i z i,0 + γz m+1,0 ) For n = 0, 1, . . . do For i = 1, . . . , m do p i,n = P C i (z i,n ) end For p m+1,n = prox J γ (z m+1,n ) cn = Q( m i=1 w i L ⊤ i p i,n + γp m+1,n ) For i = 1, . . . , m do z i,n+1 = z i,n + λn(L i (2cn -un) -p i,n ) end For z m+1,n+1 = z m+1,n + λn(2cn -un -p m+1,n ) u n+1 = un + λn(cn -un) end For
Algorithm 1: Constrained version of PPXA+.

APPLICATION TO THE DISPARITY ESTIMATION PROBLEM

Problem statement

In the class of convex optimization methods in stereo vision, an interesting approach that was shown to be competitive with other state-of-the-art methods, was proposed in [START_REF] Miled | A convex optimisation approach for depth estimation under illumination variation[END_REF]. It is based on a subgradient projection method. In what follows, we will adopt the same problem formulation and then we will solve it using the tools described in Section 2. The function J defined in (2) is nonconvex despite the convexity of the function φ. To alleviate this problem, we propose to perform a Taylor expansion of the nonlinear term I R (x-u(x, y), y) around an initial estimate ū(x, y). We have then:

I R (x -u(x, y), y) ≃ I R (x -ū(x, y), y) -(u(x, y) -ū(x, y)) I x R (x -ū(x, y), y) , (7)
where I x R is the horizontal gradient of the right image. Based on (2) and ( 7), we deduce that:

J(u) = (x,y)∈D φ(T (x, y) u(x, y) -r(x, y)) (8)
where

T (x, y) = I x R (x -ū(x, y), y), (9)
r(x, y) = I R (x -ū(x, y), y) + ū(x, y) T (x, y) -I L (x, y).

(

) 10
Since the occluded image areas (denoted by O) yield very large disparity errors, they will be discarded in the expression of the similarity criterion. This is performed by rewriting J as:

J(u) = (x,y)∈D\O φ(T (x, y) u(x, y) -r(x, y)) (11)
The proposed PPXA+ algorithm provides an efficient solution to solve this problem. Our main contribution here lies in the various choices that can be made for the error measure. As already mentioned, the algorithm is based on iterating computations of the proximity operator of J and projections onto convex sets related to constraints. We will first give the explicit forms of the proximity operators of interest in this work which correspond to error measures J taken as ℓ p -norms or Kullback-Leibler divergence.

Error measures

ℓ p -norm

Based on Table 1, it can be deduced that

u(s) =prox γ -1 φ(T (s)•-r(s)) ū(s) = 1 T (s) prox γ -1 |T (s)| 2 φ (T (s)ū(s) -r(s)) + r(s)
where s = (x, y). Then the proximity operator of

γ -1 |T (s)| 2 φ with φ = | • | p and ξ(s) = T (s)ū(s) -r(s), is [15]: prox γ -1 |T (s)| 2 φ(•) ( ξ(s)) =                    sign( ξ(s)) max( ξ(s) -γ -1 |T (s)| 2 , 0), if p = 1; ξ(s) 2γ -1 |T (s)| 2 +1 , if p = 2; sign( ξ(s)) √ 1+12γ -1 |T (s)| 2 | ξ(s)|-1 6γ -1 |T (s)| 2 , if p = 3; ( µ+ ξ(s) 8γ -1 |T (s)| 2 ) 1 3 -( µ-ξ(s) 8γ -1 |T (s)| 2 ) 1 3 , if p = 4, where µ = ξ(s) 2 + 1 27γ -1 |T (s)| 2 . ( 12
)

Kullback-Leibler divergence

One can note that ( 7) can be rewritten as

I R (x -u(s), y) = r(s) + I L (s) -u(s)T (s) = r(s) -u(s)T (s) = ζ(s). (13)
Then, the criterion to be minimized takes a form slightly different from [START_REF] Pustelnik | Parallel ProXimal algorithm for image restoration using hybrid regularization[END_REF]:

J(u) = (x,y)∈D\O Φ(I L (s), ζ(s)) (14)
where Φ(I L (s), ζ(s)) is the Kullback-Leibler divergence function given by:

Φ(I L (s), ζ(s)) =        IL(s) ln I L (s) ζ(s) + ζ(s) -IL(s), if (IL(s), ζ(s)) ∈ ]0, +∞[ 2 ζ(s), if IL(s) = 0, ζ(s) ≥ 0 +∞, otherwise. ( 15
)
The properties in Table 1 allow us to re-express the proximity operator u(s) at ū(s) of the function γ -1 Φ(I L (s), r(s) -T (s) •) as:

(∀s ∈ D \ O) u(s) = 1 T (s) r(s) -prox γ -1 |T (s)| 2 Φ(IL(s),•) ζ(s) where ζ(s) = r(s) -ū(s)T (s). If I L (s) ∈ ]0, +∞[, then prox γ -1 |T (s)| 2 Φ(IL(s),•) ( ζ(s)) = prox γ -1 |T (s)| 2 (-IL(s) ln(•)+•) ( ζ(s)) = ζ(s) -γ -1 |T (s)| 2 2 + | ζ(s) -γ -1 |T (s)| 2 | 2 + 4γ -1 |T (s)| 2 I L (s) 2 .
Furthermore, it can be shown that the above expression remains valid if I L (s) = 0. Let us now turn our attention to the convex constraints that can be introduced.

Convex constraints

The motivation behind introducing these constraints is to incorporate additional prior knowledge about the solution.

First, we will constrain the range values of the disparity u by setting a first constraint:

S 1 = {u | u min ≤ u ≤ u max } . (16)
Then, we can impose an upper bound on the total variation measure. Indeed, this constraint allows us to recover piecewise and homogeneous areas with sharp edges:

S 2 = {u | TV(u) ≤ τ u } . (17)
Finally, one can introduce some prior information on the wavelet coefficients of the disparity. Based on Section 2.3, we can represent the associated constraint in terms of Besov space B 1 1,1 semi-norm as:

S 3 = u | j≥1,k∈Z 2 ,o∈{H,V } |c B j,k,o (u)| ≤ κ . (18)
One can link these constraints (S i ) 1≤i≤3 to some convex sets

(C i ) 1≤i≤3 by noticing that S 1 = C 1 , S 2 = L -1 2 (C 2 ) where L 2 = H V , and
S 3 = L -1 3 (C 3 ) with L 3 = F .
The projection onto C 1 is straightforward. The projections onto the convex sets C 2 and C 3 were implemented using the iterative projection algorithm proposed in [START_REF] Van Den | Probing the pareto frontier for basis pursuit solutions[END_REF].

SIMULATION RESULTS

The objective of this section is to present numerical results obtained with the proposed method for a variety of synthetic stereo noisy images. In [START_REF] Gheche | Proximal splitting methods for depth estimation[END_REF], we provided some comparisons with the state-of-the-art DDE method presented in [START_REF] Miled | A convex optimisation approach for depth estimation under illumination variation[END_REF]. We briefly report the results of the application of our approach to the standard Middlebury stereo dataset. The objective of this first experiment is to demonstrate the validity of the proposed method without any perturbation. Fig. 1 shows the stereo image pairs considered in this work, namely Corridor, Saw, Teddy and Cones, along with their ground truth images. Our method allows us to compute a smooth disparity map with accurate depth discontinuities as also shown by the preliminary results in [START_REF] Gheche | Proximal splitting methods for depth estimation[END_REF]. We specify in each case the considered error measure and constraints (the choice was made according to the best performances). We evaluate our method using two error measures between the computed field u c and the ground truth field u ref :

SNR = 10 log 10 s∈D |u ref (s)| 2 s∈D |u c (s) -u ref (s)| 2 , (19) MAE = 1 N s∈D |u c (s) -u ref (s)|, ( 20
)
where N is the number of pixels.

In the presence of noise, we evaluate our method using different kinds of perturbation. To introduce a significant noise variation, we modified the left and the right images by using the imnoise function in Matlab. Fig. 2 shows the computed disparity maps from the noisy stereo pairs. As we can see, our method is less affected than DDE by noise changes and provides accurate depth maps. We display in Fig. 2 results obtained with the proposed error measures (ℓ 1 -norm for salt and paper noise, Kullbackleibler divergence for Poisson noise and quadratic norm for Gaussian noise). Different constraints were also taken into account: we chose the one that leads to the best results either for DDE or PPXA+. As shown in the last column, our approach outperforms the DDE method [START_REF] Miled | A convex optimisation approach for depth estimation under illumination variation[END_REF], which proves the effectiveness of our method even in the presence of Gaussian noise. Furthermore, as we can see in the first and second columns, numerical values and visual comparisons on the obtained disparity maps confirm that our PPXA+ algorithm produces good estimates in the presence of other kinds of noise. Finally, Fig. 3 displays the MAE with respect to the Poisson noise intensity α (a smaller α corresponding to a higher noise). PPXA+ algorithm outperforms the DDE approach especially when the noise intensity is high.

CONCLUSIONS

An efficient proximal method that deals with disparity estimation problems has been proposed in this paper. Being very flexible, it allows us to minimize various criteria (not necessarily differentiable nor strictly convex), which can be useful when perturbations are present during image acquisition. The effectiveness of the proposed approach was demonstrated on different stereo image pairs, while taking into account different kinds of perturbations. In the future, we plan to extend our approach to the case where illumination variations exist.

Fig. 1 .

 1 Fig. 1. From top to bottom: Corridor, Saw, Teddy, Cones. From left to right: left reference images, ground truth images, PPXA+ results.

Fig. 2 .

 2 Fig. 2. From left to right: Salt and paper noise, Poisson noise, Gaussian noise. From top to bottom: noisy left image, DDE results, PPXA+ results.

Fig. 3 .

 3 Fig. 3. MAE with respect to Poisson noise intensity α applied to Cones pair.

dom J is the domain of J and the relative interior of a set C is denoted by riC.


Thomas Moulard
email: thomas.moulard@gmail.com

Pablo F Alcantarilla

Florent Lamiraux

Olivier Stasse

Frank Dellaert
email: dellaert@cc.gatech.edu

Reliable Indoor Navigation on Humanoid Robots using Vision-Based Localization

Keywords: Humanoid robots, Localization, Motion Planning, Robot Control, Stereo Vision 1

Reliable localization on a humanoid robot is a sine qua non condition to succeed in realizing complex robotics scenarios. Before dealing with perturbations and online modification of the environment, one has to make sure that the planned trajectory alone will be correctly followed. This paper demonstrates that a control framework suited to humanoid robots relying on a vision-based localization system can achieve this goal. Our localization framework is based on a real-time vision-based localization system that assumes that a pre-existing 3D map of the environment exists and allows to obtain accurate results in complex robotics scenarios. By compensating for execution errors such as drifts and robot model errors, the HRP-2 robot is able to achieve high precision tasks.

I. INTRODUCTION

In this paper, we consider the problem of real-time reliable trajectory execution for humanoid robots. A precise localization on a humanoid robot is both crucial to achieve reliable trajectory execution and challenging due to humanoid robots specific features. We propose a novel control architecture for humanoid robots where the current localization is used by the robot controller and planner trajectory. The objective of this work is to localize the humanoid robot HRP-2 [START_REF] Kaneko | Humanoid robot HRP-2[END_REF] while it is navigating through an indoor environment.

Localization is crucial on a humanoid robot as its position cannot be controlled directly. Legs motion are planned under the assumption that all contacts will be perfect (i.e. no friction) during the movement. This is not the case in practice and leads to execution errors which cannot be ignored. Moreover, reactive control algorithms rely on simplified dynamics models such as the Linear Inverted Pendulum model [START_REF] Kajita | The 3D linear inverted pendulum mode: a simple modeling for a biped walking pattern generation[END_REF]. These imply a gap between the generated motion and the executed one. Therefore, localization is necessary to ensure a consistent robot behavior.

However, localization on a humanoid robot is challenging. Mobile robots odometry can be computed using wheel encoders and give a reasonable hint on the current robot motion. 2D maps can also be used to simplify the navigation through indoor environments. This is not the case on a legged robot on which 3D localization is required and where no encoder based reliable odometry exists. All these reasons make localization a cornerstone of reliable trajectory execution on humanoid robots. To obtain an accurate localization of the robot in the environment, we employ fast vision-based localization techniques [START_REF] Alcantarilla | How to localize humanoids with a single camera? Autonomous Robots[END_REF] that assume that a prior 3D map of the environment exists. By means of visibility prediction [START_REF] Alcantarilla | Visibility Learning in Large-Scale Urban Environment[END_REF] a fast and robust data association between a large map of 3D points and perceived 2D features in the image can be performed very efficiently yielding an accurate pose estimate of the robot in the 3D map. This prior 3D map of the environment can be generated by standard visual Simultaneous Localization and Mapping (SLAM) techniques [START_REF] Davison | MonoSLAM: Real-Time Single Camera SLAM[END_REF][START_REF] Konolige | FrameSLAM: from bundle adjustment to real-time visual mapping[END_REF]. In this article, the past results in vision processing on humanoids robots are presented in Section II. Section III details how the vision is linked to the robot controller via our proposed control scheme. The vision-based localization algorithm is described in Section IV. The experimental results are discussed and compared to motion capture data in Section V. Finally, conclusions and future work are described in Section VI.

II. RELATED WORK

Cameras seem to be an appealing sensor for humanoid robotics applications: they are small, cheap and light-weight compared to other sensors such as laser scanners. Despite of this, there have been only limited attempts at visionbased localization for humanoid robots. Ozawa et al. [START_REF] Ozawa | Using visual odometry to create 3D maps for online footstep planning[END_REF] proposed to use stereo visual odometry to create local 3D maps for online footstep planning. The main drawback of this approach is the drift created by the accumulation of errors that typically occur in visual odometry systems [START_REF] Kaess | Flow separation for fast and robust stereo odometry[END_REF].

Other past experiments where vision algorithms have been used to construct navigation maps for humanoid robots are [START_REF] Michel | Vision-guided humanoid f ootstep planning for dynamic environments[END_REF][START_REF] Michel | Online environment reconstruction for biped navigation[END_REF]. However, most of these works were focused on the humanoid trajectory generation replanning problem. In this paper, our approach is different: the objective is to execute a planned trajectory without any online replanning. The computation of a whole-body trajectory remains very costly and do not meet the real-time requirements of small robots such as Nao 1 . Therefore, local trajectory deformation is important to achieve high reactivity and to avoid unnecessary computations.

The work by Dune et al. [START_REF] Dune | Cancelling the sway motion of dynamic walking in visual servoing[END_REF] relies on visual servoing to compute a center of mass reference velocity to control the walking trajectory generation algorithm described in [START_REF] Herdt | Online Walking Motion Generation with Automatic Foot Step Placement[END_REF]. By computing the foot placement online, one can control directly the robot center of mass velocity on a 2D plane and get rid of most of the complexity of the walking process. This approach is extremely interesting, but unfortunately does not take into account obstacles in the environment. Also, this algorithm is not suited for trajectory tracking which is our objective. The following works [START_REF] Harada | An analytical method on real-time gait planning for a humanoid robot[END_REF][START_REF] Morisawa | Experimentation of humanoid walking allowing immediate modification of foot place based on analytical solution[END_REF] aim at allowing sudden changes in the robot trajectory. Again, the objectives pursued by these works and our work differ as they alter the initial plan to react to changes in the environment.

Davison et al. [START_REF] Davison | MonoSLAM: Real-Time Single Camera SLAM[END_REF] showed successful monocular SLAM results for small indoor environments using the HRP-2 robot. This approach, known as MonoSLAM, is a monocular Extended Kalman Filter (EKF) vision-based system, that allows building a small map of sparse 3D points. However, acceptable results were only obtained when the pattern generator, the robot odometry and inertial sensing were fused to aid the visual mapping into the EKF framework as it was shown in [START_REF] Stasse | Real-time 3D SLAM for a humanoid robot considering pattern generator information[END_REF]. The fusion of the information from different sensors can reduce considerably the uncertainty in the camera pose and the 3D map points involved in the EKF process, yielding better localization and mapping results.

Despite of this, EKF-based approaches have important drawbacks such as the limited number of 3D points that can be tracked and divergence from the true solution due to linearization errors. As shown in [START_REF] Strasdat | Realtime Monocular SLAM: Why Filter[END_REF], nonlinear optimization techniques such as bundle adjustment [START_REF] Mouragnon | Generic and Real-Time Structure from Motion using Local Bundle Adjustment[END_REF] are superior in terms of accuracy to filtering based methods, and allow tracking many hundreds of features between frames. In this paper we use the framework described in [START_REF] Alcantarilla | How to localize humanoids with a single camera? Autonomous Robots[END_REF], where a stereo visual SLAM algorithm based on local bundle adjustment is 1 http://www.aldebaran-robotics.com/en/ used to compute a 3D map of the robot's environment. Then, we use the prior 3D map for an efficient data association obtaining a real-time accurate localization of the robot within the 3D environment.

III. ROBOTICS FRAMEWORK FOR RELIABLE TRAJECTORY EXECUTION

A. Notations

Before describing the robot architecture, let us introduce the following notations: SE(2) and SE(3) the Special Euclidian groups of dimension 2 and 3 denote rigid transformations in the 2D and 3D Euclidian spaces respectively. q(t) denotes the robot configuration in its configuration space C at time t. It is divided into a controllable part q(t) ∈ C and the robot position x(t) ∈ SE(3). The robot position is defined as the position of a particular robot body. In this paper, we choose the left ankle as the reference body used to compute the robot position. q(t) describes the joint velocities. 3) is the robot trajectory as perceived by the localization system. c(t) ∈ R 2 and z(t) ∈ R 2 are respectively the projection of the robot center of mass on the floor (z = 0) and the Zero Momentum Point position at time t. γ la (t) ∈ SE(3), γ ra (t) ∈ SE(3) denote respectively left ankle and right ankle trajectories. θ(t) ∈ SE( 6)is the estimated camera pose (translation and rotation) at time t.

x ref (t) ∈ SE(3) is the robot planned trajectory. x(t) ∈ SE(
The set of trajectories γ la , γ ra and c defines a complete balanced walking movement and is denoted Γ. It can be completed by an additional trajectory which defines the upper-body configuration during the walk. Changes in the upper-body posture do not impact the walking movement as long as the center of mass trajectory is unchanged and the induced angular momenta variations remain small. Therefore, this article will not mention upper-body trajectory although the final experiment contains arms and head movement.

B. Robotics Architecture

The proposed architecture can be divided in two parts:

1) The vision processing and localization components which are responsible for acquiring data from the robot's stereo rig system, computing an estimation of the camera pose and deducing the robot localization. 2) The control system which follows the planned trajectory while closing the loop on the localization data.

These two systems are running on two different computers embedded on the robot. They are composed of several robotics components running simultaneously. Fig. 2 illustrates the complete robotics infrastructure. To communicate, we rely heavily on both OpenHRP and ROS middlewares. b) Error Estimation: The error estimation component computes δx ∈ SE(3), the transformation from the perceived robot position to the planned robot position and is defined at time t as:

δx t = x t .x -1 t (1)
However, x-1 t is not directly provided by the localization which estimates the camera pose and not the robot pose. The transformation from the camera to the left ankle (defining the robot position) can be easily deduced from the control and is not subject to significant execution errors as joint encoders allow low-level reliable servoing.

Additionally, the control is used to help the camera pose estimation. Indeed, to estimate a 3D pose, it is necessary to estimate six parameters; three for the translation: x, y, z and three for the rotation: roll, pitch, and yaw. Nevertheless, we know for sure that the camera height z and two of the three parameters of the rotation cannot drift nor have significant execution errors while walking on a flat ground. This explains why, in practice, only three parameters are estimated using the vision components whereas the z, roll and pitch parameters are known from the plan. c) Closed-Loop Trajectory Following: To achieve walking and compute motor control, the task based control framework described in [START_REF] Mansard | A versatile Generalized Inverted Kinematics implementation for collaborative working humanoid robots: The Stack Of Tasks[END_REF] is used. One interesting feature of this approach is that we can directly add tasks to follow ankles and center of mass reference trajectories. The inverse kinematics computation will then be implicitly realized by the task resolution.

Regarding trajectory generation, we use the planner described in [START_REF] Dalibard | Small-space controllability of a walking humanoid robot[END_REF] to generate the required reference trajectories.

Our control scheme adds a step which processes the reference trajectories to compensate for execution errors based on the error estimation.

Let us consider that at time t, the robot is on double support (i.e. both feet on the ground) and that an estimation of the execution error δx(t) is available. We have to alter the trajectories of the left ankle, right ankle and center of mass simultaneously while preserving robot balance. As we cannot change a foot position while it is on contact with the floor, we have to modify the trajectories during the two following steps to be able to correct both feet trajectories. To do so, we will sum each trajectory with a trajectory dependent third order polynomial denoted by

∆ la ∈ SE(3), ∆ ra ∈ SE(3), ∆ c ∈ R 2 : γ la (t) = ∆ la (t).γ la (t) γ ra (t) = ∆ ra (t).γ ra (t) c (t) = ∆ c (t) + c(t) (2)
Let t 1 be the end of the next step and t 2 the end of the second next step. If the first step is realized with the left foot (i.e. the right foot is the support foot). In Eq. 2, if ∆ provides a correction from t start to t end , we consider the following constraints:

∆(t) = 0., if t≤ t start δx, if t≥ t end ∆(t 1 ) = ∆(t 2 ) = 0 (3)
The constraints expressed in Eq. 3 ensures us that the correction is smooth: zero velocity at the beginning and end of the correction and that the error is correctly compensated. They are independent of the corrected trajectory (ankles or center of mass). These four constraints fully determine the four coefficients of the polynomial. The starting and timing correction time depends on the corrected trajectory: center of mass (∆ c (t)):

t start = t, t end = t 1 left ankle (∆ la (t)): t start = t, t end = t 1 right ankle (∆ ra (t)): t start = t 1 , t end = t 2
The remaining issue is ensuring that the updated trajectory will be balanced. A classical approach is using the Zero Momentum Point (ZMP): this virtual point acts as a criterion to determine whether a trajectory is balanced or not. If it lies in the convex hull of the robot contact points all the time, the trajectory is balanced, otherwise it is not. The ZMP is defined as:

z = c + 1 m(c z + g) 0 -1 0 1 0 0 L - c z cz + g c (4)
where c z is the height of the center of mass with respect to the ground, m is the mass of the robot, g is the gravity constant, c = (c x , c y ) is the projection of the center mass of the robot on the ground and L is the angular momentum of the robot about the center of mass. Solving this equation is time consuming as it is nonlinear and requires dynamics computation. Under the asumption that L and c z is constant in time [START_REF] Kajita | The 3D linear inverted pendulum mode: a simple modeling for a biped walking pattern generation[END_REF], we can obtain the follwing simplified linear model:

z = c - c z g .c (5)
Considering r a polynomial depending only of z, (V x , V y , W x , W y ) free parameters used to constrain the initial position and velocity of the center of mass, a general solution of Eq. 5 is:

c(t) = cosh( g z c .t).V + sinh( g z c .t).W + r(t) (6)
Given the formulation in Eq. ( 6), it is possible to continuously modify the center of mass trajectory to make it follow z (t) the corrected trajectory. This new trajectory can be expressed as the sum of two polynomials:

c (t) = cosh( g z c .t).V + sinh( g z c .t).W + r(t) + ∆(t) (7)
It appears then that moving the center of mass final position of δx is equivalent to moving the ZMP final position of δx then solving the linear differential equation. As the flying foot end position is transformed similarly, the ZMP will stay in the convex hull contact points for any correction value δx. The only limitations are feet velocities limits and assumptions realized when using the simplified linear model.

To correct properly the robot trajectory, a correction of about two to three cm every two steps is sufficient. In this case, the simplified linear model is particularly well suited as the dynamical effects of such changes are so low they can be safely ignored.

When a correction is completely applied, i.e. t > t 2 , a new correction is computed using the current estimation of the execution error.

By combining vision-based localization, error estimation and closed-loop trajectory following, our control framework can reshape walking trajectories on the fly to reliably track the planned trajectory. We will focus the discussion in the next section on how to localize the camera precisely.

IV. VISION PROCESSING

Our computer vision module comprises two different stages: first, a stereo visual SLAM module for building a persistent 3D map of the environment and second a visionbased localization framework with visibility prediction that assumes that a prior 3D map of the environment is given. Prior to any SLAM or localization processing, we correct the distortion of the images and perform stereo rectification. Stereo rectification simplifies considerably the stereo correspondences problem and disparity maps can be easily obtained with the rectified images. More details of the visual SLAM and vision-based localization algorithms can be found in [START_REF] Alcantarilla | How to localize humanoids with a single camera? Autonomous Robots[END_REF]. In this section, we will briefly describe the main components of our localization framework assuming that we already have a prior 3D map of the environment. Given a prior 3D map of the environment we can use the computed map for fast and robust localization. In this context we employ the visibility prediction technique described in [START_REF] Alcantarilla | Visibility Learning in Large-Scale Urban Environment[END_REF][START_REF] Alcantarilla | How to localize humanoids with a single camera? Autonomous Robots[END_REF] to perform an efficient data association between known 3D points and detected 2D features.

Our localization framework is composed of two different modules: initialization and a combination of vision-based localization with visibility prediction and stereo visual odometry.

A. Initialization and Re-Localization

At this stage, the robot is lost and can be located in any area of the map. Therefore, we need to find an initial camera pose to start the vision-based localization algorithm. For this purpose we detect 2D features in the new image using the Harris corner detector [START_REF] Harris | A combined corner and edge detector[END_REF] at different scale levels and compute an appearance descriptor for each detected corner. For this purpose, we compute the appearance descriptors of the detected 2D features in the new image and match this set of descriptors against the set of descriptors from the list of stored keyframes from the prior 3D reconstruction. Similar to Speeded Up Robust Features (SURF) [START_REF] Bay | SURF: Speeded Up Robust Features[END_REF], for a detected feature at a certain scale, we compute a unitary descriptor vector of dimension 16 in order to speed up the descriptor computation. We use the upright version of the descriptors (no invariance to rotation) since upright descriptors perform better in scenarios where the camera only rotates around its vertical axis, which is often the case for humanoid robots, and are also faster to compute than its rotation invariant counterpart.

In the matching process between the camera frame and the list of keyframes, we perform a RANSAC [START_REF] Bolles | A RANSAC-based approach to model fitting and its application to finding cylinders in range data[END_REF] procedure forcing epipolar geometry constraints. We recover the camera pose from the stored keyframe that obtains the highest ratio of inliers. If this ratio is lower than a certain threshold, we do not initialize the localization algorithm until the robot moves into a known area yielding a higher ratio.

B. Localization

Given a prior map of 3D points and perceived 2D features in the image, the problem to solve is the estimation of the camera pose with respect to the world coordinate frame.

Once the system has a good initialization, the vision-based localization system works through the following steps: 1) While the robot is moving, the stereo pair acquires a new set of images which are rectified. Then, from the rectified images, a disparity map is computed. 2) A set of image features Z t = {z t,1 . . . z t,n } is detected by Harris corner detector only for the left image. Then, a feature descriptor is computed for each of the detected features.

3) Then, by using the visibility prediction algorithm, a promising subset of highly visible 3D map points is chosen and re-projected onto the image plane based on the estimated previous camera pose θ t-1 and known camera parameters. 4) Afterwards, a set of putative matches C t is formed where the i-th putative match C t,i is a pair {z t,k , x j } which comprises a detected feature z k and a map element x j . A putative match is created when the Euclidean distance between the appearance descriptors of a detected feature and a re-projected map element is lower than a certain threshold. 5) Finally, we solve the pose estimation problem minimizing the following cost error function, given the set of putative matches C t :

arg min R,t m i=1 z i -K (R • x i + t) 2 (8)
where z i = (u L , v L ) is the 2D image location of a feature in the left camera, x i represents the coordinates of a 3D point in the global coordinate frame, K is the left camera calibration matrix, and R and t are respectively the rotation and the translation of the left camera with respect to the global coordinate frame. The pose estimation problem is formulated as a nonlinear least squares procedure using the Levenberg-Marquardt algorithm. The set of putative matches may contain outliers, therefore RANSAC is used in order to obtain a robust model free of outliers.

There can be some frames where the pose estimation problem cannot be solved efficiently since we may have textureless areas or slightly different viewpoints from the ones captured at the mapping sequence. For those situations, we employ stereo visual odometry to update the pose of the robot with respect to the map coordinate frame. Notice here that our system does not suffer from the typical drift of visual odometry systems, since in the next frame the system will try to localize with respect to the prior 3D map. When the number of consecutive frames where the pose estimation fails is higher than a fixed threshold (e.g. 100 frames), we declare that the tracking is lost and start a re-localization process.

V. EXPERIMENTAL RESULTS

A. Experimental Setup

The experiment demonstrates that by localizing the robot while walking, HRP-2 can reach a goal position independently of the execution errors. In the chosen scenario, the robot must drop a ball on a shelf after walking 2 m. A more precise description of the scenario is depicted by Fig. 4. Empirically, we estimated the HRP-2 mean drift to be around one cm per step. Considering that 32 steps are required to reach the goal without hitting the obstacles, the usual drift would prevent the task from being accomplished. The camera trajectory estimated by the localization algorithm is validated using motion capture data. Markers have been placed on both the robot head and left ankle to provide ground truth.

Finally, the camera position is given by the algorithm described in detail in the previous sections. The localization During the experiments, the robot sometimes failed to drop the ball properly. It was either due to a bad initialization which sometimes lead to localization failure or to robot flexibility. Indeed, the shelves is quite narrow and the robot ankles passive joints were generating hands movements which were causing collisions, even though the robot localization is precise enough. It would be interesting in a future work to not only localize the robot with respect to an absolute frame but also switch to a task-based localization at the end to correct the hand trajectory and avoid collisions.

On our platform, the vision localization module returned a pose at around 16 Hz. On the opposite, the HRP-2 robot uses a control loop at 200 Hz. To preserve real-time, the realtime component was communicating the robot configuration asynchronously at around 75 Hz. The error estimation components was running at a mean rate of 10 Hz. As a correction is added every two steps, it was not useful to evaluate the error at a higher frequency. The mean error is 0.2 m in translation. Drop in the algorithm precision can occur when the robot enters a part of the map which is not dense enough or where the environment does not incorporate enough texture to detect enough interest points. It happens in our scenario near the end when the camera is too close to the shelf to perceive its environment correctly.

B. Comparison to Motion Capture Data

The used motion capture system is a Motion Analysis system relying on six Eagles and four Hawks cameras. It provides an estimation of the camera position at 200 Hz with a precision error of less than one cm.

VI. FUTURE WORK AND CONCLUSIONS

In this article, we demonstrated that vision-based localization allows humanoid robots to achieve complex tasks. By integrating the vision algorithm into the control framework, we have been able to gradually reshape the trajectory in order to compensate for execution errors. This allowed the robot to achieve complex tasks which would be difficult to realize otherwise.

However, several intesresting additions remain to be integrated. First, the local trajectory modifications are not checked to make sure that no auto-collision occurs. Currently, a conservative maximum correction (i.e. 2 cm every two steps) is imposed for safety. By applying recent works such as fast feasibility tests described by Perrin et al. [START_REF] Perrin | Approximation of Feasibility Tests for Reactive Walk on HRP-2[END_REF], we would be able to increase the maximum correction without compromising the robot safety. Second, the HRP-2 robot embeds an Inertial Measurement Unit which could be used to estimate the robot chest attitude. It would definitely be interesting to provide an initial estimation of the robot motion to help the localization process. In addition, we are interested in improving the capabilities of our visual SLAM and visionbased localization systems towards the goal of long-term localization to deal with possible changes in the environment.

Fig. 1 .

 1 Fig. 1. HRP-2 dropping a ball on a shelf after having navigated through an indoor environment. Navigation followed by manipulation are typical scenarios where walking drift can prevent a task from being accomplished.

Fig. 2 .

 2 Fig. 2. Complete robotics framework overview. Let I left , I right be respectively the left and right camera images. Tcam, Tcam respectively the planned and estimated camera position. γ the concatenation of the left ankle, right ankle, upper body and center of mass trajectories. γ the trajectories reshaped by taking into account execution errors.

Fig. 3 .

 3 Fig. 3. The HRP-2 robot drops a ball on a shelf while localizing itself using vision based localization.

Fig. 4 .

 4 Fig. 4. Experimental setup description (top view). The dotted line displays the robot waist trajectory.

Fig. 5

 5 Fig.5contains the camera position estimated by the motion capture system and by the localization algorithm. The mean error is 0.2 m in translation. Drop in the algorithm precision can occur when the robot enters a part of the map which is not dense enough or where the environment does not incorporate enough texture to detect enough interest points. It happens in our scenario near the end when the camera is too close to the shelf to perceive its environment correctly.The used motion capture system is a Motion Analysis system relying on six Eagles and four Hawks cameras. It provides an estimation of the camera position at 200 Hz with a precision error of less than one cm.

Fig. 5 .

 5 Fig. 5. Camera position estimated by both the motion capture system and the localization algorithm.

ACKNOWLEDGMENTS


Philippe Marie Durand

Maurice Marin

François Faure

Bruno Raffin

Philippe Marin

DEM-based simulation of concrete structures on GPU

Keywords: GPU, DEM, reinforced concrete structures MOTS-CLÉS : GPU, DEM, structures en béton armé

come

Introduction

The design of some particular civil engineering structures must take into account the risk of severe dynamic loadings due to natural or anthropogenic hazards such as rock falls, aircraft or missile impacts. Often, these severe loadings lead to localized fractures and fragmentation in the concrete structure. The Discrete Element Method (DEM) [START_REF] Cundall | A discrete numerical model for granular assemblies[END_REF] is an appropriated method for modeling such discontinuities. The model uses disordered assembly of spherical and rigid elements of different sizes and masses to reproduce an isotropic and homogeneous behavior at a macroscopic scale. This method is very well adapted to dynamical problems, and does not rely upon any assumption about where and how a crack or several cracks occur and propagate as the medium is naturally discontinuous.

These DEMs were used first to model the behavior of granular materials, but they also provide very accurate results for cohesive materials like concrete (D'addeta et al., 2002). The studies of Camborde in 2D [START_REF] Camborde | Numerical study of rock and concrete behaviour by discrete element modelling[END_REF] or Rousseau in 3D [START_REF] Rousseau | Damage prediction in the vicinity of an impact on a concrete structure: a combined FEM/DEM approach[END_REF] demonstrated the efficiency of such a discrete approach to deal with impact problems on reinforced concrete structures. They also pointed out the heavy computational load of DEM, limiting its use to small structures. To reduce the computational cost, we can use a coupling between the Discrete Element Method and the Finite Element (FE) Method [START_REF] Xiao | A bridging domain method for coupling continua with molecular dynamics[END_REF], [START_REF] Dhia | The Arlequin method as a flexible engineering design tool[END_REF], [START_REF] Frangin | On the use of combined finite/discrete element method for impacted concrete structures[END_REF], [START_REF] Rousseau | Multidomain finite and discrete elements method for impact analysis of a concrete structure[END_REF]. In the vicinity of the impact, where important non-linear phenomena occur, the medium will be modelled by means of Discrete Elements (DEs). The use of the FE method far from the impacted area is a way to reduce this limitation since in most cases severe degradation phenomena are localized in the vicinity of the impact.

But for shell structure [START_REF] Rousseau | A discrete element/shell finite element coupling for simulating impacts on reinforced concrete structures[END_REF] impacted by large projectile like an aircraft, the size of the area that is represented by DEs is very important leading to an important computation time. Code parallelization is a classical approach to decrease simulation time or enable larger simulations. The goal of this paper is to evaluate the benefits of the parallelization of DEM based concrete structure simulation on Graphics Processing Unit (GPU). GPUs are coprocessors, usually having their own memory, communicating with the CPU through the PCI Express bus. GPUs were first dedicated to 3D graphics rendering, but as GPUs evolved towards more programmable architectures capable of executing user developed codes, it became possible to use them for performing generic computations. Today, using GPUs as general coprocessors is a major trend in high performance computing, usually called General Purpose GPU (GPGPU). Their high peak performance associated with their moderate cost make them good candidates to dramatically boost the computing power of a PC or the node of a super computer. Some GPUs are dedicated to GPGPU like the NVIDIA Tesla family (660 GigaFlops peak double precision floating point performance for the NVIDIA Tesla M2090), and software environments like OpenCL or CUDA are targeted at GPGPU programming. GPUs are highly parallel architectures, significantly diverging from CPUs (cf. section 3). An implementation on GPU can lead to a significant performance improvement over a CPU implementation [START_REF] Lee | Debunking the 100X GPU vs. CPU myth: an evaluation of throughput computing on CPU and GPU[END_REF].

However, such performance gains often require a significant programming effort and may be limited if the parallelism the application can exhibit does not match the GPU architecture. DEMs are good candidates for an efficient GPU implementation. A significant part of DEM related computations are data parallel (one instruction can be executed in parallel on many different data) and memory needs of today's DEM simulations usually fit the GPUs capabilities.

We first give a brief discription of the DEM we use in Section 2. We then sketch in Section 3 the principles of GPU architecture and programming, and their consequences on DE simulation. Our implementation is presented in Section 4, and we present experimental results in Section 5, showing speed-ups of an order of magnitude compared to an execution on a single core CPU. We finally conclude and discuss future work in Section 6.

Discrete Element Model

The DEM is based on the modeling of the continuum by means of rigid particles with 6 degrees of fredoom, 3 translations and 3 rotations. Interaction laws between DEs determine the macroscopic constitutive behavior. In the early developments particle interactions relied on friction laws for non-cohesive materials like sands [START_REF] Cundall | A discrete numerical model for granular assemblies[END_REF]. Interaction laws for cohesive materials were defined later on [START_REF] Hentz | Identification and Validation of a Discrete Element Model for Concrete[END_REF]. To guarantee reasonable calculation durations, a model based on the Distinct Elements Model [START_REF] Cundall | A discrete numerical model for granular assemblies[END_REF] with rigid spheres was chosen. Two types of interaction are defined. The initial interaction between two elements is generally a link interaction (the two elements are not necessarily in contact). Initially, two DEs interact if the distance between their centroids is less than a given radius of interaction. During the simulation additional interactions of contact type can be added. For concrete material, we used a modified Mohr-Coulomb model with softening [START_REF] Rousseau | Damage prediction in the vicinity of an impact on a concrete structure: a combined FEM/DEM approach[END_REF]. More sophisticated laws taking into account compaction phenomena can be elaborated [START_REF] Tran | A discrete element model of concrete under high triaxial loading[END_REF] but they are not necessary to model a thin slab mainly subject to flexion and tension effects.

We can find in [START_REF] Rousseau | Damage prediction in the vicinity of an impact on a concrete structure: a combined FEM/DEM approach[END_REF] a procedure to identify all the material parameters based on the simulation of quasi-static compression and tension tests. The point is to identify local parameters to model macroscopic values such as compressive strength σ c , tensile strength σ t and fracture energy G f . Special links are used to represent steel reinforcements and steel-concrete interface [START_REF] Potapov | Advanced simulation of damage of reinforced concrete structures under impact[END_REF]. In this paper, we rely on the interaction laws proposed by [START_REF] Rousseau | Damage prediction in the vicinity of an impact on a concrete structure: a combined FEM/DEM approach[END_REF] and [START_REF] Potapov | Advanced simulation of damage of reinforced concrete structures under impact[END_REF]. Refer to these publications for more details.

GPU Architecture and DE Parallelization Issues

Since our target parallel architecture is the GPU, we briefly introduce its architecture and its consequences on DE simulation programming. The Graphics Processing Units (GPUs) are high-performance many-core coprocessors initially designed to accelerate graphics rendering. They are increasingly used to accelerate general purpose scientific and engineering computing (GPGPU). The GPU, or device, is controlled by a CPU, also called host. The GPU architecture being significantly different form the CPU, we first introduce its most important characteristics necessary to understand the design and implementation of the GPU specific DEM algorithm we propose. We then present CUDA, the de facto standard GPGPU programming environment. This allows us to finally discuss the parallelization of the DE simulation.

GPU Architecture Overview

The GPU is built from multiple computing cores called multiprocessors. Each multiprocessor is an array of synchronous scalar processor units enabling concurrent SIMD computations. SIMD stands for Single Instruction Multiple Data streams, meaning that all scalar processors in a multiprocessor can only execute the same instruction, on different data, at a given time. It is well adapted for data-parallel applications, where a large amount of data undergoes similar computations. Its limitation is that when executing a conditional block, all scalar processors evaluating the condition to true proceed synchronously with the block execution while the other scalar processors are kept idle. Obviously, full efficiency is achieved when all scalar processors follow the same execution path.

Multiprocessors have less double precision floating point units than single precision ones. Thus, performance can significantly drop when moving floating operations from single to double precision. For example, recent GPUs from the NVIDIA Tesla family feature up to 665 GigaFlops of double precision perfomance and 1 TeraFlops of single precision performance (the ratio single/double was about 8 in previous generations).

The GPU manages its own memory, separate from the CPU. Data transfers between the CPU and GPU are explicitly controlled by the application. The GPU memory is divided between a global memory all multiprocessors can access, and a shared memory local to each multiprocessor. Data transfers between global and shared memory are also explicit. GPUs also have a read-only global memory (constant memory), one part having special access functions (texture memory), but we will not detail these memories which we do not use in this work.

When the data accessed by different threads of a multiprocessor is contiguous, these accesses are coalesced, i.e. the data is gathered to be transfered in a single pass. Otherwise they are performed sequentially, incurring a strong performance penalty.

Concurrent writings to the same address can lead to inconsistencies. It is up to the programmer to avoid this, or to use atomic operations. An operation is said atomic when it can not be interrupted by any concurrent processes before it ends up. Atomic operations are transparently sequentialized by the GPU, at the price of a loss of efficiency.

Each multiprocessor has a limited number of registers shared by scalar processors. Only scalar processors that can have their register needs filled will be able to perform computations, the other will be idle. Register use must therefore be carefully managed to optimize the performance.

GPU Programming Overview

NVIDIA is currently the leader in GPGPU, and we therefore focus on their architecture and associated programming environnement, CUDA (NVIDIA Corporation, 2011). The implementation presented in this paper could be ported to the OpenCL emerging standard for GPGPU programming without major difficulty (OpenCL encompasses the CUDA programming model), with probably lower performance, as OpenCL compilers still need to gain in maturity.

The CUDA programming model closely matchs the GPU architecture. The base parallelization unit is a thread, which is a sequence of instructions executed on a scalar processor. The CUDA programmer writes programs to be executed on a scalar processor, also called kernels. The kernels are transfered to the GPU and their parallel execution is parameterized by block sizes. A block is a group of threads executed by the same multiprocessor. Inside a block, each thread is identified by a unique index. The address of the thread data in the memory is easily computed using the block and thread indices. When the number of threads per block is greater than the actual number of scalar processors, the execution is composed of a sequence of warps. Each warp contains one thread per scalar processor. The number of threads per block should thus be a multiple of the warp size. If the active threads become idle waiting for data from the memory, they will be transparently suspended to allow an other warp to become active, enabling to overlap computations with memory accesses. All threads within the same warp are executed synchronously as imposed by the SIMD architecture of the multiprocessor. Communications between the threads belonging to the same block are enabled through memory sharing and synchronization instructions. Threads in different blocks can not directly interact.

Parallelization of the DE Simulation

In this section, we review the main phases of the DE simulation loop and we discuss the parallelization issues due to the GPU architecture, especially the concurrent writing. The simulation relies on an explicit integration scheme. The integration loop (Figure 1) can be divided in three main parts. In the first part, the positions and ve-locities of the DEs are updated based on known forces using Newton's laws of motion within a given timestep. The parallelization is straightforward, each DE being processed independently.

The second part consists in detecting new contacts between the DEs, which interact if the distance between their centroids is less than the sum of their interaction radii. This can be compute intensive, as the number of potential contact pairs grows quadratically along with the number of DEs. A common approach to accelerate the collision detection process is to define a 3D grid that divides the space filled by DEs into cubic cells, with size equal to the maximum interaction diameter of all DEs. The first collision detection step consists in identifying the cell each DE belongs to. Then for each DE we search for new contacts browsing other DEs from the cell it belongs to and from the 26 neighbor cells. This is again data-parallel, but appending a new contact to the list requires a concurrent access. To avoid duplicating contacts, we need for each element an up-to-date list of every other DE it interacts with. Managing this dynamic set of links significantly complexifies the implementation. Without altering the simulation correctness, the collision detection can be executed periodically and not at every iteration. This period depends on the simulation parameters. In our simulations, collision detection is typically scheduled every 10 time step.

The last part is force computation and link update. Each force computation is costly, as the high number of parameters required to represent the concrete structure behavior at a macroscopic scale leads to numerous operations. The interaction force can be computed in parallel on each link. However, the accumulation of forces on each DE requires a concurrent access to the force vector of each DE. We then detect the links and contacts that break due to excessive distance or interaction force. This is straightforwardly parallel, but the list of contacts needs to be updated. The initial links can be marked as broken, to avoid inefficient list compactions.

Finally, since stability is a major issue with explicit integration schemes, the time step is dynamically adjusted at the end of every step (Subsection 4.4). Using DEs' stiffnesses for translation and rotation, we compute interactions stiffnesses that are then accumulated for each DE, which requires another concurrent write operation. The time step is then set to the minimum over every coordinates of the square root of the ratio of inertia over stiffnesses.

Comparison with Molecular Dynamics

DEM simulations are close to Molecular Dynamics (MD) simulations, notably because both are dealing with moving particles, with interactions between particles leading to forces, as well as periodic collision detections (neighbor search in MD). [START_REF] Anderson | General purpose molecular dynamics simulations fully implemented on graphics processing units[END_REF] proposed a GPU implementation for MD, later improved by [START_REF] Rapaport | Enhanced molecular dynamics performance with a programmable graphics processor[END_REF], for example using an interaction cutoff range for the dimension of the detection grid during the neighbor search. Though their implementation shows several similarities with the one we propose in this paper, there remain important differences. First, interactions between DEs have a stateful history while forces corresponding to interacting MD particles observe some well defined laws. Thereby, in MD simulations when the neighbor list is updated there is no need to unselect already existing interactions, which requires storing the list of existing interactions. Moreover DEM computations are more complex, since DEs have rotational degrees of freedom, and are of different types (steel, concrete), requiring different interaction types. Consequently, DEM computations are so greedy in term of operations and registers need that it is efficient to compute one force per interaction, and then to accumulate these forces at each DE. Conversely, because MD is significantly less computationally intensive, it is more efficient on GPU to compute all interaction forces for each atom (each force is thus computed twice).

Implementation on the GPU

We now detail our GPU implementation composed of a sequence of 19 CUDA kernels (Figure 2). Many of these kernels perform data parallel operations running one thread per link, contact or DE. The implementation challenge is to manage memory accesses as well as double precision computations. For that purpose data is maintained in special structures in the global memory to ensure correct alignment in memory to favor coalesced parallel accesses.

Positions and Velocities

Update Ú ÐÓ Ø × ÖÓÑ Ð Ö Ø ÓÒ× ÔÓ× Ø ÓÒ× ÖÓÑÎ ÐÓ Ø × Collision Detection ÓÑÔÙØ Ö Ñ Ò× ÓÒ× Ø Ð Ñ ÒØÈÖÓÔ ÖØ × Ö ÓÖ Ö Ø Ò Ò ÐÐËØ ÖØ ÓÖÄ Ò × Forces, Moments and Internal Variables Computation ÓÑÔÙØ ÓÖ ÓÒÖ Ø ÓÖ ÓÒÖ Ø ÓÑÔÙØ ÓÖ ËØ Ð ÓÖ ËØ Ð ÓÑÔÙØ ÓÖ ÓÒÖ Ø ËØ Ð ÓÖ ÓÒÖ Ø ËØ Ð ÓÑÔÙØ ÓÖ ÓÒØ Ø ÓÖ ÓÒØ Ø Timestep Evaluation Ö × ØËØ Ò ×× × ÓÑÔÙØ ËØ Ò ×× × ÓÑÔÙØ Ö Ø ÐÌ Ñ Ò ÄÓ Ð × Å Ò Acceleration Computation Ð Ö Ø ÓÒ× ÖÓÑ ÓÖ ×
Figure 2: CUDA kernels related to integration steps.

Memory Layout

Let B be the number of threads per block and N the number of DEs. All state vectors (position, velocity, forces and accelerations) are stored per block of B elements. All the vectors have 3 or 6 components. For instance a force has six coordinates, three for centroids positions and three for the element rotations. The 6 × N force components are stored continuously in memory, with padding data at the end to obtain a ⌈ 6×N B ⌉ × B size array (Figure 3). Within a block, all threads fetch the force components from the GPU global memory into the shared memory by series of B contiguous, aligned, concurrent and thus coalesced memory accesses. Next, after a synchronization to ensure all memory transfers are effective, each thread loads from the shared memory the data it needs for computations. At the end of the kernel, the updated forces stored in shared memory are written back to the global memory following the same pattern. This data organization ensures efficient memory transfers while keeping data conveniently organized in the GPU global memory if the CPU needs to access them. The rest of the internal data is more directly stored in global GPU memory in arrays of B size blocks. For instance the Ð Ö Ø ÓÒ× ÖÓÑ ÓÖ × kernel computes accelerations using masses and inertia. For sake of commodity and to favor data locality, masses and inertia are stored in a single array interleaving B masses and B inertia as illustrated in figure 3. Thus the threads of a block access masses or inertia data taking advantage of the GPU coalescing capabilities (Figure 4). In this case, the data is directly loaded to registers without using the shared memory. This layout is applied to every set of internal variables stored in the GPU global memory and used only by the device.

Force Computation

A series of kernels are in charge of force computations (ÓÑÔÙØ ÓÖ ¶). Each kernel takes care of one interaction type, operating on links or contacts only. We adopted this multi-kernel organization to avoid the divergent conditional branches required to support different interaction types. The amount of data parallelism available would have been reduced. Next, forces are accumulated on each DE, a very data parallel operation the ÓÖ takes care of. In a recent paper, [START_REF] Shigeto | Parallel computing of discrete element method on multi-core processors[END_REF] directly accumulate computed force contributions using an atomic add operation. However this operation is costly. Though our approach requires to maintain the list of interactions each DE is involved in, it does not require any atomic

Collision Detection

Collision detection is performed using four CUDA kernels. ÓÑÔÙØ Ö Ñ Ò× ÓÒ× computes the bounding box of the current simulation state. Ø Ð Ñ ÒØÈÓ× Ø ÓÒÁÒ Ö identifies the cell each DE belongs to. The third one, Ö ÓÖ Ö Ø Ò Ò ÐÐËØ ÖØ, sorts DEs according to cell indices and computes the memory offset to identify for each cell the starting address of its list of DEs. Then the kernel ÓÖ ÓÒØ Ø× detects the new contacts. More precisely, the first step computes the dimension of the bounding box that includes all DEs in the three-dimensional space. This is done by an optimized parallel search of both maximum and minimum among DE positions . The cell size is set to a multiple of the size of the maximum interaction radius in every of the three dimensions. Next, a kernel computes the cell index of each DE. We store a cell index as a key and the DE index as the value in two corresponding arrays. Next, we sort the DE indices according to their cell index using the parallel radix sort operation from the CUDPP library [START_REF] Satish | Designing efficient sorting algorithms for manycore GPUs[END_REF]. The size of these structures is proportional to the number of DEs and not to the grid resolution. To avoid processing the whole cell array searching for the neighbors of a given DE, an other kernel builds an indirection table giving for each cell the beginning and end index of its DE list. To improve the local coherence when searching for neighbors, DE positions and radii are also sorted into new structures according to their particle order. This method has been implemented by NVIDIA in a simple demo called Particles in the CUDA SDK (NVIDIA, 2010), although they do not need to compute the bounding box at each step as the simulation space does not change during the simulation.

Everything being set up, the interaction search can begin ! Search is parallelized over DEs. Every thread is in charge of testing new interactions with the particles present in its cell and neighbor cells. To avoid redundant computations, only DEs with higher indices than the current DE are tested in the DE's cell, and only half of the neighbor cells are tested. We also maintain the neighbor list of each DE, used to check if DEs are already connected together (either with an initial link or a previously detected contact). When a new interaction is detected, it is pushed into a list.

Time step

To compute the new time step (Figure 5), we first reset the values computed during the previous step with a call to Ö × ØËØ Ò ×× × relying on the GPU Ñ Ñ× Ø function. Stiffnesses are computed per interaction type by ÓÑÔÙØ ËØ Ò ×× × and cumulated in the corresponding elements stiffnesses. Then the ÓÑÔÙØ Ö Ø ÐÌ Ñ × kernel computes a local critical time T c per DE from the 6 stiffnesses (three position and three rotation components) and the inertia. This kernel ends by executing a parallel minimum search to obtain the smallest critical time, i.e. the time step for the next step.

Finally, accelerations are computed based on forces by the Ð Ö Ø ÓÒ× ÖÓÑ ÓÖ × kernel, followed by an update of velocities and positions (Ú ÐÓ Ø × ÖÓÑ Ð Ö Ø ÓÒ× and ÔÓ× Ø ÓÒ× ÖÓÑÎ ÐÓ Ø ×). Attention must be paid to floating-point operations. Depending on the GPU and the CUDA version used, default operators are not IEEE-754 compliant and must be changed for IEEE-754 compliant ones. In our case kernels are implemented using IEEE-754 with round to nearest compliant operations and can be compiled to support either single or double precision floating point data.

Results

We measured the performance of our implementation for the simulation of an impact on a reinforced concrete slab (Figure 6). The scene involves 14, 274 DEs with 784 being reinforcing steel rods, 12, 718 the slab concrete and 772 the concrete representing the impactor. The slab dimensions are 1m x 1m x 0.1m. There are 79211 cohesive links at start. The GPU used is a NVIDIA Fermi C2050 with 448 cores and 3Go of memory while the CPU is a 6-core Intel Xeon X5650 @2.67GHz (only one core is used for running the simulation). are constant because the impactor is not yet colliding with the slab. The collision detection execution is 10 times slower than the other calculation stages on the CPU (6 times on GPU). Thus, when performed every 10 steps, it takes about 50% of the average execution time on the CPU (70% on the GPU). The speed-up reached by our GPU implementation is about 15 when the collision detection is not performed and about 30 otherwise compared to the execution time on the GPU.

Results during the impact are presented to show the evolution of the simulation behavior (Figure 7b). Considering steps without collision detection, on the CPU, the execution time increases during the simulation because of the creation and removal of interactions. On the GPU, the execution time appears constant. The data structures for initial interactions are not altered during the simulation, broken links being simply invalidated, which does not significantly change the execution time needed to launch and perform the kernels for these interactions. The execution time can increase when contacts are created. However, as these extra contacts are processed in parallel, their number needs to be important to have a significant impact on the overhall execution time.

During the collision detection on the CPU, the particles are sorted in cells using linked lists. Dynamic memory allocation can lead to an execution time overhead, but this is efficient regarding memory space, and execution time is not affected by the spreading of particles in space. Conversely, on the GPU, the parallelism is reduced when particles start being fired away so the collision detection becomes less efficient. As a consequence, the execution time on the GPU tends to increase for collision detections.

Conclusion

We presented how the DEM can be used to simulate concrete structures on GPU. After a reminder of the method and its parametrization, we detailed our GPU implementation and simulation results. One challenge of the DEM is its high computational cost. However, its inherent data parallel structure and limited memory footprint make it a good candidate to take advantage of the high performance of GPUs. Our experiments show a speed-up of an order of magnitude for 14274 particles compared to an execution on a single core CPU. This implementation work revealed several difficulties. A careful management of memory access patterns is essential to ensure that the GPU processing units do not stall waiting for data. The necessity of managing 2 different types of particles, for concrete and for steel, reduces the data parallelism of the application and almost duplicates the number of GPU kernels to execute. Only simulations with a high particle count amortize the associated overheads. The neighbor computation step that is required to identify the interacting particles exhibits a limited level of data parallelism that, creates a performance bottleneck, even using optimized sorting GPU implementations.

The high implementation cost may counterbalance the performance gain, especially given today's availability of CPUs with up to 12 cores. A careful parallel implementation on such processors could probably lead to a significant speed-up. Moreover, good tools for debugging and monitoring GPU implementations are still missing, making the development painful. However GPGPU is still in its early stage and is quickly evolving. The next GPU generation will provide features that may give a significant advantage to GPGPU. The Intel MIC coprocessor [START_REF] Seiler | Larrabee: a many-core x86 architecture for visual computing[END_REF] has a large number of x86 cores (32 for the current development prototype), each core having a large vector unit for data parallel operations. This architecture will offer a higher level of flexibility, enabling cores to work asynchronously, locally having to extract only a limited amount of data parallelism to fill the vector unit, while globally showing a similar amount of parallel processing power. Because the MIC relies on a classical x86 architecture, we can expect to benefit from the development tools already available for x86 processors. The AMD APU (Advanced Processing Unit) integrates on a single chip a CPU and a GPU. It will enable the CPU and the GPU to share the same memory, significantly decreasing the cost of data exchange between the CPU and the GPU.

Figure 1 :

 1 Figure 1: Global simulation integration scheme.

Figure 3 :Figure 4 :

 34 Figure 3: GPU structure layout for DE masses and inertia. An array of B elements masses is associated with an array of B elements inertia in a structure called gpuED. The whole set is constructed by an ⌈ N B ⌉ length array of these structures.

Figure 5 :

 5 Figure 5: Detailed process for checkTimestep with CUDA kernels and Cuda structures names.

Figure 6 :

 6 Figure 6: Simulation of an impactor projected against a reinforced concrete slab. The picture has been taken at step 2340, during the impact itself.

Figure

  Figure 7a presents an example of each step execution time at different stages of the simulation. This plot figures the starting of the simulation. The peaks correspond to the steps performing the collision detection. At the beginning, execution times

Figure 7 :

 7 Figure 7: Execution time (ms -logarithmic scale) per step on the GPU or the CPU. Time peaks occur when the collision detection is performed.