A Sliding-rod Variable-strain Model for Concentric Tube Robots Federico Renda 1 , Conor Messer 2 , Caleb Rucker 3 , and Frederic Boyer 4 Abstract-In this work, the Piecewise Variable-strain (PVS) approach is applied to the case of Concentric Tube Robots (CTRs) and extended to include the tubes' sliding motion. In particular, the currently accepted continuous Cosserat rod model is discretized onto a finite set of strain basis functions. At the same time, the insertion and rotation motions of the tubes are included as generalized coordinates instead of boundary kinematic conditions. Doing so, we obtain a minimum set of closed-form algebraic equations that can be solved not only for the shape variables but also for the actuation forces and torques for the first time. This new approach opens the way to torquecontrolled CTRs, which is poised to enhance elastic stability and improve interaction forces' control at the end-effector.

Index Terms-Modeling, Control, and Learning for Soft Robots; Surgical Robotics: Steerable Catheters/Needles.

I. INTRODUCTION

C ONTINUUM robots can be considered a class of soft manipulators, particularly suited for Minimally Invasive Surgery (MIS) [START_REF] Burgner-Kahrs | Continuum robots for medical applications: A survey[END_REF]. One of the most promising continuum robotic systems developed so far is the Concentric Tube Robot (CTR), a collection of nested millimeters tubes whose elastic interaction is used to control the system's overall shape and the iteration force at the end-effector for surgical intervention [START_REF] Mahoney | A review of concentric tube robots: modeling, control, design, planning, and sensing[END_REF]. CTRs have demonstrated promising results in a variety of MIS applications. For a thorough survey of the clinical applications of CTRs, see [START_REF] Alfalahi | Concentric tube robots for minimally invasive surgery: Current applications and future opportunities[END_REF] and the references therein.

The development of CTRs has been made possible thanks to rapid advancement in modeling these kinds of complex non-linear systems. These research efforts have converged on a select type of Cosserat rods model with a particular kinematic structure representing the internal tubes' additional rotational motion [START_REF] Dupont | Design and control of concentric-tube robots[END_REF], [START_REF] Rucker | A geometrically exact model for externally loaded concentric-tube continuum robots[END_REF]. Although this Cosserat approach has been improved over the years to include tubes' clearance, friction [START_REF] Ha | Modeling tube clearance and bounding the effect of friction in concentric tube robot kinematics[END_REF], and inertial dynamics [START_REF] Till | A dynamic model for concentric tube robots[END_REF], it has been essentially maintained as initially proposed [START_REF] Hunter | Concentric Tube Robots: The State of the Art and Future Directions[END_REF]. Here, we refer to it as the standard CTR model, and we briefly recall it in section II.

On the other hand, a general soft robot is composed of flexible and rigid elements arranged in a parallel or serial fashion and capable of locomoting in the environment. A novel coordinate system has been proposed to model these kinds of systems, which discretizes the continuous Cosserat model of the flexible components onto a finite set of strain basis functions [START_REF] Renda | A geometric variable-strain approach for static modeling of soft manipulators with tendon and fluidic actuation[END_REF], [START_REF] Boyer | Dynamics of continuum and soft robots: A strain parameterization based approach[END_REF]. This Piecewise Variable-strain (PVS) approach is a generalization of traditional robotics' geometric model [START_REF] Murray | A Mathematical Introduction to Robotic Manipulation[END_REF] to the case of highly flexible or soft robots [START_REF] Renda | A geometric and unified approach for modeling soft-rigid multi-body systems with lumped and distributed degrees of freedom[END_REF]. Thus, it provides the theoretical framework for applying traditional control strategies to the field of soft robots. One of the aims of the present work is to apply the PVS approach to CTR systems. In this way, we further extend the general PVS approach to this essential soft robotic technology. At the same time, we open new prospects for the control design of CTRs. In particular, the PVS model provides the equilibrium equations as a minimum set of closed-form algebraic equations, which may be easier to handle for control, design optimization, and stability assessment.

Furthermore, in this work, we tackle one of the main limitations of the standard CTR model concerned with the sliding structure problem, also known as the spaghetti problem [START_REF] Carrier | The spaghetti problem[END_REF]. In the standard model, the tube rotations and insertions dictate the boundary conditions and locations of the section domains. While the rotational actuator torques are easily computed from the model solution, the insertion forces are not present in the model equations. The recent extension to dynamics also does not model the insertion forces and assumes negligible insertion velocity and acceleration [START_REF] Till | A dynamic model for concentric tube robots[END_REF]. Here we relax this assumption, although still in a static setting, which allows calculation of the actuation insertion forces required for the equilibrium for the first time. The spaghetti problem appears in the planar models of animal locomotion [START_REF] Cicconofri | A study of snake-like locomotion through the analysis of a flexible robot model[END_REF], structural stability [START_REF] Armanini | Configurational forces and nonlinear structural dynamics[END_REF], and finite element analysis [START_REF] Humer | General sliding-beam formulation: A non-material description for analysis of sliding structures and axially moving beams[END_REF]. Thus, to the authors' knowledge, it is also the first time the spaghetti problem is tackled in three dimensions and with an additional relative rotation.

II. MODEL PRELIMINARIES

Before applying the PVS approach to solve the CTR equilibrium, we revise the standard CTR model [START_REF] Dupont | Design and control of concentric-tube robots[END_REF], [START_REF] Rucker | A geometrically exact model for externally loaded concentric-tube continuum robots[END_REF] using geometric notation. For simplicity, we will present the case of fully overlapping tubes here. Note that additional details about the CTR kinematics will be provided for the more general cases analyzed in section III. Furthermore, a new kinematic model that includes the sliding and rotation of the the bases of the tubes as generalized coordinates of the system will be introduced in section IV.

A. Kinematics

Each tube of a concentric tube system can be modeled as a Cosserat rod, a continuous set of rigid cross sections of infinitesimal thickness along a material curvilinear abscissa X ∈ [0, L] where L is the total length of the rod. Identifying each rigid cross-section with the moving frame rigidly attached to it, the configuration of the tube j is completely defined by a curve g j (•) : X → g j (X) ∈ SE(3) represented by the homogeneous matrix:

g j (X) = R j (X) r j (X) 0 T 1 , (1)
where r j (X) ∈ R 3 is the position vector of the origin of the moving frame and R j (X) ∈ SO(3) is the rotation matrix representing its orientation with respect to the spatial frame.

The high aspect ratio and material property of the conventional CTR allow assuming inextensibility and shearless deformation (Kirchhoff-Love kinematics). Furthermore, neglecting tube clearance, the rods are perfectly concentric. Thus, their centerlines must be the same at any configuration. Then, it follows that an inner tube j differs from its outer tube j -1 only by a rotation around the tangent vector to the centerline [START_REF] Rucker | A geometrically exact model for externally loaded concentric-tube continuum robots[END_REF], defined by:

g θj (X) =     1 0 0 0 0 cos(θ j (X)) -sin(θ j (X)) 0 0 sin(θ j (X)) cos(θ j (X)) 0 0 0 0 1     , (2)
where θ j (X) is the rotation angle and we have assumed a local coordinate frame as shown in Figure 1. Thus, for an ordered set of tubes we obtain:

g 1 (X) = g 1 (X) g j (X) = g j-1 (X)g θj (X), (j = 2, ..., N ) (3)
where N is the total number of tubes and the outer tube is numbered with 1. A pictorial representation of this kinematics is shown in Figure 1. Note that in this model, the concentric tubes share the same material abscissa X. Let's obtain the time (∂/∂t = ˙) and space (∂/∂X = ) derivative of the tubes configuration g with the Lie algebra se(3). Starting with tube 1 we get:

g -1 1 (X) ġ1 (X) ∨ = η 1 (X) = w T 1 , v T 1 T ∈ R 6 (4) g -1 1 (X)g 1 (X) ∨ = ξ 1 (X) = k T 1 , u T 1 T ∈ R 6 . (5)
where v(X), u(X) ∈ R 3 represent the linear strains and velocity respectively, while w(X), k(X) ∈ R 3 are the angular strains and velocity respectively. All the quantities are defined in the local frame at X. To indicate the angular strain in the reference stress-free configuration, we use the notation k * . Finally, the superscript ∨ indicates the isomorphism between the Lie algebra se(3) and R 6 [START_REF] Murray | A Mathematical Introduction to Robotic Manipulation[END_REF] (∧ will be used in the opposite direction).

Equating the mixed partial derivatives of the tube configuration, we obtain a compatibility equation between the velocity and strain twists.

∂ ∂t g 1 = ∂ ∂X ġ1 =⇒ η 1 = ξ1 -ad ξ1 η 1 . (6)
As shown in [START_REF] Renda | A geometric variable-strain approach for static modeling of soft manipulators with tendon and fluidic actuation[END_REF], integrating [START_REF] Ha | Modeling tube clearance and bounding the effect of friction in concentric tube robot kinematics[END_REF] with respect to space yields the following useful relation:

η 1 (X) = Ad -1 g1(X) X 0 Ad g1 ξ1 ds , ( 7
)
where the operator Ad is the Adjoint map in SE(3) defined in Appendix A.

For what concerns the inner tubes, first we define the derivatives of the relative rotation g θj .

g -1 θj (X) ġθj (X) ∨ = η θj (X) = θj (X) 0 0 0 0 0 T (8) g -1 θj (X)g θj (X) ∨ = ξ θj (X) = θ j (X) 0 0 0 0 0 T . ( 9
)
The same procedure that led to equation [START_REF] Till | A dynamic model for concentric tube robots[END_REF], in this case yields:

η θj (X) = X 0 ξθj ds , (10)
Finally, the chain rule together with (4), ( 5), [START_REF] Hunter | Concentric Tube Robots: The State of the Art and Future Directions[END_REF], and (9) yields, for a general tube j:

g -1 j ġj ∨ = ηj(X) = 2 i=j Ad -1 g θ i η1 (X) + j i=2 η θ i (X) (11) g -1 j g j ∨ = ξj(X) = 2 i=j Ad -1 g θ i ξ1 (X) + j i=2 ξ θ i (X). (12)

B. Statics

In this work, we focus our attention on the static equilibrium of concentric tubes with no external applied force (except for the actuation in section IV) to get more neat results. However, it should be noted that there are no particular impediments that would prevent extending the approach to the general dynamic case. We have also neglected friction forces and torques interaction between the tubes as usually done with the standard model.

The static equilibrium of a Cosserat rod in a concentric tube setting can be derived from Hamilton's principle (see [START_REF] Boyer | Poincare's equations for cosserat media: Application to shells[END_REF] for a derivation) with the addition of the constraints forces due to the concentricity constraint. Thus, for a tube j we obtain:

F ij + ad * ξj F ij + F λj = 0, (j = 1, ..., N ) , (13)
where

F ij (X) = M T ij (X), N T ij (X)
T ∈ R 6 is the vector of internal moment and force, F λj (X) ∈ R 6 is the vector of distributed constraint force, and ad * is the co-adjoint map in SE(3) defined in Appendix A. Note that the value of the internal force N ij is unknown due to the inextensibility and shearless constraints. On the other hand, the internal moment's value can be computed from a constitutive law that we assume linear for simplicity (Hooke's law).

M ij (X) = Σ j k j (X) -k * j (X) , (14)
where Σ = diag(GJ, EI, EI) ∈ R 3×3 is the elasticity matrix, E is the Young modulus, G the shear modulus, and J, I are, respectively, the polar and bending second moment of area of the circular cross-section. For what concerns the constraint contact force, although its value is also unknown, its basis can be obtained from the virtual constraints equations that enforce the concentricity constraint. These equations can be written in the form

j j • (δr j (X) -δr j-1 (X)) = 0, (j = 2, ..., N ) k j • (δr j (X) -δr j-1 (X)) = 0, (j = 2, ..., N ) , ( 15
)
where j j and k j point in the y-, and z-axis directions of tube j at X (see Fig. 1) , and δ is a virtual displacement. Equation [START_REF] Armanini | Configurational forces and nonlinear structural dynamics[END_REF] says that the distributed constraint force is equal and opposite between two consecutive tubes and that it takes the form:

F λj (X) = 0 0 0 0 λ yj (X) λ zj (X) T - Ad * g θ j+1 0 0 0 0 λ yj+1 (X) λ zj+1 (X) T , (16)
where λ y (X) and λ z (X) are Lagrange multipliers, and Ad * is the co-Adjoint map in SE(3) defined in Appendix A.

The static equilibrium of a collection of tubes imposes the following boundary conditions on equation ( 13) [START_REF] Rucker | A geometrically exact model for externally loaded concentric-tube continuum robots[END_REF].

N j=1 j k=1 Ad * g θ k M ij ( X-) = N j=1 j k=1 Ad * g θ k M ij ( X+ ).
(17) Note that ( 17) also holds at discontinuity points X where the tubes stop to overlap (in that case, the upper limit of the sum on the RHS will be smaller). Furthermore, since the concentric tubes can not transmit axial torque to one another, the following condition holds for each tube [START_REF] Rucker | A geometrically exact model for externally loaded concentric-tube continuum robots[END_REF].

i • M i ( X-) = i • M i ( X+ ) , ( 18
)
where i points in the x-axis direction, tangent to the midline (see Fig. 1).

III. VARIABLE-STRAIN CTR MODEL

In [START_REF] Renda | A geometric variable-strain approach for static modeling of soft manipulators with tendon and fluidic actuation[END_REF], [START_REF] Boyer | Dynamics of continuum and soft robots: A strain parameterization based approach[END_REF], a novel variable-strain approach has been presented to model soft manipulators driven by tendons and pneumatic chambers. This section aims to apply this technique to CTRs to obtain a minimum set of closed-form equations describing the system equilibrium.

According to the PVS approach, the configuration g j of the collection of tubes is represented by the strain fields ξ 1 and ξ θj . In particular, the infinite dimensional strain fields are discretized on a finite set of basis functions as follows.

ξ 1 (X) = B 1 (X)p 1 + ξ , (19) ξ θj (X) = B θj (X)p θj , (20)
where B 1 (X) ∈ R 6×n1 and B θj (X) ∈ R 6×nθj are matrix functions whose columns form the basis for the strain field ξ 1 and ξ θj , respectively, while p 1 ∈ R n1 and p θj ∈ R nθj are the vectors of coordinates. Note that, due to the assumptions made in section II, ξ is equal to [0 0 0 1 0 0], while the last three rows of B 1 (X) and the last five rows of B θj (X) are all equal to zero. Then, the set of generalized coordinates become q = (p 1 , p θ2 , . . . , p θ N ), and the configuration g j (q) can be reconstructed through the integration of equations ( 5), ( 9) and the recursive formula given by [START_REF] Alfalahi | Concentric tube robots for minimally invasive surgery: Current applications and future opportunities[END_REF]. The differential relation between configuration and generalized coordinates is obtained by replacing equations [START_REF] Webster | Toward active cannulas: Miniature snake-like surgical robots[END_REF], [START_REF] Webster | Mechanics of precurvedtube continuum robots[END_REF] in the velocity equations ( 7), [START_REF] Boyer | Dynamics of continuum and soft robots: A strain parameterization based approach[END_REF], which yields

η 1 (X) = Ad -1 g1 X 0 Ad g1 B 1 ds ṗ1 = S 1 (X) ṗ1 , (21) η θj (X) = X 0 B θj ds ṗθj = J θj (X) ṗθj , (22)
and, in turn, replacing the result in the velocity equation [START_REF] Murray | A Mathematical Introduction to Robotic Manipulation[END_REF]. Finally, we obtain:

η j (X) =   2 i=j Ad -1 g θ i S 1 (X)   ṗ1 + j i=2 J θi (X) ṗθi = J j (X) ṗ1 + j i=2 J θi (X) ṗθi . (23)
Note that J θi (X) is an analytical function that can be computed offline given the choice of basis B θi . The differential equation [START_REF] Ha | Optimizing tube precurvature to enhance the elastic stability of concentric tube robots[END_REF] provides the required Jacobians to project the static equilibrium ( 13) by d'Alembert's principle. In particular, equation ( 13) is projected with

L 0 J T j (•) dX and L 0 J T θi (•)
dX for all j and i. The results are then summed over all j from 1 to N to obtain as much algebraic equilibrium equation as the dimension of q. This procedure is shown in detail for the case of two fully overlapping tubes in the next section. The general case of multiple non-overlapping tubes with straight actuation, as shown in Figure 3, will be illustrated in the following.

A. Single overlapping sections

Let us consider the case of two fully overlapping tubes. Applying d'Alembert's principle as described above yields:

2 j=1 L 0 J T j F ij + ad * ξj F ij + F λj dX = 0 L 0 J T θ2 F i2 + ad * ξ2 F i2 + F λ2 dX = 0. (24)
Considering the form of the constraint force ( 16) and the identity

F ij + ad * ξj F ij = Ad * g -1 j
Ad * gj F ij , we get (note that Ad T g = Ad * g -1 ):

2 j=1 L 0 X 0 B T 1 Ad T g1 ds Ad gj F ij dX = 0 L 0 X 0 B T θ2 ds Ad * g -1 2 (Ad g2 F i2 ) dX = 0. (25)
Finally, integrating by part we obtain:

L 0 B T 1 F i1 + Ad * g θ 2 F i2 dX = 0 L 0 B T θ2 F i2 dX = L 0 J T θ2 ad * ξ1 Ad * g θ 2 F i2 dX , (26)
where we have used the strain equation ( 12), the identity

J T θ2 ad * ξ θ 2
F i2 = 0 and the boundary conditions ( 17), [START_REF] Sears | A steerable needle technology using curved concentric tubes[END_REF] at L. Note that, since the generalized coordinates q are independent, the unknown constraint forces F λj cancel out as expected.

Equation (26) can be solved numerically for the unknowns p 1 and p θ2 that appear in F ij , ξ 1 , and g θ2 through equations ( 19), ( 20), [START_REF] Cicconofri | A study of snake-like locomotion through the analysis of a flexible robot model[END_REF], and the integration of ( 9). This integration becomes simple in this case, since the rotation angle θ of (2) can be obtained by:

θ 2 (X) = α 2 + (J θ2 (X)p θ2 ) x , (27)
where () x extracts the first rotational element of a vector and α 2 is the rotation of the tube imposed at the base by the actuators. Note that the integrals of equation ( 26) can be expressed analytically with the use of hypergeometric functions, which makes their evaluation quick for numerical root finding.

Comparison with analytic torsionless case

As a first test to check the validity of the proposed model, we consider here the case of two torsionless, fully overlapping tubes with three-dimensional deformation. The torsionless approximation was used at the beginning of the CTR development, and analytic solutions are available for the constant curvature approximation [START_REF] Sears | A steerable needle technology using curved concentric tubes[END_REF], [START_REF] Webster | Toward active cannulas: Miniature snake-like surgical robots[END_REF], [START_REF] Webster | Mechanics of precurvedtube continuum robots[END_REF] as well as the general variable curvature case [START_REF] Rucker | Equilibrium conformations of concentric-tube continuum robots[END_REF].

Assuming no torsion, the equilibrium equation ( 26) becomes:

L 0 B T 1 (F i1 + F i2 ) dX = 0 , ( 28
)
with θ 2 = 0 and ξ 1 = ξ 2 = ξ. To reproduce the available analytic solution with variable curvature [START_REF] Rucker | Equilibrium conformations of concentric-tube continuum robots[END_REF], we select the following basis for the strain ξ.

B T 1 (X) =     0 (k * 1 ) y (X) 0 0 0 0 0 (k * 2 ) y (X) 0 0 0 0 0 0 (k * 1 ) z (X) 0 0 0 0 0 (k * 2 ) z (X) 0 0 0     , (29)
where k * j is the reference strain of tube j. Thus, n 1 = 4 and the generalized coordinates

p 1 = q = [q 1 q 2 q 3 q 4 ]
T represent the amount of reference curvatures making up the equilibrium configuration. Rearranging (28), we can write:

L 0 (k * 1 ) 2 y (k * 1 ) y (k * 2 ) y (k * 2 ) y (k * 1 ) y (k * 2 ) 2 y dX q 1 q 2 = 1 (EI 1 + EI 2 ) L 0 EI 1 (k * 1 ) 2 y + EI 2 (k * 1 ) y (k * 2 ) y EI 1 (k * 2 ) y (k * 1 ) y + EI 2 (k * 2 ) 2 y dX . (30)
A similar equation can be obtained for q 3 and q 4 . Solving for q, we get:

q = 1 (EI 1 + EI 2 ) [EI 1 EI 2 EI 1 EI 2 ] T , (31)
which correctly corresponds to the analytic solution k

= (Σ 1 + Σ 2 ) -1 (Σ 1 k * 1 + Σ 2 k * 2 ).

B. Multiple non-overlapping sections

Let us now move to the case of multiple non-overlapping tubes. To address this problem, we divide the system domain into piecewise variable-strain sections corresponding to the curvature discontinuity located at the end of every overlapping portion. Thus, the kinematic equation ( 3) is applicable locally at each overlapping portion and can be generalized as follows (see Figure 1):

g js (X) = j-1 i=1 g i (L i )g θ (i+1)s (L i ) g j (X) g ks (X) = g (k-1)s (X)g θ ks (X), (k = j + 1, ..., N ) , (32)
where g js indicates the configuration of tube j with respect to the spatial frame, while g j represents the configuration with respect to the local overlapping portion. L i is the length of tube i, which corresponds to the end of an overlapping portion. g θj s (X) is the total rotation between tube j and tube j -1 at X, which is obtained through the product of the partial rotation of each overlapping portion preceding X. For example, consider the total rotation between tube j and tube j -1 at X located on the third overlapping section, then g θj s (X) = g θj 1 (L 1 )g θj 2 (L 2 )g θj 3 (X).

A procedure similar to the single overlapping sections case leads to the equilibrium equations for general multiple nonoverlapping tubes in the generalized variable-strain coordinates. For example, consider the case of three non-overlapping tubes, then we obtain:

L 1 0 B T 1 F i 1 + Ad * g θ 2s F i 2 + Ad * g θ 3s F i 3 dX = 0 L 1 0 B T θ 2 1 (F i 2 + F i 3 ) dX = L 1 0 J T θ 2 1 ad * ξ 1 Ad * g θ 2s F i 2 + Ad * g θ 3s F i 3 dX L 1 0 B T θ 3 1 F i 3 = L 1 0 J T θ 3 1 ad * ξ 1 Ad * g θ 2s g θ 3s F i 3 dX+ J T θ 3 1 (L1) L 2 L 1 ad * ξ 2 Ad * g θ 3s F i 3 dX ------------------ L 2 L 1 B T 2 F i 2 + Ad * g θ 3s F i 3 dX = 0 L 2 L 1 B T θ 3 2 F i 3 dX = L 2 L 1 J T θ 3 2 ad * ξ 2 Ad * g θ 3s F i 3 dX ------------------ L 3 L 2 B T 3 F i 3 dX = 0 . (33)
Each dashed line separates an overlapping section. As with (26), equation (33) can be solved numerically for the unknown q = p 1 , p θ2 1 , p θ3 1 , p 2 , p θ3 2 , p 3 . We note that the model grows linearly with the number of tubes, and number of degrees of freedom.

Comparison with simulation in literature

The multiple overlapping tubes model is tested by comparing the equilibrium configuration obtained for a two tubes robot with a literature benchmark reported in [START_REF] Rucker | Equilibrium conformations of concentric-tube continuum robots[END_REF]. Consider two tubes of length 140 mm (outer) and 200 mm (inner) with constant pre-curvature and physical properties reported in Table I. We seek for the equilibrium configuration when the inner tube is rotated by α 2 = 180 • . The system domain is divided into two variable-strain sections. The generalized coordinates are q = (p 1 , p θ2 , p 2 ), and the equilibrium equations are a subset of (33). The resulting model takes the form of a nonlinear algebraic system of equations that can be solved numerically with the MATLAB © fsolve function. We chose a linear nonhomogeneous base for the strains involved in the first section, i.e.: k 1 , and θ 2 . Then, we have

B 1 =        
1 X 0 0 0 0 0 0 1 X 0 0 0 0 0 0 1 X 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

        , B θ2 =         1 X 0 0 0 0 0 0 0 0 0 0         , ( 34
) For the second section, since it is unloaded, we use the reference strain k * 2 as a basis for the strain k 2 , which gives p 2 = 1. As illustrated in [START_REF] Rucker | Equilibrium conformations of concentric-tube continuum robots[END_REF], the CTR has three equilibrium configurations. The solution to which the solver converges depends on the initial guess. Figure 2 illustrates the equilibrium configurations obtained with the variable-strain model. Comparing the Cartesian tip positions with respect to the standard model, we find a maximum error of 1.17mm, and a mean error of 7.8 * 10 -1 mm, which correspond to 0.58% and 0.39% of the backbone length, respectively. Furthermore, the inner tube's rotation with respect to the outer tube at the end of the overlapping section is equal to θ 2 (L 1 ) = ±85.5 • , which is only 1.3% bigger than the value reported in [START_REF] Rucker | Equilibrium conformations of concentric-tube continuum robots[END_REF] (±84.4 • ).

C. Complete CTR system

The CTR systems considered so far only include the tubes' portions that follow the insertion orifice, assuming that the rotation input α was applied there (see equation ( 27)). However, a complete CTR system is also composed of a straight portion that precedes the insertion point and terminates with an actuation system responsible for applying the input angle and insertion motion. To include this part, we place the 0 of the system domain X at the innermost tube base. Then, we divide the system domain into several sections corresponding to the different discontinuities. Those include a change in the number of overlapping tubes (as before), the insertion orifice, and, eventually, a reference strain jump. For instance, the straight portion is usually obtained by letting the tube be originally straight at the proximal end and curved toward the distal end. Finally, we constrain the straight portion by choosing only a torsional strain basis for these proximal sections. Note that this method applies to fully curved tubes as well, where external constraints enforce the straightness condition.

For example, let us consider the two tubes CTR system shown in Figure 3, where each tube has a straight proximal reference strain of length L js as indicated. There are five sections in this example indicated with an additional subscript in the following. Define D 1 , D 2 to be the distance between the insertion orifice and the outer and inner tube base, respectively, and ∆D = D 2 -D 1 . Then, the equilibrium equations for the generalized coordinates q = p θ2 1 , p 12 , p θ2 2 , p 13 , p θ2 3 , p 14 , p θ2 4 , p 25 are given by:

∆D 0 B T θ 2 1 F i 2 dX = J T θ 2 1 (∆D) L 2s D 2 ad * ξ 1 Ad * g θ 2s F i 2 dX+ ∆D+L 1 L 2s ad * ξ 1 Ad * g θ 2s F i 2 dX -------------------- D 2 ∆D B T 1 2 (F i 1 + F i 2 ) dX = 0 D 2 ∆D B T θ 2 2 F i 2 dX = J T θ 2 2 (D2) L 2s D 2 ad * ξ 1 Ad * g θ 2s F i 2 dX+ ∆D+L 1 L 2s ad * ξ 1 Ad * g θ 2s F i 2 dX -------------------- L 2s D 2 B T 1 3 F i 1 + Ad * g θ 2s F i 2 dX = 0 L 2s D 2 B T θ 2 3 F i 2 dX = L 2s D 2 J T θ 2 3 ad * ξ 1 Ad * g θ 2s F i 2 dX+ J T θ 2 3 (L2 s ) L 2 L 2s ad * ξ 1 Ad * g θ 2s F i 2 dX -------------------- ∆D+L 1 L 2s B T 1 4 F i 1 + Ad * g θ 2s F i 2 dX = 0 ∆D+L 1 L 2s B T θ 2 4 F i 2 dX = ∆D+L 1 L 2s J T θ 2 4 ad * ξ 1 Ad * g θ 2s F i 2 dX -------------------- L 2 ∆D+L 1 B T 2 5 F i 2 dX = 0.
(35) In developing the equilibrium equation ( 35), we have used the identity J T θ ad * ξ F i = 0 for the straight portion that precedes the insertion orifice. Furthermore, note that B 12 has only torsional components, and B θ2 2 (X), B θ2 3 (X) are equal to zero at X < ∆D and X < D 2 , respectively.

Considering the CTR systems of equation ( 33) and ( 35), we notice that each section follows the same structure of the single overlapping case (26).

Comparison with experiments in literature

The results of the complete system model ( 35) are compared with experimentally validated simulations available in the literature [START_REF] Rucker | Equilibrium conformations of concentric-tube continuum robots[END_REF]. Consider two tubes with a straight portion followed by a constant pre-curvature, as shown in Figure 3, and physical properties reported in Table I. We seek the equilibrium configuration in two sets of insertion and rotation conditions of the inner tube. In the first set, the full overlap case, D 2 = 208.5 mm, and α 2 goes from 0 • to 280 • in 20 • increments. In the second set, the partial overlap case, D 2 = 170.5 mm, and α 2 goes from 0 • to 200 • in 20 • increments. D 1 and α 1 are always equal to 93.5 mm and 0, respectively. We chose a quadratic basis for all the strains involved, except for the last section, where we use k * 2 . Figure 4 illustrates the equilibrium configurations obtained with the variable-strain model. Comparing the Cartesian tip positions with respect to the standard model, we find a maximum error of 4.8 * 10 -1 mm, 8.9 * 10 -1 mm and a mean error of 4.4 * 10 -1 mm, 7.2 * 10 -1 mm, respectively, for the full overlap and partial overlap case. The mean errors correspond to 0.14% and 0.23% of the backbone length.

D. Benefits of the PVS model for CTR

The proposed piecewise variable-strain approach for concentric tube robots presents several benefits. It provides the equilibrium equations as a minimum set of closed-form algebraic equations. The system is easy to handle for control and design optimization since all the boundary and continuity conditions are already intrinsically embedded in the equilibrium equations. Furthermore, extension to dynamics can be done without any change to the number and type of generalized coordinates. Finally, the PVS model is well-posed to accommodate additional DOFs and incorporate design variation of CTRs as well as generalization to other rod-driven soft robots. In the next section, we include the tubes' insertion and rotation motion as generalized coordinates of the system for the first time.

IV. SLIDING CTR MODEL

In the standard model approach, the tube rotations and insertions are actuation variables that dictate the boundary conditions of the governing ODEs. Actuator motion is thus simulated by simply varying the boundary conditions and solving the static model at sampled points in time [START_REF] Dupont | Design and control of concentric-tube robots[END_REF], [START_REF] Rucker | A geometrically exact model for externally loaded concentric-tube continuum robots[END_REF], thus producing what is known as a "quasi-static" simulation. While the Jacobian can be numerically computed along the robot [START_REF] Rucker | Computing jacobians and compliance matrices for externally loaded continuum robots[END_REF], the tube velocities are not automatically included in the standard model. The recent extension of the standard model to dynamics does describe velocities, but still assumes negligible actuator insertion velocity and acceleration [START_REF] Till | A dynamic model for concentric tube robots[END_REF]. Further, while the torque required to rotate the tubes is computable from the solution boundary values and the constitutive law, the standard model does not include information about actuator insertion forces. This section will extend the proposed PVS model for CTR to include the input motions as generalized coordinates and overcome the standard model's above limitations.

Before that, we establish a general result for the differential kinematics of rod's sections with a variable domain. Consider a rod section with material abscissa X ∈ [a(t), b(t)] and g(a) = I. The tip of the section g(b) can be formally represented by the integration of the linear time-varying differential equation ( 4) through its state-transition matrix g(b) = Φ(b, a). Then, the derivative of g(b) with respect to time has to include the variation of the domain boundaries in addition to the variation of the shape. This can be calculated as follows.

where the first can be considered a "growing from the base" term and the second a "growing from the tip" term. Note that internal cross-sections g(X) will vary only due to the proximal boundary a, while they are indifferent to the distal boundary b. Finally, the total velocity twist becomes:

η(b) = Ad -1 g(b) b a Ad g ξdX + ξ(b) ḃ -Ad -1 g(b) ξ(a) ȧ , ( 37
)
Observe the differences with equation ( 7) and ( 10) used so far.

To keep track of the potentially varying boundaries, in the following, we replace the notation g(b) with g([a, b]).

A. Complete CTR system with actuation forces and torques

Let us consider, without loss of generality, two nonoverlapping concentric tubes constrained to be straight before the insertion orifice, as shown in Figure 3. We removed the reference strain discontinuity here for simplicity, reducing the number of required sections to four. First, we define two material abscissas respectively for the outer tube X 1 = [0, L 1 ], and the inner tube X 2 = [0, L 2 ]. Overlapping cross-sections are then related by:

X 1 (∆D) = X 2 -∆D , (38)
where D 1 and D 2 are generalized coordinates now. In this example, the full set of generalized coordinates is q = D 1 , D 2 , α 1 , α 2 , p θ2 1 , p 12 , p θ2 2 , p 13 , p θ2 3 , p 24 . Contrary to section III-C, we fix the spatial frame on the insertion orifice. Then, the tubes' sections' kinematics can be expressed with equation (32) after being pre-multiplied by a translational homogeneous matrix g t , which depends on the insertions D 1 and D 2 . Let us focus on the overlapping section immediately after the insertion orifice. The kinematics there can be written as:

g 1s (X 1 ) = g t ([D 1 , 0])g 12 ([0, D 1 ])g 13 ([D 1 , X 1 ]) g 2s (X 2 ) = g t ([D 1 , 0])g 12 ([0, D 1 ])g 13 ([D 1 , X 1 (∆D)]) g θ2 1 ([0, ∆D])g θ2 2 ([∆D, D 2 ])g θ2 3 ([D 2 , X 2 ], X 1 (∆D)) , (39)
where, for the inner tube, X 1 is given by (38). The last term g θ2 3 ([D 2 , X 2 ], X 1 (∆D)) indicates the relative rotation between X 1 (∆D) of the outer tube and X 2 of the inner tube accumulated in the last section starting from D 2 . Using the general formula (37) and the discretizations [START_REF] Webster | Toward active cannulas: Miniature snake-like surgical robots[END_REF], and (20), we obtain the differential kinematics of the section.

η 1s (X 1 ) = Ad -1 g1 3 (X1) B α α1 -ξk1 (D + 1 ) Ḋ1 + S 12 (D + 1 ) ṗ12 + J 13 (X 1 ) ṗ13 η 2s (X 2 ) = Ad -1 g θ 2s (X2) η 1s (X 1 ) + B α α2 + J θ2 1 (∆D) ṗθ2 1 + J θ2 2 (D 2 ) ṗθ2 2 + J θ2 3 (X 2 ) ṗθ2 3 -ξk2 (X 2 )∆ Ḋ , (40)
where B α = [1 0 0 0 0 0] T , and ξk = [0 (k) y (k) z 1 0 0] T . Note that the variable upper boundary of g θ2 3 is due to X 1 (∆D). Thus, its variation has to be taken with X 2 fixed. Equation (40) provides the additional Jacobians corresponding to the input motions. Similar equations can be obtained for the other three sections. Projecting the differential equation ( 13) by d'Alembert's principle using these additional Jacobians yields the equilibrium equations for the insertion and rotation input forces τ D1 , τ D2 , and torques τ α1 , τ α2 .

τD 1 = EJ1 ξ T k 1 (ξ k 1 -ξ * k 1 ) D + 1 + EJ2 1 2 ξ T k 1 ξ k 1 D + 2 + L 1 +∆D D 2 ξ k 1 + ad ξ 1 ξk 1 T Adg θ 2s ξ * k 2 dX+ ξ T k 1 1 2 ξ k 1 -Adg θ 2s ξ * k 2 L1 + ∆D -- ξ T k 2 (ξ k 2 -ξ * k 2 ) L1 + ∆D + τD 2 = EJ2 ξ T k 1 1 2 ξ k 1 -Adg θ 2s ξ * k 2 D + 2 - L 1 +∆D D 2 ξ k 1 + ad ξ 1 ξk 1 T Adg θ 2s ξ * k 2 dX- ξ T k 1 1 2 ξ k 1 -Adg θ 2s ξ * k 2 L1 + ∆D -+ ξ T k 2 (ξ k 2 -ξ * k 2 ) L1 + ∆D + τα 1 = -τα 2 τα 2 = -B T α L 1 +∆D D 2 ad * ξ 1 Ad * g θ 2s F i 2 dX
(41) where ξ k = [0 (k) y (k) z 0 0 0] T . To obtain (41) we have also used the following boundary conditions.

ξT

F i1 (0) = -τ D1 , B T α F i1 (0) = -τ α1 ξT F i2 (0) = -τ D2 , B T α F i2 (0) = -τ α2 (42)
According to the Newton law of conservation of momentum, the input forces' and torques' sums have to be zero. It is an excellent verification to check if equation (41) satisfies these conditions. For the input torques, the third equation of (41) ensures that the actuation torques are equal and opposite. For what concerns the input forces, we have:

τ D1 + τ D2 = EJ 1 ξ T k1 ξ k1 -ξ * k1 D + 1 + EJ 2 ξ T k1 ξ k1 -Ad g θ 2s ξ * k2 D + 2 , τ D1 + τ D2 = ξ T k1 F i1 + Ad * g θ 2s F i2 D + 2 = 0 (43)
where the last identity, which states that the internal stress is null on the cross-section immediately after the orifice, is justified by the fact that no external forces act on the system from that point onward. Thus, also the input forces are equal and opposite as required.

B. Simulation tests

Some simulation tests are presented in this section to explore the behavior of the actuation forces and torques. We consider six configurations (three planar and three out-of-plane) of two concentric tubes with physical properties as for Table I. In all the cases, we consider linear non-homogeneous (or constant) pre-curvatures. The linear and non-homogeneous pre-curvature coefficients are varied to study the effect on the actuation forces and torques. In the out-of-plane cases, the inner tube is rotated by 90 • with respect to the outer tube. Table II reports the exact values of the non-homogeneous and linear coefficients and the calculated values of the actuation inputs.

The following observations can be made. For case a), the insertion motions D 1 and D 2 do not influence the actuation inputs value. In case b), while the insertion of the inner tube D 2 , and it is symmetric with respect to clockwise or counterclockwise rotations. In particular, the linear coefficients of e) and f) have an opposite contribution to the actuation force τ D2 . In case e), τ D2 decreases with the increment of the linear coefficient until it becomes negative, as reported in Table II. In case f), the opposite applies.

C. Benefits of the sliding-rod PVS model for CTR

The sliding-rod PVS model presented here proposes a growing (non-material) approach, which allows extending the PVS model to include the tubes' sliding motion without the need of calculating the unknown interaction forces. This new model may be used to control the CTR motion through the actuation force and torques instead of the insertion and rotation kinematics. Torque-controlled CTR can provide a new way to enhance elastic stability and improve interaction forces' control at the end-effector, currently a major concern in the design of CTR [START_REF] Ha | Optimizing tube precurvature to enhance the elastic stability of concentric tube robots[END_REF]. In the present form, the base force equations (41) can also be used to improve the fast kinematic controllers of [START_REF] Xu | Position control of concentric-tube continuum robots using a modified jacobian-based approach[END_REF], which is based on actuation load sensing.

V. CONCLUSIONS

In this paper, the recently proposed piecewise variablestrain approach for modeling highly deformable rods has been adjusted and applied to the case of concentric tube robots. The performances of the new approach have been compared with analytic, simulated, and experimental data available in the literature. Furthermore, the PVS approach has been further extended to include the tubes' insertion motion for the first time, which opens new unexplored possibilities for controlling these kinds of systems.

Future works include experimental validations of the slidingrod PVS model, the addition of external forces, and the extension to dynamics.

APPENDIX A ADJOINT REPRESENTATIONS

Ad g = R 0 3×3 rR R , Ad * g = R rR 0 3×3 R , ad ξ,η = k, w 0 3×3 u, v k, w ,ad * ξ,η = k, w u, v 0 3×3 k, w

Fig. 1 :

 1 Fig. 1: Kinematics of multiple non-overlapping sections. For each overlapping portion the kinematic equation (3) applies.

Fig. 2 :

 2 Fig. 2: Equilibrium configurations computed with the variable-strain model for two concentric tubes with constant pre-curvature. The inner tube is rotated by 180 • . The result of the variable-strain model match well with the one off the standard CTR model [21].

Fig. 3 :

 3 Fig. 3: Sketch of the simulated complete CTR system. Dimensions of the tubes are as follows: straight length L 1s = 93.5mm (outer) L 2s = 218.5mm (inner), curved length L 1c = 92.3mm (outer) L 2c = 85mm (inner)

Fig. 4 :

 4 Fig. 4: Equilibrium configurations computed with the variable-strain model for two concentric tubes with straight portion and constant pre-curvature. The inner tube is rotated from 0 • to 280 • in 20 • increments (right) and from 0 • to 200 • in 20 • increments (left). The result of the variable-strain model match well with the one of the standard CTR model [21].

  ξ(a)g(b) ȧ, and ∂Φ(b, a) ∂b db dt = g(b) ξ(b) ḃ,

TABLE I :

 I Physical Properties of Tubes

		Outer tube Inner tube
	Young's Modulus (GP a)	58	58
	Shear Modulus (GP a)	21.5	21.5
	Inner diameter (mm)	2.01	0
	Outer diameter (mm)	2.39	1.6
	Reference curvature (mm -1 )	0.0099	0.0138

TABLE II :

 II Sliding-rod test results , the force input absolute value |τ D | increase with the outer tube's insertion and curvature. In c), the force input absolute value |τ D | increase with the insertion of the outer tube and the retraction of the inner tube. For the out-of-plane cases, in general, the absolute value |τ D | increases with the rotation angle α 2 up to 90

			Outer tube	Inner tube	Actuation
		α 2	constant [rad/m]	linear [rad/m 2 ]	constant [rad/m]	linear [rad/m 2 ]	τ D 2 [N]	τα 2 [Nm]
	a)	0 •	10	0	30	0	0	0
	b)	0 •	30	-80	30	0	3.65	0
	c)	0 •	30	0	50	80	6.71	0
	d) 90 •	20	0	10	0	0.28	0.118
	e)	90 •	20	40	10	0	-1.45 0.122
	f)	90 •	20	0	10	40	1.89	0.129
	is still irrelevant

This paper was recommended for publication by Editor Cecilia Laschi upon evaluation of the Associate Editor and Reviewers' comments. This work was supported by Khalifa University of Science and Technology under Award No. RC2-2018-HEIC, CIRA-2020-074, RC1-2018-KUCARS and ADEK Award for Research Excellence (AARE-2018-105).


Sylvain Tourancheau
email: sylvain.tourancheau@univ-nantes.fr

Borje Andrén

Kjell Brunnström

Patrick Le Callet

Börje Andrén

Visual Annoyance and User Acceptance of LCD Motion-Blur

It has been recognized for some time now that LCD displays will introduce blur when showing moving objects or moving images. Common motion-blur measurement methods permit to picture the blurred profile of an edge moving with a constant velocity. A normalized blurred edge width is then measured for several grayto-gray transitions to give a motion-blur score of the display under test. However, these objective measurements are partly based on the behavior of the human visual system and it is an open question how well they correlate with subjective experience of observers. In this study, we develop a subjective experiment in order to assess the annoyance and the acceptance of motion-blur. Results are given and compare with measurements data.

Introduction

Motion blur is still an important issue for liquid crystal displays (LCD). Recently, efforts have been done in the measurement and the characterization of this artifact, in order to obtain the blurred profile of a sharp edge (defined by its gray-to-gray transition) moving on the display under test. The profile of this edge is referred to as the blurred edge profile and its width is referred to as the blurred edge width (BEW). This is generally measured between 10% and 90% of the amplitude of blurred edge profile. It has been shown from measurements as well as from analysis that, for a given gray-to-gray transition from gray level N i to gray level N j , BEW varies linearly with the velocity V [START_REF] Feng | LCD motion-blur analysis, perception, and reduction using synchronized backlight flashing[END_REF][2] [START_REF] Watson | The Spatial Standard Observer: A human vision model for display inspection[END_REF]. As BEW is expressed in space units and V in space over time units, the slope has a unity of time and is generally referred as the blurred edge time BET:

, , for each

i j i j i j BEW BET V N N = ⋅ → (Eq. 1)
These measurements will give some indication of the blurred edge profile that reaches the retina (if the effects of the optics of the eye are neglected). However, the human visual system is more complex and it needs to be further investigated when the perception of the edges will be different from the measurements or when they will be agree.

The goal of this study is to know how the blur is subjectively perceived in terms of acceptance and annoyance. As it has been shown [START_REF] Pan | Quantitative analysis of LCD motion blur and performance of existing approaches[END_REF][3] [START_REF] Tourancheau | LCD motion-blur estimation using different measurement methods[END_REF], the quantity of blur varies as a function of the gray-to-gray transition considered and as a function of the velocity of the moving object. Furthermore, it has been shown that the perception of LCD motion-blur varies with the contrast of the edge [START_REF] Tourancheau | Psychophysical study of LCD motion blur perception[END_REF]. The blur of low contrasted edges (i.e. gray-to-gray transitions) is less perceived than the blur of high contrasted edges. Regarding these results, subjective assessment of LCD motion-blur is needed in order to determine a threshold beyond which the quantity of perceived motion-blur is not acceptable.

Only few subjective studies [6][7] have been conducted to determine such acceptability levels. They consisted in the assessment of sequences on different LCDs and the results have been compared with global motion-blur index of the displays such as MPRT. In our study, we decided to assess the motion-blur perception on three liquid crystal displays. Motion-blur has been characterized on these LCDs, for 20 gray-to-gray transitions, through temporal measurements as described by Tourancheau et al. [START_REF] Tourancheau | LCD motion-blur estimation using different measurement methods[END_REF]. Moreover, a psychophysical study has been conducted on these three displays in order to determine precisely the quantity of perceived blur for each transition [START_REF] Tourancheau | Psychophysical study of LCD motion blur perception[END_REF].

As we measure the blurred edge width of 20 transitions, we choose to assess specific stimuli in which gray-to-gray transitions were controlled. Thus, we designed some synthetic sequences with gray objects moving on gray backgrounds, along with some sequences with scrolling texts on both uniform backgrounds and natural pictures. Some natural panoramic pictures have been assessed too. For each sequence, observers have to decide if the blur in the sequence is acceptable or not and are then asked to give a quality score (on a five-grade impairment scale). This would permit to explore the motion-blur perception, and to link visual annoyance and user acceptance with both physical measurements and psychophysical evaluation that have been done before. Experiments have been led at IRCCyN (France) for two displays, and at Acreo AB (Sweden) for the third one, with the same conditions and procedure. One of the goals of this study has been to help TCO Development to find a suitable requirement for motion blur.

In the following, the subjective assessment is described in Section 2. Results are given in Section 3 and discussed in Section 4.

Description of the Experiment

Displays under test

Stimuli were presented on three liquid crystal displays:

• DUT1: 26", 1920 x 1200@60Hz

• DUT2: 30", 2560 x 1600@60Hz

• DUT3: 37", 1920 x 1080@60Hz Blurred edge time (BET) has been measured on each display for the 20 gray-to-gray transitions between gray levels 0, 63, 127, 191, and 255. Measurements have been done according to the temporal method described in a previous work [START_REF] Tourancheau | LCD motion-blur estimation using different measurement methods[END_REF]. The minimum, maximum, and average values of BET are given in Table 1.

Viewing conditions

Experiments were performed in a psychophysics test room, with a surrounding illumination of 20 lux and a D65 chromaticity on the walls. Stimuli were generated in Matlab using the PsychToolbox.

All stimuli were presented in an area of 1920 x 1080 pixels at a refresh rate of 60 Hz. On displays with a larger resolution, the sequence was displayed on the center of the screen and the borders of the screen were filled with a mean luminance L mean corresponding to the geometric mean of the maximum and minimum luminance of the display, L mean = sqrt(L min L max ).

The viewing distance was adjusted in such a way that one pixel subtended 1.5 minute of arc of visual angle. In other words, the video sequence subtended an area of 48 x 27 visual degrees in the observer's visual field. This corresponded to a viewing distance of 2.08H, with H is the height of the picture (1080 pixels). Thus, stimuli had the same size on the retina whatever display which was tested. Examples of stimuli are given in Figure 1. In these 24 stimuli, the object of interest (train, text, or the whole picture) was moving at a constant speed.

Stimuli

Procedure

In each session, the observers assessed the quality of 96 sequences: 40 sequences with stimuli Train, 40 sequences with stimuli Text and 16 sequences with stimuli Picture. These 96 sequences were presented in a random order which was different for each observer.

After a sequence had been viewed by an observer, he/she was asked to decide if the blur in the sequence was acceptable or not. Then, he/she was also asked to rate the sequence according to a five-grade impairment scale: DUT2 have been tested at Acreo. 32 observers participated in the experiment.

Subjects

All subjects possessed normal or corrected-to-normal vision (visual acuity of 0.9 or better on both eyes). They were familiar with the procedure after the training sessions, the results of which have not been considered. All observers have participated in two sessions. Half of them started with one session whereas the other half started with the second session.

A rejection process has been done independently for each stimuli.

The linear correlation between each observer's results and the mean results was computed and the observer who obtain the lower correlation coefficient was rejected if the coefficient was lower than 0.8. Then the mean results were recalculated and the process was repeated until all observers' results are correlated to the mean results with a coefficient higher than or equal to 0.8.

In average, 3.3 observers have been rejected for stimuli with the synthetic train, 1.5 for stimuli with scrolling text, and 1.7 for panoramic pictures. In all cases, more than 25 observers were kept for each stimulus.

Results

For each of the 24 stimuli, we obtained the acceptance frequency and the mean opinion score (MOS) as a function of the speed.

Examples of results are given in Figure 2 with 95% confidence intervals for the MOS. For each stimulus, three thresholds can be computed:

• ThPerc: threshold of perception, i.e. the value for which MOS=4.5 (the stimulus goes from not perceptible to perceptible) • ThAnnoy: threshold of annoyance, i.e. the value for which MOS=3.5 (the stimulus goes from not annoying to slightly annoying) • ThAccept: threshold of acceptance, i.e. the value for which the acceptance frequency is 0.5

These three thresholds are represented by red dots in Figure 2. In a previous work [START_REF] Tourancheau | Psychophysical study of LCD motion blur perception[END_REF], it has been shown that observers are less sensitive to motion-blur for low-contrast transitions than for highcontrast transitions. Here, no particular tendency has been found

regarding the confidence intervals of the thresholds. Even for the Train Experiment, in which the gray-to-gray transitions were exactly known, the threshold of acceptance was not higher on low contrast as it could have been expected.

For this reason, we decided in this paper to only consider the average thresholds of each experiment. For each display, these speed values have been converted to blur values according to the measurements performed on the display. For synthetic stimuli, we took the mean between the BET of the rising transition and the BET of the falling transition considered, and for natural pictures we took the mean of all BET values. The blur value has been computed using Equation 1; results are presented in Table 2. (av. of 4) ThAccept. 12.6 10.3 12.4

Results of the train experiment were those with the most important variability intra-observers, particularly for the threshold of acceptance. When asked, some observers explained that for this synthetic sequence motion-blur added naturalness and was acceptable even for high speed, whereas some other observers decided it was not acceptable anymore once it was perceptible.

For the two other experiments, results were really more stable.

The most obvious evidence of this is the very close thresholds for DUT1 and DUT3 for which the pools of observers were the same at more than 65%.

For the three experiments, it can be noticed than threshold values were always lower on DUT2. This might be due to the difference of observers, and also to the fact that DUT2 were tested at Acreo (Sweden), when DUT1 and DUT3 were tested at IRCCyN-IVC (France). Despite of a lot of care, some small differences in viewing conditions could have happened. It might also be a display related difference, but this was not statistically significant.

Discussion

The experiment shows, see Table 2, that on average the blur was not perceived if its width was less than about 4.2 arc min = 0.07 deg. All three experiments were in agreement with this threshold, which was expected since it is a level where the blur becomes perceivable.

The annoyance threshold and the acceptance threshold were different for the three experiments. The train experiment had the highest threshold, then the picture experiment, and the text experiment showed the lowest values. It is interesting to notice that despite of variability between results, this tendency was always respected. For the text experiment, observers may have based their judgment of acceptance on the readability of the text, this could explain why this experiment obtained the lowest acceptance thresholds.

The threshold of acceptance was always higher than the threshold of annoyance. This means that observers can accept a slight annoyance due to motion-blur. On average, the threshold of acceptance corresponded to a MOS of 2.6. This value was found to be stable over the three experiments and the three displays.

As depicted in Finally, it has been reported that the maximum speed present in TV signal is around 30deg/s (except few cases). In order to be sure to not introduce blur at all, a perfect display would have a BET value < 6.66 ms. However, these high speeds are not very common and generally brief in time.

According to these typical speeds values, this would give the following design guidelines concerning LCD motion-blur: [START_REF] Vesa | Flat Panel Display Measurements (FPDM2)[END_REF] recommends to take the average BET value of 42 transitions between 7 gray levels perceptually equidistant from one to another. However, according to our previous work [START_REF] Tourancheau | Psychophysical study of LCD motion blur perception[END_REF], we would recommend to weight BET measurements according to the contrast/dynamic of the tested transitions, since motion-blur is less perceived on low-contrast edges than on high-contrast edges. In accordance with these latter results, a sharpness metric of moving edge has been described by Teunissen and Heynderickx [START_REF] Teunissen | Measuring motion blur in liquid crystal displays[END_REF] taking in account the contrast of the edge.

Conclusion

In this study, the visual annoyance and the user acceptance of LCD motion-blur has been assessed through subjective experiments. 44 sequences have been tested, with 8 different motion speeds. Stimuli were of three types: synthetic stimuli using measured gray-to-gray transitions, scrolling text sequences and natural scenes. These experiment have permitted to define three thresholds corresponding to the speeds for which blur becomes respectively perceptible, annoying and not acceptable anymore.

From the motion-blur measurements performed on the three DUT, it has been possible to convert these thresholds in blur quantities expressed in visual angles. These values have permitted to give some design guidelines concerning LCD motion-blur, based on the statistical characteristics of TV signal in terms of motion speeds.

  24 stimuli were presented to the observers; they can be classified in three different experiment: • Train experiment: 10 stimuli representing a synthetic train moving from left to right. Gray levels of train and background were varied from one stimulus to another. The height of the train was set to 2°. • Text experiment: 10 stimuli representing a text scrolling from right to left. Text was white or black and background was either a uniform area or a natural picture. The text was upper-case and its size has been set to 0.4° which is the average size of scrolling text in news channels. • Picture experiment: 4 natural panoramic pictures consisting in scenes shot in a 360° angle, these pictures were made scrolling continuously.

  annoying Observers had the possibility to replay the sequence as many times as they wanted.

Figure 1 :

 1 Figure 1: Examples of stimuli used in the experiment. First row: Train experiment (synthetic train), second row: Text experiment (scrolling text), two last rows: Picture experiment (panoramic scenes).

  DUT1 and DUT3 have been tested at IRCCyN. 21 observers have assessed motion-blur on both, 11 observers only on DUT1 and 9 only on DUT3. This led to a total of 32 observers on DUT1 and 30 on DUT3 with an overlap of around two-thirds. Among the 21 observers that have assessed motion-blur on both DUT1 and DUT3, 11 have seen DUT1 first and 10 have seen DUT3 first.

  They have been computed by the use of a fitting model on the subjective results. A logistic function has been used for acceptance and a spline has been used for the mean opinion score. Confidence interval of ThAccept is obtained from the fitting since this value is a parameter of the logistic function. For ThPerc and ThAnnoy, we made a fitting on the values of the confidence intervals of subjective data to obtain some bounds on the thresholds values. These bounds are represented by red triangles on Figure 2 (bottom row).

Figure 2 :

 2 Figure 2: Examples of results for two stimuli. Acceptance frequency (top) and Mean Opinion Score (bottom) for Picture Experiment with Seine picture (left) and for Train Experiment with black train on a white background (right) as a function of motion speed.

Table 1 : Minimum, maximum, and average values of BET (in ms) for the three displays under test.

 1

	Display min BET average BET max BET
	DUT1	13.6	17.8	21.6
	DUT2	12.2	14.1	15.8
	DUT3	13.7	15.3	17.5

Table 2 : Average threshold values for each display and each experiment. These values are blurred edge width (BEW) values, expressed in minutes of arc.

 2

			DUT1 DUT2 DUT3
	Train Exp.	ThPerc	5.3	4.0	4.4
	(av. of 10)	ThAnnoy. 10.9 ThAccept. 19.3	8.4 12.2	9.3 16.5
	Text Exp.	ThPerc.	4.1	3.6	4.1
	(av. of 10)	ThAnnoy. ThAccept.	7.1 9.8	6.6 8.7	7.1 9.6
	Picture Exp.	ThPerc.	4.3	3.7	4.4
		ThAnnoy.	7.9	7.2	8.1

Table 2

 2 , the average threshold of acceptance for the blur width 12' = 0.20 deg. In other words, blur was globally not acceptable if it was higher than 0.20 deg. According to Yoshida et al.[START_REF] Fujine | Real-life inhome viewing conditions for FPDs and statistical characteristics of broadcast video signal[END_REF] 10 deg/s is a typical speed limit for movies where there is no camera blur introduced: i.e. for motion speed less than 10 deg/s, observers are not expecting to see blur. According to Equation 1, the BET value of the display must be less than 20 ms for the display to not introduce blur.Concerning scrolling text on TV, the mean speed is 13.8 deg/s according to Yoshida et al.[START_REF] Fujine | Real-life inhome viewing conditions for FPDs and statistical characteristics of broadcast video signal[END_REF]. The results of the text experiment shows that the average acceptance threshold for this type of content is BEW = 9.4' = 0.157 deg. This means that BET must be less than 11 ms in order to display scrolling text with an acceptable quantity of blur. Subjects are more sensitive to blur in the text case according the experiment.

Acknowledgements

The authors would like to thank Romuald Pépion and Romain Cousseau who helped to conduct the experiments, as well as all the subjects who participated to the experiments.

The support from TCO Development, Intertek Semko and VINNOVA (The Swedish Governmental Agency for Innovation Systems) is hereby also gratefully acknowledged.