Zhiguang Liu

Antonio Mucherino
email: antonio.mucherino@irisa.fr

Ludovic Hoyet
email: ludovic.hoyet@inria.fr

Franck Multon
email: franck.multon@irisa.fr

Surface based Motion Retargeting by Preserving Spatial Relationship

Keywords: Computing methodologies → Animation, Optimization algorithms, Character Animation, Motion Retargeting

Retargeting motion from one character to another is a key process in computer animation. It enables to reuse animations designed for a character to animate another one, or to make performance-driven be faithful to what has been performed by the user. Previous work mainly focused on retargeting skeleton animations whereas the contextual meaning of the motion is mainly linked to the relationship between body surfaces, such as the contact of the palm with the belly. In this paper we propose a new context-aware motion retargeting framework, based on deforming a target character to mimic a source character poses using harmonic mapping. We also introduce the idea of Context Graph: modeling local interactions between surfaces of the source character, to be preserved in the target character, in order to ensure fidelity of the pose. In this approach, no rigging is required as we directly manipulate the surfaces, which makes the process totally automatic. Our results demonstrate the relevance of this automatic rigging-less approach on motions with complex contacts and interactions between the character's surface.

INTRODUCTION

Performance-driven animation of human-like characters has been widely explored, especially with the dissemination of cheap and easy-to-use motion capture devices. Users generally wish to control characters with various shapes, ranging from realistic to cartoon characters. Moreover, animations studios have gathered Terabytes of rigged animations accurately adapted to a given character in specific conditions. Being able to automatically reuse this data to make new characters move in the same way is essential to save time and associated costs. The problem is: how to accurately transfer a motion perfectly adapted to character A (either an actor/user, or a preprocessed virtual character) to character B, whereas A and B have different shapes and environments to interact with.

To answer this question, motion retargeting has been widely explored [START_REF] Gleicher | Retargetting Motion to New Characters[END_REF]]. However, most of previous approaches rely on a rigged skeleton, whereas the contextual meaning of the motion is mainly linked to relationship between skin surfaces. Hence, motions such as applauding, touching the belly, or crossing the legs, need manual editing to correct the errors introduced by skeletonbased motion controllers driven by inverse kinematics. Even though skeleton animation looks fine, manual post-processing is needed to correct artifacts introduced when animating the body surface. Moreover the rigging process is tedious and has to be done for each new character, leading to a reprocessing of all the required motions.

In order to tackle these limitations, this paper proposes a new approach to get rid of the rigging process and skeleton control, as suggested in previous works [START_REF] Baran | Semantic Deformation Transfer[END_REF][START_REF] Rhodin | Interactive motion mapping for real-time character control[END_REF]. We aim at pushing the idea further to directly retarget a database of motions to new characters while ensuring that the contextual meaning linked to the body surface is preserved, without manipulation of rigged skeletons. As shown in Figure 1, we input the source calibration pose in Figure 1(a), a sequence of deformed source poses in Figure 1(c) extracted from the prerecorded database, and a target calibration pose in Figure 1(f). Then our method outputs a set of newly synthesized target poses for the target character, preserving the relationship between body surfaces of the deformed source poses in Figure 1(i).

Because of the versatility of computer graphics applications, poses such as those explored in this paper can typically be represented either by manifold mesh or, most of the time, by nonmanifold mesh such as multiple-part objects, polygon soup or tetrahedral volumetric meshes. In order to reuse existing shape deformations and deal with various representations of shapes, we therefore propose a new surface-based motion retargeting framework. Our method is general and does not require the source and target calibration meshes to be compatible, namely, to have the same number of vertices, triangles, and connectivity. Our method needs to set a sparse set of joints inside the source and target calibration shapes. These joints can be provided by skeleton joint extraction algorithms automatically or obtained from user specification. Then we generate the retargeted animation based on the joints and available source mesh sequences. In contrast, the rigging process requires the users to perform operations such as setting up the skeleton, painting body part weight, and skeleton adjustment, which is much more tedious and time-consuming. In this sense, our method simplifies the work of the animation process. However, we assume there are clear semantic correspondences between the surface of the two 3D domains.

We use harmonic function as the mapping function, as previous work showed that it results in a smooth and detail-preserving mapping between two 3D domains [Ben-Chen et al. 2009b;[START_REF] Ferrari | Harmonic deformation of Delaunay triangulations[END_REF][START_REF] Joshi | Harmonic coordinates for character articulation[END_REF][START_REF] Zayer | Harmonic guidance for surface deformation[END_REF]]. However, the Laplacian operator and its derivative can be quite complicated functions, which makes it hard to analytically solve these integrals. Inspired by the Method of Fundamental Solutions (MFS) [START_REF] Fairweather | The method of fundamental solutions for elliptic boundary value problems[END_REF][START_REF] Li | Meshless Harmonic Volumetric Mapping Using Fundamental Solution Methods[END_REF]] that discretized the problem by using a linear combination of component functions, we construct an approximated harmonic mapping using MFS. Results presented in Section 8 show that our method using MFS produces very small reconstruction error.

The major contributions of this paper are summarized as follows:

• We propose a new unified framework for motion retargeting consisting in two steps: 1) extracting deformation descriptors and 2) synthesizing new deformations based on these deformation descriptors. During the deformation descriptor extraction process, we compute the affine transformation by solving a linear system. During the deformation synthesis process, we generate new poses by minimizing the difference between affine transformation descriptors and other constraints such as position, spatial relationship, and volume preservation. • We propose a new structure to capture the spatial relationship among different body parts on the surface of a mesh, which we call Context Graph. By minimizing the deformation of the Context Graph between a source animated pose and a synthesized pose, we are able to preserve the contextual meaning for the newly synthesized pose.

The rest of the paper is organized as below. We first review the related work on motion adaptation and shape deformation in Section 2. Section 3 gives an overview of our method. Algorithm details on deformation analysis and synthesis will be discussed in Section 4 and 5, respectively. In Section 6, we describe the construction of our Context Graph and how it is deformed to preserve spatial relationships and body parts volume, followed by Section 7 with some implementation details. Experimental analyses and evaluations are conducted in Section 8. Finally, in Section 9, we conclude the paper, as well as discuss the limitations and future research directions.

RELATED WORK 2.1 Motion Retargeting

Motion retargeting techniques tackle the problem of adapting an animated motion from one character to another, and they have been widely explored to reuse existing motion capture data in computer graphics [START_REF] Gleicher | Retargetting Motion to New Characters[END_REF][START_REF] Kulpa | Morphology-independent representation of motions for interactive human-like animation[END_REF].

Ho et al. [2010] proposed a new structure called interaction mesh to preserve the spatial relationships between body parts/objects in the scene by iteratively minimizing a Laplacian deformation energy. Their method generates sensible motions for close interactions as well as those without any contacts. However, their method is not suitable for real-time animation application as the algorithm performs the optimization over the whole animation sequence. Similarly, Al-Asqhar et al. [2013] propose to use relationship descriptors to describe the kinematics of the body parts by the weighted sum of translation vectors with respect to points sampled on the surface of the object in close proximity. Their method adapt the motion of a character to large updates of character morphologies and scenarios on-the-fly. However, re-position of contacts may be invalid due to large deformations. This work was then extended by [START_REF] Tonneau | Character contact re-positioning under large environment deformation[END_REF][START_REF] Ali Al-Asqhar | Relationship descriptors for interactive motion adaptation[END_REF]], who presented a scheme that enhances the detection and re-position of unnatural motion adaptations. [START_REF] Kim | Retargeting human-object interaction to virtual avatars[END_REF] develop a spatial map, which builds a bijection between two domains bounding the source and target objects. During the online motion generation process, the synthesized motion tries to preserve the spatial relationship between the source pose and an object. However, the algorithm needs to tetrahedralize two domains and requires the shape of source and target object to be very similar. [START_REF] Molla | Egocentric Mapping of Body Surface Constraints[END_REF] propose an egocentric normalization of body surfaces to preserve the spatial relationships between the limbs and the other body parts. They test the proposed method a variety of characters with and without self-contacts. There are important differences between the method of [START_REF] Molla | Egocentric Mapping of Body Surface Constraints[END_REF], as they start from the joints with IK solver to reconstruct new postures, while our method works on the surface of the source and target meshes directly.

We are interested in generating realistic animations transfering existing motions of a character onto other characters with very different body morphologies. All of the above-discussed motion retargeting methods work on sparse joints. Thus, the speed and quality may decrease greatly when they are extended to a surface mesh, as the dimension of a surface mesh model is much higher than the sparse joints, e.g. a typical surface mesh model contains 20K vertices compared to approximately 50 sparse joints in a skeleton.

Surface Deformation

In this Section, we focus on direct surface manipulation deformation methods. Differential coordinates have been widely used to represent local details of a mesh. [START_REF] Lipman | Differential coordinates for interactive mesh editing[END_REF] reconstruct the surface by minimizing the difference of Laplacian coordinate before and after the deformation. This method requires to explicitly find a correct rotation for the Laplacian coordinate of a local vertex, otherwise, distortion would have resulted. Similarly, [START_REF] Sorkine | Laplacian Surface Editing[END_REF] propose a Laplacian representation by encoding each vertex relative to its neighborhood. However, unlike Lipman et al. [2004], the transformation of the differential coordinate for each vertex is implicitly computed from the transformation of a control handle. However, for transformations that involve large rotation, the deformation will suffer from anisotropic scales and need a postprocess to refine it. [START_REF] Zhou | Large mesh deformation using the volumetric graph laplacian[END_REF] enhance the Laplacian surface editing scheme by introducing a volumetric graph. Specifically, they embed a cubic lattice inside the shape, as well as an offset surface wrapping the exterior of the shape. The Laplacian editing framework is then employed to preserve the volume while keeping the local surface details. [START_REF] Matthew | MoSh: Motion and Shape Capture from Sparse Markers[END_REF] estimate body shape and pose by using a set of standard motion capture marker data. They successfully approximate the soft-tissue deformation as shape variations within the space of static human body shape. [START_REF] Wang | Linear subspace design for real-time shape deformation[END_REF] design linear subspaces for controlling shape deformation, which cuts down the time complexity of variational shape deformation methods. Methods such as [Boukhayma and Boyer 2017;Boukhayma et al. 2017] transfer or generate variations of 3D models by building a probabilistic low dimensional embedding of shape poses with Gaussian Process Models. [START_REF] Casas | Interactive Animation of 4D Performance Capture[END_REF] propose the parametric motion graphs that enables real-time animation editing based on a database of temporally aligned mesh sequence reconstructions.

Most of direct surface manipulation techniques only work on the surface of a manifold mesh, which are hard to be applied to other mesh structures such as polygon soups and multi-component shape. In this paper, we apply the Laplacian deformation framework [START_REF] Sorkine | Laplacian Surface Editing[END_REF] to our Context Graph that preserves the volume and spatial relationship among different body parts. We approximate the harmonic mapping by using a linear combination of the fundamental solution of harmonic equation [START_REF] Li | Meshless Harmonic Volumetric Mapping Using Fundamental Solution Methods[END_REF][START_REF] Poullikkas | Methods of fundamental solutions for harmonic and biharmonic boundary value problems[END_REF]], which enables us to process both manifold and non-manifold shapes.

Space Deformation

A space deformation maps a source domain to a target domain inside the Euclidean space while preserving user-defined constraints. Cage-based deformation methods [Ben-Chen et al. 2009b;[START_REF] Joshi | Harmonic coordinates for character articulation[END_REF][START_REF] Lipman | Green Coordinates[END_REF]] are such kinds of approaches that can be used to deform any type of shapes inside the cage, such as polygon soup, point cloud, and manifold mesh. Compared with direct surface based deformation that manipulates a set of dense points, controlling a sparse polyhedral cage simplifies the deformation process.

Shapes inside a cage can be represented as a linear combination of the vertices of the cage by barycentric coordinates [START_REF] Michael S Floater | Mean value coordinates in 3D[END_REF][START_REF] Huang | Subspace gradient domain mesh deformation[END_REF]]. However, barycentric coordinates are not positive for a concave domain, which may result in stretching in the resulting deformation. Harmonic coordinates [START_REF] Joshi | Harmonic coordinates for character articulation[END_REF]] are then explored as they are always positive inside the cage. Unfortunately, computing a closed-form solution for harmonic Dirichlet is hard. Boundary element method (BEM) [START_REF] Sauter | Boundary element methods[END_REF] offers a way to only discretize the boundary of the domain and derive the solution by solving a linear system. Later, in order to avoid computing the integral over the boundary, method of fundamental solution (MFS) [START_REF] Li | Meshless Harmonic Volumetric Mapping Using Fundamental Solution Methods[END_REF][START_REF] Poullikkas | Methods of fundamental solutions for harmonic and biharmonic boundary value problems[END_REF]] has been proposed to use only the fundamental solution in computing the solution. [START_REF] Lipman | Green Coordinates[END_REF] introduce Green Coordinates for polyhedral cage in both 2D and 3D to extend barycentric coordinates, which leads to space deformation with a shape-preserving property. Ben-Chen et al. [2009b] propose a space deformation method that manipulates a sparse set of position and orientation constraints by finding an appropriate harmonic mapping of the space. They solve a non-linear minimization problem on the cage to find such a map that satisfies the user constraints.

In this paper, we employ harmonic mapping as our deformation function and compute the solution by a linear combination of component functions introduced in [START_REF] Li | Meshless Harmonic Volumetric Mapping Using Fundamental Solution Methods[END_REF][START_REF] Poullikkas | Methods of fundamental solutions for harmonic and biharmonic boundary value problems[END_REF]]. Thanks to MFS, our method is a meshless boundary method which does not require to enforce the connectivity of the mesh to function.

Deformation Transfer

Deformation transfer methods reuse existing animation sequences from the source shape to produce new animations. Given a source calibration shape, a target calibration shape, and a deformed source shape, the algorithm produces a deformed target shape that mimics the source shape deformation. Sumner and Popovic [2004] present an efficient method to transfer animations from a source shape to a target shape, establishing a bijection mapping between the triangles of the source and target shapes. The deformation gradient for each triangle in the source shape is then computed, which transforms a source triangle into its corresponding target triangle. One limitation is that only manifold mesh can be processed. [START_REF] Zhou | Deformation Transfer to Multi-Component Objects[END_REF] enhance this framework by connecting all adjacent components with a minimal spanning tree. In contrary, Ben-Chen et al. [2009a] project the deformation into a linear space by transfering a sparse set of deformation gradients inside a cage. However, their method requires constructing two polyhedral cages for both the source and target shape. [START_REF] Jin | Aura Mesh: Motion Retargeting to Preserve the Spatial Relationships between Skinned Characters[END_REF] propose the aura mesh, which is a volumetric mesh that surrounds a character's skin, to preserve the spatial relationships between two characters. The interaction between two characters is expressed with respect to the collision between the skin mesh of one character and the aura mesh of the other. The method is fast and produces interaction motions for characters with different sizes and skin shapes. However, their method requires a rigged skeleton and define the weights of the skin surface whereas our method just needs to specify several joints inside the characters to transfer the deformation descriptors.

Most of these deformation transfer methods just transfer the deformation descriptor of the source shape to the target shape without considering the spatial relationship between different body parts of the source shape. To solve this problem, we define a Context Graph on the surface of source/target shape, used to describe the spatial relationship among different body parts.

Compared to the Green coordinates model [START_REF] Lipman | Green Coordinates[END_REF]] that needs to move the cage to deform the bounding space, our method directly builds volumetric mappings between different 3D domains. Our method is somewhat similar to spatial deformations such as [Ben-Chen et al. 2009a,b], in the sense that we also use the Jacobian of the mapping function as the affine transformation descriptor. However, there are important differences: (1) Our method does not need a cage and we just need to compute some offset points close to the boundary surface of the pose using methods such as marching cubes and Silhouette Clipping [START_REF] Lorensen | Marching Cubes: A High Resolution 3D Surface Construction Algorithm[END_REF][START_REF] Sander | Silhouette clipping[END_REF]] and this is completely hidden to the users. In addition to the shape to be deformed, Ben-Chen et al. [2009a;2009b]

ALGORITHM OVERVIEW

As shown in Figure 1, our method consists of two phases: the deformation analysis and synthesis. During the deformation analysis process, if we do not have the skeleton together with mesh data, we first compute the joints of both the source and target calibration poses M 1 and Ω as shown in Figure 1(a) and (f), using skeleton extraction methods such as [START_REF] Oscar | Skeleton Extraction by Mesh Contraction[END_REF][START_REF] Cao | Point cloud skeletons via laplacian based contraction[END_REF]. Secondly, as shown in Figure 1(b), we generate an offset surface enclosing the source calibration pose M 1 . The points uniformly sampled on this offset surface are auxiliary points used for approximating the mapping function, and we call them singularity point. Then we approximate the harmonic function using the singularity points to map the source calibration shape to the deformed source shape using MFS [START_REF] Fairweather | The method of fundamental solutions for elliptic boundary value problems[END_REF][START_REF] Li | Meshless Harmonic Volumetric Mapping Using Fundamental Solution Methods[END_REF]]. We use the Jacobian matrices of the harmonic function at these joints as the affine transformation.

During the deformation synthesis process, similar to the deformation analysis step, we first place some offset points close to the target calibration shape as seen in Figure 1(g). Secondly, we assume that the joints of the target calibration pose Ω experience similar deformations than the joints of the source calibration pose M 1 such that we can calculate the coefficients for approximating the harmonic mapping function of the target calibration pose. Thirdly, in order to keep the spatial relationships similar to the deformed source pose M 2 and reduce the loss of volume, we introduce and construct a Context Graph, and employ the Laplacian deformation framework to further constrain the system. Finally, we synthesize the new posture by solving a linear system.

EXTRACTING DEFORMATION FROM SOURCE SHAPE HARMONIC MAP 4.1 Solving Harmonic Mapping Using MFS

Recent work [START_REF] Ferrari | Harmonic deformation of Delaunay triangulations[END_REF][START_REF] Joshi | Harmonic coordinates for character articulation[END_REF][START_REF] Zayer | Harmonic guidance for surface deformation[END_REF] shows that harmonic functions generate pleasing deformations. We therefore choose them as our deformation mapping function, and formulate our harmonic mapping problem as follows.

Let M 1 be an open region of R 3 with a smooth boundary ∂M 1 , and M 2 be another open region of R 3 with a smooth boundary ∂M 2 . We aim at constructing a mapping f : M 1 → M 2 between two boundaries with the same topology, that is, the objects have pairs of corresponding boundary surfaces, such that

∆f (P) = 0, P ∈ M 1 s.t . f (P) = P ′ , P ∈ ∂M 1 , P ′ ∈ ∂M 2 (1)
where the ∆ is the Laplacian operator defined as

∂f ∂x 2 + ∂f ∂y 2 + ∂f ∂z 2 .
Since the Laplacian operator ∆ is linear, we can construct a linear subspace of functions from R 3 to R 3 using harmonic mappings. Inside this linear space, We would like to find a mapping function that both preserves the local details and the spatial relationships. However, the Laplacian operator ∆f and its derivative are too complex to be solved analytically. We follow previous methods, such as [START_REF] Li | Meshless Harmonic Volumetric Mapping Using Fundamental Solution Methods[END_REF], and discretize the boundary using the Method of Fundamental Solutions (MFS) [START_REF] Fairweather | The method of fundamental solutions for elliptic boundary value problems[END_REF]]. Then the solution to Equation 1 can be approximated by a linear combination of component functions of harmonic equations:

f (P) = K s (P, Q)W s (2)
where

• P = [p 1 , p 2 , ..., p N s c ] is a matrix of 3 × N s c concatenating all N s
c surface points of source calibration shape ∂M 1 . • Q = [q 1 , q 2 , ..., q N s a ] is a 3 × N s a matrix consisting of all N s a three dimensional singularity points outside ∂M 1 . More details on placing q can be found in Section 4.3. • k(p, q) = 1/(4π ∥p -q∥ ) is the component function of the harmonic function f , which depends on the distance between a point p ∈ P and singular point q ∈ Q. To simplify the notation, we use K s = K s (P, Q) to represent the matrix of components function of size N s c × N s a and each row K s (p, Q) consists of N s a components for a point p ∈ P.

• W s = [w T 1 , w T 2 , ..., w T N s a
] is a N s a × 3 matrix of harmonic coefficients associated with these N s a singularity points, which offers the degree of the freedom to control the boundary fitting. The length of row vector w T i is 3 and each of its element corresponds to one dimension of point p ∈ P.

The vanishing Laplacian operator ∆ on f is enforced by the component functions k(p, q), we only need to ensure that the function satisfies the boundary condition. To perform boundary fitting on each surface point p ∈ ∂M 1 , we evaluate f (p) using Equation 2. The boundary constraints are then f (p) = p where p is the corresponding point of p laid on the boundary surface of deformed source pose. Enforcing these constraints on all collocation points reduces to a linear system:

A s W s = b s (3)
where A s is a N s c × N s a coefficient matrix with element A s i j = k(p i , q j ), N s c is the number of surface points p while N s a is the number of singularity points q, b is the image of p on the boundary of deformed source pose ∂M 2 with size N s c × 3. Therefore, the above fitting process reduces to a linear system on three axis directions. This is an over-determined linear least-squares problem as usually N s a ≪ N s c , which results in a closed-form solution [START_REF] Sorkine | Laplacian Surface Editing[END_REF]].

Calculating Descriptor of Affine Transformation

The Jacobian matrix of the mapping function can be used to prescribe an affine transformation for any point p ∈ (M 1 ∪ ∂M 1 ). Once we have approximated the source deformation f as a linear combination of component functions of Laplacian equation, the Jacobian of the deformation for any point p inside the domain M 1 can be computed using the gradients of the Laplacian equation. The transposed of the Jacobian of the source deformation of f at the point p is:

J s (p) T = D s W s , p ∈ ∂M 1 ∪ M 1 (4)
where D s is a matrix whose i th column is the derivative of the kernel and W s are linear coefficients of f computed in Equation 3. D s is a 3 × N s a matrix and W s is a N s a × 3 matrix. The Jacobian of the source deformation on each point p prescribes the affine transformation which we will transfer to the target calibration pose. There are many choices for the selection of Jacobian, e.g. we can choose any surface point's Jacobian to be transferred. However, as suggested by Ben-Chen et al. [2009b], points on the medial axis of the domain undergo the smallest distortion. On the other hand, densely sampling points inside source domain M 1 would result in high computation and even artifacts [START_REF] Zhou | Deformation Transfer to Multi-Component Objects[END_REF]. Therefore, similar to Ben-Chen et al. [2009a;2009b], we choose the Jacobian of each joint to represent the affine transformation.

Placing Singularity Points

We apply the MFS to solve the partial differential equation defined in Equation 1. The kernel function defined in Equation 2 is not well defined if a point p inside the shape equals to the singularity point q. On the other hand, the vanishing Laplacian operator cannot always be guaranteed if the singular points are sampled inside the calibration shape. Therefore, we have to place singularity points outside the surface of the calibration shape. The locations of the singularities are either preassigned or determined along with the coefficients w i of the component functions so that the approximate solution satisfies the boundary conditions as well as possible. If the locations of the singularities are to be determined, the resulting minimization problem is nonlinear and time-consuming. Therefore, in this paper, we preassigned the positions of singularity points and find the approximate solution by the least square fit of the boundary surface using Equation 1. In the current implementation, we use 260 singularity points on average.

We follow the method of [START_REF] Li | Meshless Harmonic Volumetric Mapping Using Fundamental Solution Methods[END_REF] to place singularity points uniformly on an offset surface ∂ M1 outside the boundary surface ∂M 1 as shown in Figure 1(b). Here is the process of how to sample singularity points: (1) We collapse edges of ∂M 1 until the Z. Liu et al.

desired number of faces is achieved but ensures that new vertices are placed outside all previous meshes [START_REF] Sander | Silhouette clipping[END_REF]]. (2) We uniformly sample some points Q = [q 1 , q 2 , ..., q N s a ] on the boundary points of ∂ M1 by furthest point relaxation as described in Fast Automatic Skinning Transformations [START_REF] Jacobson | Fast automatic skinning transformations[END_REF]]. Similarly, the collocation points are uniformly sampled from the boundary points of ∂M 1 .

PRODUCING NEW POSTURES BY FINDING APPROPRIATE HARMONIC MAP

Once the deformation descriptor has been extracted from the source deformation, we can apply it to the target calibration pose. We deform the target calibration pose with the same affine transformation (Jacobian constraints) extracted from Equation 4. First, we calculate our deformation descriptor of the harmonic functions of the target calibration pose, with a set of unknown coefficients w t . Then, we compute the deformed target calibration pose as a linear combination of component functions of Laplacian equation by Equation 6.

Computing a Raw Pose

Let the target calibration shape be given as a set of N t c points V = {v 1 , v 2 , ..., v N t c }, and corresponding singularity points U = {u 1 , u 2 , ..., u N t a }. Given the Jacobians J s at the m joints J s inside the deformed source pose, we want to minimize the Jacobian differences between the deformed source pose and newly synthesized target pose by:

E Af f ine = D t W t -J s 2 F (5)
where W t is the unknown linear coefficient of target calibration shape with size N t a × 3. D t = ∇K(J t , U ) is the gradient of target kernel function on m joints inside target calibration shape with size 3m × N t a and J t is a m × 3 matrix that concatenating positions of m joint inside target calibration shape. J s is a 3m × 3 dimensional matrix concatenating matrices of all m affine transformations, computed by Equation 4, of the source shape deformation.

This optimization problem for computing the coefficients W t can be solved by a linear system. Equipped with the coefficients W t , the deformed target pose is given as a linear combination of the component function of the harmonic function for the target calibration pose as:

V = K t W t (6)
where

K t = K t (V , U ) is a matrix of size N t c × N t
a that each row computes the inverse of the distance between the point v ∈ V and all singularity points u i .

Constraints

The problem defined in Equation 5 is under-determined because the number of joints m is much smaller than the number of singularity points N t a . On the other hand, we need to predetermine a position for any vertex to compute a unique solution. Lastly, we introduce the bone length constraint such that the resulting poses are physically plausible. More details on solving the system with constraints are discussed in Section 7.2.

Surface Smoothness Constraint

Here, we follow [Ben-Chen et al. 2009b] to assume that adjacent points on the surface will be transformed similarly. Thus, the second derivative of the mapping function tends to be very small. In order to make the surface as smooth as possible, we constrain the Hessian matrix H (f t ) of the target mapping function f t = K t W t for a newly synthesized surface to be as small as possible by:

E Smooth = H (K t W t ) 2 F (7)
where E Smooth is a linear function of W t with the coefficients of size 3N t c × N t a .

Positional Constraint

The gradient of a deformation function can only specify the rotation and scaling transformations. Thus, to determine all the absolute positions, we predetermine the position of any point, e.g., v 1 , as v posit ion by:

E Posit ion = K t (v 1 , U )W t -v posit ion 2 (8)
where K t (v 1 , U ) is a row vector of length N t a and U are singularity points of target calibration shape.

Bone Length Constraint

In an ideal character animation, the bone length between two joints is fixed. We introduce the bone length constraint to avoid stretching the synthesized shape.

Specifically, the bone length l i j between joint J t i and J t j for target calibration shape is prescribed by ( J t i -J t j -l i j ) 2 . We compute the Jacobian Linearization of this non-linear system [START_REF] Cao | Point cloud skeletons via laplacian based contraction[END_REF] by:

( K t i W t -K t j W t -l i j ) 2 (
E Bone = J Bone W t + C (10)
where J Bone is the Jacobian matrix of size 3B × N t a , B is the number of bones, and C is a constant matrix.

CONSTRUCTING CONTEXT GRAPH

In this Section, we compute the Context Graph to preserve the volume and spatial relationship between body parts. We first construct two compatible Context Graphs and then simplify the graphs to speed up the computation. Lastly, we deform the graph constructed on the target calibration pose into the one built on the deformed source pose such that we preserve the spatial relationships.

We project some singular points onto the boundary of source calibration pose M 1 to construct a graph such that spatial relationship of a particular pose could be preserved. Here, we project the N s a singular points of M 1 onto the surface P ∈ ∂M 1 , denoted as P a = {p 1 , p 2 , ..., p N s a }. We connect all the points P a on each deformed source pose to build a fully connected Context Graph G s , which enable us to capture the spatial relationships among points on the boundary surface of deformed source pose.

We use the non-rigid iterative closest point (ICP) algorithm [START_REF] Li | Global Correspondence Optimization for Non-Rigid Registration of Depth Scans[END_REF]] to find N s a corresponding points V a ∈ ∂Ω on the boundary of target calibration poses Ω. Connecting all the points V a found on ∂Ω would result in a compatible graph G t as the one G s built on the deformed source pose.

Pruning the Context Graph

The fully connected Context Graph G s contains too many redundant edges that will greatly slow down the deformation. Thus, we need to prune the Context Graph, which we do according to the following rules:

(1) As we can consider human poses to be close to rigid, the distance between any two given points on the same body part rarely changes during the deformation. For this reason, edges connecting two points laid on the same body-part are discarded, and we therefore only consider points distributed on different body parts. To do so, we apply K-means to cluster the body parts by the joints extracted either from the users or the skeleton joint extraction algorithms such as [START_REF] Oscar | Skeleton Extraction by Mesh Contraction[END_REF]].

(2) As we focus on close interactions and assume that nodes which are far away from each other on the graph have little influence on the others, we remove edges longer than a given threshold. Currently, we set the threshold to 50% of the calibration shape height.

(3) As we only consider interactions on the surface of body parts, edges passing through any body part would be invalid.

Edges intersecting with the interior of the mesh are therefore discarded. Figure 2 shows an example of constructing and pruning Context Graph.

Deformation of Context Graph

There are many choices to combine the existing differential deformation frameworks (reviewed in Section 2). For a clear presentation, we follow in this paper the method of [START_REF] Sorkine | Laplacian Surface Editing[END_REF]] to compute Laplacian coordinate, which implicitly computes the affine transformations from the translation of a point.

The Laplacian coordinates of points P a are defined as L(P a ) = LP a = (I -Diaд -1 Adjacent)P a with the uniform weights, where Adjacent is the adjacent matrix of Context Graph, Diaд is the degree matrix, and L is the N s a × N s a Laplacian operator matrix constructed from the Context Graph of deformed source mesh. The deformation of target Context Graph that preserves spatial relationship is computed by:

E S pat i al = T Spat i al L Spat i al (P a ) -L Spat i al ( Va ) 2 F (11)
where points Va ∈ ∂Ω are the node of spatial Context Graph that correspond to P a ∈ ∂M 1 , which are extracted from the surface of target calibration shape Ω, and L Spat i al (P a ) are the Laplacian coordinates computed from the Context Graph of the source deformed shape. The coefficients of T Spat i al are the linear combinations of unknown W t [START_REF] Sorkine | Laplacian Surface Editing[END_REF]], which implies E Spat i al is a linear function of the unknown W t with size N s a × N t a . Similarly, we define the deformation of volume Context Graph that preserves the volume of the target calibration shape by: where V is a N t c × 3 matrix that concatenating the position of point v ∈ ∂Ω. Similar to T S pat i al , the coefficients of T V ol ume are linear functions of unknown W t with size N t c × N t a . L V ol ume ( V ) are the Laplacian coordinates computed from the volume Context Graph of the target calibration shape

E V ol ume = T V olume L V olume (V ) -L V olume ( V ) 2 F (12)

IMPLEMENTATIONS 7.1 Solving the Linear System for Extracting Source Deformation

As discussed in Section 4, the deformation analysis process can be optimized by solving a linear system A s W s = b s . Each row in the coefficient matrix A s is the value of the kernel functions evaluated on a collocation point and all the other singularity points. However, the resulting matrix A s is not sparse because the value of kernel functions is almost never zero. Thus, we cannot apply a sparse solver to speed up the computation. As studied in [START_REF] Li | Meshless Harmonic Volumetric Mapping Using Fundamental Solution Methods[END_REF][START_REF] Palghat | Method of fundamental solutions: singular value decomposition analysis[END_REF], such a matrix might be ill-conditioned and will not produce a stable solution using regular linear system solvers such as Gaussian elimination. We follow previous methods of [START_REF] Li | Meshless Harmonic Volumetric Mapping Using Fundamental Solution Methods[END_REF][START_REF] Palghat | Method of fundamental solutions: singular value decomposition analysis[END_REF] to use Singular Value Decomposition (SVD) because it approaches accurate and stable results even when the coefficient matrix is highly ill-conditioned. The matrix A s is then decomposed as A s = USV T and the inverse of A s is computed as A s -1 = VS -1 U T . Once we have precomputed the inverse of A s , we are able to perform a fast fitting of new boundary condition b s * (new pose in the deformed source pose sequence) by:

W s * = A s -1 b s * = VS -1 U T b s * (13)
Then the transpose of the deformation descriptor of joints J s ∈ M 1 that fits a new boundary b s * is computed by: J s (J s ) T = D s W s * (14) 7.2 Solving the Linear System for Synthesizing a New Pose

We define the retargeting problem as the sum of Equation 5, 7, 8, 11, and 12 by:

E = E Af f ine + λ Smoot h E Smoot h + λ P osi t ion E P osi t ion +λ Bone E Bone + λ Spat i al E Spat i al + λ V olume E V olume (15)
where λ P os i t ion , λ Smoot h , λ Bone , λ Spat i al , and λ V olume are the weights for the energy terms. In our implementation, they are empirically set to be 0.01, 10, 1, 0.8, and 1, respectively. There are some principles to tune the values of the weights. The bone length and volume constraints are hard because we do not want to change the length of bone or lose the volume of the shape during the deformation. The position constraint makes the linear system solvable and it should be with a large weight. The smooth constraint should be small, otherwise, the shape would be flattened out. The spatial constraint should be large such that the spatial relationship between body parts can be preserved We rewrite Equation 15 as a linear system by:

A t W t = b t (16
) where W t is the unknown coefficient matrix of size N t a × 3, A t and b t represent a concatenation of matrices from Equation 5, 7, 8, 11, and 12 as:

A t =            A Af f ine λ P os i t ion A P osi t ion λ Bone J Bone λ Smoot h A Smoot h λ S pat i al A Spat i al λ V ol ume A V olume             , b t =         J s λ P osi t ion v posi t ion -λ Bone C ì 0         (17)
where b t is the right hand side matrix containing the affine transformation J s , position v posit ion , and bone length C constraints. The remaining elements are all 0. A Af f ine = D t (J t ) is a coefficient matrix of size 3m × N t a defined in Equation 5, where m is the number of joints J t and N t a is the number of singularity points of the target calibration pose. A P osit ion = K t (v 1 , U ) is a row vector of length N t a defined in Equation 8. A Smoot h is a Hessian matrix of size 3N t c × N t a defined in Equation 7, where N t c is the number of boundary points on the target calibration pose. A S pat i al and A V ol ume are the linear coefficients defined in Equation 11 and 12 that preserve the spatial relationship and volume for the synthesized pose. J s is the Jacobian matrix of size 3m × 3 from m joints of deformed source pose defined in Equation 5 and v posit ion is a row vector of length 3.

The analytical least-squares solution is defined by the following normal equation:

W t = ((A t ) T A t ) -1 (A t ) T b t (18)
Once we calculate the unknown coefficient W t , we produce a new pose by using Equation 6.

EXPERIMENTAL RESULTS

Harmonic Mapping Error

In order to quantitatively analyze the correctness of the proposed method, we use the source calibration shape as the target calibration shape. Thus, we treat the source animation as our ground truth data. We compute the average approximation error for each pose from all the animation sequences by:

Error Mappinд (P) = 1/m m i=1 f i (P) -Pi ( 19
)
where function f i maps each vertex P ∈ ∂M 1 of i th frame onto deformed source pose Pi and m is total number of frames.

To simplify the calculation of reconstruction error, we normalize each calibration shape to a height of 1. Figure 3 shows the resulting error for a set of punching motions. Each pose contains around 18K vertices. As shown in Figure 3b-3d, the reconstruction error is coded with color from blue (0) to red (0.04). We can observe that our method generates meshes close to the ground truth meshes. Note that, as illustrated in Table 1, the average error per pose over the source animation sequence is very small, which shows that the approximated harmonic function is able to map from one domain onto another one precisely.

Retargeting Results

All experiments were performed on a PC with an Intel core i7-7700HQ CPU and with 16GB of RAM. Our method takes around 4 seconds to synthesize each frame. The deformation analysis takes around 1 second. We first demonstrate the influence of different terms on producing the resultant poses. As shown in Fig. 5, bone length constraint defined in Equation 10 ensures that the deformation of the target calibration shape is physically correct, while preserving the volume during the deformation prevents unnatural morphological changes. As illustrated in Fig. 4(b), the belly in blue color is collapsed if we do not add the volume constraint. In addition, by preserving the Context Graph, the spatial relationship between the legs is more similar to the original source shape as shown in Fig. 4(d).

We test the proposed method on various morphologies. As shown in Fig. 6, we compare the performance of our method to the Linear Blend Skinning using the same skin weights as the source animation, and the results created by an artist. We can see that our method can capture the contextual meaning of the source animation, and the resulting poses are similar to the artist's manual corrections. Serious artifacts can result if we apply the skin weights of source animation to a new target calibration shape directly. As illustrated in Fig. 6(f), the whole body is twisted and interpenetrated due to the change of morphology and the lack of contextual information. Thus, it is inefficient to produce retargeting animation based on the naive Linear Blend Skinning.

Although our method generates sensible retargeting poses, we can not totally avoid interpenetrations because we use a sparse set of nodes to construct the Context Graph. As shown in Fig. 6(d), the right hand may penetrate the left hand due to insufficient nodes of Context Graph to capture the spatial relationship. On the other hand, due to the limitation of the bone length and the avoidance of interpenetration, the resulting pose can be slightly different from the deformed source pose. As shown in the fourth example of Fig. 6(d), the position of the two hands relative to the head is a little different from the one in the deformed source shape. Such a kind of deformation is common if the morphology of the target calibration is extremely different from the source calibration shape.

CONCLUSIONS

We proposed a new context-aware framework for surface based motion retargeting, which can be applied to both manifold and non-manifold meshes. Our algorithm works in two stages: (1) We learn the affine transformation descriptor from the source animation sequences. (2) We produce a new deformation of the target calibration shape that respects the affine transformation descriptor and spatial relationships among different body parts. The contextual meaning is preserved by a new structure called the Context Graph. We solve a linear system to compute the deformation of the target shape such that the resulting motions mimic the source deformation.

Currently, in order to preserve the spatial relationship, our method requires the source and target calibration shapes to be in the same pose. While it makes sense to assume the source and target calibration poses are similar in motion retargeting area, a typical target calibration shape provided by the user could have arbitrary pose. Given an arbitrary target calibration shape, we may have two possible solutions: (1) We deform the given target calibration pose such that it has a similar pose as source using shape deformation methods such as [Ben-Chen et al. 2009a;[START_REF] Lipman | Green Coordinates[END_REF]; Wang

  (b) Place Singularity Points q in Red (g) Place Singularity Points u in Red Deformation Analysis Deformation Synthesis (c) Deformed Source Pose M 2 (a) Source Reference Shape M 1 with Joints l S in Green (f) Target Reference Shape Ω with Joints l T in Green (d) Approximate Harmonic Mapping using MFS Compute Deformation Descriptor for Joints l S (e) Compute Spatial Relationship (h) Newly Synthesized Target Pose Synthesize Deformed Target Pose by Solving a Linear System Preserve Spatial Relationship (i) Output

Figure 1 :

 1 Figure1: Overview of synthesizing new poses using harmonic mapping. During the deformation analysis process, we calculate the deformation descriptor of the source calibration pose. During the deformation synthesis process, we compute the coefficients of approximating harmonic mapping by assuming that the target calibration pose perform similar deformation to the source calibration pose. We further constrain the synthesized motion to have an as similar spatial relationship as that of the source deformed pose by adding Laplacian framework.

  Figure 2: Pruning the Context Graph. (a) Source calibration pose with joints. (b) Clustering of body-parts. (c) Projecting singular points onto the boundary surface and extracting them. (d) Removing edges on the same body part. (e) Removing long edges. (f) Removing edges interpenetrating the body parts. (g) The simplified Context Graph.

  Figure 3: Reconstruction error of harmonic mapping computed as per pose for 371 frames of punching motion. We map the source calibration pose (a) with 18K vertices onto the other punching motion (b-d). The coded color indicates the scale of approximation error from 0 (blue) to 0.04 (red). The reconstruction error curve shows that most of reconstructed vertices are very close to ground truth points (e).

Figure 4 :

 4 Figure 4: Influence of constraints preserving the volume and using the Context Graph: Source shape (a) and Target shape deformed using our method without volume preservation (b), with volume preservation (c), and with both volume and Context Graph preservation (d).

Figure 5 :

 5 Figure 5: Our method with (a) and without (b) bone length constraints. Overlay of both (c).

Figure 6 :

 6 Figure 6: The 'tree pose' retargeting comparison: source poses (a), retargeting results produced by our method (b-d), by an artist (e), and using traditional Linear Blend Skinning (f).

Table 1 :

 1 Average reconstruction errors over different motion sequences

	Motion Name	Number of Frames	Number of Vertices	Mean Error Std
	Punching	371	17996	0.021	0.14
	Sitting on a chair	490	19296	0.125	0.26
	Tree pose	329	19296	0.193	0.21


 Abstract-Pathological gait assessment and assistive control based on functional electrical stimulation (FES) in post-stroke individuals, brings out a common need to robustly quantify kinematics facing multiple constraints.

This study proposes a novel approach using inertial sensors to compute dorsiflexion angles and spatio-temporal parameters, in order to be later used as inputs for online close-loop control of FES. 26 post-stroke subjects were asked to walk on a pressure mat equipped with inertial measurement units (IMU) and passive reflective markers. A total of 930 strides were individually analyzed and results between IMU-based algorithms and reference systems compared. Mean absolute (MA) errors of dorsiflexion angles were found to be less than 4°, while stride lengths were robustly segmented and estimated with a MA error less than 10 cm. These results open new doors to rehabilitation using adaptive FES closed-loop control strategies in "foot drop" syndrome correction.

Being able to observe and assess a pathological motion is a recurrent need for practitioners. After a stroke, individuals are often hampered by walking difficulties [START_REF] Woolley | Characteristics of Gait in Hemiplegia[END_REF].

Step-by-step kinematic and spatio-temporal parameters have been shown to be clinically relevant markers of impaired walking performances in post-stroke subjects [START_REF] Balasubramanian | Variability in spatiotemporal step characteristics and its relationship to walking performance post-stroke[END_REF]. Partly related to a deficit of voluntary control, individuals with hemiparesis are often unable to perform complete ankle dorsiflexion during swing phase, resulting in poor foot clearance, slow gait speed and compensatory strategies while walking [START_REF] Chen | Gait differences between individuals with post-stroke hemiparesis and non-disabled controls at matched speeds[END_REF]. Between 20% and 50% of individuals with post-stroke hemiplegia would be affected by this commonly called "foot drop" syndrome [START_REF] Kralj | Enhancement of hemiplegic patient rehabilitation by means of functional electrical stimulation[END_REF]. Among different rehabilitation approaches, one conventional treatment consists in using functional electrical stimulation (FES) applied to the peroneal nerve in order to elicit ankle dorsiflexion [START_REF] Everaert | Effect of a foot-drop stimulator and ankle-foot orthosis on walking performance after stroke: A multicenter randomized controlled trial[END_REF]. Classically, the stimulation is triggered by a heel switch inserted in the paretic side shoe (e.g. Odstock © ODFS III) or by a tilt angle sensor (e.g. WalkAide © , Bioness L300 © ). Novel strategies to real-time modulate FES intensity or to regulate dorsiflexion in the presence of muscle state changes, such as fatigue and spasticity, could increase treatment efficiency and the number of potential users [START_REF] Coste | Continuous gait cycle index estimation for electrical stimulation assisted foot drop correction[END_REF]. This brings out a common need of accurately quantifying kinematics in post-stroke F.Feuvrier and J.Froger are with Nîmes University Hospital (CHU), Nîmes, 30000 France (e-mail: jerome.froger@chu-nimes.fr). subjects, whether for assessing gait performance or for artificially controlling motion. Nevertheless, few systems are able to perform such ambulatory observations with a sufficient accuracy. Indeed, wearable sensors have been increasingly explored over the past few years as a mobile gait analysis solution. Numerous authors have suggested methods involving Inertial Measurement Units (IMU) sensors in order to estimate gait temporal or spatio-temporal parameters [START_REF] Caldas | A systematic review of gait analysis methods based on inertial sensors and adaptive algorithms[END_REF] as well as to compute joint angles during walking [START_REF] Seel | IMU-Based Joint Angle Measurement for Gait Analysis[END_REF]. Inertial sensors basically measure acceleration, magnetic field and angular rate in their local coordinate system. Some devices directly provide the estimates of the sensor's orientation (in the form of quaternions, Euler angles or rotation matrices) with respect to a global fixed frame (e.g. Earth frame) by fusing sensors data. This is usually done by strap-down integrating the gyrometer values to obtain a first orientation, used then to remove Earth gravity from acceleration and combined with magnetometer measurements to provide azimuth [START_REF] Savage | Strapdown Inertial Navigation Integration Algorithm Design Part 1: Attitude Algorithms[END_REF]. Position can be afterward determined by double integrating acceleration without gravity. However, numerical integration introduces a drift error, while magnetic measurement can be disturbed by ferromagnetic materials, resulting in significant errors in the estimation of IMU attitude and position. During walking, a common solution is to update drift integration by simply zeroing the velocity during each detected stance phase, also called Zero Velocity Update (ZUPT) [START_REF] Mariani | 3D gait assessment in young and elderly subjects using foot-worn inertial sensors[END_REF] or by using Attitude Heading and Reference System (AHRS) fusion methods proposed in the literature [START_REF] Cavallo | Experimental Comparison of Sensor Fusion Algorithms for Attitude Estimation[END_REF]. This requires to accurately segment gait cycle by detecting gait events (initial contact, toe off, etc…) from IMU data [START_REF] Rueterbories | Methods for gait event detection and analysis in ambulatory systems[END_REF]. Only a few of the above mentioned methods have been compared to the gold standards (i.e. optical motion capture system or instrumented mats). In addition, they usually segment gait phases using angular velocity patterns, which is satisfying for normal gait [START_REF] Sijobert | Implementation and Validation of a Stride Length Estimation Algorithm , Using a Single Basic Inertial Sensor on Healthy Subjects and Patients Suffering from Parkinson ' s Disease[END_REF] but tends to fail when gait is impaired and subjects develop compensatory strategies [START_REF] Yang | Estimation of spatio-temporal parameters for post-stroke hemiparetic gait using inertial sensors[END_REF]. In this study, we propose a novel approach dedicated to pathological motion assessment from IMUs in 26 participants with post-stroke hemiplegia, integrating two aspects: 1) robust stance phase detection based on acceleration and angular rate combination and 2) estimation of joint angles based on an AHRS algorithm and gravity cancellation for reconstructing 3D trajectory of individual steps. The method has been validated by estimating 4 parameters initially defined as relevant by the practitioner to monitor rehabilitation progresses and to later adapt stimulation parameters: dorsiflexion angle at initial contact and at mid-swing instants, stride length and speed velocity.

A sensor fusion approach for inertial sensors based 3D kinematics and pathological gait assessments: toward an adaptive control of stimulation in post-stroke subjects*

B. Sijobert, F.Feuvrier, J.Froger, D.Guiraud, Member, IEEE EMBS and C.Azevedo Coste

II. MATERIAL AND METHODS

A. Experimental Protocol 1) Instrumentation Subjects were equipped with 4 IMUs (Fox HikoB © Villeurbanne, France, L45mm x W36mm x H17mm, weight: 22g) featuring a 3-axis accelerometer, a 3-axis magnetometer and a 3-axis gyrometer respectively mounted on the feet and shanks. Each IMU was strapped on a rigid support together with 4 reflective markers (Fig. 1a) tracked by an optical motion capture system (OMCS, Vicon © Bonita MX) which cameras were installed along a Gaitrite © (CIR System Inc) walkway system (Fig. 1b). 2) Subjects This study included 29 subjects (mean 58.5 ± 10.4 years old; 9 females) after supratentorial ischaemic or hemorrhagic stroke, presenting a foot-drop, able to walk 10 meters without human help, with or without a walking stick. The protocol was approved by a national ethical committee (CPP) and by the local ethical committee of the University Hospital (CHU Nimes, France, RCB 2015-A00572-47), all subjects provided informed consent prior to the experiment.

3) Stimulation In order to facilitate walking and elicit measurable dorsiflexion angles, subjects were also equipped with a classical "foot drop" stimulator (Odstock © ODFS III) configured at the beginning of the experiment. Stimulation was triggered at heel off using a heel switch inside the shoe. Two electrodes of 23 cm² delivered the stimulation either to the peroneal nerve or directly to the tibialis muscles of the affected side.

4) Protocol conduct In order to best match Emory functional Ambulation Profile (E-FAP) test, participants were asked to walk five meters at a comfortable self-selected speed on the gait mat, to turn at the end of the carpet and walk back to their initial position. In case of technical issues or data loss, each trial was repeated three times to record at least one set of data by subject. IMUs, OMCS and gait carpet were synchronized at a hardware level via a trigger sent by the Gaitrite © to all the acquisition systems.

B. Algorithms

In the following formulas we use the quaternion notation defined in [START_REF] Hamilton | Lectures on Quaternions : containing a systematic statement of a new mathematical method[END_REF] as four scalar numbers, one real dimension and an imaginary (or also called vector) part . The norm of a quaternion is defined as:

. A quaternion with a unity norm is called unit quaternion. In the following algorithms, quaternions were systematically normalized as: to be able to use unit quaternion properties. Let us note the quaternion describing orientations of B relative to F, for a united quaternion, q 𝐵 𝐹 equals q ̅ 𝐹 𝐵 , with q ̅ defined as the quaternion conjugate : q ̅ = (w, -x i, , -y j , -z k ).

To compute quaternions from IMUs, AHRS fusion methods are usually based on nonlinear observers [START_REF] Mahony | Nonlinear Complementary Filters on the Special Orthogonal Group[END_REF], [START_REF] Martin | Design and implementation of a lowcost observer-based attitude and heading reference system[END_REF], on a gradient descent method [START_REF] Madgwick | Estimation of IMU and MARG orientation using a gradient descent algorithm[END_REF] or on estimators based on Extended Kalman Filters (EKF). In a previous work [START_REF] Sijobert | Validation of inertial sensor fusion algorithms in the context of hemiparetic gait assessment[END_REF] we compared three different fusion algorithms to compute attitude from inertial measurement units and found out Martin et al. [START_REF] Martin | Design and implementation of a lowcost observer-based attitude and heading reference system[END_REF] was offering the best compromise between computational cost and reliability over time and magnetometer disturbances. Therefore, each quaternion representing the 3D orientation in global coordinate system was computed from 200Hz sampled magneto-inertial data using Martin et al. nonlinear observer. To reduce computational error, sensors were individually calibrated on each day of experimentation using Frosio et al. calibration model [START_REF] Frosio | Autocalibration of MEMS Accelerometers[END_REF] for fitting in static conditions Earth magnetic field and gravity, incorporating not only bias and scale factor for each axis but also cross-axis symmetrical actors calculated through Gauss-Newton nonlinear optimization.

1) Dorsiflexion angles Let us defined and the orientation of the shank and foot sensors in the global (i.e. ground) reference frame. The dorsiflexion angle corresponds to , the quaternion expressing rotations from shank to foot relative to the shank frame 𝑆 𝐹 and computed using the following formula (Equ. 1), based on quaternion unit properties (Fig. 2):

q 𝑠ℎ𝑎𝑛𝑘-𝑓𝑜𝑜𝑡 𝑆 𝐹 = q 𝑑𝑜𝑟𝑠 = q ̅ 𝑠ℎ𝑎𝑛𝑘 𝐺 * * q 𝑓𝑜𝑜𝑡 𝐺 [START_REF] Woolley | Characteristics of Gait in Hemiplegia[END_REF] with * * the specific quaternion multiplication, known as Hamilton product [START_REF] Hamilton | Lectures on Quaternions : containing a systematic statement of a new mathematical method[END_REF]. At the beginning of each trial, subjects were asked to stand still during five seconds for getting rid of potential biases while defining zeros between IMU system and VICON system and for enabling Martin et al. algorithm to converge to an initial attitude.

2) Trajectory Let us define and , the position and velocity of the feet sensors, respectively calculated by integrating twice and once the recorded acceleration without gravity . Given , computed from [START_REF] Martin | Design and implementation of a lowcost observer-based attitude and heading reference system[END_REF] and unit quaternion properties , we define the following formula (Equ. 2) to rotate acceleration vector and remove gravity 𝑔: q 𝑎𝑐𝑐 𝐺 = (0, 𝑎 𝑥 𝐺 , 𝑎 𝑦 𝐺 , 𝑎 𝑧 𝐺 ) = q 𝐵 𝐺 * * q 𝑎𝑐𝑐 𝐵 * * q ̅ 𝐵 𝐺

Therefore, to compute foot trajectory on each stride (Equ.

3):

where The method described in [START_REF] Sijobert | Effects of Sensitive Electrical Stimulation-Based Somatosensory Cueing in Parkinson's Disease Gait and Freezing of Gait Assessment[END_REF] has been used to detect gait events, not only based on sagittal angular rate but also taking into account acceleration measurements.

Stride length was then defined as the distance between two consecutive initial contacts, computed in the transverse plane (Equ. 4):

III. RESULTS

29 subjects took part to the protocol. 3 subjects were finally too weak to perform any trial. Therefore the analysis and results refer to 26 subjects. A total of 930 strides (457 on paretic side and 473 on healthy side) were recorded. For each stride, the spatio-temporal parameters were computed from IMUs and compared to parameters estimated from gold standard devices. Dorsiflexion angles were compared between IMUs and Vicon © at mid-swing and heel-on, based on gait events instants extracted from Gaitrite © . Table 1 shows the results for all the analyzed strides: a Root Mean Square (RMS) error between OMCS and IMUs estimations of dorsiflexion of 5.51° at initial contact and 5.01° at midswing and a Mean Absolute (MA) error of respectively 3.39° and 3.74°. The average dorsiflexion angle in all participants was 8.76° at initial contact and 9.32° at mid-swing. We observed a RMS error of 12.64 cm and a MA error of 9.84 cm regarding stride length estimations between Gaitrite © and IMUs, and a RMS error of 6.17 cm/s and MA error of 5,06 cm/s for gait speed computation. The average stride length in all participants was 57.49 cm and the average gait speed 30.36 cm/s. IV. DISCUSSION -CONCLUSION

In this study, we compared gait variables estimated from IMUs data and gold standard systems. The initial aim was to analyze the reliability and accuracy of using IMUs to analyze motion, for later considering them as FES-based assistive control inputs. In our case, assessing pathological motion has been challenging and has required multiple non optimal choices from the data collection to the algorithm design. OMCS data acquisition suffered from many data losses and artefacts due to marker occultation by cane or markers breakage. As a result, we had to use rigid objects instead of a complete set of markers. Using a simple single strap, they were approximately positioned in sagittal plane, in order to improve visibility and ensure an eased installation time. An important advantage of our approach is the absence of specific procedure and preparation to locate sensors and use them. Except of a 5s static posture at the beginning of each trial, no calibration motions [START_REF] Seel | IMU-Based Joint Angle Measurement for Gait Analysis[END_REF] were requested from the participants and no manual measurements of body dimensions had to be done nor individual adaptation. To estimate stride length and velocity for each step in each trial, feet trajectories have to be computed in a global frame. Usually, the solution consists in double integrating 2D linear acceleration in sagittal plane with an angular rate based gait cycle segmentation [START_REF] Sijobert | Implementation and Validation of a Stride Length Estimation Algorithm , Using a Single Basic Inertial Sensor on Healthy Subjects and Patients Suffering from Parkinson ' s Disease[END_REF], [START_REF] Benoussaad | Robust Foot Clearance Estimation based on Acceleration Integration of Mounted Inertial Measurement Unit[END_REF] and an angular rate integration to estimate orientation needed to gravity removal. In our study, pathological gait was often associated to compensatory strategies (e.g. circumduction walk) and slow motions. These existing methods were not applicable to assess the gait of post-stroke subjects with a complex forward swing. Therefore, we had to adapt an algorithm we developed in an earlier work for Parkinson's disease [START_REF] Sijobert | Effects of Sensitive Electrical Stimulation-Based Somatosensory Cueing in Parkinson's Disease Gait and Freezing of Gait Assessment[END_REF], to segment impaired gait cycle and we took advantage of Martin et al. [START_REF] Martin | Design and implementation of a lowcost observer-based attitude and heading reference system[END_REF] quaternion computation to accurately remove gravity. Only a few studies have been conducted in hemiplegic participants in the literature [START_REF] Trojaniello | Estimation of step-by-step spatio-temporal parameters of normal and impaired gait using shank-mounted magneto-inertial sensors : application to elderly , hemiparetic , parkinsonian and choreic gait[END_REF], thus we mainly compared our results with publications on healthy individuals [START_REF] Yang | Inertial Sensor-Based Methods in Walking Speed Estimation: A Systematic Review[END_REF]. They seem in accordance with a MAE less than 4° for the dorsiflexion angle and less than 10 cm for stride length estimation.

One challenge of our study was to estimate very small dorsiflexion angles (about 10° in all participants). To be able to compute dorsiflexion angle with an error of 4° is a sufficient accuracy to detect the dorsiflexion tendency to decrease in the presence of muscle fatigue and to adapt FES parameters for counteracting fatigue effects. Furthermore, we have designed the algorithms to be straight forward implemented for online use, as they intend to be used for FES control. To shorten calculation time and be able to track any orientations without singularities, we chose a quaternion representation computed from a low-cost observer-based attitude and heading reference system [START_REF] Martin | Design and implementation of a lowcost observer-based attitude and heading reference system[END_REF].

Euler angles well known gimbal lock may not be an issue when dealing with 3D rotations, but in our protocol it could have been a serious limitation when tracking lower limbs movements.

Finally, as we compared three different systems outputs (IMU, Gaitrite © and Vicon © ) at different sampling rates but at a similar time mark (e.g. mid-swing), synchronization had to be accurately done not to introduce additional error. These results break new grounds toward adaptive online control of the dorsiflexion. The complete knowledge of 3D trajectory and attitude of each stride could enable new closeloop strategies, such as adapting stimulation to obstacles (Fig. 3), stairs, or fatigue based on joint angle estimation.

A

  sensor fusion approach for inertial sensors based 3D kinematics and pathological gait assessments: toward an adaptive control of stimulation in post-stroke subjects Benoît Sijobert, François Feuvrier, Jérôme Froger, David Guiraud, Christine Azevedo Coste To cite this version: Benoît Sijobert, François Feuvrier, Jérôme Froger, David Guiraud, Christine Azevedo Coste. A sensor fusion approach for inertial sensors based 3D kinematics and pathological gait assessments: toward an adaptive control of stimulation in post-stroke subjects. EMBC 2018 -40th International Conference of the IEEE Engineering in Medicine and Biology Society, Jul 2018, Honolulu, HI, United States. pp.3497-3500, 10.1109/EMBC.2018.8512985. hal-01924841 I. INTRODUCTION

Figure 1 .

 1 Figure 1. Each subject is equipped with IMUs strapped on rigid support together with 4 reflective markers (a) tracked by an optical motion capture system (Vicon © ). Subject is instructed to stand still for 5 seconds and then walk five meters on an instrumented walkway mat (Gaitrite © ), turn at the end of the carpet and walk back to initial position (b).

Figure 2 .

 2 Figure 2. Dorsiflexion angle computation. The Hamilton product of the shank unit conjugate quaternion in global (Earth) frame by the quaternion of the foot also expressed in global frame gives the quaternion corresponding to the rotation of the foot relative to the shank.

Figure 3 .

 3 Figure 3. 3D trajectory of the foot IMU from a post-stroke subject obtained by integrating 3D accelerations with gravity removed from quaternion AHRS estimation. A Zero Velocity Update (ZUPT) is performed at each beginning and ending of the strides. In this example, the fourth stride corresponds to an obstacle crossing put on the gait mat. (x-IMU based Matlab 3D plotting : http://x-io.co.uk/gait-tracking-with-x-imu/)

Table 1 .

 1 Dorsiflexion angles RMS (root mean square) and MA (mean absolute) error at heel on and mid-swing, between IMU computation and Vicon. Spatio-temporal parameters compared between Gaitrite and IMUs computation.

		RMSE	MAE (±SD)	RMSE	MAE (±SD)
	Dorsiflexion at heel on (initial contact)	Stride length (cm)
	All (S=930)	5.51°	3.39° (±3.37)	12.64	9.84 (±7.94)
	Paretic side (S=457)	4.89°	3.12° (±3.17)	12.61	9.86 (±7.87)
	Healthy side (S=473)	6.06°	3.68° (±3.56)	12.67	9.81 (±8.02)
	Dorsiflexion at mid-swing		Speed (cm/s)
	All (S=930)	5.01°	3.74° (±3.83)	6.17	5.06 (±3.45)
	Paretic side (S=457)	4.89°	3.73° (±3.72)	6.22	5.08 (±3.47)
	Healthy side (S=473)	6.10°	3.75° (±3.95)	6.11	5.06 (±3.43)


Imen Zomaya

HAL is

Parallel Optimization using/for Multi and Many-core High Performance Computing N. Melab a, , A. Zomaya b, , I. Chakroun c, a Université Lille 1 / INRIA Lille Nord Europe / CNRS CRIStAL -Lille, France.

b The University of Sydney -Sydney, Australia. c Exascience Life Lab, IMEC Research Center -Leuven, Belgium.

On the path to exascale which is fast approaching, multi-and many-core (co)processors are increasingly becoming key-building blocks of many computing platforms, as evidenced in the Top500 international ranking. On the other hand, plenty of hard problems in a wide range of areas (biology, energy, etc.) are often modeled and tackled using optimization approaches. These approaches include mainly greedy algorithms, exact methods (dynamic programming, Branch-and-X, constraint programming, A * , etc.) and meta-heuristics (evolutionary algorithms, particle swarm, ant or bee colonies, simulated annealing, Tabu search, etc.). In many research works, on the one hand, optimization techniques are used to face high performance computing (HPC) challenges including hardware design, compiling, scheduling, auto-tuning, etc. On the other hand, optimization problems become increasingly large and complex, forcing the use of parallel computing for their efficient and effective resolution. The design and implementation of parallel optimization methods raise several issues such as load balancing, data locality and placement, fault tolerance, scalability, thread divergence, etc.

This special issue invited high-quality contributions on the joint use of advanced (discrete, continuous or mixed, single or multi/many-objective, static or dynamic, deterministic or stochastic, hybrid) optimization methods and distributed and/or parallel multi/many-core computing, and any related issues. While the submissions were on the topics of the POMCO'2016 workshop organized in conjunction with HPC&S'2016, the call for papers was open to all contributors. We received 17 submissions from several different countries. The submissions were rigorously reviewed by at least three expert external reviewers, and further carefully evaluated by the guest editors. After the review process, only 5 were finally accepted for publication. We thank all the reviewers for their critical and expert help. Below, we provide an overview of the papers appearing in this volume.

In [START_REF] Aparicio | Parallel algorithms for computing the smallest binary tree size in unit simplex refinement[END_REF], the authors address the irregular combinatorial problem of finding the smallest binary tree size generated by the iterative partitioning of the regular simplex by a longest edge. Solving this problem requires an increasing computational burden in size (dimension and stopping criterion). Several multi-core algorithms are proposed and compared using TBB and PThread for threaded implementation. The TBB-based approach uses a static number of threads and features task stealing while the Pthread-based one uses a dynamic number of threads and the creation-destruction policy implicitly implementing load balancing. According to the reported experimental results several interesting observations related mainly to memory management and dynamic load balancing can be highlighted. First, the TBB version outperforms its Pthreaded counterpart for the small instances of the problem, but fails to solve the largest one. The Pthread version could solve more instances than TBB more efficiently, especially for the hardest instances. This can be explained by the way in which threads get work units from the search tree. Indeed, Breadth-First search is activated at deeper levels of the tree which results in less memory requirement. In addition, the LLAlloc memory allocator of LockLess Inc has been demonstrated to be the best among multiple experimentally evaluated allocators.

Another irregular combinatorial problem is tackled using parallel exact optimization in [START_REF] Sundfeld | PA-Star: a Disk-Assisted Parallel A-Star Strategy with Locality-Sensitive Hash for Multiple Sequence Alignment[END_REF]. The problem, called Multiple Sequence Assignement (MSA), is a basic operation in bioinformatics highly useful to find out the similarities among a set of sequences. MSA is modeled as a search for the path with minimum cost in a graph and solved using the state-of-the-art A * algorithm revisited in the paper using in single-core as well as multi-threaded computing for multi-core platforms. The design of a multi-threaded A * (PA-Star) for solving MSA raises two major challenges: irregular data access patterns and substantial memory requirements. The first challenge is raised using a one-level hash function with locality preserving characteristics allowing to reduce the overhead induced by the communications between threads, thereby increasing the scalability of the algorithm. The memory issue is addressed using out-of-memory computing. Indeed, PA-Star uses two primary lists of nodes, called OpenList and ClosedList, to explore the search with a best-first strategy. OpenList (resp. ClosedList) is used to store the nodes to be (resp. already) analyzed. The disk-assisted strategy consists in saving the nodes of ClosedList to the disk when the RAM is being saturated. The PA-Star algorithm has been experimented on multi-core systems using the BALiBASE reference set. The reported results demonstrate the efficiency of the proposed mechanisms and its superiority over the P F A * -DDD state-of-the-art algorithm.

The paper [START_REF] Menouer | Solving Combinatorial Problems using a Parallel Framework[END_REF] is focused on parallel software framework for solving large problems on multi-core systems. The author proposes an extension, called IBobpp, of the Bobpp framework.

Bobpp supports the node-based parallel model which is widely used for parallel tree search. In this model, a global shared pool is used to store the nodes generated by the parallel exploration process.

The tree nodes are distributed to threads according to the work stealing paradigm. The shared pool becomes a bottleneck when it comes to solve large problem instances limiting thereby the scalability of the framework. To deal with this issue the author proposes in IBobpp another approach using multiple pools. In addition, he proposes an implementation of the treebased coarse-grained parallel model. In this model, each thread explores locally and serially a sub-tree limiting the amount of memory accesses to the shared pool. Extensive experiments reported in the paper show that the combination of the node-based and tree-based parallel models together with the use of multiple pools allows to improve the performance and scalability of the framework.

The work presented in [START_REF] Atashpendar | A Scalable Parallel Cooperative Coevolutionary PSO Algorithm for Multi-objective Optimization[END_REF] deals with the parallelization of metaheuristics, which are near-optimal algorithms. More exactly, the focus is put on the Particle Swarm Optimization (PSO) metaheuristic for solving multi-objective optimization problems on multi-core systems. The authors propose a parallel variant of the SMPSO, which stands for Speed-constrained Multi-objective PSO, algorithm. The new algorithm called CCSMPSO, which stands for Cooperative Co-evolutionary SMPSO, is based on the cooperative co-evolution paradigm. This later consists in breaking down the population of the algorithm into multiple subpopulations that evolve specific parts of the global solution by cooperating during the evolution. Extensive experiments have shown that CCSMPSO outperforms the CC versions of state-of-the-art multi-objective metaheuristics, namely NSGA-II, SPEA2 and MOCell as well as their serial counterparts. In addition, the reported results highlight that CCSMPSO scales better when increasing the size of the problem (number of decision variables) as well as the number of processing cores. In particular, the algorithm improves the quality of obtained Pareto fronts (solutions) while the others decrease it.

In the four previous papers, fast execution is targeted but not energy efficiency. In [START_REF] Oxley | Rate-based Thermal, Power, and Co-location Aware Resource Management for Heterogeneous Data Centers[END_REF], the focus is put on this later criterion within the context of heterogeneous data centers. More exactly, the authors propose a new detailed model of a data center considering the power consumption, the performance of the compute nodes and cooling system, thermal constraints, DVFS and co-location interference. This later is crucial in modern data centers including massively multi-core compute nodes. Co-located tasks run on the processing cores within the same processor leading to co-location interference effect because of the use of shared memory and last-level cache. In addition, the authors propose, experiment and rigorously analyze three resource management techniques: a greedy heuristic, a genetic algorithm combined with a local search and a revisited non-linear programming approach. Co-location is considered together with power and thermal constraints.

Finally, we warmly thank all the authors for their high quality contributions. We hope that this volume will be of timely value to the readers.

Preprint submitted to ElsevierSeptember 13, 2017

GUEST EDITORS

• Nouredine Melab

Université Lille