R Chellal

L Cuvillon

E Laroche
email: laroche@unistra.fr

Model identification and Vision-Based H ∞ Position Control of 6-DoF Cable-Driven Parallel Robots

Keywords: Cable-driven parallel robots, model identification, vision-based control, H∞ control, tension distribution

This paper presents methodologies for the identification and control of 6-Degrees of Freedom (DoF) Cable-Driven Parallel Robots (CDPRs). First a two-step identification methodology is proposed to accurately estimate the kinematic parameters independently and prior to the dynamic parameters of a physics-based model of CDPRs. Second, an original control scheme is developed, including a vision-based position controller tuned with the H∞ methodology and a cable tension distribution algorithm. The position is controlled in the operational space, making use of the end-effector pose measured by a motion-tracking system. A four-block H∞ design scheme with adjusted weighting filters ensures good trajectory tracking and disturbance rejection properties for the CDPR system, which is a non-linear coupled MIMO system with constrained states. The tension management algorithm generates control signals that maintain the cables under feasible tensions. The paper makes an extensive review of the available methods and presents an extension of one of them. The presented methodologies are evaluated by simulations and experimentally on a redundant 6-DoF INCA 6D CDPR with 8 cables, equipped with a motion-tracking system.

Introduction

Cable-Driven Parallel Robots (CDPRs) are special parallel robot manipulators where the endeffector or platform is connected to the base via cables only, its movement resulting from the winding and unwinding of the cables around pulleys actuated by motors. This kind of robot manipulator appeared in the late 80s by replacing the segments of the Gouth-Stewart platform manipulators by cables in order to solve the major problem of the reduced workspace of rigid parallel robots [START_REF] Dagalakis | Stiffness study of a parallel link robot crane for shipbuilding applications[END_REF].

Compared to serial and parallel manipulator robots actuated by rigid links, cable-actuated robots benefit from interesting features like: very large workspace, high speed motion due to the low moving mass, modular geometry, portability and adaptability to multi-scales. However, their control is a more complex issue as the cables must remain under tension at any time [START_REF] Oh | Cable suspended planar robots with redundant cables: controllers with positive tensions[END_REF]. The approaches available in the literature to solve this issue can be classified in two main categories. In the off-line solutions, a path planning step is used prior to motion in order to design a reference trajectory which guarantees that the cables remain under tension during the predefined motion [START_REF] Gosselin | Global planning of dynamically feasible trajectories for three-DoF spatial cable-suspended parallel robots[END_REF][START_REF] Trevisani | Experimental validation of a trajectory planning method avoiding cable slackness and excessive tension in under-constrained translational planar cable-driven robots[END_REF]. In the on-line solutions, an algorithm of tension distribution (also known as force calculation or redundancy resolution) is used to ensure that the control signals maintain the cable tensions inside a predefined feasible workspace during motion (Lafourcade, 2004;[START_REF] Ming | Study on multiple degree-of-freedom positioning mechanism using wire-concept, design and control (part 1)[END_REF]. This is a common solution for redundant manipulators, where the number of cables exceeds the number of Degrees of Freedom (DoF) of the end-effector, and it is the solution considered in this work.

Dealing with position control of CDPR manipulators, most of the proposed methods rely on joint position measurements. According to the coordinate space chosen to solve this control problem, there are two alternatives. In the first one, the controllers are designed in the joint space coordinates. Using the Inverse Position Kinematic Model (IPKM), the reference end-effector pose is converted into reference joint positions which are then controlled by a feedback loop. Some related works are the joint space PD controller proposed by [START_REF] Kawamura | High-speed manipulation by using parallel wire-driven robots[END_REF] applied to the SEGESTA robot and later to the KNTU robot by [START_REF] Gholami | On the control of the KNTU CDRPM: A cable driven parallel manipulator[END_REF], and the joint space PID controller for the redundant suspended ReelAx8 prototype presented by Lamaury et al. (2012a). In the second one, the controllers are designed in the task space coordinates. Assuming that a Direct Position Kinematic Model (DPKM) is available, the end-effector pose is calculated from the joint position measurements and a feedback control allows to track a reference pose. [START_REF] Gholami | On the control of the KNTU CDRPM: A cable driven parallel manipulator[END_REF] evaluated such a task space PD controller and compared it to the previous approach. However, for parallel manipulators, the DPKM is difficult to obtain (see for instance [START_REF] Carricato | Direct geometrico-static problem of under-constrained cable-driven parallel robots with three cables[END_REF] for an intensive study on the matter). In the previously mentioned control schemes, the modelling errors and the deformations of the cables result directly in errors on the end-effector pose. One solution for improving the accuracy is then to use some exteroceptive sensors in order to obtain a direct measurement of the end-effector pose. Some preliminary works using cameras have been proposed by Dallej et al. for controlling the redundant suspended ReelAx8 robot (2011) or the large-dimension CoGiRo robot (2012) .

CDPRs are non-linear coupled MIMO systems with constrained states. For such systems, stability issues arise with linear LTI or model inversion based controller as soon as the model is not well known or when the parameters evolve during operation. Multi-variable control methodologies that have been developed since the 90s now allow to manage the trade-off between performances and robustness and are then good candidates for CDPRs. Some preliminary results were proposed by [START_REF] Laroche | A preliminary study for H ∞ control of parallel cabledriven manipulators[END_REF] for the design of a LTI robust controller that manages both position and tension, including simulation results of a 3-DoF cable-robot with 4 cables.

The purpose of this paper is twofold. First, a two-step methodology is proposed to accurately estimate the kinematic parameters prior and independently from the dynamic parameters of a physics-based CDPRs model. It requires both the measurement of the motors angular positions (measured by optical encoders) and the end-effector pose (provided by a system of cameras). After a presentation of the available identification methods, experimental results are provided, showing that the identified model fits the system behaviour with good accuracy, and can then be used for control. Second, an original control scheme is developed in two parts: a vision-based position control scheme, for which an H ∞ methodology is proposed, in addition to a cable tension distribution management algorithm. The position control is achieved in the operational space, making use of an end-effector pose directly measured by a vision-tracking system. A four-blocks H ∞ design scheme with adjusted weighting filters ensures good trajectory tracking and disturbance rejection for CDPR systems which are Non-Linear Time Invariant (NLTI) coupled MIMO systems with constrained states. In conjunction with this position control loop, a tension management algorithm aims at ensuring that the control signals maintain the cables under feasible tensions. The paper makes an extensive review of the available methods and presents an algorithm inspired from one of them, to account for the motor dynamics. Results from simulations and experiments are then reported using the redundant 6-DoF CDPR INCA 6D prototype with 8 cables manufactured and sold by Haption, equipped with a motion-tracking system Bonita developed by Vicon.

The paper is organised as follows: Section 2 describes the setup composed of the INCA robot and the Bonita motion-capture system. In Section 3, a physics-based model of the 6-DoF CDPRs is developed. In Section 4, the identification methodology is described and implemented. Section 5 is dedicated to the control strategies and includes both simulation and experimental results.

System description

INCA robot

The INCA robot developed by Haption is a haptic device with force feedback that is used in this work as a manipulator to address the control issues raised by CDPRs.

On its 6D version, the INCA 6D has a cubic configuration of 3 m by side, and uses eight driving cables to move the end-effector and 8 balancing cables to ensure pretension in the driving cables when the motors are non-powered (Fig. 1(b)). Each actuator located at one of the eight vertices of the workspace (Fig. 1(a)) is composed of a DC motor with a current driver, coupled to both the driving and balancing winches (Fig. 1(b)) to store the cables, and is associated with a pulley to guide the cables, thus avoiding shearing. The contact points between the driving cables and the end-effector are called attachment points, and the contact points between the driving cables and the pulleys are called output points.

A measurement of the motor positions and currents are respectively achieved by incremental optical encoders and current sensors.

Bonita motion-capture system

The Bonita motion-capture system used to measure the pose of the INCA end-effector is composed of six infrared (IR) cameras (Fig. 1(a)) and a tracker software running on a Windows PC, both from Vicon company. Each camera has its own emitting source and delivers a grayscale image with VGA resolution up to a 250 Hz frame-rate. Assuming that this stereo system has been previously calibrated, the pose of the INCA end-effector fitted with five retro-reflective fixed markers can be tracked by the software.

The temporal and spatial performances of the pose reconstruction are critical for the robot control and have been evaluated as follows:

• the delay between the start of the image acquisition and the availability of the pose measurement has been evaluated at 10.7 ± 0.7 ms with a 200 Hz camera frame rate. This latency of roughly twice the acquisition period is the sum of one period of image acquisition and one period for the pose reconstruction. • the accuracy of the pose is of 1.7 ± 0.4 mm. It has been estimated by the RMS error in the IR camera images between the current position of the visible markers and their expected positions given by the reconstructed pose of the end-effector.

Real-time control architecture

To allow the controller to operate at a higher frequency than the frequency of the Bonita motioncapture system, the control architecture summed up in Fig. 2 is designed as the following two asynchronous real-time tasks:

• the main periodic task is run by a PC under RT Linux at 1 kHz. It consists in reading the motor currents and the motor positions, checking the consistency of the robot state and finally computing the reference for the motor currents. • a secondary task is run by a PC under Windows asynchronously at 200 Hz. It waits for the availability of the end-effector pose provided by the vision-tracking system.

Modelling of the 6-DoF CDPRs

The physics-based model considered in this work is derived from the general model for the m-DoF CDPR manipulators with n cables, which is augmented with the pretension system (balancing cables, winches and springs) of the INCA prototype. It is assumed that the cables are of negligible mass (straight, no sagging) and of infinite stiffness (inextensible).

Kinematics modelling

Position kinematics

Denoting R o the base reference frame and R e the end-effector frame centred at its center of mass, the pose of the end-effector can be represented by the vector X e = [P T e Φ T e ] T , with P e the position of the origin and Φ e a representation of the orientation of R e with respect to the base reference frame R o (Fig. 3(a)). Considering Φ e = [φ r φ p φ y ] T , the rotation matrix o R e from R o to R e has been chosen as a composition of three successive rotations Roll-Pitch-Yaw [START_REF] Khalil | Modélisation, identification et commande des robots[END_REF] of angles (φ r , φ p , φ y ) respectively around the principal axes

(X o , Y o , Z o ) of R o : o R e (Φ e ) =   c y c p c y s p s r -s y c r c y s p c r + s y s r s y c p s y s p s r + c y c r s y s p c r -c y s r -s p c p s r c p c r   (1)
with s k = sin φ k and c k = cos φ k , k standing for r, p or y.

The attachment point E i with coordinates e r i fixed in R e is located at the position P Ei = P e + o R e (Φ e ) e r i in R o . With P Oi denoting the position of an output point O i fixed in R o , the length of cable #i writes L i (X e ) = l i (X e ) 2 where l i (X e ) = P e + o R e (Φ e ) e r i -P Oi . We then stack the lengths of the cables in L = [L 1 . . . L n ] T .

By convention, the angular position is considered as increasing during winding. With a null reference of the motor positions θ o = O 8×1 for the initial pose of the end-effector

X eo = [P e T o Φ e T o ]
T , the Inverse Position Kinematic Model (IPKM) is then given by:

θ(X e , α K ) = -R -1 pm (L(X e , α K ) -L(X eo , α K )) (2)
where α K is the vector of the parameters involved in the kinematic model and R pm = diag(r pm1 1 , ..., r pm1 n ) contains the radii of the driving winches.

Velocity kinematics

Differentiating Eqn.

(2) with respect to time yields the Inverse Velocity Kinematic Model (IVKM):

θ = -R -1 pm J(X e ) V e (3)
where the i th row of the inverse kinematics Jacobian matrix J is the same as for a rigid parallel manipulator [START_REF] Merlet | Les robots parallèles[END_REF] and is given by J i (X e ) = [u i T (X e ) (r i (Φ e ) × u i (X e )) T ], in which u i (X e ) = l i (X e ) / L i (X e ) is the unit direction vector of the i th driving cable (Fig. 3(a)). The endeffector velocity V e = [v T e w T e ] T (including the linear v e and angular w e velocities) can then be converted into the time derivative of end-effector pose Ẋe by:

V e = A v (Φ e ) Ẋe (4)
with A v (Φ e ) = diag(I 3×3 , J rpy (Φ e )), in which the matrix J rpy maps the angular velocity w e to the time derivative of the chosen orientation representation Φe such as w e = J rpy (Φ e ) Φe . For a Roll-Pitch-Yaw representation of the rotation, J rpy has the following form:

J rpy (Φ e ) =   c y c p -s y 0 s y c p c y 0 -s p 0 1   (5)
Denoting J(X e ) = R -1 pm J(X e ) A v (Φ e ) and recalling the dependance with respect to the kinematic parameter vector α K , the IVKM is finally rewritten as:

θ = -J(X e , α K ) Ẋe (6)
3.2 Dynamics modelling

End-effector dynamics

The Newton-Euler equations [START_REF] Khalil | Modélisation, identification et commande des robots[END_REF] applied to the end-effector rigid body of mass M e and moment of inertia tensor I e , written at its center of mass O e gives:

F e = W (X e ) T + G (8)
where W (X e ) = -J T (X e ) is the wrench matrix and G = [0 0 -M e g 0 0 0] T with g the modulus of the gravity acceleration vector.

Actuators dynamics

With the diagonal matrix of the inertia moments J eq = j eq I n , viscous F veq = f v I n and Coulomb F ceq = f c I n friction coefficients, rotational stiffness K eq of the balancing spring, and neglecting the tension loss due to the output roll (T a = T e = T in Fig. 3(b)), the equation of the actuator dynamics writes:

J eq θ + F veq θ + F ceq sign( θ) + K eq θ = K em I m -R pm (T a -T ao ) (9)
with: pm L(X eo ) the pretension of the driving cables for an initial pose X eo of the end-effector, due to a change of length of the balancing springs by R pe R -1 pm L(X eo ).

• K eq = K sp R pe 2 , R pe =

System dynamics

Combining ( 7), ( 8), ( 9) with the use of Eqn. ( 2), (4), Eqn. ( 6) and their time derivative, the Direct Dynamic Model (DDM) of the system written in the operational space is given by:

M x (X e ) Ẍe + C x (X e , Ẋe ) + K x (X e ) + G x = A T v (Φ e ) F ev (10)
under the n constraints expressing that the cables are kept inside a feasible tensions workspace [T min T max ]:

T min ≤ T (I m , θ, θ, θ) ≤ T max (11)
The matrices involved in the model write:

                     M x (X e ) = A T v (Φ e ) A e (X e ) A v (Φ e ) + JT (X e ) J eq J(X e ) K x (X e ) = JT (X e ) K eq R -1 pm L(X e ) G x = -A T v (Φ e ) G = -G C x (X e , Ẋe ) = A T v (Φ e ) A e (X e ) Ȧv (Φ e , Φe ) + JT (X e ) F veq J(X e ) + JT (X e ) J eq J(X e , Ẋe ) Ẋe +A T v (Φ e ) B e (Φ e , Φe ) + JT (X e ) f ceq sign( J(X e ) Ẋe ) (12
) and the wrench F ev resulting from the motor-current vector I m is:

F ev = W I (X e ) I m (13)
where

W I (X e ) = W (X e ) R -1 pm K em .
This model relies on the kinematic parameters stacked in α K and also on the dynamic parameters (inertias, torque constant, friction ratios...) stacked in α D .

Nominal linear model

As can be seen in the system dynamics given by the ODEs system (10), CDPRs are NLTI coupled MIMO systems. So the linearisation of ( 10) is required for the use of the H ∞ methodology presented in Section 5. It leads to the linear coupled ODEs system:

M o δ Ẍe + C o δ Ẋe + K o δX e = A T v (Φ eo ) δF ev (14)
where the linearised dynamic matrices could be evaluated analytically by differentiating F ev from (10) at the operating point op = (X eo , Ẋeo , Ẍeo ):

M o = ∂(A T v F ev ) ∂ Ẍe op C o = ∂(A T v F ev ) ∂ Ẋe op K o = ∂(A T v F ev ) ∂X e op ( 15
)
and are expressed at the operating or equilibrium point op = (X eo , 0, 0) by:

             M o = A T v (Φ e0 )A e (X e0 )A v (Φ e0 ) + JT (X eo ) J eq J(X eo ) C o = JT (X eo ) F veq J(X eo ) K o = JT (X eo ) K eq J(X eo ) +   ∂ JT ∂p x Xe o K eq R -1 pm L(X eo ) . . . ∂ JT ∂φ y Xe o K eq R -1 pm L(X eo )   (16)
Coulomb frictions are naturally removed in the linear model. They can be considered as disturbances acting on the wrench that will be rejected by the controller. Hence, the state-space representation of the linearised system based on state R 12 , control u = δF ev ∈ R 6 and measurement y = δX e ∈ R 6 vectors writes:

x = [δX T e δ ẊT e ] T ∈
ẋ = A o x + B o u y = C o x + D o u (17)
with the state-space matrices expressed as

A o = O 6×6 I 6×6 -M -1 o K o -M -1 o C o , B o = O 6×6 M -1 o A T v (Φ eo ) , C o = I 6×6 O 6×6 and D o = O 6×6 .

Model identification

A number of approaches are available in the literature for identification of robot models. [START_REF] Yu | Kinematic parameter identification of parallel robots for semi-physical simulation platform of space docking mechanism[END_REF] have considered the estimation of the kinematic parameters of a rigid parallel manipulator.

For estimation of dynamic parameters, the reader can refer to [START_REF] Gautier | Identification en boucle fermée par modèle inverse des paramètres physiques de systèmes mécatroniques[END_REF]; [START_REF] Gautier | A new closed-loop output error method for parameter identification of robot dynamics[END_REF] for the electromechanics systems in general and for the serial robots particularly and to [START_REF] Poignet | Robust estimation of parallel robot dynamic parameters with interval analysis[END_REF] for the case of parallel robots.

Identification methodology

Following [START_REF] Renaud | Kinematic and dynamic identification of parallel mechanisms[END_REF], we propose to estimate the parameters sequentially in two steps, thus reducing the number of parameters to be considered at one time:

•

Step 1: The kinematic parameters are first estimated using the IPKM (see [START_REF] Vischer | Kinematic calibration of the parallel Delta robot[END_REF] for an example of determination of the kinematic parameters of a parallel robot). • Step 2: The dynamic parameters are then estimated from the DDM, using the identified kinematic parameters.

Step 1: Kinematic parameters identification

Considering a set of N s measurements, the estimate α K of the vector α K of the kinematic parameters can be determined by minimising the following criterion E(α K ) on the motor positions (see Fig. 4.a):

E(α K ) = n i=1 Ns k=1 (θ ik -θ ik (α K )) 2 (18)
with :

θ ik = θ i (t = k T s ) motor position measurements. θ ik (α K ) = θ i (t = k T s , α K )
motor position estimations using the IPKM.

Minimising E(α K ) is a non-linear least squares optimisation problem that can be solved iteratively using numerical algorithms such as: Gradient-Descendant, Gauss-Newton or Levenberg-Marquardt. The identification of the dynamic parameters of a robot is a general complex problem, due to the non-linear nature of both the input-output behaviour and the parameter dependence of these MIMO systems. Identification methodologies for robots can be classified according to the minimization criterion [START_REF] Gautier | A new closed-loop output error method for parameter identification of robot dynamics[END_REF]):

• Inverse Dynamic Identification Model (IDIM): the Inverse Dynamic Model (IDM) is written in a linear form according to the parameters to be estimated, thus allowing to compute their estimate in one shot by minimizing the quadratic error between the measured input and its estimation computed on the same output trajectory [START_REF] Gautier | Identification en boucle fermée par modèle inverse des paramètres physiques de systèmes mécatroniques[END_REF]. This method requires the estimation of the speed and the acceleration, through bandpass filtering of the joint position at high sampling rate. • Output Error (OE): the quadratic error between measured outputs and their estimate is minimised, where the output estimate is computed with the DDM, the system being fed with the same input signal. This technique can be implemented in an open-loop (OLOE) or closed-loop (CLOE) situations [START_REF] Janot | Comparison between the CLOE method and the DIDIM method for robots identification[END_REF]. This very popular method often yields to a non-convex problem with possible local minima and is then sensitive to initialisation.

Despite of its drawbacks, the OE method has been chosen as it is more realistic (the errors are mainly on the measurement signals) and thus leads to a more accurate model for control synthesis. Beneficing from the stabilization effect of the balancing springs, the experiments were done in open loop. For a robot without balancing springs, a closed-loop framework would better suit to ensure that the platform remains within a limited workspace. An excitation is first designed, leading to a set of input and output data collected from N s samples at a sampling rate f e . Then, the estimate α D of the vector α D of the dynamic parameters can be determined by minimising the following dynamic identification error criterion E( α K , α D ), defined in the operational space in term of end-effector posture (see Fig. 4.b with c = 0):

E( α K , α D ) = 6 i=1 Ns k=1 (X eik -X eik ( α K , α D )) 2 (20)
or also in the joint space in term of motors positioning (see Fig. 4.b with c = 1):

E( α K , α D ) = n i=1 Ns k=1 (θ ik -θ ik ( α K , α D )) 2 (21)
with :

X eik = X ei (t = k T s ) end-effector pose measurement, X eik ( α K , α D ) = X ei (t = k T s , α K , α D ) end-effector pose estimation using the DDM, θ ik = θ i (t = k T s ) motor position measurements, θ ik ( α K , α D ) = θ i (t = k T s , α K , α D )
motor position estimations using the DDM and IPKM.

As for the kinematic parameters, minimising E( α K , α D ) is a non-linear least-square optimisation problem that can be solved iteratively using the same algorithms mentioned for the kinematic parameters estimation. However, the evaluation of the objective function requires the simulation of the DDM, which is very demanding in computation time.

Implementation and results

The proposed identification methodology has been implemented on the INCA 6D CDPR. For the i th output of a MIMO system of input vector u (dim(u) = N u ) and output vector y (dim(y) = N y ) among a set of N s samples, the fit is defined by:

F IT i (%) = 1 - Ny i=1 Ns k=1 (y ik -y ik ) 2 Ny i=1 Ns k=1 (y ik -y i ) 2 × 100 (22)
where y is the estimation of the output vector y and y i the average of the i th output y i .

Initial values of the parameters

The reference frame R o is located at the center of the workspace (Fig. 5(a)), and the initial position for the frame R e attached to the end-effector is chosen equal to R o (Fig. 5(b)).

Each motor is controlled in current and the current loop, i.e. the transfer between the current reference I * mi and the effective motor current I mi , is modelled by a continuous-time first-order transfer function

G i (s) = 1 / (1 + τ i s) with a time constant τ i .
The initial values of the kinematic (Fig. 5) and dynamic parameters of the INCA 6D robot, known by design or calculated, are given in the Table 1. The different actuators are assumed to have the same parameters.

Coulomb friction parameter f ceq has been evaluated by powering each motor separately with a ramp signal for the reference current slowly varying from I min = 0 A to I max = 3 A with a slope of 1 A/s, whereas the other motors are controlled to 0 A. The Coulomb friction coefficient is then calculated when detecting the minimum value of the current that produces a motion of the end-effector. The worse case over the 8 winders has been considered for f ceq given in Tab. 1.

Estimation of the kinematic parameters

The test consists in moving the effector in eight different directions by applying successively a current reference of 1 A on one motor and 0 A on the other ones. The positions of the motors and the pose of the end-effector are acquired at steady-state.

The identification scheme presented in subsection 4.1.1 and summarised in Fig. 4(a) has been implemented using the Levenberg-Marquardt optimiser. The estimates of the kinematic parameters are given in Table 1. Notice that parameters a and r pe are not involved in the kinematics modelling and then were not involved in the estimation procedure. The comparison between the positions of the motors estimations issued from the IPKM and the experimental identification data is shown in Fig. 6. The error on the motor position is centered (less than 0.4 • , i.e. about 1 mm on the cable length) and the standard deviation is 45 • , corresponding to 13.7 mm on the cable length. The standard deviation of the errors on the estimates provided in Table 1 are computed from the diagonal terms of Σ 2 αK . The corresponding numerical values of the fit (22) on the motors positions are provided in the first column of the Table 2. One can notice that the model reproduces accurately the measured positions.

Estimation of the dynamic parameters

4.2.3.1 Experiment. The choice of the excitation trajectories is an important point for the identification procedure. Indeed, the model needs not only to be structurally identifiable, but the excitation trajectory must also be rich enough. The excitation trajectories must contain both slow (for friction and stiffness) and fast (for inertia) dynamics.

For the current case of a CDPR, the trajectories should be also large enough to excite the nonlinear behavior, the cable tension must be sufficiently high to avoid slackness, and sufficiently slow not to excite the vibrations of the cables.

Motor current signals have been chosen as decoupled Pseudo-Random Binary Sequence (PRBS) that are usually used in system identification as they can both excite the different inputs of a multivariable system in a decoupled fashion and are able to excite over a large bandwidth. For the The validation trajectories Val(e) and Val(f) with much higher frequency f max allow to evaluate the effects of the cable vibrations that were observed during these experiments and that are not accounted for in the model.

Estimated dynamic parameters.

Applying the dynamic identification scheme presented previously in Fig. 4(b) (considering the joint positions as outputs, i.e. with c = 1), the dynamic parameters are estimated iteratively using the Levenberg-Marquardt optimiser for solving the nonlinear least-square optimisation problem of Eqn. ( 21).

The estimates of the dynamic parameters are given in Table 1. Parameters k em and τ i have been maintained at their initial values, that have been confirmed by previous experiments. The nondiagonal components of the moment of inertia are set to: e I exy = e I eyz = e I exz = 0 kg•m 2 , and the diagonal components e I eyy and e I ezz are set as equal ( e I ezz = e I eyy ), because of the symmetry of the end-effector.

An evaluation of the standard deviation of the estimates are given in the last column of table 1. This evaluation has been done by considering the sum of the parameter variances when considering all the errors on the kinematic parameters. More precisely, the dynamic parameter identification has been run seven more times, considering an error of 1 % on each kinematic parameter, thus resulting in a sensitivity matrix expressing the effect of an error on each kinematic parameter on the dynamic parameter estimate. The variances of all the kinematic parameters, weighted by their sensitivities, are then summed up to produce the variances of the dynamic parameter estimates and finally their relative standard deviations given in table 1. One can notice that the highest deviations are for the inertias and are of 13 %. In Fig. 7 are provided for the first three motors: the reference current (on the top), the estimated and measured motor angular positions (2nd row) and a zoom on the motor position error (lower row).

) 0 0* / e I eyz (kg•m 2 ) 0 0* / e I exz (kg•m 2 ) 0 0* / j eq (kg•m 2 ) 2.91 • 10 -5 3.88 • 10 -5 0.4 f v eq ((N•m)/(rad/s)) 3.1 • 10 -3 6.68 • 10 -5 3.8 f ceq ((N•m)/(rad/s)) 1.8 • 10 -3 4.4 • 10 -3 3.3 k sp (N/m) 16
Table 2 provides the numerical values of the fit ( 22) on the motor positions. The fit obtained with the identification trajectory (see the column labeled "Dynamic ident") is complemented with the fit obtained for six different validation trajectories with variable amplitudes and frequencies. One can see that the model reproduces accurately the measurement, even on the validation data that were not used for identification. These results validate the assumptions used for deriving the model such as the negligible mass of the cables and their infinite stiffness.

Frequency-domain behaviour

The linearisation has been performed for the end-effector velocity and acceleration equal to zero. The nominal model corresponds to a position of the platform at the center of the workspace, with the nominal values of the kinematic and dynamic parameters previously estimated.

The frequency behaviour of the nominal linear model G n (s) is reported in Fig. 8. The singular values for both position and orientation models behave as a second order under-damped LTI system with a resonance at the frequencies close to 5 rad/s for position and of 22 rad/s for orientation.

I min (A) 0 0 0 0.5 0.5 0 0.5 Imax (A) 1 1 1 1.5 1.5 1 1.5 fmax (Hz) 1 1 2 1 2

Position control and tension distribution

The system can be considered as redundant as the number of actuators exceeds the number of DoF of the platform (refer to [START_REF] Kanoun | Kinematic Control of Redundant Manipulators: Generalizing the Task-Priority Framework to Inequality Task[END_REF] for a more general reflexion on redundant systems). The extra DoF are used to handle the cable tension. The considered control scheme presented in Fig. 9 is composed of two parts detailed in this section:

• Position control that aims to control the pose of the end-effector at a reference value. It is typically a dynamic position-based visual control scheme (also known as Position-Based Visual Servoing (PBVS) or 3D Visual Servoing (3DVS)). • Tension distribution that aims to maintain the tension of the cables inside a feasible workspace.

H ∞ synthesis

The H ∞ synthesis provides a framework for tuning dynamic output-feedback MIMO controllers in the frequency domain, allowing to handle easily the robustness versus performance trade-off [START_REF] Duc | Commande H ∞ et µ-analyse[END_REF].

In the standard synthesis scheme described in Fig. 10, the controller K(s) to be designed closes the loop of the extended plant G e (s), thus modifying the performance channel from v to z. Weighting filters located at the input (W i (s)) and output (W o (s)) complete the so-called augmented scheme. A number of problems can be recast into the following H ∞ synthesis problem that consists in computing K(s) that:

• stabilises the weighted closed-loop system in Fig. 10; • minimises γ := min K(s) T zṽ (s) ∞ , i.e. the H ∞ norm1 of the weighted closed-loop system T zṽ (s) = W o (s) T zv (s) W i (s).

Assuming that the weighting filters W i (s) and W o (s) are diagonal, and denoting W i (s) = diag(W i1 (s), . . . , W inv (s)) and W o (s) = diag(W o1 (s), . . . , W onz (s)), where n v and n z being the dimensions of v and z respectively, then any SISO transfer T zkvl (s) between input v k and output z l satisfies:

|T zkvl (jω)| ≤ γ |W ol (jω)| • |W ik (jω)| (∀ ω ∈ R + ) (23)
Therefore, 1 Wol(s)Wik(s) is a template for the performance channel T zkvl (s) and the template is satisfied as soon as γ ≤ 1.

The resolution of the standard problem of the H ∞ synthesis can be done by solving the resulted Riccati [START_REF] Dolye | State-space solutions to standard H 2 and H ∞ control problems[END_REF] or LMI equations [START_REF] Gahinet | A linear matrix inequality approach to H ∞ control[END_REF] resulting in full-order controllers, i.e. controllers with the same order as the augmented plant. Solvers are also available for the structured H ∞ synthesis: let us mention the HIFOO package [START_REF] Burke | HIFOO -a Matlab package for fixed-order controller design and H ∞ optimisation[END_REF] and hinfstruct available in the the Robust Control Toolbox [START_REF] Apkarian | Nonsmooth H ∞ synthesis[END_REF].

To summarise, the method consists in the following steps:

• select the performance channels,

• design the weighting functions according to the requirements,

• compute the augmented plant and design the controller,

• iterate the two previous steps until γ is close to one.

Controller design

Fig. 11 depicts the synthesis scheme. The control scheme is deduced by ignoring the weighting filters W 1 (s), W 2 (s) and W 3 (s). It allows to control the measurement vector y = X e to the reference r = X * e , while rejecting the disturbance d that is added on the control vector u and considering the measurement noise b m . As b m and r have the same effects on the signals of interest e and u, only r is considered for building the extended plant G e (s), leading to a four-block design scheme equivalent to the scheme of Fig. 10 T and z = [e T u T ] T , and the weighting filters

W i (s) = diag(I 6×6 , W 3 (s)) and W o (s) = diag(W 1 (s), W 2 (s)).
The weighing filters are chosen as following:

• W 1 (s) allows to tune the sensitivity function S(s) = T er (s), i.e. the transfer from r to e. It is chosen as constant, to impose the modulus margin

∆ M . As ∆ M = 1/||S(s)|| ∞ , we select W 1 (s) = ∆ * M
where ∆ * M denotes the required modulus margin. • W 2 (s) allows to tune T ur (s), the transfer between r to u, which is equal to T ubm (s). So W 2 (s)

has the charge of reducing the effect of the measurement noise on the control signal. Thus, W 2 (s) must amplify in the high frequencies and 1/W 2 (s) is chosen as a low-pass filter. • W 3 (s) is used to penalise the disturbance in the low frequency thus enhancing the disturbance rejection properties.

The desired performances are: a bandwidth of 5 rad/s, a modulus margin of 6 dB and an accuracy of 10 -3 due to d for both position and orientation, with a reasonable amplitude of the control inputs IJC2015_v7_withNames Figure 11. Design scheme with the weighting filters W 1 (s), W 2 (s) and W 3 (s) for the H∞ synthesis for both forces and moments. The weighting filters W 1 (s), W 2 (s) and W 3 (s) are chosen under the shape W k (s) = diag(w k1 (s) I 3×3 , w k2 (s) I 3×3 ), where the SISO weighting filters w k1 (s) correspond to the positions or forces and w k2 (s) corresponds to the orientations or moments. The weighting filters were tuned to achieve the requested closed-loop performances, resulting in: Notice that the more simple two-bloc scheme, based on weighting filters W 1 (s) and W 2 (s) only, was not sufficient to reject the disturbances induced by the non-linear behaviour of the system. The frequency response of the controller (singular values) is given in Fig. 12. For both position and orientation, the controller exhibits a high gain in low frequencies that favors the disturbance rejection, then the gain increases in the medium frequencies, bringing phase to allow a high bandwidth and finally decreases in the high frequencies, thus reducing the effects of the measurement noise. Figure 13 shows the singular values of the closed-loop transfer functions of interest involved in the H ∞ synthesis and their corresponding templates. The resolution has been done with the Riccati method. The obtained value of γ = 0.92 being slightly smaller than 1, the templates are fully satisfied and the tracking trajectories performances slightly exceed those intended with a bandwidth of 6 rad/s for position and of 2 rad/s for orientation, a modulus margin of 0.57 for position and of 0.61 for orientation and a static error less than 1 % on position and orientation (see Fig. 13(a)).

          

Tension distribution

Problem statement and state of the art

The wrench F ev computed by the H ∞ controller must be converted into motor current vector I * m . This requires to solve the system of linear equations (13) under the inequality constraints (11), in order to ensure that the cable tension vector T remains inside the interval [T min T max ]. Due to the actuation redundancy, the system of equations ( 13) is under-determined and then has an infinity of solutions (assuming that W I has full rank r = nm). This can be exploited to minimise the distance towards some objective value I obj corresponding to an objective tension vector T obj . This is done by solving the quadratic optimisation problem of the objective function E(I * m , λ) given by:

E(I * m , λ) = 1 2 (I * m -I obj ) T (I * m -I obj ) + (F ev -W I I * m ) T λ (25)
where λ ∈ R m is the Lagrange multiplier associated to the equality constraints F ev -W I I * m = 0. The solution writes:

I * m = I * mp + I * mn ( 26
)
where:

• I * mp , the particular solution that minimises the norm of I * m -I obj , acts on the end-effector pose:

I * mp = W + I (X e ) F ev (27)
with

W + I = W T I (W I W T I )
-1 the Moore-Penrose pseudo-inverse of the wrench matrix W I .

• I * mn , the homogeneous solution of W I null-space, acts on the tensions, without effect on the end-effector pose:

I * mn = [I n×n -W + I (X e ) W I (X e )] I obj (28)
A number of approaches are available in the literature the determination of the objective tension vector T obj that maintains the cable tension vector T inside a feasible tensions workspace. The available works can be classified in two categories. The first approach opted for iterative algorithms, so that efficient constrained optimisation methods can be used such as Linear Programming Methods (LPM) [START_REF] Gosselin | On the determination of the force distribution in overconstrained cable-driven parallel mechanisms[END_REF][START_REF] Ming | Study on multiple degree-of-freedom positioning mechanism using wire-concept, design and control (part 1)[END_REF][START_REF] Oh | Cable suspended planar robots with redundant cables: controllers with positive tensions[END_REF], but the cable tension continuity is not guaranteed. Other optimisation methods are also used such as Non-Linear Programming Methods (NLPM) in the particular case of Quadratic Programming Methods (QPM) [START_REF] Oh | Cable suspended planar robots with redundant cables: controllers with positive tensions[END_REF][START_REF] Vafaei | Integrated controller for an over constrained cable driven parallel manipulator: KNTU CDRPM[END_REF], and the general NLPM with the gradient descent method to solve the problem in a quadratic formulation [START_REF] Gholami | On the control of the KNTU CDRPM: A cable driven parallel manipulator[END_REF]. These quadratic methods guarantee the tensions continuity but have a non-predictable runtime, which is a drawback for real-time implementation. The second approach relies on non-iterative algorithms to handle the real-time control constraints. For instance, [START_REF] Mikelsons | A real-time capable force algorithm for redundant tendon-based parallel manipulators[END_REF] proved that the Center of Gravity (CoG) of the feasible tension distribution workspace (the set of solutions of equations (13) satisfying the tension constraints ( 11)) is a solution that ensures the tensions continuity. Recently, [START_REF] Gouttefarde | A Versatile Tension Distribution Algorithm for n -DOF Parallel Robots Driven by n+2 Cables[END_REF] have proposed contributions to improve the CoG method to the case of two degrees of redundancy.

Considered algorithm

The proposed approach is inspired from the algorithm proposed by Lafourcade ( 2004) appropriate to satisfy real-time constraints (currently less than 1 ms). It has been improved to account for the motor dynamics (9) by considering a variable I obj and accounting for the current limits of the motors. It consists in:

(1) selecting T obj inside the feasible tension workspace and calculating I obj from Eqn. (9) by considering T = T obj and solving in I obj = I m ; (2) computing the current reference I * m that satisfies constraints (13) while minimizing the meansquare error on I * m -I obj without considering the tension inequalities constraints (11) as given by ( 27);

(3) if some inequality constraints on the tensions are violated, then q concerned inequalities selected among them are transformed into equality constraints on the currents (not more than r tension inequality constraints can be saturated simultaneously) and included into the optimisation problem. All combinations of 1 to r violated tension inequality constraints are considered until one solution is met that satisfies all the non-saturated tension inequalities constraints. The process is repeated until a feasible solution is met. If no solution is found, the vector I obj can be scaled by a scalar factor without impact on the trajectory. If the problem has no solution, the reference trajectory must be decelerated.

More precisely, for a combination of q tension constraints, the current vector I * m is obtained by solving the new quadratic optimisation problem of the objective function E sat (I * m , λ, µ) given by:

E sat (I * m , λ, µ) = 1 2 (I * m -I obj ) T (I * m -I obj ) + (F ev -W I I * m ) T λ + (S T I * m -I sat ) T µ (29)
where:

• the selection matrix S = [s 1 . . . s q ] ∈ R n×q concatenates the vectors s k of the canonical basis of R n to select the combinations of the violated tension constraints to be saturated, • the vector µ ∈ R q is the Lagrange multiplier associated to the current equality constraints S T I * m -I sat = 0, such as the vector I sat ∈ R q is the current vector that corresponds to the saturated tension vector T sat and is calculated from Eqn. ( 9).

The solution of the optimisation problem (29) writes:

I * m = I * mp + I * mt (30)
where:

I * mt = I * mn + W + I (X e ) W sat I (X e ) -S S T W + I (X e ) W sat I (X e ) -I q×q -1 sat (31)
in which the resulting saturated wrench matrix is W sat I (X e ) = W I (X e ) S, and the vector of the excessive motor currents sat = I sat -S T (I * mp + I * mn ) is an image of the excess in cable tensions evaluated from Eqn. (9).

Experiments

Evaluation tasks

A set of trajectories have been chosen in order to evaluate the performance of the controlled system. Starting from the center of the workspace identified by a null position and orientation, the trajectory is composed of steps: the end-effector is requested to follow the edges of a square centred into the workspace, belonging to the plane (X o , Y o ) and of width 0.2 m.

The feasible tension workspace is defined by the boundaries T min = 1.48 N and T max = 18.5 N (the maximum tension supported by the cables being of 247 N), which have been calculated based on the static model of the actuators, considering the limits of the motor currents I min = 0 A and I max = 3 N, and the limits of the unwinding cable lengths L min = 0 m and L max = 4.82 m. The objective tension vector is chosen as T obj = 10 N.

Results

The experimental results obtained with the square trajectory are reported in Fig. 14. One can notice that the position is properly controlled in a decoupled fashion (Fig. Simulation results provided in Fig. 14 include Coulomb frictions that were found as a major source of misfit between the simulated and the measured signals. As shown in Fig. 15, the prediction provided by the model is much more accurate as soon as Coulomb friction is included.

) X -A x is Y -A x is Z-A x is (

Conclusion

This paper is a methodological contribution for the identification and control of 6-DoF CDPRs. First, a two-step identification methodology has been proposed to estimate sequentially the combined kinematic and dynamic parameters of 6-DoF CDPRs. The method has been successfully implemented on the INCA 6D robot, allowing to improve the input-output behaviour of the model, as shown by the good fit on the motor positions. Moreover, the estimated values of both kinematic and dynamic parameters remain close from their initial guessed values, except for the viscous friction term. The validity of the identified model remains satisfying even when the cables are subject to non-negligible vibrations, as could be observed on the validations with faster trajectories. Then, an H ∞ methodology has been proposed for the vision-based position control of 6-DoF CDPRs, which shown that the H ∞ position controller combined with the tension distribution algorithm allows good bandwidth, accuracy and disturbance rejection, while maintaining the cables under feasible tensions. Notice that strong disturbance rejection properties had to be conferred to the controller in order to dominate the non-linear behaviour of the system. One strong limitation of the H ∞ methodology is that it leads to a LTI controller that is restricted to a neighborhood of the nominal position. However, the robustness of the designed controller was sufficient for an evolution in a reasonably large domain. For a larger workspace, gain-scheduling techniques would be fruitful to adapt the controller behaviour.

  Figure 1. The INCA robot

F

  e = A e (Φ e ) Ve + B e (Φ e , Φe ) (7) with A e (Φ e ) = M e I 3×3 O 3×3 O 3×3 o I e (Φ e ) and B e ( Φe , Φ e ) = O 3×1 w e (Φ e , Φe ) × ( o I e (Φ e ) w e (Φ e , Φe )) , in which the moment of inertia tensor o I e expressed in R o is calculated from its constant tensor e I e expressed in R e such as: o I e (Φ e ) = o R e (Φ e ) e I e o R T e (Φ e ). The resulting wrench F e = [f T e m T e ] T (including the forces f e and moments m e ) applied on the end-effector center of mass O e includes contributions of the cables tensions T = [T 1 . . . T n ] T applied on the end-effector side (Fig. 3(a)) and of the weight G:

  r pe I n being the diagonal matrix of the balancing winches radius and K sp , the diagonal matrix of the balancing spring stiffness, • K em = k em I n , a diagonal matrix with the torque constants, • I m = [I m1 , . . . , I mn ] T the motor current vector related to the motor torque vector τ m = [τ m1 . . . τ mn ] T by the linear equation τ m = K em I m , • and T ao = K eq R -2

Figure 3 .

 3 Figure 3. Diagram of the 6-DoF CDPR with the dynamic parameters

  Denoting σ 2 θ the variance of the measurement error on the joint positions and A the sensitivity matrix of the model of dimension nN s × n αK of the error on the estimates writes Σ 2 αK = σ 2 θ (A T A) -T (Walter and Pronzato, 1997).

  Figure 4. Parameter estimation schemes

  IJC2015_v7_withNames

Figure 5 .

 5 Figure 5. Kinematic parameters of the INCA 6D robot

Figure 6 .

 6 Figure 6. Estimation of the kinematic parameters. Left-hand side: experimental (plain line) and estimated (dashed line) angular position for the four first motors. Right-hand side: corresponding estimation errors

Figure 7 .

 7 Figure 7. Identification and validation results of the dynamic model, for the first three motors (left to right). On top: reference (plain line) and measured (dashed) values of the motor current; on the second row: predicted value (plain) and measured values (dashed) of the motor positions; lower row: motor prediction error. Before 5 s: data used for identification; after 5 s: validation data PRBS Dynamic Dynamic Dynamic Dynamic Dynamic Dynamic Dynamic parameters Ident Val(a) Val(b) Val(c) Val(d) Val(e) Val(f ) I min (A) 0 0 0 0.5 0.5 0 0.5 Imax (A) 1 1 1 1.5 1.5 1 1.5 fmax (Hz) 1 1 2 1 2 5 5

Figure 8 .

 8 Figure 8. Singular values of the nominal model Gn(s) : position and orientation parts

Figure 10 .

 10 Figure 10. Standard H∞ synthesis scheme with the extended plant Ge(s), the weighting filters W i (s) and Wo(s) and the controller K(s)

  Fig.11depicts the synthesis scheme. The control scheme is deduced by ignoring the weighting filters W 1 (s), W 2 (s) and W 3 (s). It allows to control the measurement vector y = X e to the reference r = X * e , while rejecting the disturbance d that is added on the control vector u and considering the measurement noise b m . As b m and r have the same effects on the signals of interest e and u, only r is considered for building the extended plant G e (s), leading to a four-block design scheme equivalent to the scheme of Fig.10by selecting the signals v = [r T d T ]T and z = [e T u T ] T , and the

Figure 12 .

 12 Figure 12. Singular values of the H∞ position controller K(s)

  Figure 13. Singular values of the closed-loop transfer functions of interest Tz k v l (s) and their corresponding templates γ/W ik (s) W ol (s)

  14.a), while the currents (Fig. 14.b) and tensions (Fig. 14.c) remain inside the prescribed workspace and the tensions are close to T obj = 10 N. The cables tensions have been estimated using the dynamic model of the actuators (9).

  Platform position (x and y coordinates): reference and measurement

  Figure 15. Effect of Coulomb friction on H∞ position controller: the comparison of step responses on models with and without Coulomb friction and on the experimental system validate the sensitivity of the closed-loop system with respect to Coulomb friction

Table 1 .

 1 Initial and estimated values of the kinematic and dynamic parameters of the INCA 6D robot ( * : parameters that have been kept equal to their initial value)

	14.4	0.05

Table 2 .

 2 Identification results: fit on the motors positions

	5	5

  b) Control signals: force fe and moment me

		2.6							13
		2.2 2.4						Mot or 1 Mot or 2 Mot or 3 Mot or 4 Mot or 5 Mot or 6 Mot or 7 Mot or 8	11 12						C able 1 C able 2 C able 3 C able 4 C able 5 C able 6 C able 7 C able 8
		2
	( A)	1.8						( N )	10
	I	1.6						T	9
		1.4
									8
		1.2
									7
		1
		0.8							6
		39	40	41	42	43	44	45	39	40	41	42	43	44	45
					t ( s)							t ( s)
				(c) Measured currents

(d) Estimated tensions Figure 14. Experimental results with the H∞ position controller

The H∞ norm of a linear system is the maximum singular value of its complex gain over frequency. It is also the maximum amplification when considering the L

norm on signals. For a system G(s) of input u(t) and output y(t), ||G(s)||∞ = max ω∈R σ(G(jω)) = max u(t) ||y(t)|| 2 ||u(t)|| 2 .


Nader Mahmoud Elshahat

Oscar Elsayed Ali

Garcia Grasa

Stéphane A Nicolau

Christophe Doignon

Luc Soler

Jacques Marescaux

Jose Maria Martinez Montiel

Nader Mahmoud
email: nader-mahmoud.ali@etu.unistra.fr

Óscar G Grasa

J M M Montiel

On patient see-through Augmented Reality based on visual SLAM

Keywords: Augmented Reality, Visual SLAM, Registration, Operating room, Surface meshes

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

Introduction

Patient pre-operative 3D model (P3DM) is readily available through various imaging modalities such as computed tomography (CT) or magnetic resonance (MRI). The P3DM is typically displayed on a desktop monitor, laptop or tablet-PC. However, the practitioner has to mentally project that information onto the patient. An augmented reality (AR) superimposition of the P3DM onto the patient can provide the practitioner with a kind of "X-ray vision", easing the information transfer from the P3DM to the actual patient for a surgical procedure. Such AR technology can overcome some of minimally invasive surgery (MIS) limitations such as trocar/instrument placement in thoracic surgery [START_REF] Hallet | Trans-Thoracic Minimally Invasive Liver Resection Guided by Augmented Reality[END_REF]. Indeed, the trocar placement can then be decided before the surgery on the P3DM (cf. Fig. 1[a]) and this location is superimposed intra-operatively on a static view provided by a fixed camera. This AR view allows fast, safe and optimal trocar set-up to provide the best surgical approach to reach the target (cf. Fig. 1[b-c]). However, this technique suffers from two important drawbacks. First, the P3DM registration is performed manually and needs to be recomputed after every change of the relative position of the camera with respect to the patient. Second, this kind of relative motion is difficult to avoid even using bulky fixing methods for both the camera and the patient. The accuracy of the image augmentation with the P3DM depends on two factors: 3D camera tracking and P3DM registration. 3D camera tracking implies the computation, in real-time, of the 3D camera position with respect to the patient body. P3DM registration implies anchoring, also in real-time, of the P3DM to the patient body. In this paper, we propose a Visual Simultaneous localization And Mapping (VSLAM)-based approach for on-patient AR visualization, while the tablet is moved by the practitioner around the patient. Our system satisfies the following constraints: seamless integration in the operating room (OR), realtime performance, minimal interaction from the medical staff and robustness to occlusion and failure.

The reminder of the paper is organized as follows. Section 2 provides a review of the related work on on-patient AR visualization and VSLAM for camera tracking and mapping. Section 3 gives an overview of the proposed on-patient AR visualization system, followed by a detailed system description in Section 4. In section 5, the experimental results and evaluations are discussed. Finally, a conclusion summarizes our achievements and future work are presented in section 6.

Related work

On-Patient AR visualization

Various approaches for on-patient AR have been proposed in recent years, based on two different techniques to track 3D camera location: surface-based registration and 2D/3D point correspondences. In surface-based registration techniques [START_REF] Kilgus | Mobile markerless augmented reality and its application in forensic medicine[END_REF][START_REF] Dos Santos | Pose-independent surface matching for intra-operative soft-tissue marker-less registration[END_REF][START_REF] Macedo | High-Quality On-Patient Medical Data Visualization in a Markerless Augmented Reality Environment[END_REF][START_REF] Lee | Medical augment reality using a markerless registration framework[END_REF][START_REF] Chen | Development of a surgical navigation system based on augmented reality using an optical see-through head-mounted display[END_REF], a tablet-PC is mounted with a range camera, RGB-D sensor or stereo-vision to continuously capture the depth and color information, from which the skin surface is automatically extracted. This surface is then registered with the P3DM acquired from CT images, typically ICP (Iterative Closest Point) is used. This process is repeated for every frame at 5-10Hz [START_REF] Macedo | High-Quality On-Patient Medical Data Visualization in a Markerless Augmented Reality Environment[END_REF]. The major drawbacks of this type of techniques are the computation cost of depth image segmentation and ICP. Secondly, a good initialization for the ICP registration is required and has to be provided manually. Depending on the interface, it is not clear whether practitioners/surgeons can accept this task. Thirdly, this kind of methods are not robust to partial occlusions of the skin. To achieve real-time performance, either parallel processing [START_REF] Macedo | High-Quality On-Patient Medical Data Visualization in a Markerless Augmented Reality Environment[END_REF] or client/server architecture [START_REF] Chen | Development of a surgical navigation system based on augmented reality using an optical see-through head-mounted display[END_REF] or both [START_REF] Kilgus | Mobile markerless augmented reality and its application in forensic medicine[END_REF] are used: a powerful server PC is necessary to process data and the tablet-PC is used as a display tool only.

In 2D/3D point correspondences techniques [START_REF] Rassweiler | iPad-Assisted Percutaneous Access to the Kidney Using Marker-Based Navigation: Initial Clinical Experience[END_REF][START_REF] Muller | Mobile augmented reality for computer-assisted percutaneous nephrolithotomy[END_REF][START_REF] Sun | Validation of anatomical landmarks-based registration for image-guided surgery: an in-vitro study[END_REF][START_REF] Schneider | Landmark-Based Surgical Navigation. Clinical Image-Based Procedures[END_REF], markers need to be visible in the CT image and can be either natural landmarks or artificial ones placed on the patient before scanning. Their 2D positions in each frame are used to solve a 2D/3D geometrical relationship. In this case, the obvious drawback is the use of markers which should be visible in CT/MRI. Moreover, a minimum number of markers must be visible in every frame to register the P3DM [START_REF] Rassweiler | iPad-Assisted Percutaneous Access to the Kidney Using Marker-Based Navigation: Initial Clinical Experience[END_REF][START_REF] Schneider | Landmark-Based Surgical Navigation. Clinical Image-Based Procedures[END_REF], which impedes the surgeon movements. Indeed, markers are likely to be occluded by surgeon hand or a surgical instruments. Moreover, current on-patient visualization techniques typically evaluate their accuracy by measuring the registration errors of skin fiducials [START_REF] Kilgus | Mobile markerless augmented reality and its application in forensic medicine[END_REF][START_REF] Dos Santos | Pose-independent surface matching for intra-operative soft-tissue marker-less registration[END_REF][START_REF] Lee | Medical augment reality using a markerless registration framework[END_REF][START_REF] Chen | Development of a surgical navigation system based on augmented reality using an optical see-through head-mounted display[END_REF], the discrepancy in pixels in the image [START_REF] Muller | Mobile augmented reality for computer-assisted percutaneous nephrolithotomy[END_REF][START_REF] Schneider | Landmark-Based Surgical Navigation. Clinical Image-Based Procedures[END_REF], and/or the processing time [START_REF] Macedo | High-Quality On-Patient Medical Data Visualization in a Markerless Augmented Reality Environment[END_REF]. To address these drawbacks, a VSLAM-based method for on-patient visualization is proposed. Our system is rigorously evaluated in terms of: processing time and robustness on human data, registration accuracy on pigs during both breath-hold and respiration phases using fiducials and registration accuracy on a liver phantom using fiducials.

VSLAM-based camera tracking

VSLAM is a popular topic in robotics and computer vision, as it aims at building a 3D map of an unknown environment while simultaneously tracking camera location. Davison [11] proposed the first real-time VSLAM, based on an Extended Kalman Filter (EKF). Further improvements over the EKF SLAM have been proposed [START_REF] Civera | Inverse Depth Parametrization for Monocular SLAM[END_REF][START_REF] Civera | 1-Point RANSAC for extended Kalman filtering. Application to realtime structure from motion and visual odometry[END_REF]. EKF approaches initialize robustly, but have a poor scaling given its limitation of the map to a few hundreds points. EKF approaches have been successfully applied to MIS notably for abdominal surgery to track endoscope motion, to provide 3D scene structure [START_REF] Grasa | Visual SLAM for Handheld Monocular Endoscope[END_REF] and to live AR annotations [START_REF] Grasa | EKF Monocular SLAM 3D Modeling, Measuring and Augmented Reality from Endoscope Image Sequences[END_REF].

A significant leap with respect to EKF VSLAM was Klein and Murray's PTAM (Parallel Tracking and Mapping) [START_REF] Klein | Parallel Tracking and Mapping for small AR Workspace[END_REF]. PTAM algorithm performs in real-time all the steps of the classic photogrammetric 3D reconstruction [START_REF] Mikhail | Introduction to modern photogrammetry[END_REF]: matching, view selection, initialization and non-linear optimization termed Bundle Adjustment (BA). Recently, ORB-SLAM [START_REF] Mur-Artal | ORB-SLAM: A Versatile and Accurate Monocular SLAM System[END_REF] builds on PTAM and extends its performance to large scale mapping including loop closure, large scale relocation and position graph optimization. In this system, ORB image features play a key role in data association and relocation. ORB [START_REF] Rublee | ORB: An efficient alternative to SIFT or SURF[END_REF] is a newly developed image point descriptor that is able to handle rotation and scale changes with a comparable performance to SIFT [START_REF] Lowe | Object recognition from local scale-invariant features[END_REF] but with only a fraction of the computational cost. Our VSLAM combines some of these systems in order to balance good performance and reasonable hardware requirements in the medical environment. We use EKF for automatic initialization, PTAM approach for operation in room size environments where loop closure is not needed, and ORB features are used for relocation.

Our contributions are: 1) A VSLAM for on-patient AR visualization, which only requires a tablet-PC. 2) A usage strategy that fits the clinical constraints and is easy to setup and use inside the OR. 3) Interactions from the surgeon are reduced to the identification of 4 to 6 anatomical references at the beginning of the procedure. 4) The system is validated providing geometrical accuracy and computing cycle time.

System overview

The workflow of our system is shown in Fig. 2. It firstly consists of an offline stage. A CT volume of the patient is acquired and segmented to generate the P3DM. The P3DM is composed of surface meshes corresponding to the skin surface and to the selected body structure surfaces. The practitioner/surgeon selects at least 4 (typically between 4 and 6) anatomical landmarks (L i , i ∈ {1..6}) on the skin mesh (called anchor points). The anchors should be easily identifiable on the skin of the patient during the procedure.

Fig. 2: Overall system workflow

In the OR, the practitioner directs the tablet-PC camera at the patient and performs a translational motion to bootstrap the VSLAM, then he/she identifies the anchor points (L i ) by clicking on the tablet-PC live video stream. Those anchor points provide the 3D/3D registration of the P3DM within the VSLAM map. Afterwards, a synthetic image of the P3DM can be overlaid on live video stream. This interactive procedure of anchor points identification is required only once at the beginning of the procedure, then, the practitioner can move the tablet-PC around the patient and experience the AR visualization (cf. Sec. 4.3), even if none of the anchors remains visible in the tablet-PC camera field of view.

System description

VSLAM architecture

Our VSLAM aims at running on mobile devices and is intended for small scenes. It is a PTAM-like algorithm that automatically and sequentially computes a 3D photogrammetric reconstruction from the live video stream (cf. Fig. 3). To do so, a set of interest points has to be matched along the sequence; we use sparse features detected in the image by the Features from Accelerated Segment Test (FAST) detector [START_REF] Rosten | Machine learning for high-speed corner detection[END_REF]. In order to reduce number of outliers, we keep only the most salient features, whose Shi-Tomasi score [START_REF] Shi | Good features to track[END_REF] is over 100. The ORB descriptor [START_REF] Rublee | ORB: An efficient alternative to SIFT or SURF[END_REF] is then used to describe the detected features (block A in Fig. 3). A set of frames -named keyframes-has to be selected because the complexity scales cubically with the number of frames. Then a non-linear BA optimization yields 3D location for the map points and positions for the keyframes accurately. The BA is run in a thread termed Mapping described below. In parallel, another thread -termed trackingestimates the position of the frames that are not keyframes.

Camera tracking This task operates sequentially on all frames of the live video (block C in Fig. 3). The 3D locations of the map points are assumed to be available. Then the position of each frame is computed by non-linear optimization of the reprojection error for the matched points. To avoid the influence of spurious matches, a two-stage optimization is applied. In the first stage the Huber influence function is used as it is less sensitive to outliers. While in the second stage the optimization is switched to the Tukey kernel to achieve a robust optimization. An initial guess is needed for the camera position and is computed from the camera position and velocity estimated for the previous frame.

The difference among the various VSLAM methods is how the matches between the map points and the current frame are computed. We estimate a region where the map points are expected to be found by reprojecting the 3D map points onto the predicted camera position. The ORB descriptor of each map point is compared with those of all the features detected inside the predicted region, using the ratio between closest to second-closest neighbors as a score [START_REF] Lowe | Distinctive Image Features from Scale-Invariant Keypoints[END_REF]. If no matching is found, a correspondence is searched by patch correlation in the prediction region. If the number of matches is below 20 (empirically defined), the camera is assumed to be lost and the relocation process is started. Additionally, the tracking thread chooses a keyframe among the processed frames using the standard VSLAM criteria of minimal parallax distance with respect to all map keyframes. For each of the keyframes, we compute the median parallax with respect to the current frame. If the smallest parallax angle is over 2 degrees, the current frame becomes a new keyframe. This parallax angle threshold is chosen to be small to increase the number of keyframes, map points and to avoid rapid tracking loss.

The median parallax is estimated using the median of the XYZ coordinates of the map points detected in the current frame. Mapping The mapping thread runs in parallel with the tracking thread but at a lower frequency, continuously improving the map points estimation (block D in Fig. 3). The BA minimizes total reprojection error with respect to the keyframe positions X W C i and the 3D map point locations X W j : arg min

X W j ,X W C i i,j ρ ||u ij -CamProj(X W j , X W C i )| | 2 (1)
where u ij is the matched observation of the j -th map point by the i-th keyframe.

The CamProj codes the projection function including perspective and radial distortion. ρ denotes the robust Huber influence function. The BA non-linear optimization has been implemented using Google Ceres [24]. After each BA iteration, spurious map points are removed if: 1) the reprojection error of the point on the keyframe used for its creation is over a threshold (the median of reprojection errors of all measured map points). 2) the point is detected only in two keyframes, although it is visible in more than two.

3) The ratio between number of times point is measured and number of times point is predicted is smaller than 0.3 (empirically defined). The mapping thread is also responsible for the initialization of new map points. Once a new keyframe is added to the map, matches between the new keyframe ORB points and all other keyframes in the map are sought. We use standard patch correlation guided by epipolar geometry.

Bootstrapping Previous mapping and tracking processes assume that there is a map. Next we describe how the map is initialized from scratch (block B in Fig. 3). For bootstrapping, the system has to select two keyframes that render enough parallax, this selection has proven an issue in VSLAM. We use a simple EKF VSLAM with all features encoded in inverse depth [START_REF] Civera | Inverse Depth Parametrization for Monocular SLAM[END_REF]. This approach can handle low parallax geometries, being able to exploit every single image in the sequence to estimate the map and the camera position. We process images until most of the map points are detected with enough parallax. Then we consider the first and the last processed images as the two initial keyframes for BA. Given these two keyframes and their relative locations, robust new point matches are computed by epipolar search, as in the mapping thread. Afterwards, an initial guess for the map points and two keyframe positions is fed to BA. The proposed EKF VLSAM only initializes points in the first frame. If those points fail to be tracked, or go out of the field of view before rendering enough parallax, the map is discarded and a new initialization is launched automatically.

Camera relocation Tacking can be lost because of camera occlusion, feature deletion due to fast camera motion or failure to track enough map points. Then the camera has to be located with respect to the map from scratch. Our system detects all the ORB points in the current image (block E in Fig. 3). They are matched with respect to the ORB descriptors of all the map points using as score the ratio between closest to second-closest neighbors to compute the putative matches. Then a perspective-three-point (P3P) [START_REF] Gao | Complete solution classification for the perspective-three-point problem[END_REF] from random samples of size 4 is executed.

The number of Random Sample Consensus (RANSAC) iterations are limited by the frame rate. To validate the relocation, the tracking algorithm has to produce a coherent position for the next frame in the sequence, otherwise relocation is re-attempted with the new frame.

Registration of pre-operative model with VSLAM map

Registration is initialized interactively by the practitioner once the VSLAM has been bootstrapped. The practitioner selects the 2D anchor points over the live video stream by tapping on the tactile screen of the tablet-PC. The 3D coordinates of the 2D anchor points are computed and appended to the map following the procedure described in Algorithm 1 (a variant of [START_REF] Gálvez-López | Real-time Monocular Object SLAM[END_REF]). The P3DM is then translated, rotated and scaled to align the landmarks in the model with the anchors in the map. Initialized anchor points store a correlation patch to match in other keyframes. The correlation matching is guided by the epipolar geometry. Once two keyframes with proper matches are found, the 3D coordinates of the anchor point are triangulated, then the matches are propagated among previous keyframes. The BA is iterated to refine the 3D geometry of the map anchors. These anchor points are never removed from the map. On the arrival of a new keyframe, geometry-guided search for correlation matches is performed and T pm (rotation, translation and scaling) is recomputed using [START_REF] Horn | Closed-form solution of absolute orientation using unit quaternions[END_REF] by minimizing the alignment error err = n i=1 l i -T pm • L i (cf. Fig. 4 [b]). After each mapping iteration of the mapping thread, T pm is refined with the newly available l i estimations.

Input : List of keyframes with their estimated positions Input : Query image in which anchor points are selected Input : 3D coordinates of landmarks from P3DM (L i ) Output: Transformation Tpm from P3DM to map foreach selected 2D anchor point do Extract a square patch around the point Perform epipolar-guided cross correlation with all the map keyframes Estimate 3D coordinates (l i ) using triangulation Find more matches in previous keyframes Append l i to the map end Apply BA to refine (l i ) estimations Compute Tpm from L i and l i [START_REF] Horn | Closed-form solution of absolute orientation using unit quaternions[END_REF] Algorithm 1: P3DM registration using anchor points from a real image

See-through AR

To provide AR overlay on the live video: 1) The VSLAM tracking thread provides a position estimate for each frame of the live stream, then a virtual camera with the same intrinsic parameters of the tablet-PC camera is located at the estimated position in the virtual scene (cf. Fig. 4[b]). 2) The image acquired by the virtual camera, taking into account the tablet-PC camera distortion, is rendered (cf. Fig. 4[c]). 3) The fusion is performed (cf. Fig. 4[d]) between the real camera image (cf. Fig. 4[a]) and the rendered one (cf. Fig. 4[c]). 5 Results

The proposed system has been implemented in C++ with OpenCV and VTK libraries and executed on a Sony VAIO Duo 13 tablet-PC with Intel(R) Core i7 (1.8 GHz), 8 GB RAM, with a camera of 640x480 at 30fps. Firstly, the system performance was evaluated on in-vivo data in terms of computation time with two volunteers, each of them laying on the table while the practitioner holds the tablet-PC and moved around them. The CT scans of these volunteers were performed several years ago. Secondly, the system accuracy was assessed by means of several experiments with fiducials; the first experiments were on four in-vivo pigs, the second on a phantom. All computations were performed exclusively on the tablet-PC. Beforehand, the focal and distortion of the tablet-PC camera were calibrated using [START_REF] Zhang | A Flexible New Technique for Camera Calibration[END_REF]. Every P3DM in our experiments was segmented using our own software but can be obtained, for clinical applications, using a commercial service like Visible Patient [29]. More details can be appreciated in the accompanying video.

Volunteers experiments: computation time evaluation.

In this experiment, the time required for each step of the system was evaluated. For both volunteers, the five anatomical landmarks chosen as anchors for registering the P3DM were: right nipple, left nipple, umbilicus, right iliac crest and left iliac crest. The left and right iliac crests were marked with a pen on the skin of both volunteers, to easily identify them in the 2D images (cf. Fig. 5[a] and Fig. 6[a]). Fig. 5 and Fig. 6 show AR annotated frames for both volunteers from different points of view. VSLAM initialization: the VSLAM bootstrapping did not fail in any of the experiments. It was initialized using on average less than 20 frames. If failed, the initialization was automatically relaunched, and eventually succeeded.

Camera tracking and VTK rendering: Average tracking time was approximately 32 ms per frame for a map size that ranged between 180-200 points and 30-40 keyframes for video sequences composed of 750-900 frames. The average VTK rendering time was approximately 33 ms per frame, including the ideal projective imaging, the distortion and the fusion with the real frame. After each mapping step, the anchors 3D locations were updated, hence the AR insertion location in the map had to be recomputed, which took less than 1.2 ms. The time of initial insertion of the P3DM into the map can take up to 3 seconds depending on the sequence, due to searching for the anchor matches in all the keyframes. Therefore, total average time was 66.2 ms per frame. Loss of tracking and relocalization performance: In case of lateral (cf. Fig. 5[f] and Fig. 6[c].) or close up (cf. Fig. 6[f]) tablet-PC movements, the P3DM can still be registered even if most of the anchor points are not visible. Camera tracking is robust to partial scene occlusion since few map points are needed for VSLAM to estimate the tablet-PC position (cf. Fig. 9[b] and [c]). In case of full scene occlusion or severe camera motion, the relocalization module in VSLAM always relocated the tablet-PC position once a few map points were visible again, which required approximately 15 ms. As a result of this module, re-initialization of the whole system in the case of tracking loss is not necessary.

(a) (b) (c) (d) (e) (f)

Accuracy evaluation

To assess the registration accuracy of the proposed system, experiments on four pigs were performed and the surface fiducial registration error (FRE) as well as the target registration error (TRE) were reported. Additionally, a plastic phantom was used to evaluate the registration accuracy on internal body structures that are far from the anchor points used for registration.

Data acquisition

Each pig was placed on the CT table, and nine radio-opaque markers were stuck on its skin before acquisition (cf. Fig. 7[a]). The CT scan was performed with breath-hold via a mechanical ventilation system. For each pig, two videos were recorded, one with breath-hold and another during respiration. The P3DM and the 3D coordinates of the markers were extracted from the CT images.

Registration accuracy on pigs

The 3D coordinates of the markers extracted from CT were considered as ground truth. All markers were clicked on 2D image and their 3D coordinates computed and appended to the map. Only five markers were used as anchors to compute the 3D/3D registration those are displayed in blue in Fig. 7[a]. Fig. 7[b] and [c] shows the registration results on two pigs from different directions. The keyframe locations during camera motion around one of the pigs are displayed in Fig. 7[c] and represented by axes. The averaged F RE over all the frames in the sequence was calculated from the five markers used in the registration. The averaged T RE over the sequence was computed from the remaining four markers. F RE and T RE are defined in eq. ( 2):

F RE = 1 F F f =1 1 5 5 i=1 ||l i -T pm • L i || T RE = 1 F F f =1 1 4 9 i=6 ||l i -T pm • L i || (2)
where F refers to the number of processed frames. As defined in eq. ( 2), the distance between the two point clouds was computed for every frame. In the Inner summation of eq. ( 2), the average distances of the five markers used for the registration and average distances of the remaining four markers were computed. Then F RE and T RE over all frames in the sequence were defined from the outer summation in eq. ( 2). The length of all video sequences ranged between 600 and 800 frames with 30 to 40 keyframes and map sizes between 176 and 279 points.

Each video was processed five times, each time the same frame was used to select the anchors. For the five registration trials on each pig sequence, the minimum, maximum and mean values of F RE and T RE are reported in Table 1 during breath-hold. Table 2 shows the influence of the breathing. After the initial insertion into the map, the P3DM is affected by a small jittering, due to the low number of keyframes and the poor geometrical conditioning. On the arrival of keyframes with a wider baseline, thus rendering bigger parallax, this jittering disappears within a few seconds, according to our experiments. Afterwards, the estimation of the anchor points 3D coordinates becomes accurate and so does the 3D/3D registration (cf. Fig. 8).

Registration accuracy on phantom

A phantom with a plastic liver was used to evaluate the system accuracy for points far from the body surface. 13 markers were attached on the external surface, 2 markers on the plastic liver and 4 markers on the phantom base (cf. Fig. 9[a]). The phantom sequences and CT were obtained following the same steps as those of the pig experiments in Sec. 5.2.1. Five markers on phantom surface were used to compute the registration (cf. Fig. 9[a]). Table 3 shows F RE of the five markers used for the registration, T RE of the two liver markers and T RE of the four markers at the phantom base. All were computed after processing the full sequence and averaging the results of five registration trials. As shown in Table 3, the closer the target to the skin, the better the registration accuracy. The four markers on the phantom base represent the worst target position, i.e. close to the skin of the back. Therefore, 11.8mm can be considered the worst system accuracy. For the sake of completeness, the registration accuracy using all the 13 markers stuck on the surface has also been computed and provide a reduction between 1.5-2.0 mm on both F RE and T RE. It is worth noting that the F RE is larger than in case of pigs due to utilization of different markers. The markers used with pigs were covered and a pen was used to mark their center of mass to be easily identifiable in the images (cf. Fig. 7[a]). In the phantom the markers were not covered and hence there were some inaccuracies in clicking the center of the cross shape of each markers.

Conclusion and Future work

A VSLAM-based on-patient AR visualization system is presented, which can be seamlessly integrated into the OR as the only external device is a commercial tablet-PC computer. The proposed system provides real-time performance, robustness to occlusion and detection failure. It requires minimal interaction with medical staff, i.e. the definition of the anchors by clicking on the live video. This is interaction is considered non-disruptive by most surgeons. In contrast to markerbased AR, our system is able to provide AR overlays even if none of the anchors used for the registration remains visible. Experimental results show the applicability of the proposed system, both in terms of computation time and accuracy. Although the system can already provide a great assistance, it can be further improved. Extension from rigid to non-rigid registration would allow taking breathing motion into account.

Fig. 1 :

 1 Fig. 1: Port positioning with AR guidance during trans-thoracic minimally invasive hepatectomy [1]. (a) Pre-operative trocar placement planning. (b) and (c) Marking of the chosen port site.

Fig. 3 :

 3 Fig. 3: VSLAM architecture.

Fig. 4 :

 4 Fig. 4: AR insertion. (a) Tablet-PC camera frame with projected (red) and matched (green) map points. (b) Virtual 3D scene including the registered P3DM. (c) Virtual camera image. (d) Fused AR image.

Fig. 5 :

 5 Fig. 5: Experiment on first volunteer. (a) Real tablet-PC image. (b) Skin registration with anchor points (in blue) and automatic map points. (c-f) Skin AR overlays over frames from different points of view.

Fig. 6 :

 6 Fig. 6: Registration of transparent skin, liver, left kidney and right kidney on the body of the second volunteer from different points of view.

Fig. 7 :

 7 Fig. 7: Experiments on pigs. (a) Nine radio-opaque markers were attached to the surface of the pig. (b) P3DM composed of skin, bones, liver, left kidney and right kidney overlaid on an image of the first pig. (c) Lateral view of skin registration on the second pig. (d) Keyframe locations during camera motion around a pig.

Fig. 8 :

 8 Fig. 8: Evolution of the average distances between l i and L i (in mm) over all frames of one video sequence in case of breath-hold. (a) The average distances between the 5 markers used in the registration. (b) The average distances of the remaining 4 markers.

Fig. 9 :

 9 Fig. 9: Experiments on phantom. (a) Markers in blue were used for the registration. (b) and (c) Liver AR overlay with partial scene occlusion.

Table 1 :

 1 F RE and T RE (in mm) of the four pigs sequences recorded during breath-hold.

			Pig 1			Pig 2			Pig 3			Pig 4
		min	max mean	min	max mean	min	max mean	min	max mean
	F RE	2.72	3.07	2.94	2.41	2.61	2.52	3.42	4.27	3.82	1.01	2.66	1.55
	T RE	3.36	3.99	3.74	3.38	3.98	3.69	3.75	4.28	4.07	1.5	2.88	2.2

Table 2 :

 2 F RE and T RE (in mm) of the four pigs sequences recorded during breathing

			Pig 1			Pig 2			Pig 3			Pig 4
		min	max mean	min	max mean	min	max mean	min	max mean
	F RE	2.53	3.64	3.11	2.76	3.7	3.1	4.62	6.09	5.32	2.62	4.22	3.1
	T RE	3.49	4.56	3.92	3.65	4.08	3.88	4.95	6.24	5.6	2.62	4.62	3.44

Table 3 :

 3 F RE and T RE (in mm) after processing the whole phantom sequence.

		F RE		T RE (two liver markers)	T RE (four markers at the base)
	min	max mean	min	max	mean	min max	mean
	2.28	7.74	5.1	5.16	9.7	6.61	9.9	13.7	11.8

Acknowledgment: This work is supported by the Direccíon General de Investigacíon Centífica y Téctinca of Spain under Project DPI2015-67275-P

Conflict of Interest:

The authors declare that they have no conflict of interest.

Ethical approval: All procedures performed in studies involving human participants were in accordance with the ethical standards of the institutional and/or national research committee and with the 1964 Helsinki declaration and its later amendments or comparable ethical standards. In addition to All applicable international, national, and/or institutional guidelines for the care and use of animals were followed.

Informed consent: Informed consent was obtained from all individual participants included in the study.