Randa Mallat

Vincent Bonnet
email: vincent.bonnet@u-pec.fr

Raphaël Dumas

Mohamed Adjel

Gentiane Venture

Mohamad Khalil

Samer Mohammed

Raphael Dumas

R Dumas

Sparse Visual-Inertial Measurement Units Placement for Gait Kinematics Assessment

Keywords: Extended Kalman Filter, Gait rehabilitation, Inertial Measurement Unit, Augmented Reality Markers. I

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

are often used. When accurately estimated, those parameters can be combined to build assessment indicators to assist the therapists in their diagnosis. For example, the Gait Deviation Index (GDI) introduced by Schwartz et al. is increasingly used by the clinical community to estimate the level of gait pathology among disabled persons [START_REF] Schwartz | The Gait Deviation Index: a new comprehensive index of gait pathology[END_REF]. The GDI compares gait kinematics data of a given patient to a healthy group, by estimating the distance between healthy strides and patient's strides. This distance (with respect to the healthy group), which reflects the pathology level, requires an estimation of the kinematics parameters to be assessed properly. Stereophotogrammetric Systems (SS) are considered as the reference tools in kinematics motion analysis. However, their use require a large financial investment and a complex experimental protocol. An alternative solution proposed by Cappozzo was called the minimum measured input model approach, that is aiming at maximizing the functional information extracted from simplified and affordable experimental protocols [START_REF] Cappozzo | Minimum measured-input models for assessment of motor ability[END_REF].

A. Related works

In this context, mainly two technologies were proposed: Inertial Measurement Unit (IMU) and RGB/RGB-Depth camerabased systems. Small, low-cost IMUs measure three accelerations and three angular velocities. Thus, in principle, it is possible to estimate their position and orientation through double and single integration, respectively. Unfortunately, drift jeopardizes the time integration of raw signals [START_REF] Ferrari | A Mobile Kalman-Filter Based Solution for the Real-Time Estimation of Spatio-Temporal Gait Parameters[END_REF]. To overcome this problem, Kalman adaptive filters are often proposed in the literature to ensure real-time estimates of the IMU orientation by fusing accelerometer and gyroscope data. Since only two orientation angles can be obtained with this approach [START_REF] Sabatini | Quaternion-based extended Kalman filter for determining orientation by inertial and magnetic sensing[END_REF], several studies have proposed to use at least one IMU for each segment of interest together with a multi-body model representing the human skeleton as a mean to reduce driftinfluence. For example, El Gohary et al. proposed to use a zero velocity update approach with an Unscented Kalman Filter relying on a kinematics model for arm motions tracking [START_REF] El-Gohary | Human joint angle estimation with inertial sensors and validation with a robot arm[END_REF]. Remarkable results were obtained when operating with a rigid robotic arm so that human soft-tissue-artifact and sensor-tosegment calibration errors were not considered in this study. To reduce the sensors count, Sy et al. proposed recently a method to estimate lower limbs 3D joint angles during gait [START_REF] Sy | Estimating Lower Limb Kinematics using a Reduced Wearable Sensor Count[END_REF]. Despite a very interesting approach based on the use of an Extended Kalman Filter (EKF) including assumptions and constraints such as rigid body constraints, zero velocity update at the foot contact, assuming flat floor or making use of the total pelvis height, the presented results displayed a large Root Mean Square Difference (RMSD) superior to 10deg when compared to the same estimates using SS. These relatively poor results might be due to the fact that a proper sensor-to-segment calibration was missing and that a constant acceleration model was used. Indeed, IMUs require a calibration phase to align each sensor with respect to its segment. Static and functional calibration techniques based on a set of pre-defined postures and/or motions are commonly used in literature. However, these methods are time-consuming, might be limited by the subject's ability to perform pre-defined postures or motions, and thus can be hardly used with some population or in clinical settings for example [START_REF] Seel | IMU-Based Joint Angle Measurement for Gait Analysis[END_REF]. Large errors might also be due to sub-optimal weights tuning leading to an emphasized influence of IMU drift on joint angle estimates. Joukov et al. showed the influence of parameter's tuning when using an adaptive filter such as the EKF [10]. Moreover, during gait, they showed the superiority of using Fourier series over using standard constant acceleration model of joint temporal evolution. The assumption that the joint trajectories can be modelled by using low-order Fourier series has already been made in the literature [10]- [START_REF] Ziegler | Simultaneous identification of human body model parameters and gait trajectory from 3D motion capture data[END_REF]. Fourier representation was also largely used for spatio-temporal parameters assessment using IMUs [START_REF] Sabatini | Fourier-based integration of quasi-periodic gait accelerations for drift-free displacement estimation using inertial sensors[END_REF]. Fourier series coefficients were identified using a feedback adaptive frequency phase oscillator. As Joukov et al. discussed, the learning rate and parameters tuning of the oscillator was decisive in the performance of their method [10]. Nevertheless, by performing a grid search, it was possible to fine-tune the filter and obtain a low RMSD of 2.4deg for the hip and knee joint angles. However, it should be noted that the proposed method could not estimate ankle angles, required a cumbersome sensor-to-segment calibration process and was based on one IMU per investigated segment. To cope with the IMU drawbacks, some authors have proposed to use camera based systems. Recently, markerless visual motion capture based on a RGB camera and machine learning algorithms were developed for estimating human motion in real-life scenarios [START_REF] Chakraborty | Accuracy of image data stream of a markerless motion capture system in determining the local dynamic stability and joint kinematics of human gait[END_REF], [START_REF] Zongmian | Estimating 3D Motion and Forces of Person-Object Interactions From Monocular Video[END_REF]. However, such approaches are not yet real-time or accurate enough to reliably estimate dependent subject's joint angles and thus are not ready to be used yet in clinical settings for example. Moreover such methods need often to be trained on large scale datasets requiring a large computational power. Devices, based on RGB-Depth sensors such as the Kinect camera (Microsoft), with embedded skeleton tracking algorithms are available. However, it was shown that the Kinect camera is not accurate enough for rehabilitation assessment mainly due to segment length variations [START_REF] Colombel | Physically Consistent Whole-Body Kinematics Assessment Based on an RGB-D Sensor. Application to Simple Rehabilitation Exercises[END_REF]. For gait analysis, the Kinect sensor can only estimate timing and spatial characteristics [START_REF] Bonnechère | 3D Analysis of Upper Limbs Motion during Rehabilitation Exercises Using the KinectTM Sensor: Development, Laboratory Validation and Clinical Application[END_REF] as joint angles are not reliable enough for a consistent analysis. Numerous methods have been developed in the literature to improve kinematic estimates of RGB-D like sensors. Most of them add the use of a multi-body model to obtain a more robust joint angle estimate [START_REF] Zongmian | Estimating 3D Motion and Forces of Person-Object Interactions From Monocular Video[END_REF], [START_REF] Colombel | Physically Consistent Whole-Body Kinematics Assessment Based on an RGB-D Sensor. Application to Simple Rehabilitation Exercises[END_REF]. However, acceptable differences can only be obtained for motions of very large amplitude [START_REF] Colombel | Physically Consistent Whole-Body Kinematics Assessment Based on an RGB-D Sensor. Application to Simple Rehabilitation Exercises[END_REF].

Another camera based approach is homography, which is used in Augmented Reality (AR) to accurately estimate 3D rotations and translations of a given AR marker after a proper, yet simple, camera calibration process.

Nagymate et al. proposed the use of AR markers for gait analysis [START_REF] Nagymate | Affordable gait analysis using augmented reality markers[END_REF]. They have used a commercial Multi-body Kinematics Optimization (MKO) approach that was fed with segments 3D pose obtained from very large AR markers that were located at each investigated segments. They displayed a RMSD of 2.3cm for the step length, and a joint angle estimate varying from 2.5deg to 6.7deg. The authors state that they have chosen to use such large AR markers to avoid occlusion and undetected markers. In a recent survey [START_REF] Chen | A survey of depth and inertial sensor fusion for human action recognition[END_REF], it was suggested that merging RGB-Depth and IMU data could prevent occlusions and improve joint angle estimate accuracy [START_REF] Feng | Fusing Kinect sensor and inertial sensors with multi-rate Kalman filter[END_REF]. Ahmed et al. for example were able to estimate head and feet trajectories during overground walking, with a RMSD that was less than 2% of the walked distance, based on the fusion of only three sensors: a head mounted IMU-camera and two feet mounted IMU [START_REF] Ahmed | A Visual-Inertial Approach to Human Gait Estimation[END_REF]. However, the joint kinematics were not assessed in this study. Moreover, it is worthy noting that MKO and EKF were alternatively used with IMU [10] and AR markers [START_REF] Nagymate | Affordable gait analysis using augmented reality markers[END_REF], [START_REF] Mallat | Upper Limbs Kinematics Estimation Using Affordable Visual-Inertial Sensors[END_REF] (as well as with SS markers [START_REF] De Groote | Kalman smoothing improves the estimation of joint kinematics and kinetics in marker-based human gait analysis[END_REF]). Both methods seem to provide similar results however the MKO can hardly be implemented in real-time [START_REF] Begon | Multibody kinematics optimization for the estimation of upper and lower limb human joint kinematics: a systematized methodological review[END_REF].

In this context, it seems that there is a consensus in the community on the fact that a multi-body model based approach should be used to perform inverse kinematics with a sparse placement of sensors. Fusing camera based information, obtained from AR markers, and IMU data could lead to reliable and robust 3D sensor pose estimate. Proper sensor-to-segment calibration and EKF parameters tuning are also essential to provide accurate joint kinematics estimation. Our group has already proposed such an approach in a recent study [START_REF] Mallat | Upper Limbs Kinematics Estimation Using Affordable Visual-Inertial Sensors[END_REF] but it was dealing with simple arm tasks. Moreover, it required one sensor module per investigated segment. Nevertheless it allowed to show that it is possible to reduce IMU drift influence and to cope with visual occlusion problems.

B. Contribution

The objective of the present study is to demonstrate the validity of using sparse calibrated newly designed Visual-Inertial Measurement Units (VIMU) and adaptive filtering, during treadmill gait, to estimate relevant kinematics quantities that may be of great benefit to the clinical community. Moreover, besides new constraints implemented in the EKF, two formulations of the proposed state-vector are studied and compared: a classical constant acceleration model associated to an optimal tuning process of the EKF's parameters, and a model based on Fourier series that allows a simple tuning of the EKF's parameters based on a priori knowledge. Thus one of the sub-objective of this study is to assess if by having a more realistic state variable evolution model the tuning of the EKF parameters will be simplified.

II. METHOD

The overall principle of the proposed approach is depicted in Fig. tracked with two affordable RGB cameras. The raw data collected from both IMUs and AR markers require a prior and unique calibration procedure that is based on a process detailed in our previous study [START_REF] Mallat | Upper Limbs Kinematics Estimation Using Affordable Visual-Inertial Sensors[END_REF] and briefly described in section II.B. A set of three AR markers is mounted on the top of each IMU sensor to form a so-called VIMU (see Fig. 3). Then, a practical subject-specific sensor-to-segment static calibration (section II.B.2) is done to define each VIMU pose in the segment coordinate systems (subject's segment axes and joint centres are therefore defined at this step). Finally, an EKF based on an eighteen Degrees-of-Freedom (DoF) multi-body kinematics model that is described in section II.A, fuses all measured quantities to estimate lower-limbs joint angles. The tuning of the EKF's parameters along with its state-vector expressed in two different versions with respect to the evolution of the joint trajectories are also compared and discussed (section II.C and section II.D).

A. Lower-limbs Multi-Body Model

A multi-body model of the lower-limbs is devised to relate the joint kinematics with the three VIMU measurements. It consists of = 7 segments articulated with = 18 DoF. As represented in Fig. 2.a, it is composed of two spherical joints (hip), two hinge joints (knee) and two universal joints (ankle). Moreover, the model floating base is located with respect to a global coordinate system ( 0 ) using three prismatic and three rotational virtual joints. Three VIMUs are assumed to be rigidly attached to the sacrum, and to the left and right heels, respectively (Fig 2 .b). A VIMU is composed of three AR markers and a single IMU (see Fig. 3 and Fig. 4). The pose (position and orientation) of each VIMU marker is directly estimated by the off-the-shelf ArUco library [START_REF] Romero-Ramirez | Speeded Up Detection of Squared Fiducial Markers[END_REF]. The orientation provided by ArUco is expressed using the quaternion formalism as presented by Diebel [START_REF] Diebel | Representing attitude: Euler angles, unit quaternions, and rotation vectors[END_REF]. Thus, each VIMU allows measuring the absolute 3D position p 0 and orientation (quaternion) q 0 with respect to 0 of each AR marker as well as the 3D linear acceleration a and angular velocity with respect to the IMU coordinate system. Three coordinate systems are therefore used to represent VIMU different elements: a coordinate system corresponding to both IMU and central AR marker and two additional coordinate systems for the left and right AR markers, respectively. The forward kinematics model (FKM) and its first and second derivative (J and J) relating joint kinematics to the VIMU measurements are given as follows:

p 0 = FKM( ) q 0 = f (R 0 v ) = (R 0 ) J ( ) a = (R 0 ) (J ( ) + J ( ) ) + b (1)
where , , are the ( × 1) joints positions, velocities and accelerations vectors, respectively. R 0 (3 × 3) is the rotation matrix expressing the VIMU's central AR marker orientation in 0 . It is converted into a (4×1) quaternion vector q 0 through a custom mapping function f (). Details on the implementation of f () can be found in section 6.5 of the review paper proposed by Diebel [START_REF] Diebel | Representing attitude: Euler angles, unit quaternions, and rotation vectors[END_REF]. The Jacobian matrix

(6 × 3) is composed of J = [J J ],
and b is a (3 × 1) acceleration bias vector.

B. Calibration

Two easy-to-use calibration procedures are required to use the raw measurements effectively:

• a calibration for each VIMU, independent of the subject,

• a subject-specific sensor-to-segment calibration. On the one hand, the VIMU calibration is done without participation of the subject by positioning and moving the unit in the field of view of the RGB camera. On the other hand, thanks to the use of AR markers, an anatomical calibration of the VIMU can simply be done [START_REF] Cappozzo | Position and orientation in space of bones during movement: anatomical frame definition and determination[END_REF].

1) VIMU calibration:

The gyroscope and accelerometer were calibrated, without the use of any external device, as proposed by Tedaldi et al. to take into account inaccuracies in scaling factors, misalignment and non-orthogonality of axes, and non-zero biases [START_REF] Tedaldi | A robust and easy to implement method for IMU calibration without external equipments[END_REF]. Moreover, as shown in our previous study [START_REF] Mallat | Upper Limbs Kinematics Estimation Using Affordable Visual-Inertial Sensors[END_REF], the coordinate systems of the IMU and of the central AR markers were aligned thanks to a dynamic calibration step making use of centripetal acceleration, to calculate the rotation matrix R (3 × 3). Two rigid transformations also exist between the side coordinate systems of the AR markers and of the central one:

R = (R 0 ) R 0 q = (R ) r = (R 0 ) (r 0 -r 0 ) (2)
where R (3 × 3), q (4 × 1), r (3 × 1) represent the rotation matrix, its corresponding quaternion formulation, and position vector of = (left) or (right) AR marker expressed in the central one and obtained after calibration, respectively. These relative transformations were easily determined using 50 different static postures of the VIMU and solving an overdetermined system [START_REF] Soderkvist | Determining the movements of the skeleton using well-configured markers[END_REF]. These rigid transformations between the coordinate systems of the AR markers will be considered as rigid body constraints in the EKF (section II.C). Finally, a VIMU is fully calibrated given that all its measurements are expressed with respect to the coordinate system attached to the central AR marker frame . 2) Sensor-to-Segment Calibration: As exemplify in Fig. 4 the information from an additional AR marker mounted on a large calibration wand, making its detection accurate and robust, can be used to pin-point anatomical landmarks similarly to what is done when using a SS [START_REF] Cappozzo | Position and orientation in space of bones during movement: anatomical frame definition and determination[END_REF]. The calibration wand is used to pin-point eleven lower-limbs anatomical landmarks in 0 . These are the midpoint between the postero-superior illiac spines, the right and left antero-superior illiac spines, the lateral and medial femoral epicondyles, and the lateral and medial malleolus. This operation can be done in less than a minute. Joint center positions are then calculated using classical linear regression methods [START_REF] Bell | Prediction of hip joint centre location from external landmarks[END_REF]. The norm between consecutive joint center positions is used to determine the segments lengths in the multi-body model (section II.A, Fig. 2). Note that this sensor-to-segment calibration is valid as long as the VIMUs location on their corresponding segments are not modified. Otherwise, it has to be repeated.

C. Extended Kalman Filter

An EKF is proposed to estimate lower-limbs joint kinematics, gathered in the state-vector x, that minimizes at each time step the least-square difference between the VIMU measurements, gathered in y, and their estimate from the measurement model h, gathered in h, including the estimate of the VIMUs 3D positions p 0 , orientations q 0 , linear accelerations a and angular velocities (section II.A).

h = [p 0 q 0 a h 1 h 2 ] h 1 = [( - ) 2 ] with = 7, ..., h 2 = [(r -r ) 2 (q -q ) 2 (r -r ) 2 (q -q ) 2 ]
(3) The vector of estimated VIMUs measurements h is based on the forward kinematics of the multi-body model but also incorporates two type of constraints terms (c1 and c2). The elements of h 1 ( × 1) are used to solve kinematics indetermination. Since the number of VIMU is reduced, i.e. a single VIMU per kinematics chain, there exist multiple kinematics solutions leading to the same VIMU measurements in the Cartesian space. However, there is not an infinite number of solutions leading to the same VIMU pose as the proposed leg model has 6DoF similarly to the number of estimated quantities from a VIMU. This kinematics indetermination is well known as the elbow-up/elbow-down problem in robotics. It has been solved in different ways. In this study, we propose to push the joint solution toward the average joint angle value, , reported for walking in the literature [START_REF] Rose | Human Walking[END_REF], [START_REF] Ramakrishnan | On the estimation of joint kinematics during gait[END_REF]: 13deg, 8deg, 11deg, 16deg, 1deg, 6deg for 7 , 8 , 9 , 10 , 11 , 12 angles of the left leg and, symmetrically for the right leg, respectively. In the context of the EKF, this is done by adding -6 virtual measurements aiming to minimise the square difference with at each sample of time by setting the corresponding measurements in y to zero. This is equivalent to minimizing a cost function in an optimization process [START_REF] Aravkin | Optimization viewpoint on Kalman smoothing with applications to robust and sparse estimation[END_REF].

The elements in h 2 (42 × 1) additionally ensure rigid body constraints between the three AR markers in the VIMU. The position and orientation of the right and left AR markers with respect to the central AR marker should match their values obtained during the VIMU calibration (section II.B).

the measurement model was established, the joint kinematics estimation was conducted, first, through the prediction of the a priori state-vector xthanks to a state evolution model , that is approximated by a linear form F, and a process covariance matrix Q:

x-= Fx -1 P -= FP -1 F + Q (4)
where P -is the a priori estimation of the error covariance matrix. Two cases were investigated for the state-vector formulation:

• case 1 EKF based on a Constant Acceleration (EKF-CA) model: A first classical state-vector formulation proposes to estimate the joint angles, velocities and accelerations and assumed a constant acceleration model [START_REF] De Groote | Kalman smoothing improves the estimation of joint kinematics and kinetics in marker-based human gait analysis[END_REF]:

x = [ b ] with = 1, ..., =               I × Δ I × T 2 s 2 I × 0 × I × Δ I × 0 × 0 × I ×       0 ×9 0 9 3 I 9 9         (5)
where , 0 and I are the sampling time, and the null and identity matrices, respectively. To reduce the drift effect, a random acceleration bias b (9 × 1) was added to better track the 3D measured accelerations of each VIMU (see Eq. 1) [START_REF] El-Gohary | Human joint angle estimation with inertial sensors and validation with a robot arm[END_REF].

• case 2 EKF based on a Fourier Series (EKF-FS) model: Since gait can be assumed to be a quasi-periodic task, the angle of the ℎ joint can be represented using a Fourier series:

= a 0 + =1 (a ( ) + b ( )) (6)
where is the number of harmonics, and is the main motion pulsation.

As described in the experimental validation (see section III), the gait activity takes place on a treadmill at constant velocity. Consequently, Fourier series coefficients are expected to converge towards near constant values and identity matrix can be used as state-vector evolution model. The state-vector and its evolution model were defined as follows:

x = [a 0 a b b ]
with = 1, ..., and = 1, ...,

x +1 = I × x (7) with = (1 + 2 ) + 1 + 9.
Secondly, depending on the measurement error, v, and the measurement covariance matrix, R, the Kalman gain K was calculated to update the state-vector:

v = y -ℎ( x-) S = H P -H + R K = P -H S -1 x = x-+ K v P = K R K + (I -K H )P -(I -K H ) (8)
where H is the Jacobian matrix ℎ

x and S is the measurement covariance matrix.

D. EKF Parameters Tuning

The proposed EKF requires the tuning of the following elements: the initial state-vector value x0 , the initial estimation of the error covariance matrix P 0 , the process noise covariance matrix Q and the measurement noise covariance matrix R.

For the initial state-vector and error covariance matrices, it has been shown that their values mainly affect the initial part of the estimation [START_REF] Zarchan | Fundamentals of Kalman Filtering: A Practical Approach[END_REF]. Thus, x0 was set with the result of the inverse geometric model calculated from the AR markers data collected at the first sample of time. P 0 was set equal to the identity matrix giving the same influence to all joints.

The tuning of the elements of Q and R is more sensitive as they impact the filter convergence rate and stability. The elements of R corresponding to the IMU data can be experimentally determined. However, the accuracy of the visual measurements is variable as a function of different aspects such as marker occlusion, motion blur, or distance and orientation with respect to the camera and thus is difficult to be assessed. The elements of R corresponding to the constraints c1 and c2 (section II.C) must also be tuned leading to a large number of parameters. Moreover, the ratio between the tuning of the elements of R and Q, that is more sensitive and task and joint dependent [START_REF] De Groote | Kalman smoothing improves the estimation of joint kinematics and kinetics in marker-based human gait analysis[END_REF], affects the filter performance [START_REF] Zarchan | Fundamentals of Kalman Filtering: A Practical Approach[END_REF]. The tuning of Q reflects the trust level that could be given to the state-vector evolution model . For example, the constant acceleration model is obviously incorrect for human motion leading to inconsistencies in the state-vector [START_REF] Bonnet | A constrained extended Kalman filter for the optimal estimate of kinematics and kinetics of a sagittal symmetric exercise[END_REF] and thus this should be reflected in Q. Similarly, despite the bracelet form, AR markers have shown some sensitivity to occlusions. Fortunately, these are easily detectable. Thus, if the loss of all three markers data was detected, the corresponding elements of the EKF noise covariance matrix R will be automatically updated to reflect the new state of input data. Consequently, two R matrices were considered: the one corresponding to the whole markers trajectories when no occlusion of at least one marker has occurred and another one when all AR markers were occluded. In the latter, the missing markers were replaced by their previous values and their associated R elements were given a large value (1 3 ) indicating a lack of trust in these measurements. Considering the assumption that the parameters tuning is task-dependent [START_REF] De Groote | Kalman smoothing improves the estimation of joint kinematics and kinetics in marker-based human gait analysis[END_REF], which is commonly used in EKF for motion analysis, the data of one randomly chosen subject was used to identify the elements of Q. The same identified values were then used for the remaining subjects. Then, depending on the state-vector formulation, two methods to tune the parameters of R and Q were proposed in this study.

In the case 1, corresponding to the EKF-CA, an optimisation process aiming at tuning the parameters of the EKF-CA was used. The optimization process consists of finding V = [

Ω R ] (13 × 1) that minimizes the square difference between the estimated joint angles and the reference ones, obtained from the SS and MKO, respectively (see section IV). The optimization process is computationally costly. Thus, it is proposed to group some parameters and to set the elements corresponding to the constraints c2 to 1 -3 . The choice of this value reflects the fact that the reference values are almost exact. The thirteen parameters that need to be identified are: three parameters , , and corresponding to the joint angles, velocities and accelerations elements of the process covariance matrix Q, and ten parameters corresponding to the measured positions , quaternions , linear accelerations , angular velocities Ω , and to the six joints R = in constraint c1 with = 1...6 (assuming left and right symmetry). The problem of the optimal tuning of V then boils down to:

Find V * ∈ min V∈R 13 =1 =1 ( - ) 2 s.t. V -≤ V ≤ V + (9)
where is the number of samples of the considered trial. V -and V + represent the lower and upper boundaries, respectively.

(18 × 1) is the vector of the reference joint angles calculated using the SS (see Section IV). It should be noted that the gradient of this formulation is difficult to calculate since the equation in the EKF needs to be integrated. Thus, Eq. 9 was solved using a local derivative-free solver proposed in the popular IPOPT library [START_REF] Powell | Direct search algorithms for optimization calculations[END_REF], leading to V = [165.23 3.18 0.17 0.13 0.17 1.09 0.46 [0.09 0.04 0.03 2.32 0. 1 1]] .

In the case 2 corresponding to EKF-FS, it is crucial to determine firstly the number of harmonics required to accurately estimate the joint angles, velocities and accelerations. To do so, it is proposed to a priori solve the following fitting problem for different values of :

Find a * , b * , * ∈ min a * ,b * , * ∈R 2 + +1 =1 =1 ( - ) 2
s.t. 0 ≤ (10) The Fourier representation is also used to obtain an analytical estimate of the joint velocities and accelerations. Thereafter, the results of this fitting process were used to determine a good trade-off between the number of parameters and the accuracy. Then, the parameter tuning is only based on a priori knowledge of the biomechanics of the task. The diagonal elements in Q were set to a small value (0.01) to enforce their convergence to pseudo-constant values [START_REF] Bonnet | A constrained extended Kalman filter for the optimal estimate of kinematics and kinetics of a sagittal symmetric exercise[END_REF]. This guarantees the convergence of the corresponding error covariance P to a value close to zero. The elements of R related to the measurements were all empirically set to 0.1. As it regards the constraints (c1) related to the difference between estimated and average joint angles, the elements were set accordingly to the joint angle amplitude such that R = [1 1 1 10 1 10] while the corresponding elements of the rigid body constraints (c2) were set to 1 -3 .

III. EXPERIMENTAL VALIDATION

A. VIMU Prototype

The proposed motion capture system was developed to be affordable, user-friendly, compact and light-weight. VIMU prototypes are shown in Fig. 3.a. Each VIMU consists of one IMU (MPU-6050, Invensense) embedding a 3D gyroscope and a 3D accelerometer, and communicating at 100Hz using a Bluetooth wireless module. On the top of the IMU enclosure, measuring 4x4x2cm, three AR markers of 3.6x3.6cm are mounted in a 3D printed bracelet form. The bracelet weights 40g and has different forms depending on the targeted segment. Its form is convenient to improve the AR markers detection rate during motions involving complex rotations. VIMUs were attached to human body segments using a velcro strap. Two standard 60Hz high-definition RGB USB cameras (USBFHD01, ELP) were used to track the AR markers poses. A custom program is developed in C++ to read all VIMUs data synchronously at a frequency rate of 60 Hz. The total cost of the proposed system was below 100 e.

B. Experimental protocol

Ten healthy male volunteers (age 25.4 ± 4.1 years old, weight 71.2 ± 16.0 Kg, height 1.7 ± 0.08m) participated in the experiments. They gave their informed consent following ethical regulation for non-invasive experiments of the University of Tokyo Agriculture and Technology (Tokyo, Japan) and of the University of Paris-Est Créteil (Créteil, France) where experiments took place. Once the VIMUs and retroreflective markers were worn by the subject, at the different locations shown in Fig2.b, every subject was told to stand statically, for less than one minute, in any comfortable posture that allows detection of VIMUs AR markers using the RGB cameras. Two RGB cameras were setup on tripods to be able to visualize the whole scene and detect all AR markers. The two RGB cameras parameters were obtained thanks to a standard calibration process based on a chArUcoBoard [START_REF] Bradski | The OpenCV Library[END_REF]. The data collected by the two cameras were thereafter expressed in the same coordinate system using a large AR marker located on the floor in the field of view of the two cameras (Fig. 4). Meanwhile, anatomical landmarks of interest were pin-pointed using the affordable calibration wand and also captured using the SS (see Fig. 2.b). As a result, joint centres' positions, segments' lengths, as well as the local pose of each VIMU sensor relative to its corresponding segment anatomical frame, could be accurately calculated (Section. II.B). Hereafter, all subjects performed normal gait on a treadmill for one minute while wearing their running shoes. Each participant was free to choose his natural and comfortable gait speed.

Accuracy Assessment

The ability of the two proposed EKF state-vector formulations to accurately estimate the lower-limbs joint kinematics was assessed by calculating the RMSD and the Pearson Correlation Coefficient (CC) between the estimated and reference joint angles obtained by means of the SS. This choice is due to the fact that we were interested in estimating the temporal joints evolution, independently of the alignment and calibration offsets between the global reference frames of the different motion capture systems. Note that the sensor data are presented in international system of units while the outputs of the proposed system are given in degrees and centimeters since these units are more easily interpreted in a clinical context. The reference joint angles were obtained using a standard reference SS (6 Prime cameras, Optitrack) and a state of the art MKO tracking 22 reflective markers [START_REF] Begon | Multibody kinematics optimization for the estimation of upper and lower limb human joint kinematics: a systematized methodological review[END_REF] with the multi-body model shown in Fig. 2. The joint angles computed with the SS follow the Conventional Gait Model [START_REF] Leboeuf | The conventional gait model, an open-source implementation that reproduces the past but prepares for the future[END_REF]. Moreover, the knee flexion axis was defined orthogonal to the long axis of the femur and, at the hip, the joint coordinate system corresponded to 2-by-2 mutually orthogonal axes, similarly to the multibody model as shown in Fig. 2.a. The data of the three VIMUs, located at the pelvis and both heels of the subjects as shown in Fig. 2.b, were synchronously collected. The gait deviation index (GDI) [START_REF] Schwartz | The Gait Deviation Index: a new comprehensive index of gait pathology[END_REF] was calculated using the estimated joint angles obtained from EKF-CA and from EKF-FS. Both results were compared with the mean GDI obtained from reference SS data, with a leave-one-subject-out cross validation. Estimates of nine subjects were used to build the basis required to estimate GDI of the 10 ℎ subject. A total of approximately 260 healthy strides were used to build the GDI basis matrix on 9 subjects. The strides of the 10 ℎ subject were labelled as test strides, and used to calculate the GDI relatively to healthy strides. This operation was performed for each subject.

Moreover, based on the joint angles estimates and the forward kinematic model (see Eq. 1), the sacrum and heels positions were estimated and used to compute the gait stride lengths. The occurrences of heel strikes on the treadmill were first determined using the method developed in [START_REF] Zeni | Two simple methods for gait events during treadmill and overground walking using kinematic data[END_REF] and then used to estimate the corresponding stride lengths as in [START_REF] Supakkul | Using positional heel-marker data to more accurately calculate stride length for treadmill walking: a step length approach[END_REF].

Finally, we have studied the differences between EKF-CA and EKF-FS results by verifying that they were statistically nonsignificant. This was done using the corresponding RMSD and CC of all estimated gait kinematics parameters and of the GDI through a paired t-test with an alpha level set a priori to 0.05. The Pearson Chi-square goodness-of-fit test was conducted to check for normal distribution in the data. Fig. 5. Mean RMSD between the joint angles, velocities and accelerations estimated using Fourier series expansion and those obtained using the SS as a function of the number of harmonics .

IV. RESULTS

The mean RMSD of the joint angles, velocities and accelerations calculated using Fourier series representation and those estimated using the SS as a function of the number of harmonics for a randomly selected subject are shown in Fig. 5. To minimise the RMSD for the joint trajectories the number of harmonics was set to = 4 accordingly to Fig. 5. Fig. 6 shows the tracking performance using the proposed EKF-CA of the data collected by the VIMU sensor attached to the left heel while walking on the treadmill. Overall the input data were correctly tracked. The RMSD on the IMU data tracking were on average of 0.017 m.s -2 and of 0.029rad.s -1 for the linear accelerations and for the angular velocities, respectively. For the AR marker positions and quaternions, the RMSD were on average of 3.8cm and 2.2 -2 , respectively. These relatively large tracking errors were due to the fact that the AR markers were not detected just after toe-off as exemplify by the zoomed windows of Fig. 6. In fact at this moment the AR markers were positioned almost perpendicularly to the cameras. Despite these tracking errors, the EKF was able to ]

Quaternions Estimated Measured 𝑟 !" # 𝑟 !$ # 𝑎 %& % 𝑎 %' % 𝑎 %( % 𝜔 !" ! 𝜔 !) ! 𝜔 !$ ! 𝑞 !* # 𝑞 !" # 𝑞 !) # 𝑞 !$ # 𝑟 !) # 1 N 1 N 1 N 1 N 1 N 1 N 1 N 1 N 1 N N 1 N 1 N 1 N 1
Samples Fig. 6. Comparison between the input data of the VIMU attached to the left heel (red) and their tracking with the EKF-CA (dashed black). These are the 3D position and quaternion, 3D linear acceleration and 3D angular velocity. Note that the vertical grey lines are used to separate data that are presented for all N samples on the same subplot.

accurately estimate the eighteen joint angles of the lower-limbs model described in Fig. 2. This is valid for both the CA and FS state-vector formulations using the optimal parameters tuning. Fig. 7 exemplifies the accurate estimates of the 6 DoF that were used to calculate the pose of the pelvis with respect to 0 for both approaches. In this figure, the corresponding average RMSD of the 3D pelvis position was equal to 0.4±0.0cm (CC 0.86 ± 0.20) and 1.6 ± 0.9deg (CC 0.75 ± 0.16) in the case of EKF-CA and 0.4 ± 0.0cm (CC 0.80 ± 0.30) and 2.1 ± 0.7deg (CC 0.76 ± 0.21) in the case of EKF-FS. Fig. 8 shows a representative comparison of the left leg joint angles estimate obtained using both EKF-CA and EKF-FS approaches. The corresponding joint angles were satisfactorily estimated with an average RMSDs and CCs of 3.3 ± 1.6deg and 0.76 ± 0.24 using the EKF-CA and of 3.3 ± 1.2deg and 0.74 ± 0.20 using the EKF-FS, respectively. Table I summarizes the results of the comparison between the lower-limb joint trajectories obtained with the proposed system and the SS, for each considered measure, based on the data from all subjects. Note that the results for the leg joints (i.e., hip, knee and ankle) were obtained over the average of left and right leg values. In general, results were similar between both EKF-CA and EKF-FS approaches. For the pelvis the average RMSDs were of 0.6±0.1cm (CC 0.84±0.21) and 2.0±0.25deg (CC 0.60 ± 0.09) when using the EKF-CA and 0.7 ± 0.0cm (CC 0.82±0.22) and 2.2±0.4deg (CC 0.63±0.18) when using the EKF-FS. For the leg joint angles, i.e., the hip, knee, and Sterephotogrammetric system EKF-CA EKF-FS Fig. 7. Representative comparison of the pelvis joint trajectories obtained for a random selected subject while walking on a treadmill using the stereophotogrammetric system (black) and the proposed approach based either on a constant acceleration model (red) or on Fourier series expansion (dashed blue).

ankle joint angles, the average RMSDs for all subjects were equal to 3.9 ± 1.0deg (CC 0.74 ± 0.19) and 3.9 ± 0.8deg (CC 0.71 ± 0.19) when using EKF-CA and EKF-FS, respectively. Table . II summarizes the mean and standard deviation RMSDs for right and left stride lengths estimates. These were calculated first using the VIMUs raw data (without processing) obtained from the AR markers pose estimation, then using both EKF-CA and EKF-FS estimated data, and were compared with respect to the SS. In both cases, the use of the EKF reduced the RMSD in stride length estimates. Specifically, the EKF-FS showed an improvement by 50% approximately passing from 2.6cm to 1.3cm. Moreover, using the joint angles estimates for the detected strides, the GDI was also calculated. Results presented in Table II show that all the mean of GDIs (the reference and the estimated ones), are very close to 100 indicating a normal walking pattern.

The paired t-test results for RMSD and CC of each joint trajectory estimated using both the EKF-CA and EKF-FS implementations are reported in Table I. For all joints, results showed nonsignificant statistical differences (p > 0.05) between EKF-CA and EKF-FS with the exception of the pelvic tilt angle (p = 0.001, p = 0.039). Nevertheless, considering the above results, the gait kinematics estimates do not seem to be deteriorated by the poor estimate of the pelvic tilt angle. Errors in pelvis tilt angle also mechanically affects hip flexion-extension. However, while pelvis tilt is of Sterephotogrammetric system EKF-CA EKF-FS Fig. 8. Representative comparison of the left leg joint angles obtained for a random selected subject while walking on a treadmill using the stereophotogrammetric system (black) and the proposed approach based either on a constant acceleration model (red) or on Fourier series expansion (dashed blue) limited amplitude during gait, hip flexion-extension is of large amplitude. Therefore, the segment angle errors do not translate into significant joint angle errors between EKF-CA and EKF-FS (p = 0.200 for RMSD and p = 0.071 for CC). A similar statistical test was performed in Table III for the RMSD of the stride lengths and for the GDI calculation. Both stride lengths' RMSD and GDI showed nonsignificant differences between EKF-CA and EKF-FS. The difference in GDI values resulting from each EKF implementation and SS data was also nonsignificant. On the other hand, the raw stride length values compared to their estimates based on the two EKFs (see Table III) revealed, as expected, significant differences. This emphasizes the improvement provided by the EKFs in either implementation. Fig. 9 shows the evolution of the nine Fourier series coefficients corresponding to the ankle flexion angle estimated using EKF-FS. The coefficients converge toward pseudo-constant values showing, as expected, the gait periodicity.

V. DISCUSSION

The proposed system proposes to fuse visual and inertial data of three VIMUs attached to the pelvis, left and right heels into an EKF and a lower-limb kinematic model. This allowed handling the inaccuracies that result from using affordable sensors. Visual occlusion was greatly limited thanks to the bracelet form of the VIMUs. Rigid body constraints were integrated in the EKF formulation allowing to reconstruct marker loss based on its non-occluded partners within the same VIMU. This was possible as the rigid transformation between each two VIMU markers was priory calibrated (section II.B.1). Nevertheless, depending on the motion phases, the complete occlusion of the three AR markers could not be entirely avoided. Yet, the EKF was able to cope with such occlusion by relying more on IMU data thanks to the online update of the measurement covariance matrix R. This is shown in the zoomed capture of Fig. 6 where one can note that the EKF is not tracking inconsistent visual data. However, it is important to note that the occlusion was never longer than few seconds.

We have shown in a recent study that EKF approach was able to cope efficiently with such short occlusion [START_REF] Mallat | Upper Limbs Kinematics Estimation Using Affordable Visual-Inertial Sensors[END_REF]. Thus we can conclude that the proposed method is not very sensitive to occlusion with the current setup. Moreover, constraints having the effects of joint springs pushing the joint angles toward mean values were added to the EKF as soft constraints in the measurement model h. This was motivated by the fact that a sparse placement of the sensors was used. Thus, for the same VIMUs tracking, different joint configurations, possibly within the joint limits, are possible. Basically one could imagine to have springs of different stiffness attached to each joint. This concept is similar to the popular potential field approach [START_REF] Khatib | The potential field approach and operational space formulation in robot control[END_REF] commonly used within the robotics community.

The proposed approach has been validated with a direct comparison to the joint kinematics estimated using a reference SS. For consistency in the joint trajectories comparison, the SS data were also processed with MKO using the same multibody model. When tracking the same data, EKF and MKO are expected to provide very similar results [START_REF] Begon | Multibody kinematics optimization for the estimation of upper and lower limb human joint kinematics: a systematized methodological review[END_REF]. Over the fifteen lower-limbs joint angles the mean RMSD was of 3.5deg. These results show the ability of the proposed approach to handle the 3D IMU drift without magnetometers involvement while compensating a sparse sensors placement.

For both EKF implementations, the knee and hip flexionextension angles estimates presented a relatively high CC of 0.927 on average and a low RMSD ( 5deg). In contrast, the CCs of the ankle dorsi/plantar flexion (0.435) and the pelvic tilt angle (0.485) estimates are the poorest but with a satisfactory RMSD of 3.1deg and 2.0deg, respectively. Moreover the p-value of the pelvic tilt angle is the only angle displaying significant different results between the two EKF implementations. This might be explained by the small angle amplitude of this joint while walking on the treadmill. Moreover, the tracking error is not equally distributed among all joint trajectories due to constraints c1 (similar to joint springs) and to different parameters values in the covariance matrix R.

When compared to reference kinematics (bone pins or fluoroscopy), a systematic review reported errors between 1 and 22deg at the knee joint when tracking retro-reflective markers placed on both thigh and shank segments [START_REF] Begon | Multibody kinematics optimization for the estimation of upper and lower limb human joint kinematics: a systematized methodological review[END_REF]. The proposed method outperforms related studies dealing with a sparse sensor placement [START_REF] Sy | Estimating Lower Limb Kinematics using a Reduced Wearable Sensor Count[END_REF] (RSMD of 10deg) or when using a simple visual observation (RMSD of 9deg [10]) that is often still used by physiotherapists (see section I.A for further details on these studies). When using a single IMU per investigated segment and a quasi-periodic assumption it is possible to reach a better accuracy without monitoring the ankle joints (RMSD of 2.4deg) [10]. Of course this increases the complexity and the price of the overall system. While similar results were observed by comparing the use of a constant acceleration model and a Fourier series expansion in the EKF (see Table . I), it is important to note that preprocessing to obtain these results is very different with the two methods. In the former, the model assumption is inherently inaccurate in human motion. The literature related to the use of EKF for motion analysis rarely proposes methods to model these inaccuracies and to reflect them in the parameters tuning. Cerveri et al. [START_REF] Cerveri | Real-time human motion estimation using biomechanical models and non-linear state-space filters[END_REF] were the sole authors to propose a frequency based approach but this requires to have a prior knowledge of the joint angle evolution, applies to cyclic motions and was not tested with low-cost sensors data. We have shown in a previous study [START_REF] Colombel | Physically Consistent Whole-Body Kinematics Assessment Based on an RGB-D Sensor. Application to Simple Rehabilitation Exercises[END_REF], using a RGB-D camera and a markerless joint center estimate algorithm, that an optimization process outperforms frequency based method and is absolutely required to tune the EKF parameters. Nevertheless, when using markerless data as input and despite an optimal tuning of the EKF parameters and a simplified experimental setup [START_REF] Colombel | Physically Consistent Whole-Body Kinematics Assessment Based on an RGB-D Sensor. Application to Simple Rehabilitation Exercises[END_REF], it was not possible to achieve a joint estimate RMSD lower than 9.7deg on average when compared to a reference SS. Specifically, for tasks closely related to walking, a RMSD up to 20deg was observed for the hip joints. This emphasizes that markerless approaches are not yet accurate enough to estimate human motion, and the superiority of our proposed system being a good and affordable compromise between reference SS and markerless systems.

Unfortunately, an optimal parameters tuning approach is time consuming, since it requires to run the EKF several thousands of times leading to at least 10h of calculation time. Moreover, it can not be performed independently from a SS. Although this optimisation is task dependant and can be performed once for different subjects, it requires that one subject performs the task wearing simultaneously VIMUs and retro-reflective markers at least once. In the latter, when using a Fourier series model, the EKF parameters can be directly set based on a priori knowledge. As exemplified in Fig. 9, the Fourier series converge towards pseudo-constant values, confirming, as expected, the quasi-periodicity of the joint trajectories. The quasi-periodicity of gait was already used for motion analysis in the literature [10]- [START_REF] Ziegler | Simultaneous identification of human body model parameters and gait trajectory from 3D motion capture data[END_REF] and seems to be a valid assumption even for patients.

Regarding values obtained by the optimization process for the tuning of the EKF-CA, they are difficult to be interpreted as the elements of Q and R are interdependent. However we can observe that for the noise covariance matrix they are of the same order of magnitude than the ones set empirically for the EKF-FS and that for the joint soft constraints the ratio of one for ten between some angles is also respected. Despite the differences in the two EKF implementations, the paired t-test results obtained for GDI and most gait kinematics (see Tables I andIII) were not statistically significant. This suggests that both EKF-CA and EKF-FS can be used interchangeably for gait kinematics estimation while emphasizing the relevance of using an optimal tuning of the EKF parameters and/or a realistic process model that describes the joints evolution during movement. The position of the joint centers were estimated thanks to an anatomical calibration even if functional calibration movements could have been used [START_REF] Seel | IMU-Based Joint Angle Measurement for Gait Analysis[END_REF]. However, the proposed system targets patients that might not be able to perform exciting postures or motions [START_REF] Seel | IMU-Based Joint Angle Measurement for Gait Analysis[END_REF]. Nevertheless, for patients having a certain level of recovery, such as patients at home, the fact that with the proposed system position and orientation data could be directly provided by the AR markers would make the estimate of joint centers straightforward [START_REF] Seel | IMU-Based Joint Angle Measurement for Gait Analysis[END_REF]. Note that in both EKF formulations, the difference in joint angles estimates might be explained by the constant body segment lengths, computed statically as the distance between subsequent joint centre positions, as well as the VIMU-to-body-segment translational and rotational offsets. The latter may occur due to inaccuracies in AR markers measurements at initial static calibration, soft tissue artifacts, and the assumption of sensors-tobody-segments rigid connection. In such case, both segments lengths and VIMUs local poses could be modelled in the EKF but at the cost of increasing the state-vector dimensions.

The NRMSD of the strides lengths estimated with the SS and estimated when using the corrected EKF data was on average 1.5cm or 2% for both implementation cases (see Table II). This result is better than the 3cm reported in a recent study when using IMU data collected at both heel levels [START_REF] Hao | Smoother-Based 3-D Foot Trajectory Estimation Using Inertial Sensors[END_REF]. Many studies have intended to estimate similar spatial gait parameters. For instance, Kose et al. reported a maximal error up to 3% for step length estimate [START_REF] Kose | Bilateral step length estimation using a single inertial measurement unit attached to the pelvis[END_REF]. However, their approach, despite its ease of use, suffers from signal features dependency limiting its scope to normal walking. It relies directly on the IMU signal for heel strike detection and assumes a support leg in full extension. The approach proposed by Ahmed et al. uses a camera but only at the head level and thus it suffers from the same problems [START_REF] Ahmed | A Visual-Inertial Approach to Human Gait Estimation[END_REF]. These authors do not directly compute stride or step length, however their foot trajectory estimates, even if acquired over long travelled distance, showed a relatively large RMSD varying from 5 to 13cm. A recent study compared different spatio-temporal gait parameters between healthy and amputee patients [START_REF] Kim | Wearable Sensors Quantify Mobility in People with Lower Limb Amputation during Daily Life[END_REF]. The patient's average stride length was about 1.05 m, while the normal walking was on average of 1.24 m. This difference is much larger than the accuracy provided by the proposed system. Thus, it would be easy to distinguish between the different walking modes.

Having correctly estimated joint trajectories and stride lengths allow to estimate other indexes such as the GDI. The GDI difference between each EKF and SS data was not statistically significant. Moreover, GDI results were all very close to 100. By definition, a GDI greater or equal to 100 means that the test strides are close enough to the healthy strides to be classified as normal ones. The obtained GDI illustrates that the joint trajectories, as well as stride lengths, estimated using the proposed approach have the potential to provide clinically useful indicators. In addition, from the EKF results, the evolution of the Fourier coefficients could be further used to evaluate the gait variability of patients [START_REF] Tsang | Characterization of statistical persistence in joint angle variation during walking[END_REF].

The proposed system is composed of two RGB cameras. This number could be easily increased to match a regular motion capture system and target more complex tasks where the subject will not be facing the same set of cameras. Note that unlike classical stereophotogrammetric approaches, the proposed one is less sensitive to visual occlusion thanks to the simultaneous use of IMU data. However, as previously demonstrated [START_REF] Mallat | Upper Limbs Kinematics Estimation Using Affordable Visual-Inertial Sensors[END_REF], the subject should not be too long away from the camera(s) field of view as the IMUs drift will impact the accuracy of the joint angle estimate. This is inherent to any purely IMU based system [10]. Another limitation of purely IMU based system is that the absolute position of the body with respect to the ground, and for instance to a forceplate, remains undetermined which impedes joint kinetics to be computed. The addition of externally fixed RGB camera(s) solves this issue. This way, the proposed affordable system shall be suitable for both kinematics and kinetics. The proposed system could be used outside of the laboratory but it will still require trained staff to perform the anatomical calibration. It could be used to monitor gait in small clinical centers or at doctor or physiotherapist practice for example. This is already the case of much more expensive devices that can monitor solely the knee joint angle [START_REF] Yu | The KneeKG system: a review of the literature[END_REF].

VI. CONCLUSION

A new approach for motion analysis based on the use of VIMU sensors and an EKF has been proposed and validated in this paper. VIMU sensors can be considered a good compromise between, on the one hand, SS systems which are costly and cumbersome and, on the other hand, purely IMU based methods which requires one IMU per segment and advanced anatomical calibrations to reach an accuracy below 5deg. On the computational side, the propose EKF-FS represent one of the most advance implementations of this approach (including constraints and a clear rationale for the filter parameter tuning).

Fig. 1 .

 1 Fig. 1. Overview of the proposed approach based on Visual-Inertial Measurement Units (VIMU) and an Extended Kalman Filter (EKF) taking advantage of the human body and task biomechanics.

Fig. 2 .

 2 Fig. 2. (a) Lower-limbs multi-body model composed of segments, of lengths , = 1, .., , and Degrees-of-Freedoms (DoF) denoted , = 1, .., . (b) Experimental setup, three VIMU attached to the pelvis, left, and right heels in a bracelet form are collecting data during gait thanks to two RGB cameras located behind the treadmill.

Fig. 3 .

 3 Fig. 3. (a) The VIMU prototype with three different bracelet arrangements to fit different subjects and/or segment sections. (b) Description of the coordinate systems associated to a VIMU. The coordinate system of the IMU sensor, left, and right AR markers are all linked to the central marker coordinate system through rotation matrices R and a translation vectors r.

Fig. 4 .

 4 Fig. 4. Example of lower-limbs anatomical landmarks pin-pointed using a calibration wand base on an AR marker.

Fig. 9 .

 9 Fig. 9. Representative evolution of the Fourier series coefficients corresponding to the ankle rotation joint ( 12 ) estimated by the EKF-FS.

-specific calibration Subject-specific calibration

  1, and is based on the simultaneous use of affordable IMUs, incorporating only 3D accelerometers and gyroscopes sensors, i.e., no magnetometers, and a set of AR markers

		Pre-processing	Extended Kalman Filter
	RGB Cameras IMU subject New	Raw 3D linear acc. and ang. vel. Cameras calibration (OpenCV) Cameras intrinsic and extrinsic parameters static wand IMU calib. VIMU calib. calib. using a VIMU-to-segment Initial joint configuration VIMU local poses SensorsJoint centres estimate Calibrated linear acc. and ang. vel.	Measurement update kinematic model -VIMU rigid contraints -Joint soft constraints -Lower-limbs	and acc. velocities positions, Prediction accel. model + optimal tuning Case 2: Fourier serie model Predicted state: joint Case 1: Constant
		OUTPUT
			-Pelvis pose and lower-limbs joint kinematics
			-Stride length estimate
			-Gait deviation index

TABLE I RESULTS

 I OF THE COMPARISON BETWEEN EIGHTEEN LOWER-LIMBS JOINT TRAJECTORIES OBTAINED USING THE PROPOSED AFFORDABLE SYSTEM AND THE SS. RESULTS HAVE BEEN REPORTED AS MEAN ± SD OVER ALL SUBJECTS AND RESULTS OF HIP, KNEE, AND ANKLE JOINTS HAVE BEEN OBTAINED OVER THE AVERAGE OF LEFT AND RIGHT LEG VALUES. RESULTS OF THE PAIRED T-TEST BETWEEN EKF-CA AND EKF-FS FOR RMSD AND CC OF EVERY JOINT TRAJECTORY HAVE ALSO BEEN REPORTED.

		Constant acceleration model assumption	Fourier series model assumption	p-Values
	Pelvic Translations	RMSD [cm]	CC	RMSD [cm]	CC	RMSD	CC
	Translation Y 1	0.7±0.2	0.96±0.03	0.7±0.2 0.94±0.05	0.444	0.608
	Translation X 2	0.7±0.1	0.96±0.02	0.8±0.1 0.96±0.03	0.472	0.273
	Translation Z 3	0.5±0.2	0.59±0.19	0.7±0.2 0.56±0.25	0.182	0.796
	Mean	0.6±0.1	0.84±0.21	0.7±0.0 0.82±0.22
	Joint Angles	RMSD [deg]	CC	RMSD [deg]	CC	RMSD	CC
	Pelvis tilt	1.7±0.8	0.54±0.26	2.4±0.7 0.43±0.22	0.001	0.039
	Pelvis obliquity	2.0±0.8	0.56±0.4	1.7±0.6 0.77±0.13	0.372	0.136
	Pelvis rotation	2.2±0.8	0.70±0.19	2.4±0.7 0.69±0.29	0.610	0.382
	Hip Flex./Ext.	3.6±0.9	0.95±0.02	3.9±1.1 0.93±0.02	0.200	0.071
	Hip Abd./Add.	2.6±1.3	0.64±0.18	2.7±1.2 0.57±0.25	0.583	0.205
	Hip Int./Ext. Rotation	4.1±1.5	0.79±0.20	4.2±1.1 0.73±0.25	0.872	0.421
	Knee Flex./Ext.	5.1±1.1	0.91±0.06	4.8±0.9 0.92±0.04	0.152	0.127
	Ankle Flex./Ext.	3.1±0.8	0.42±0.27	3.1±1.1 0.45±0.25	0.841	0.941
	Ankle Inv/Eve.	5.0±1.6	0.72±0.20	4.7±1.5 0.64±0.25	0.516	0.085
	Mean	3.5±1.2	0.71±0.17	3.5±1.0 0.69±0.18

TABLE II COMPARISON

 II BETWEEN USING THE AFFORDABLE RAW/SS DATA AND THE DATA ESTIMATED BY EKF, FOR OBTAINING THE STRIDE LENGTHS AND GAIT DEVIATION INDEX (GDI), OVER ALL SUBJECTS AND STRIDES.

		Raw data	EKF-CA	EKF-FS
	Stride lengths	RMSD [cm]	NRMSD [%]	RMSD [cm]	NRMSD [%]	RMSD [cm]	NRMSD [%]
	Right	2.6±1.1	3.8±1.6	1.7±0.0	2.3±0.5	1.4±0.0	1.9±0.2
	Left	2.6±1.0	3.8±1.5	1.6±0.0	2.2±0.8	1.3±0.0	1.7±0.5
		SS data	EKF-CA	EKF-FS
	GDI	99.66±3.87 98.44±7.82 98.38±7.97

Markers-to-Markers &Markers-to-IMU 3D rigid pose

ACKNOWLEDGMENT

We thank all subjects of the GVLAB for their participation in the experiments as well as Patrick Karasinski and Florian Huet for their help in developing the VIMU prototype.


Khaled Metwally

Chiara Bastiancich

Florian Correard

Anthony Novell

Samantha Fernandez

Benjamin Guillet

Benoit Larrat

Serge Mensah

Marie-Anne Estève

Anabela Da Silva
email: anabela.dasilva@fresnel.fr

A multipurpose integrated preclinical device for the treatment of glioblastoma

A multipurpose integrated preclinical device that combines Blood-Brain Barrier permeabilization, Photo-Thermal Therapy (PTT) with photoacoustic thermometry, was developed as a novel preclinical versatile research instrument for the development of new PTT for Glioblastoma multiforme.

Introduction

Glioblastoma multiforme (GBM) is one of the most common and aggressive malignant brain tumors in adults. As current therapies are ineffective, there is an urgent need to find novel and efficient therapeutic strategies.

One major limiting factor for brain tumors chemotherapies is the Blood-Brain Barrier (BBB), which prevents drugs from spreading from the vascular network to the central nervous system (CNS) limiting their penetration at the tumor site. To permeabilize the BBB and enhance drug delivery into the CNS, we use Focused Ultrasound Sonication (FUS) combined with intravascular microbubbles. Different in vitro and in vivo studies have demonstrated the interest of Photo-Thermal Therapy (PTT) in GBM treatment, chosen here as the main therapy. To improve heating performances and localization of the therapy, metallic nanoparticles (NPs) have emerged as potential photoabsorbing agents for PTT. Recently, laser-synthesized Titanium Nitride NPs (TiN NPs) have proven as new candidates [START_REF] Popov | Laser-synthesized TiN nanoparticles as promising plasmonic alternative for biomedical applications[END_REF][START_REF] Bastiancich | Photothermal therapy for the treatment of Glioblastoma: potential and preclinical challenges[END_REF], for which we have demonstrated in vitro (U87-MG GBM cell cultures) a low cytotoxicity combined with a good cellular uptake and PTT effect. To monitor the treatment, a photoacoustic (PA) thermoetry system was implemented for monitoring the temperature in real-time, in situ and in depth [START_REF] Larina | Real-time optoacoustic monitoring of temperature in tissues[END_REF][START_REF] Shah | Photoacoustic imaging and temperature measurement for photothermal cancer therapy[END_REF][START_REF] Wang | Photoacoustic temperature measurements for monitoring of thermal therapy[END_REF][START_REF] Nikitin | Temperature dependence of the optoacoustic transformation efficiency in ex vivo tissues for application in monitoring thermal therapies[END_REF]. In the present study, we present this novel preclinical integrated device that allows performing on the same bench BBB opening with confocal FUS to improve local accumulation of NPs at the tumor site, highly localized PTT and non-invasive in situ PA thermometry.

The next section presents the fully integrated device developed. To demonstrate that the different features of this device are fully operational, we performed preliminary proof-of-concept in vivo studies that will be presented at the conference.

Description of the integrated device

A schema of the developed multi-functional setup is presented in Fig 1.

The device can be decomposed into three parts: -The FUS/PTT/PA probe (Fig 1, Bottom) comprises a single-element focused transducer, presenting a hole in its center. The ultrasound transducer was then coupled to the head of the animal via the water balloon filled with deionized and degassed water to reach 7 to 10 mm thickness, depending on the targeted depth, and covered with degassed echographic gel. Electrical power sent to the transducer was monitored during the BBB opening session. The PA signal is collected through a digital oscilloscope. A hermetic glass window has been mounted and positioned on the top of the transducer mount in order to receive an optical fiber through which the light sources are delivered to the sample.

-The light delivery part includes two types of lasers, whose beams were piped through the same optical fiber: a CW laser diode for the PTT (830 nm, 2 W), controlled in current and temperature and a pulsed laser for the PA (Nd:Yag, 1064 nm, 330 mJ, 5 ns, 10 Hz). The pulsed laser was frequency doubled with a KTP crystal at 532 nm.

-A platform comprising a small animal holder, a transducer holder mounted on fully programmable 3D scanning stages and a programmable ultrasound single channel amplifier. Hot water is conveyed to the animal bed via plastic tubes to keep the body temperature of the animal stable during the therapy. The temperature of the circulating water is stabilized.

The different acoustic and optical emissions and data acquisition were controlled with a computer. The proposed FUS/PTT protocol includes three consecutive steps: i) systemic administration efficient photoabsorbers (NPs or others); ii) injection of intravascular microbubbles followed by the application of FUS, by using the transducer in its active mode, with calibrated sonication parameters allowing BBB opening; iii) application of the PTT with PA temperature monitoring with a calibrated optical irradiation dose inducing tumor cell death.

PhotoThermal Therapy and temperature monitoring with PhotoAcoustics: Preliminary tests on subcutaneous GBM model

The laser parameters selected for the in vivo studies (laser wavelength, power and heating duration) were based on values taken from literature [START_REF] Bastiancich | Photothermal therapy for the treatment of Glioblastoma: potential and preclinical challenges[END_REF] and the efficacy results obtained in vitro on U87-MG spheroids using TiN NP [START_REF] Popov | Laser-synthesized TiN nanoparticles as promising plasmonic alternative for biomedical applications[END_REF]. Within the linear approximation, for water-based and fatty tissues, the photoacoustic signal can be directly related to temperature [START_REF] Shah | Photoacoustic imaging and temperature measurement for photothermal cancer therapy[END_REF][START_REF] Wang | Photoacoustic temperature measurements for monitoring of thermal therapy[END_REF][START_REF] Ke | Photoacoustic thermography of tissue[END_REF]: ∆P/P o = α∆T , where P o is the amplitude of the PA signal at a reference temperature T o and α is a calibration factor to be evaluated for the probed sample.

Here is reported a preliminary illustrative experiment of PTT and PA thermometry, performed on mice bearing subcutaneous tumors (n=6). The group was devided into two groups, one (n=3) without NPs, the second (n=3) with NPs. The temperature variations during PTT (∆T), estimated by considering a previously determined calibration factor are plotted in Fig. 2.

Conclusion

An original multipurpose integrated setup has been developed for efficient PTT of the glioblastoma in small animals [START_REF] Metwally | Development of a multi-functional preclinical device for the treatment of glioblastoma[END_REF]. It allows combining for the first time i) BBB opening by FUS for an increased efficiency in the targeted drug delivery to the brain; ii) PTT and iii) in situ temperature monitoring with PA. Sub-cutaneous and transcranial preliminary experiments were performed and demonstrated that the whole device is functional. Through skull PTT has also been shown feasible.

Fig. 1 :

 1 Fig. 1: Schema of the integrated experimental setup (PBS: Polarizing beam splitter; BS: Beam sampler; KTP: Frequency doubling KTP crystal; KG3: Filter glass (Bandpass: 315 -710 nm); PD: Si Photodetector (Detector, 200 -1100 nm); λ /2: Achromatic half-wave plate).

Fig. 2 :

 2 Fig. 2: Evolution of the estimated temperature rise during the PTT on mice from group without NPs (Left) and with NPs (Right).


Julien Signoles
email: julien.signoles@cea.fr

The E-ACSL Perspective on Runtime Assertion Checking

Keywords: Software and its engineering → Specification languages, Dynamic analysis, Source code generation, Software verification, Automated static analysis Runtime Assertion Checking, Behavioral Interface Specification Language, Source Code Generation, Combination of Static and Dynamic Analyses

Runtime Assertion Checking (RAC) is the discipline of verifying program assertions at runtime, i.e. when executing the code. Nowadays, RAC usually relies on Behavioral Interface Specification Languages (BISL) à la Eiffel for writing powerful code specifications. Since now more than 20 years, several works have studied RAC. Most of them have focused on BISL. Some others have also considered combinations of RAC with others techniques, e.g. deductive verification (DV). Very few tackle RAC as a verification technique that soundly generates efficient code from formal annotations. Here, we revisit these three RAC's research areas by emphasizing the works done in E-ACSL, which is both a BISL and a RAC tool for C code. We also compare it to others languages and tools.

includes others kinds of formal annotations, such as function contracts and invariants. Yet, at that time, RAC approaches remained very practical and were not studied from a theoretical perspective.

Since now more than 20 years, several works have studied RAC. Most of them have focused on BISL. Some others have also considered combinations of RAC with other techniques, e.g. deductive verification (DV). Very few tackle RAC as a verification technique that soundly generates efficient code from formal annotations. Here, we revisit these three RAC's research areas by emphasizing the works done in E-ACSL, which is both a BISL [START_REF] Delahaye | Common Specification Language for Static and Dynamic Analysis of C Programs[END_REF] and a RAC tool [START_REF] Signoles | E-ACSL, a Runtime Verification Tool for Safety and Security of C Programs[END_REF] for C code. We also compare it to others languages and tools.

BEHAVIORAL INTERFACE SPECIFICATION LANGUAGES

After the success of Eiffel, new BISL were designed for mainstream programming languages [START_REF] Hatcliff | Behavioral Interface Specification Languages[END_REF] at the turn of the millennium: first JML [START_REF] Leavens | JML: A Notation for Detailed Design[END_REF] for Java, then Spec# [START_REF] Barnett | [END_REF] for C#, VCC [START_REF] Cohen | VCC: A Practical System for Verifying Concurrent C[END_REF], ACSL [START_REF] Baudin | ACSL: ANSI/ISO C Specification Language[END_REF] and E-ACSL [START_REF] Delahaye | Common Specification Language for Static and Dynamic Analysis of C Programs[END_REF] for C, CodeContract [START_REF] Fähndrich | Embedded Contract Languages[END_REF] for .NET, Spark2014 [START_REF] Hoang | SPARK 2014 and GNATprove -A competition report from builders of an industrial-strength verifying compiler[END_REF] for Ada and, more recently, GOSPEL [START_REF] Charguéraud | GOSPEL -Providing OCaml with a Formal Specification Language[END_REF] for OCaml. Dedicated BISL such as Boogie [START_REF] Leino | This is Boogie 2[END_REF] and WhyML [START_REF] Filliâtre | Why3 -Where Programs Meet Provers[END_REF] were also designed in the meantime: they do not target a mainstream specification language but, instead, propose their own general-purpose specification and programming framework.

Even if formal specifications provided by BISL are useful in their own rights as precise code documentation, their main usage consists in checking them with verification tools. RAC is particularly suitable and lightweight for this purpose. Yet, DV à la Dijkstra [START_REF] Hähnle | Deductive Software Verification: From Pen-and-Paper Proofs to Industrial Tools[END_REF] is the tool of choice for proving them for all possible executions. Abstract interpretation à la Cousot [START_REF] Rival | Introduction to Static Analysis: An Abstract Interpretation Perspective[END_REF] is also possible [START_REF] Fähndrich | Static Contract Checking with Abstract Interpretation[END_REF][START_REF] Ouadjaout | A Library Modeling Language for the Static Analysis of C Programs[END_REF], as well as model checking [START_REF] Beckert | Modular Verification of JML Contracts Using Bounded Model Checking[END_REF], even if several logical properties cannot be easily verified with these techniques. Nowadays, most tools based on BISL support both RAC and DV, e.g. OpenJML [START_REF] Cok | OpenJML: JML for Java 7 by Extending OpenJDK[END_REF] for JML, Spec#, Spark2014, GOSPEL and, since recently, Why3 for WhyML. Frama-C [START_REF] Kirchner | Frama-C: A Software Analysis Perspective[END_REF] relies on ACSL for DV and its subset E-ACSL for both RAC and abstract interpretation, E-ACSL being compatible with ACSL.

Core Features

Beyond code assertions, the core feature of BISL are function contracts that allow the user to specify what a function assumes from its callers (its preconditions), and what its implementation guarantees (its post-conditions), assuming the preconditions are satisfied.

Figure 1 shows a function contract for a function binary_search written in the E-ACSL language. E-ACSL, being defined as a syntactic subset of ACSL, it is also an ACSL specification. Preconditions (resp. post-conditions) are introduced with the keyword requires (resp. ensures). The second precondition shows the use of the built-in predicate \valid, stating that a range of memory locations (here, the memory cells of a between 0 and len-1) has been properly allocated. Such memory predicates are specific to ACSL and E-ACSL that handles C code: the other BISL tackle higher-level programming languages that ensure memory safety by design. Yet, such properties are critical for C code. The third precondition calls the user-defined predicate sorted: most BISL provide the ability to define such predicates and logic functions. This one relies on an universal quantification: BISL are usually based on first-order logic. This specification also defines two behaviors (exists and not_exists) guarded by assumes clauses that specify when each of them are activated: behaviors provide a convenient way to distinguish different cases in a specification. complete and disjoint behaviors specifies that the behaviors covers all possible cases and do not overlap, respectively. Most BISL, including JML, Spark2014, ACSL and E-ACSL, provide such a notion of behavior or, at least, provide a way to specify exceptional cases, e.g. Spec# and WhyML. The assigns clause specifies a frame condition, i.e. the modifications in the program memory made by the function. Here, \nothing means that the function is effect-free (from an observational point of view). Common features of most BISL also include data (or object) invariants, loop invariants, loop variants (i.e., decreasing measures inside loops), and ghost code. Ghost code is a powerful way to write code as specification without interfering with the original code [START_REF] Filliâtre | The Spirit of Ghost Code[END_REF]. Another common feature consists in referring to the old value of some memory location, i.e. its value in the the initial state of the contract from a post-condition. For instance, writing ensures G == \old(G)+1; for a function f means that the global variable G is incremented by one after having executed f . ACSL and E-ACSL generalize this feature to any program point thanks to the \at keyword: \at(x,L) refers to the value of x at program point L. They are the only BISL providing this generalized form. In the following, we will refer to such properties as multi-state properties.

Main Semantic Differences between BISL

BISL are usually based on typed first order logic. Their semantic differences mainly come from whether the languages target more RAC (e.g., JML, Spec#, Spark2014 and E-ACSL) or DV (e.g., WhyML and ACSL).

The first difference is about numbers. All BISL support bounded numbers (machine integers and floating-point numbers), with two possible flavors when they overflow: either raising an error, or reducing into the bounds with the appropriate modulo operation. For instance, JML provides both mode, while ACSL and E-ACSL use the latter. Yet, even if convenient for formal reasoning, providing mathematical numbers (integers over Z and real numbers) were historically not provided for BISL targeting RAC because they are hard or even impossible (for real numbers) to be computed in finite time. Nowadays, JML and Spark2014 provide a mode that supports them, while they are natively provided in E-ACSL. They are also provided in all BISL that target DV.

A related difference is about quantifiers that must be bounded for RAC in order to be executable in finite time. Indeed, iterating over all natural numbers is not possible at runtime. Even if bounded quantifiers restrict the BISL's expressive power, it is not an issue in practice since they affect only very few mathematical properties that rarely occur in real code.

A key difference is related to partial functions (e.g., a division): in classical first order logic implemented in theorem provers used by DV tools, such functions are represented by under-specified total functions, meaning that terms such as 1/0 are well defined, even if their values are not specified. However, such terms and predicates cannot be evaluated without error at runtime and so are problematic for RAC. This issue is known as the undefinedness problem [START_REF] Cheon | A runtime assertion checker for the Java Modeling Language[END_REF]. For solving it, Chalin have proposed [START_REF] Chalin | A Sound Assertion Semantics for the Dependable Systems Evolution Verifying Compiler[END_REF] to rely on the strong validity principle [START_REF] Konikowska | A three-valued logic for software specification and validation[END_REF], which considers that a property is valid if and only if it is valid and defined. This semantics is consistent with the classical semantics, since any well-defined predicate has the same truth value in both. Nowadays, RAC-oriented BISL that also support deductive verification (e.g. JML, Spark2014, and E-ACSL) follow this semantics.

Another difference is pure functions, which are functions from the underlying programming language that can be safely used in specifications. For that purpose, they must always terminate and be side-effect free. Such functions are usually supported by RACoriented BISL (JML, Spark2014, and Spec#, whilst the E-ACSL tool has a mode that supports them), while they are not supported by DV-oriented BISL (ACSL, Why3, and Krakatoa [START_REF] Filliâtre | The Why/Krakatoa/Caduceus Platform for Deductive Program Verification[END_REF], a variant of JML based on Why, the ancester of Why3). Indeed, pure methods are easy to execute for runtime assertion checkers, but lead to several issues for consistently verify them for deductive verifiers [START_REF] Kosmatov | Static vs Dynamic Verification in Why3, Frama-C and SPARK[END_REF].

Last, one more semantic difference between BISL is actually independent from being RAC oriented or not. Indeed, frame conditions may have two flavors: assigns (or writes), that make explicit the memory locations that may be written to, and modifies, that make explicit the memory locations whose values may be changed. The former is stricter than the latter since, in the latter case, the other memory locations may change if (and only if) their value is restored before the end of the scope of the contract. There is currently no consensus between both semantics. For instance, JML relies on the assigns flavor, while ACSL and E-ACSL relies on the modifies flavor, even though they use the keyword assigns.

USING RUNTIME ASSERTION CHECKING IN PRACTICE

Checking annotations at runtime has various applications. First, when directly applied, RAC helps to strengthen code [START_REF] Casalnuovo | Assert Use in GitHub Projects[END_REF][START_REF] Kochhar | Revisiting Assert Use in GitHub Projects[END_REF] by discovering bugs upstream, when a well-identified assertion is violated, as opposed to downstream, by observing a consequence of the bug in the form of an unexpected outcome, typically a crash or a weird output. Second, combining it with testing techniques, e.g. fuzzing, allows the user to detect more issues than using these techniques alone, since formal annotations make more invalid behaviors (e.g., broken invariants or invalid API uses) observable. For instance, a buffer overflow that would overwrite data in a supposedly-unchanged memory area A would usually remain undetected with standard testing techniques if the test verdict would not depend on A.

Third, RAC may be combined with other verification techniques, notably DV, in various ways [START_REF] Kosmatov | Static vs Dynamic Verification in Why3, Frama-C and SPARK[END_REF][START_REF] Maurica | Runtime Assertion Checking and Static Verification: Collaborative Partners[END_REF], in particular when they are provided within a shared framework, such as Frama-C, OpenJML, Spark2014, Spec#, Why3, or GOSPEL. An usual combination with any static technique (typically, abstract interpretation, model checking, or DV) consists in discharging as properties as possible statically, and relying on RAC to verify the remaining ones at runtime: it lowers the amount of effort and the level of expertise required by the static verifications, while reducing the runtime overheads induced by RAC. When combined with DV, RAC also helps debug formal annotations before proving them, find counterexamples in case of proof failure and, more generally, help to understand why a property is not automatically proven [START_REF] Becker | Explaining Counterexamples with Giant-Step Assertion Checking[END_REF][START_REF] Petiot | How testing helps to diagnose proof failures[END_REF].

Besides, even if not new, a recent trend consists in seeing BISL as low-level code-oriented specification languages which high-level properties may be generated to. It leads to enhance the area of applications of tools based on such languages, including RAC. Such highlevel properties include liveness properties [START_REF] Giorgetti | Verification of class liveness properties with Java Modeling Language[END_REF], non-interference properties [START_REF] Barany | Hybrid Information Flow Analysis for Real-World C Code[END_REF], security automata [START_REF] Huisman | A Formal Connection between Security Automata and JML Annotations[END_REF], relational properties [START_REF] Blatter | Static and Dynamic Verification of Relational Properties on Self-composed C Code[END_REF], and system-wide properties [START_REF] Robles | Tame Your Annotations with MetAcsl: Specifying, Testing and Proving High-Level Properties[END_REF]. In the same spirit, implicit properties can be made explicit thanks to BISL. An example is security weaknesses [START_REF] Pariente | Static Analysis and Runtime Assertion Checking: Contribution to Security Counter-Measures[END_REF] for which the implicit semantics depends on the underlying programming language. It worth noting that these approaches usually combine RAC with other verification techniques.

COMPILING FORMAL ASSERTIONS

Checking assertions at runtime requires to compile them into executable code. This process has only received little attention [START_REF] Cheon | A runtime assertion checker for the Java Modeling Language[END_REF][START_REF] Lehner | A Formal Definition of JML in Coq and its Application to Runtime Assertion Checking[END_REF] up to now. It might be explained by two reasons. First, it looks like most BISL's designers and tool providers do not consider it as challenging enough as expressed by the authors of Spec#: who state that "The run-time checker is straightforward" [START_REF] Barnett | [END_REF], while dedicating the rest of their paper to DV. Second, even if designing procedure for checking properties at runtime is at the heart of Runtime Verification, the relevant community is more appealed by checking temporal properties [START_REF] Falcone | A Tutorial on Runtime Verification[END_REF] than behavioral properties.

Yet, generating sound and efficient code from formal properties as powerful as the ones allowed in BISL is not straightforward. The rest of this section considers a few concrete examples for supporting this claim. All of them are related to research works that have been carried out within E-ACSL [START_REF] Signoles | E-ACSL, a Runtime Verification Tool for Safety and Security of C Programs[END_REF], the runtime assertion checker of Frama-C [START_REF] Kirchner | Frama-C: A Software Analysis Perspective[END_REF]. While Frama-C supports ACSL as specification language, its E-ACSL plug-in focuses on E-ACSL, its executable subset whose semantics is based on strong validity.

Three other properties are usually important for RAC tools: trustworthyness, transparency and isolation [START_REF] Kandziora | Runtime assertion checking of JML annotations in multithreaded applications with e-OpenJML[END_REF]. Trustworthyness means that, if a failure is detected, then there indeed should be one. Transparency means that the monitored program should functionally behave as the original one if no failure is detected. Isolation means that a failure, if any, must be reported at the point where the annotation is written, and not later.

Mathematical Properties

Checking mathematical properties efficiently and soundly at runtime is hard: relying on machine numbers is unsound if BISL supports mathematical numbers, while relying on dedicated libraries that encode mathematical numbers in the target programming language, such as GMP1 for C code, is not efficient.

For being both efficient and sound, E-ACSL relies on a dedicated static type system [START_REF] Kosmatov | Efficient Runtime Assertion Checking for Properties over Mathematical Numbers[END_REF] that automatically infers a precise type in which it is sound to compute a mathematical operation over mathematical integers in Z or rationals in Q. For instance, on a 64bit architecture, if x is a variable of type int, it infers that x+1 may be soundly computed in type long: x+1 might overflow if computed over int, but it would be unnecessary inefficient to use a GMP addition in such a context. Spark2014 has adapted the E-ACSL's type system to its context, but does not support rational numbers yet: as far as we know, E-ACSL is the only RAC tool that support properties over rational numbers.

Even if the code generated thanks to this type system has been proven to be very efficient in practice, it is worth noting that additional research should be carried out for handling definitions of recursive logic functions and inductive predicates in a precise and sound way. Also, no RAC tool currently supports irrational numbers. Actually, it is not possible to get a verdict at runtime in finite time in general in that context (for instance, for equalities over irrational numbers). Yet practical solutions based on incomplete verdicts should be possible.

Multi-State Properties

Checking multi-state properties, such as P: x == \at(x,L)+1, at runtime might seem straightforward: it just requires to save the value of x at program point L and use it when evaluating P. The first versions of E-ACSL followed this way. Yet, it is unsound in the general case. Consider for instance the predicate Q: \forall integer i,j; 0<=i<LEN ==> 0<=j<i ==> t[i][j] == \at(t[i][j],L)+1: the variables i and j are local to predicate Q and so cannot be used at L. Soundness requires to copy the relevant memory cells of t. Yet, knowing at compile time the exact memory cells that must be copied out is undecidable in the general case. For instance, if LEN is a variable, it might be unbound at L, or its value might have changed between L and the definition of Q. To circumvent this issue, existing approaches copy the whole memory block that contains t [START_REF] Kosiuczenko | An Abstract Machine for the Old Value Retrieval[END_REF][START_REF] Petiot | Instrumentation of Annotated C Programs for Test Generation[END_REF]. Yet, it is not always efficient since it often copies too many memory cells. Therefore, the current version of E-ACSL implements an hybrid approach that statically computes an over-approximation of the memory cells that need to be copied out. In the best case, it is similar to the former E-ACSL implementation (but sound), whilst, in the worst case, it is similar to the other approaches. Therefore, it reconciles soundness and efficiency. Yet, it is worth noting, that no approach solves efficiently the problem in the general case: there is still room for improvements here.

Memory Properties

As already explained, memory properties, such as the predicate valid(p+(0..len-1)) used in Figure 1, are specific to ACSL and E-ACSL. Yet, checking such properties at runtime is the purpose of memory debuggers, which focus in detecting memory errors at runtime, such as buffer and heap overflows, and accesses to uninitialized data in the program memory.

In this context, the most efficient technique is memory shadowing that stores pieces of information about the program memory in order to detect memory issues. Memory shadowing usually associates addresses from program memory to values stored in a disjoint memory space, called shadow memory, and accessible in constant time. The way of structuring these shadow values is refered to as shadow state encoding and may vary from one tool to another. One of the most efficient memory debuggers is currently AddressSanitizer [START_REF] Serebryany | AddressSanitizer: A Fast Address Sanity Checker[END_REF]. Despite their efficiency, classical memory shadowing techniques are not able to detect all memory defects, in particular block-level memory errors (typically, accessing to an allocated block from another block) [START_REF] Jakobsson | Fast as a Shadow, Expressive as a Tree: Optimized Memory Monitoring for C[END_REF] or temporal memory errors (typically, accessing to a pointed-to object that is not the same as when the pointer was created) [START_REF] Simpson | MemSafe: Ensuring the Spatial and Temporal Memory Safety of C at Runtime[END_REF][START_REF] Vorobyov | Runtime Detection of Temporal Memory Errors[END_REF]. Even if cleverly using the space left between aligned allocated blocks help circumvent partially this issue [START_REF] Jones | Backwards-compatible bounds checking for arrays and pointers in C programs[END_REF][START_REF] Serebryany | AddressSanitizer: A Fast Address Sanity Checker[END_REF], it is not enough for all memory properties.

Beyond memory shadowing, others techniques do exist, such as fat pointers [START_REF] Austin | Efficient Detection of All Pointer and Array Access Errors[END_REF] or dictionary-based approaches [START_REF] Jones | Backwards-compatible bounds checking for arrays and pointers in C programs[END_REF]. The former extends the pointer representation with information about memory block bounds. Yet, it suffers from modifying the size of pointed data, which makes this technique very hard to use safely in practice. The latter associates to each allocated memory address the necessary pieces of information and storing them in a separated dictionary. While powerful, this technique has efficiency issue since the dictionnary is not accessed in constant time.

E-ACSL has first tried a dictionary based approach [START_REF] Kosmatov | An Optimized Memory Monitoring for Runtime Assertion Checking of C Programs[END_REF]. Yet, it was not efficient enough in practice. Then, an hybrid approach mixing a dictionary and a shadow memory has been experimented with better results [START_REF] Jakobsson | Fast as a Shadow, Expressive as a Tree: Optimized Memory Monitoring for C[END_REF]. However, its implementation was hardly maintainable. Finally, we design a shadow-based technique relying on two new shadow state encodings (for the heap and the stack). It is efficient enough in practice [START_REF] Vorobyov | Shadow State Encoding for Efficient Monitoring of Block-level Properties[END_REF], while more expressive than the existing shadow-based techniques in most cases [START_REF] Vorobyov | Detection of Security Vulnerabilities in C Code using Runtime Verification[END_REF]. In particular, it is able to detect that the program incorrectly accesses to a memory block from another memory block, even if both have been properly allocated earlier. It is worth noting that the instrumentation is independent from the underlying representation [START_REF] Ly | Verified Runtime Assertion Checking for Memory Properties[END_REF].

Static analysis can also optimize the instrumentation by preventing monitoring irrelevant memory blocks [START_REF] Ly | Soundness of a Dataflow Analysis for Memory Monitoring[END_REF]. Thanks to the collaborative facilities provided by Frama-C, it is also possible to discharge most memory properties before running E-ACSL [START_REF] Pariente | Static Analysis and Runtime Assertion Checking: Contribution to Security Counter-Measures[END_REF], typically by using Eva [START_REF] Blazy | Structuring Abstract Interpreters through State and Value Abstractions[END_REF], the abstract interpretation-based static analyzer of Frama-C. Both techniques makes the code generated by E-ACSL much more efficient.

CONCLUSION AND FUTURE WORK

BISL are now well understood and more and more used in frameworks that combine RAC and DV, such as Frama-C, OpenJML, Spark2014 and Why3. Yet, the process of generating executable code from formal assertions has not been extensively studied. The works done on E-ACSL, the RAC tool of Frama-C, show that being both sound and efficient is often challenging and several open questions are still to be addressed, including efficient monitoring of properties over irrational numbers, monitoring of advanced multistate properties, and monitoring of advanced memory properties.

Figure 1 :

 1 Figure 1: Example of an E-ACSL function contract.

https://gmplib.org/

ACKNOWLEDGMENTS

We would like to thank Thibaut Benjamin who proofreads an early version of this paper.

We would also like to thank all the persons who have contributed to E-ACSL: Nikolai Kosmatov, Kostyantyn Vorobyov, Fonenantsoa Maurica, Basile Desloges, Thibaut Benjamin, Guillaume Petiot, Arvid Jakobsson, Dara Ly, Romain Maliach-Auguste, and Félix Ridoux.

This work has received funding from the European Union's Horizon 2020 research and innovation programme under grant agreement N • 883242, project ENSURESEC.