Riemannian geometric approach to optimal binocular rotation, pyramid based interpolation and bio-mimetic pan-tilt movement

Over the past several years, we have been studying the problem of optimally rotating a rigid sphere about its center, where the rotation is actuated by a triplet of external torques acting on the body. The control objective is to repeatedly direct a suitable radial vector, called the gaze vector, towards a stationary point target in IR3. The orientation of the sphere is constrained to lie in a suitable submanifold of SO(3). Historically, the constrained rotational movements were studied by physiologists in the nineteenth century, interested in eye and head movements. In this paper we revisit the gaze control problem, where two visual sensors, are tasked to simultaneously stare at a point target in the visual space. The target position changes discretely and the problem we consider is how to reorient the gaze directions of the sensors, along the optimal pathway of the human eyes, to the new location of the target. This is done by first solving an optimal control problem on the human binocular system. Next, we use these optimal control and show that a pan-tilt system can be controlled to follow the gaze trajectory of the human eye requiring a nonlinear static feedback of the pan and tilt angles and their derivatives. Our problem formulation uses a new Riemannian geometric description of the orientation space. The paper also introduces a new, pyramid based interpolation method, to implement the optimal controller.


Introduction
In this paper we consider the Binocular Sensory Control problem, where each sensor is tasked to mimic the movement dynamics of the human eye. Typically we assume that the centers of the sensors are fixed in space and the sensor gaze directions rotate to inspect point targets that are located in 3D. The gaze directions are always constrained to pass through a point and the goal of the sensing mechanism is to initially start with a target fixed in its *Correspondence: bijoy.ghosh@ttu.edu 1 Department of Mathematics and Statistics, Texas Tech University, Lubbock, Texas 79409-1042, USA 2 School of Automation Engineering, University of Electronic Science and Technology of China, Chengdu, Sichuan 610054, China view and to switch to an alternate target in the visual field in a unit interval of time, assumed to be [ 0, 1]. The paper analyses rotation for a binary pair of sensors and the goal is to compute optimal control over the chosen fixed time interval, extending two of our prior papers [1,2]. In our earlier research, binocular eye rotation has also been studied as a cascade of version and vergence eye movement applied to the two eyes separately (see [3] and [4]). This paper also extends [5], a recently published paper by the authors on Riemannian geometric formulation for the optimal control of binocular eye motion. The optimal control, we show, can be implemented by a pyramid based linear interpolation introduced in this paper. Ghosh  Anatomy of the eye is such that it is only able to rotate with three degrees of freedom [6,7], but is unable to translate [8]. The eye movement system is a relatively simple mechanical control system compared to other complex human movement systems [9]. Modeling the dynamics of monocular eye rotation has been an important goal in Neuroscience (see [10] for a short review article on how brain controls the eye movement) and Biomechanics [11]. Since the early half of the 19 th century [12,13], scientists have tried to create dynamic models in order to understand various eye movement trajectories (see [14] for some historical details). Starting from some of the initial papers of [15], the principles from geometry, for example as in [16][17][18], are central to many of the key questions in nonlinear systems theory (see [19]), applied to rotational dynamics. Specific to the eye movement control system, we would also like to refer to [20][21][22] and many references therein. For a single eye, optimization problems associated with gaze control [23], have been extended to optimal control problems studied by [24][25][26].
In the last few decades, there has been considerable research on exactly how the eye rotations are controlled. It was well known since the 19 th century by physiologists such as Helmholtz, Listing and Donders that when head is kept fixed, the axes of rotation is confined to a plane called the Listing's plane (see [27][28][29]). Questions arise as to how the Listing's constraint is satisfied by the 'motion controller' . Current literature seems to support the view that the constraint is met by active neural control in the brain (see [23,[29][30][31][32][33][34]) as opposed to an alternative view that the constraints are forced by mechanical properties of the eye plant using muscle pulleys [35,36]. In this paper and many of our earlier papers [5,25] geometric methods were used to study optimal control problems within the constrained space (called the configuration space) using a Riemannian formulation. We demonstrate that the synthesized control, to be viewed as actively generated by the neural circuit in the brain, continues to enforce for example the Listing's constraint (and an additional co-planarity constraint to be introduced later in this paper), even when implemented on the ambient space SO(3) 1 .
Geometric methods have a long history in the study of eye movement rotation (see [27,28,37,38]). Riemannian geometry (see [39,40]) has been introduced for monocular optimal control problems on the configuration spaces LIST in [6] and DOND in [24]. We have recently (see [5]) extended the Riemannian geometric formulation to binocular control problems, wherein a configuration space LBIN for the binocular eye pair is described as a subset of SO(3)×SO(3). In this paper, we first recall from [5] the construction of the optimal eye rotation controller 1 When states of the dynamical system satisfy the Listing's constraint at the initial time.
for the binocular eye pair. We then propose that these controllers can be implemented using a straightforward pyramid [41] based interpolation scheme, where the plan is to synthesize the control function as a convex combination of four corner points of a pyramid. Controlling the binocular system to each of the four corner points can be learnt and kept in the memory to be used subsequently as a lookup table. Finally we replace the eye-pair with a pair of mechanical visual sensors, capable of rotating, only along pan and tilt. We describe a simple control strategy, to track the optimal (human) gaze pathway, using the mechanical pan/tilt system (see [26] and Fig. 3 for a figure of generalized gimbal. For a pan/tilt system the axial rotation angle φ 3 is constrained to zero) 2 . Hence a mechanical pan/tilt system is able to follow human eye movement only up to its gaze direction but cannot follow the changing roll of the eye.
Next we briefly talk about the importance of a mechanical system following the human eye. In one line of research [42], a mechanical system, such as a humanoid robot, tries to generate human like behaviors in robots by recording and mimicking human motions in the execution of a particular task. For example, a human instructor can explain the process of assembling a device with several parts distributed over a table [43], and his associated gaze behaviors are learnt and encoded in dynamic gaze controllers implemented in a robot [44]. Human gaze behavior can also be replicated in a robot by directly programming the control laws into the robot. For example in [45] a neurophysiological model of human saccadic motion was implemented in a robot head where each eye had pan/tilt mobility. Although the emphasis of this paper is not on humanoid robots, we show how a pan/tilt system can imitate the human ocular movements upto gaze directions only, without an explicit reprogramming of its controllers, requiring only a nonlinear static feedback.
Finally, we outline the sections of this paper. In Section 2, we introduce notations that describe the axisangle parametrization [6,24,25] of unit quaternion. Section 3 describes a recently introduced Riemannian metric on the configuration space LBIN for the binocular pair of human eyes [5]. Using the Riemannian metric from Section 3, we write down the controlled Euler-Lagrange equation in Section 4. The control variables are the external torque vector functions applied to each of the two eyes and the problem we propose is to optimally control (using a suitably chosen cost function), the rotation of the eye-pair, in time interval [ 0, 1], so that the fixated target point in the visual space switches between two points. To solve for the optimal control, an associated Two Point Boundary Value Problem needs to be solved and in our prior papers [5,26], such a boundary value problem has been solved using a program called COMSOL [46]. In [5], many examples are showcased where the eye centers are kept horizontal, i.e. the head is kept straight relative to the torso. In Section 5, an example is introduced where the head is tilted sideways i.e., the separation vector joining the two eyes contained in the Listing's Plane 3 . In Section 6, we propose a pyramid based interpolation to approximate the control torque function required to transfer the binocular gaze point from an a priori chosen initial point to a final point. This final point is assumed to be contained within a cubic interval in the visual space.
The main point of the interpolation algorithm is that the optimal control does not have to be computed for every new final point. It is enough to have precomputed the optimal control for each of the eight corner points of the cube, eventually suitably recombining four of the corner points on a pyramid by convex combination. This procedure needs to be scaled up by discretizing the visual space into tiles of cubes and the associated optimal control functions, with final gaze points at the corner points of every cube, stored in the form of a look-up table. In Section 7, we introduce Tait-Bryan parametrization with the explicit goal of talking about pan/tilt/roll type of rotation [26] on SO(3). We compute the controlled Euler-Lagrange's equation with external torque as the associated control vector. We also repeat the same calculation when the roll component is forced to zero and we have a mechanical pan/tilt system with only two degrees of freedom for rotation.
In Section 8, we mainly show how an optimal control computed for the binocular system can be used as an input to a pair of pan-tilt system, where the goal is to match the gaze trajectories of the two systems point by point. Typically this step of biomimetic matching requires a dynamic feedback. Under the specific hypothesis that the axes of the binocular system of eyes satisfy Listing's condition, the feedback structure can be reduced to a nonlinear static feedback. We also show via simulation that re-parameterizing SO(3) from axis-angle to Tait-Bryan parametrization does not change the optimal gaze trajectories of the binocular system on the human eyes. In Section 9 we briefly discuss three main topics introduced in this paper. In the first topic we point out why it is impractical to solve a two point boundary value problem in real time, every time a control function is required to be computed. It is desirable, specifically for the purpose of inspecting targets in the visual space, to have the control functions pre-computed over a fixed initial and a discrete set of final target points. The required control is synthesized using interpolation based approximation. 3 Other tilts of the head including head roll has not been considered.
The second topic we point out is that for a certain class of binocular optimal control, the control function has a linear structure. This simplifies parameterizing and tabulating the control. The third topic we introduce is biomimetic control of a pair of mechanical visual sensors, and how pan/tilt rotation can be used to track the gaze direction of human eyes, using a feedback control. The feedback structure is particularly simple (nonlinear static) when eye axes are constrained to a Listing's plane. Finally, Section 10 concludes the paper emphasizing the importance of pyramid based interpolation and bio-mimetic pan/tilt rotation.

Notations and terminology
We start this section by introducing the axis-angle parametrization (see Fig. 1) of quaternions, where the notations are borrowed from [6] and [24,25]. This parameterization is particularly used in our study of human eye rotation control. In later part of this paper, we have also introduced Tait-Bryan parametrization [26] in order to talk about Bio-Mimetic pan/tilt movement.
Let us begin with the space of quaternions denoted by Q, see [47] and write each q ∈ Q as q 0 1 + q 1 i + q 2 j + q 3 k. Space of unit quaternions will be identified with the unit sphere S 3 , and can be written as where φ ∈[ 0, 2π] is an angle variable and n = (n 1 , n 2 , n 3 ) is a unit axis vector in R 3 . We denote by rot, the standard map from S 3 into SO(3) which maps the quaternion q to an orthogonal matrix that rotates a vector in R 3 around the axis of rotation n by a counterclockwise angle φ (see Fig. 1). It is easy to verify (see (3) in [25]) that The orthogonal matrix (2) can be associated with the orientation of a rotating rigid body as follows: Each column of (2) is a mutually orthogonal unit vector. We can associate the three column vectors to three body coordinates that describe the orientation.
The rotating rigid body, viz. the human eye, has a specific 'gaze direction' , a vector whose direction is what we propose to control. We use the convention that the gaze direction is given by the third column of the rotation matrix (2). We therefore have the following projection map, projecting the orientation matrix (2) to a gaze direction vector Typically our interest is to control the gaze vector (3) so that it is pointing towards a suitable point target (see Fig. 2 where a pair of eyes are pointing towards a target). As has been commented in Section 1, additional constraints on the quaternion q need to be imposed (such as the Listing's constraint given by q 3 = 0 or perhaps a more general Donders' constraint [26]), so that the constrained orientation matrix with a specific gaze direction is unique (see [48]). For pan/tilt rotation (see [26]), considered later in the paper, the Listing's constraint is replaced by Fick Gimbal constraint q 0 q 3 = −q 1 q 2 .
The pan/tilt system is a part of a generalized gimbal system (see Fig. 3) already introduced in [26], where the parametrization of the associated quaternion uses the Tait-Bryan angles φ 1 , φ 2 and φ 3 4 . In this paper, two pan/tilt systems (see Fig. 11) are controlled simultaneously, so that it bio-mimetically follows the gaze directions 4 For additional details on Tait-Bryan parameterization, please refer to [26]. The centers of the left and the right eyes are located respectively at (0, 0, 0) and (0, 1, 0). Note that the center of the inertial coordinate is assumed to coincide with the left eye center, upward-pointing axis is the positive x-axis while the forward-pointing axis is the positive z-axis. The y-axis joins the centers of the two eyes and the direction from left to right is chosen positive. Curly arrows (in blue) indicate the positive (clockwise, viewed from the axis center) rotation about each axis. The blue dot is the target, and the two red arrows are the eye-gaze directions (but not the roll) of the binocular human eyes, optimally saccading between two target points in the visual space.

Riemannian metric on the space LBIN
This section is essentially replicated from [5]. We begin by considering parametrization of a point in S 3 , as introduced in (1) and further describe the unit vector n, the axis of rotation, as n = (cos θ cos α sin θ cos α sin α), where θ ∈[ 0, π] and α ∈ − π 2 , π 2 . The parameterization (1), (4) together describes, what is known in [25], as the axis-angle parameterization of S 3 and SO(3) using the mapping 'rot' (essentially the three axis-angle parameters are θ, φ and α). In order for the orientation of a single eye to satisfy Listing's constraint, q 3 = 0, we impose α = 0, forcing the axis of rotation to always lie on the Listing's plane z = 0. This reduces the quaternion in (1) to the form (see [6]) We now introduce two such quaternions, one for the left eye q L and other for the right eye q R , described as and The pair q L , q R is thus an element of S 3 ×S 3 . Note that the gaze directions corresponding to the left and the right eye are given by Let us now assume that the left and the right eye has centers separated by the vector e = (x, y, z) = (0, 1, 0) as shown in Fig. 2. The figure shows that the centers of the left and the right eyes are located respectively at (0, 0, 0) and (0, 1, 0). The configuration space of the binocular system is now described by imposing that the vectors g L , g R and e are coplanar. Such a co-planarity condition will impose that the gaze directions of the left and the right eyes always meet at a point. Thus we have We denote by LBIN (L stands for Listing, BIN stands for Binocular), the subset of SO(3) × SO(3) where the orientation matrices separately obey Listing's Law and together the corresponding gaze directions satisfy the coplanarity condition (10). Equivalently, LBIN is a subset of LIST × LIST (see [6]) where the gaze directions of each component satisfy (10). Let ρ be the mapping described as where from (10) we have Note that the co-planarity condition (10) or (13) can change when the separation vector e changes. This will be the case when the two sensors are fixed but their centers are located at a different set of points.
A Riemannian metric on LBIN is easily induced from SO(3) × SO(3). We define elements g ij of the symmetric Riemannian matrix 5 G LB as 6 and compute 7 the Riemannian metric g given by where G is the corresponding Riemannian matrix of inner products (see [6] for a single eye and [5] for a binocular system of two eyes).

Euler-Lagrangian formulation of binocular eye movement
Since a Riemannian metric defines kinetic energy on the manifold, we use g in (15) to define the Lagrangian L = 1 2 g of the binocular system 8 . The controlled Euler-Lagrange equations are given by where μ ∈ θ L , φ L , θ R . It follows from [25] that (16) can be written as, where G LB is the Riemannian matrix, ∇ is the gradient of generalized torques τ μ . Further as [25] describes, we define the external torque vector T (a 6-vector), in the inertial coordinate to be, 5 The acronym LB stands for Listing and Binocular. 6 We define · to be the standard Euclidean inner product on the product space R 4 × R 4 . The inner product of ((x 1 , x 2 , x 3 , x 4 ), (x 5 , x 6 , x 7 , x 8 )) and ((y 1 , y 2 , y 3 , y 4 ), (y 5 , y 6 , y 7 , y 8 )) is given by x 1 y 1 + x 2 y 2 + · · · + x 8 y 8 . 7 The details of this computation is omitted, see [6]. 8 We assume that the potential energy function V is zero. 9 For a detailed computation of the M matrix and proof of the statement (19), we refer to [5].

Remark 1
The columns of the matrix M are called the Euler basis vectors, see [49], where T has been described as the resultant moment relative to the center of mass on the body. Now we setup our dynamical system for the binocular eye rotation by defining We require that the states go from some a priori agreed Z(0) to Z(1) while minimizing the control energy in a fixed interval of time where T is the vector of external torques given by Eq. (22). 10 We denote the costate variables by, and define the Hamiltonian as, Using the Hamilton's equations [17], the system is now obtained. Using Eqs. (17), (18), and (20), one can recast (24) as Finally using the Pontryagin's Maximum Principle, the expressions for optimal external torques (see [6]) are obtained: which we write symbolically as The control torques can now be eliminated from the state space system (25) and we obtain the following dynamical system 10 The superscript L and R are for left eye and right eye respectively. The subscripts refer to the x-, y-and z-axes.
Since we know only the initial and the final value of Z, we have a two-point boundary value problem (BVP). The resulting problem is solved using COMSOL Multiphysics program (see [46]) 11 . The computed Z and variables are plugged in (28) to obtain the optimal vector T, which is denoted by T BVP .
The optimal vector has been computed for a large number of examples in [5] where the gaze of the binocular eye pair moves between target points going from 'left to right' , 'bottom to top' and 'near to far' in the visual field. The corresponding visual trajectories and the optimal torque trajectories have been plotted. The eye centers are assumed fixed and located on the (horizontal) y−axis. This would be the case when the vector separating the two eyes is parallel to the ground. In the next section we discuss one simulation when the head is fixed but tilted.

Optimal eye movement when head is fixed and tilted
The main message of this section is the following. Although in Section 4 we had said that the optimal external torque controls are computed by solving a BVP using COMSOL, it is often possible to approximate the external torque function using a linear function. This had been the case in various examples discussed in [5] 12 , where separation vector between the two eyes are horizontal. In this section we consider a case when the head is tilted (see Fig. 4) and the eye separation vector is not horizontal. We demonstrate via simulation that 'a linear approximation is still good' .

Example 1: Eye separation vector is not horizontal, i.e., the head is tilted
Let us assume a binocular system as in Fig. 4 where the left eye is centered at (0, 0, 0) and the right eye is centered at (1, 1, 0). The separation vector e is given by (1, 1, 0) and the co-planarity of g L , g R , and e is appropriately defined. The angle variable φ R needs to be redefined in comparison to (13) and this changes the Riemannian matrix G LB in (15). We now consider the two point boundary value problem sketched in Section 4 assuming the initial and final values of the angles θ L , φ L , θ R are respectively π 4 , π 8 , π 6 and π 3 , π 5 , π 4 . The time derivatives of the angle variables are assumed to be 0 at the initial and final times. As indicated in Section 4, we use COMSOL to calculate the optimal trajectories for Fig. 5) and using the co-planarity condition (not explicitly described The centers of the two eyes are not horizontal, because the head is tilted. This changes the co-planarity condition (10) and subsequently the Riemannian matrix G LB in (14) in this section), the variable φ R (t),φ R (t). The optimal trajectories of states and costates are plugged in (28) to solve for the optimal external torque vector Fig. 6). As was noted in [5], the graph of T BVP (t) appears linear from the figure and we approximate this function by T LIN (t) = T BVP (0) [ 1 − 2t] where t ∈[ 0, 1]. The approximate linear function is plotted in Fig. 6 (using the symbol +) and T LIN (t) and T BVP (t) appear indistinguishable.
Let us now consider an initial value problem by combining Eqs. (16) and (18) to obtain d dt Initial conditions for (30) are chosen by setting θ The initial value problem is now solved and the the results are plotted in Fig. 5 using the symbol +. Once again, these trajectories mimic very close to the optimal trajectories. 6 Approximating the optimal external torque function using pyramid based interpolation In the last Section 5 we observed via simulation that for the binocular human eye movement, the optimal external torque functions have a linear graph. We have demonstrated this 'linearity' when the eye separation vector lies on the frontal plane, which in our construction also happens to be the Listing's plane of the two eyes. In this section we go back to Fig. 2 and consider the binocular system where the eye centers are at (0, 0, 0) and (0, 1, 0). Assume that initially the eyes are focused at the point C (see Fig. 7). The task is to move the gaze optimally to another final point D. However the exact location of point D is uncertain and we shall assume that it could be any where on or inside a cube CU with corners at D 1 , D 2 , · · · , D 8 . The main point of this section is to illustrate the following: If we denote by T CD the optimal control of transferring the gaze from C to D, then T CD can be approximated by a linear pyramid based interpolation using four of the optimal control functions from T CD i , i = 1, .., 8. First of all we write CU as a union of six pyramids PY i , i = 1, · · · , 6 where any two pyramids may intersect only at their surfaces. We now ascertain the membership of D in one of the six pyramids and call this pyramid PY. The corner points of the pyramid PY are four of the eight points D 1 , D 2 , · · · , D 8 . Assume without any loss of generality that these corner points are D 1 , D 2 , D 3 , D 4 . We claim that T CD is approximated by a convex combination of T CD i , i = 1, 2, 3, 4. Via simulation Fig. 5 The right eye is located on the Listing's plane of the left eye centered at (0, 0, 0). The right eye is centered at (1, 1, 0). The separation vector between the two eyes are on the coronal or the frontal plane but not parallel to the y-axis (horizontal axis). The angles θ L , φ L , θ R shifts from π 4 , π 8 , π 6 to ( π 3 , π 5 , π 4 ). The figure shows the generalized angles and the generalized velocities when the input external torques are computed by solving the Boundary Value Problem and compared with the corresponding trajectories when the input external torques are approximated by a linear function we show the effectiveness of the proposed approximation. We now digress a little and describe how a cube can be written as a union of six pyramids.

Pyramid construction
We start our discussion with a cube CU of edge distance 1 unit (see Fig. 8). Assume that one corner of the cube is at (a, b, c) chosen in such a way that D lies in the interior of the unit cube. The other corners are at the following 7 points: (a+1, b, c), (a, b+1, c), (a, b, c+1), (a+1, b+1, c), We claim that CU can be decomposed into 6 pyramids PY 1 , . . . , PY 6 such that 6. Each pyramid that we construct in I R 3 has 4 vertices. Let us now write down the following unordered   are (a, b, c), (a + 1, b, c), (a, b + 1, c), (a, b, c + 1), (a + 1, b + 1, c), (a + 1, b, c + 1),  (a, b + 1, c + 1), (a + 1, b + 1, c + 1) triplets of elements, using tokensā,b,c, a + 1, b + 1 and c + 1 13 . Let us now generate the following array of 6 rows of 4 triplets as follows: Array I: 1. (ā,b,c), (a + 1,b,c), (a + 1, b + 1,c), ā,c,b), (a + 1,c,b), (a + 1, c + 1,b), Note that the above array of triplets follow a pattern. In the first column of the array, the tokensā,b,c are ordered in each of its possible 6 combinations. The second column of the array is same as the first column of the array except that the first element of a triplet in the second column is obtained by incrementing the first element of the corresponding triplet in the first column.
Likewise, the third column of the array is same as the second column of the array except that the second element of a triplet in the third column is obtained by incrementing the second element of the corresponding triplet in the second column.
Finally, the fourth column of the array is same as the third column of the array except that the third element of a triplet in the fourth column is obtained by incrementing the third element of the corresponding triplet in the third column.
For every triplet in Array I, the unordered triplets of tokens are now reordered by removing the overbar and writing a or a + 1 to the left of b or b + 1 which is written to the left of c or c + 1. We perform this operation for each triplet of tokens in the above Array I term by term. We now get the following Array II of 6 rows of 4 ordered triplets of elements.
For every triplet in Array II, we treat each element of a triplet as coordinates of a point in I R 3 . We define PY i , i = 1, · · · , 6 as the i th pyramid generated by the four points of the i th row, to be treated as corner points. The 6 pyramids are generated one for each row of Array II. We now state and prove the following theorem. Proof First of all we show that for (x, y, z) to lie in the pyramid PY i , i = 1, · · · , 6 the coordinates (x − a, y − b, z − c) needs to satisfy certain inequality. This is now described for each of the six pyramids.
All points in PY 2 has the property that For PY 3 the convex combination is All points in PY 3 has the property that For PY 4 the convex combination is All points in PY 4 has the property that For PY 5 the convex combination is All points in PY 5 has the property that For PY 6 the convex combination is All points in PY 6 has the property that Since the coordinates x − a, y − b and z − c are distinct, by assumption, it would follow that one and only one of the six inequalities described in (31), (32), (33), (34), (35), (36) would be satisfied. − (a, b, c) are not distinct, then D is at the surface of more than one of the six pyramids. The details are quite evident and is omitted.

Approximating external torque via interpolation
Recall that T CD (t) is the external torque of transferring the gaze point of the binocular system from C to D in I R 3 . Typically T CD (t) is calculated by solving a boundary value problem (using COMSOL) as outlined in Sections 4 and 5. In this section we argue that T CD (t) need not be calculated once for every D in I R 3 . Using pyramid based linear interpolation, we can approximate T CD (t) from T CD i (t), i = 1, · · · , 8 by writing where μ i , i = 1, 2, 3, 4 are unique coefficients obtained as follows. If D belongs to one of the six pyramids with corner points at D 1 , D 2 , D 3 and D 4 , (assumed without any loss of generality), then μ i -s are obtained uniquely by solving Our final step in the interpolation based approximation is valid when the graph of the functions T CD i , i = 1, · · · , 8 are linear (see Section 5). Using the 'linearity' of the external torque function we modify (37) as follows.

Example 2: (Final gaze point is in the interior of the corresponding pyramid)
In this example we consider the binocular system where the left eye is at (0, 0, 0) and the right eye is at (0, 1, 0). The two eyes are initially gazing at a target located at (1, 2, 1). The goal is to optimally move the gaze to a new target located at (5.7, −0.5, 2.3). By solving the boundary value problem, we compute and plot the optimal trajectories of the generalized angles ( Fig. 9(a), solid lines) and generalized velocities ( Fig. 9(b), solid lines). The optimal gaze trajectory is shown in (Fig. 9(c), blue, dotted line). Finally the optimal external torques are graphed in (Fig. 9(d), solid lines). Note that these graphs are straight lines.
We now use the approximation formula (38) to obtain a linear interpolation of the external torque function T L INT and T R INT . In (Fig. 9(d), dotted lines) the graphs of T L INT and T R INT are plotted. Finally we use (30) and the linear interpolation T L INT and T R INT to solve for the generalized angles and velocities (see (Fig. 9(a), dotted lines) and ( Fig. 9(b), dotted lines)). We plot the gaze trajectory in (Fig. 9(c), red, dotted lines) for the interpolated external torque function.

Example 3: (Final gaze point is on the surface of the corresponding pyramid)
In this example, we have the binocular system with sensors located as in Example 2. The goal is to optimally move the gaze from an initial gaze at (1, 2, 1) to a final gaze point at (6, 0, 4). We proceed by solving the boundary value problem, compute and plot the optimal trajectories of the generalized angles ( Fig. 10(a), solid lines) and generalized velocities (Fig. 10(b), solid lines). The optimal gaze trajectory is shown in (Fig. 10(c), blue, dotted line). Finally the optimal external torques are graphed in (Fig. 10(d), solid lines). Once again, as in Example 2, these graphs are straight lines.
Pyramid based interpolation is now used as in Example 2. It turns out that the final point (6, 0, 4) is on the surface on the cube (39) and specifically on the surface of the pyramid with corner points at  (Fig. 10(a), dotted lines) and ( Fig. 10(b), dotted lines)). We plot the gaze trajectory in (Fig. 10(c), red, dotted lines) for the interpolated external torque function.

Remark 3 The main difference between examples 2 and 3 is the location of the final target point. In the former, the target point is in the interior of the associated pyramid and in the latter, it is on the surface. The main point to observe, and it is the essence of the two simulations, is evident from Figs. 9(c) and 10(c). 'The trajectory of the gaze point does not deviate appreciably even when the input control in the form of external torque is computed using a linear interpolation from four corner points of an associated pyramid' .
Thus control can be synthesized using a lookup table and one does not have to execute a more demanding COMSOL, in real time. (38). Recall from (28) that

Remark 4 Our final remark of this section is about the approximation formula
Let us define the symbols . We can rewrite (40) symbolically as where the matrix W is a 3×6 matrix whose entries depend only on the initial gaze point C and not on the final gaze point D (see [5]). The vector λ, on the other hand, is a 1 × 3 vector whose entries depend both on C and D. We can rewrite (38), using (41), as follows.
Note in (42) that when C is fixed, then W is fixed but λ i depends on D i , the corner points of the pyramid. The points λ i , i = 1, · · · , 8 in I R 3 form the corner points of a cuboid, which can be calculated apriori. We make the following statemnt about the interpolation algorithm. For every cube CU (with corner points D i ) there is a cuboid CB (with corner points λ i ) that can be precalculated and stored. This cuboid is dependent on the initial point C.

Rotation dynamics with Tait-Bryan parameterization
In this section we introduce Tait-Bryan (TB) parameterization [50,51] and make connection with Fick Gimbals and Pan-Tilt rotation. Our introduction will be brief and we will refer the readers to a previous paper [26], see also Figs. 3 and 11. In the TB parameterization there are three angle variables φ 1 , φ 2 and φ 3 where φ 1 is the angle of rotation about axis 1 (see Fig. 3), φ 2 is the angle of rotation about axis 2 rotated by φ 1 -rotation about axis 1. Rotations about axis 1 and axis 2 are respectively called Pan and Tilt. Finally φ 3 is the angle of Axial rotation about axis 3 rotated by the previous two rotations. The three angle variables φ 1 , φ 2 and φ 3 completely parameterizes the orientation space  1) and (4), for the axis-angle parameterization, we now have the following unit quaternion for the TB parameterization (See (9) in [26]).
where M TB is the M-matrix for the TB parameterization described as follows (already reported in page 323, [25]) The Eqs. (46)- (48), describe rotation dynamics on SO(3) with the external torque vector T as the control. As in Eqs. (20)-(29), we can setup a dynamical system for the monocular unconstrained eye rotation on SO(3). We can also require the states to go from a priori agreed initial state to a final state in a unit interval of time minimizing the control energy (21). Repeating the steps in Section 4, an optimal external torque T can now be computed using COMSOL to solve the associated two point boundary value problem BVP (see also [24][25][26]).
In order to describe human eye-rotation, the orientation space is not the unrestricted SO(3) but a submanifold LIST (see [25] and Fig. 15) of SO(3). We now parameterize LIST using TB parametrization.

Example 4: LIST using Tait-Bryan parametrization
Using the quaternion parametrization (43), it follows that the Listing's constraint q 3 = 0 is described by Let us now define Substituting the Listing's constraint (50) into the quaternion (43), we obtain the following parametrization of LIST in the TB parametrization ⎡ Note that in the parametrization (52) of LIST, only the pan (φ 1 ) and tilt (φ 2 ) are used. Construction of a left invariant Riemannian metric on LIST in the TB parametrization would now be standard [6] and the details are not elaborated here.

Example 5: Euler-Lagrange equation for the Pan-Tilt system
We now impose the Fick Gimbal constraint φ 3 = 0 into the TB quaternion (43). Recall that for a Fick Gimbal, the axial rotation angle φ 3 is permanently frozen to zero. The Fick Gimbal quaternion is described as follows.
We now write and Computing the inner products we write and From (56) and (57)  We choose the potential energy V to be zero and the kinetic energy and we write the Lagrangian as The Euler-Lagrange's Eq. (16) for the Pan-Tilt system is now written as 1 4φ Finally we would like to relate the external torque vector T = (T 1 , T 2 , T 3 ) T to the generalized torque τ = (τ φ 1 , τ φ 2 ) T . This is done by first calculating the M-matrix M FG .
where • denotes quaternion multiplication. Substituting the parameters of the Fick Gimbal quaternion (53) into (59) and carrying out the algebraic computation we write ⎛ Solving (60) It would now follow from a relation similar to (18) that τ φ 1 = T 2 and τ φ 2 = T 1 . The Euler-Lagrange's equation for the Pan-Tilt system (58) now reduces to For the purpose of the next section, we shall define from (61) and write the following Pan-Tilt dynamics 8 Biomimetic Pan-Tilt movement following the optimal gaze trajectories of the human binocular system Our goal in this section is to control the Pan-Tilt dynamics (62) so that the gaze of the pan tilt system matches the gaze of the human ocular system. Since the ocular dynamics on SO(3) satisfies (46)- (48), roughly speaking, we need to match the angle variables φ 1 and φ 2 in (62) with the same variables in (46) and (47). Note that the angle φ 3 is frozen permanently to zero for the pan-tilt system, and therefore cannot match the axial movements of the human eye i.e., φ 3 in (46). Let us rewrite part of the dynamical system (46)-(48) as follows where M TB is defined in (49) and we define and H TB = sec 2 φ 2 0 sec 2 φ 2 sin φ 2 0 1 0 .
In order to match the Fick Gimbal dynamics (62) with the SO(3) dynamics (63) the right hand sides have to match and we get the following where T is any external torque input to the ocular dynamics on SO(3). The Eq. (64) can be viewed as a torquetransformer, which transforms the external torque T for the SO(3)-system to T FG , the torque input to the Pan-Tilt-system. The equality (64) would ensure that the angle variables φ 1 and φ 2 of the two systems match, provided they have the same initial conditions. In order to implement the torque-transformer, we will needφ 3 which unfortunately is not available to the Pan-Tilt system. Hence, we will need aφ 3 -generator to be implemented with the following equation: + 4 sec 2 φ 2 sin φ 2 , 0, 4 sec 2 φ 2 M T TB T. Implementation of the torque-transformer and theφ 3 generator has been sketched in Fig. 12.
If T = T * is chosen in such a way that the angle variables φ 1 , φ 2 and φ 3 satisfy the Listing's constraint (50) one can computeφ 3 explicitly aṡ As shown in Fig. 13 that in this case theφ 3 generator is implemented by a pan-tilt feedback requiring variables φ 1 , φ 2 and their derivatives.

Remark 5
In this remark we would like to put into perspective how we have synthesized and implemented the Fig. 12 The torque-transformer is a static device. Theφ 3 -generator is a dynamical system requiring T and ρ as input. ρ = φ 1 ,φ 1 , φ 2 ,φ 2 T , consists of Pan, Tilt angles and their derivatives. T is the external torque input to the unconstrained monocular plant. T FG is the external torque input to the mechanical pan-tilt device that follows the human eye gaze movement Fig. 13 The torque-transformer is a static device. Theφ 3 -generator is a static system requiring ρ as input. T * is chosen such that the states φ 1 , φ 2 , φ 3 evolve on LIST. ρ and T FG are as described in Fig. 12 external torque optimal control. First of all note that in Section 4 we have proposed to synthesize the optimal control for the binocular system where the axes of human eye rotation satisfy Listing's law and the eye gazes always remain coplanar together with its separation vector. The associated space LBIN is parametrized using axis-angle parameters as a subset of SO(3) × SO(3) (see Fig. 14 We consider an example 16 from [5] wherein we have a binocular system with eye centers located as in Fig. 2. The goal is to shift the gaze point from (7,2,4) to (3,2,8)  . LBIN is to be viewed as a parametrization of the binocular system using a dynamical system in the form of (17), (18) while minimizing a cost function (21). The optimal control was obtained in [5] and its graph was plotted (see Fig. 13 in [5]). The problem was solved on LBIN using axis-angle parametrization. The optimal control vector T, from (22) is now separated between the left and the right eye and applied separately to (46)-(48), a dynamical system on SO(3) described using Tait-Bryan parametrization. Solving the initial value problem (46)- (48), with the input torque obtained from T as indicated, we show that the gaze direction vector of each eye follow a trajectory identical to what was obtained as solution to the optimal control problem in [5]. In Fig. 16(a) we show that the Fig. 15 LIST is shown as a subset of the ambient space SO(3). LIST is to be viewed as a parametrization of the monocular system gaze points of the two binocular systems are identical. We also show that the solution to the initial value problem actually evolved in LIST (see Fig. 16(b)) although it was solved on SO(3). It shows that the optimal controller is able to maintain the LBIN constraints (specifically the Listing's constraint on any of the two eyes) although these constraints are not mechanically imposed on the eye pair.

Discussion
Synthesis of the optimal control function, encountered in eye and head rotation problems, as introduced in many of our earlier papers [5,6,[24][25][26]53] requires solving a Two Point Boundary Value Problem. These boundary value problems are typically computationally intensive (see text books [54] and [55]), and require the framework of a recursive solver, for example one based on the algorithm of [56]. In our research, we had started with MATLAB Toolbox [6] to solve boundary value problems, Pseudospectral Methods [4,[57][58][59] and subsequently used the COMSOL Program already indicated in section 4. In any application, wherein a binocular gaze directing robot would be used for target localization, inspection and reaching 17 , it is not very convenient to have to solve a boundary value problem for every chosen pair of initial and final points. This is because the iterative solution to the two point boundary value problem, required to solve for the optimal control, suffers from singularity issues and do not always converge. Often parameters in the COMSOL Program, and in the KNITRO solver, that we have used for pseudospectral methods, require tweaking. The configuration space LBIN may need to be re-parameterized to avoid singularities. The point we would like to put forward in this paper is that it is convenient to precompute the control function for a fixed initial point and over a set of discretely chosen final points, viz. the corner points of a cube. In real time, the gaze directing controller uses the precomputed control functions as a lookup table and computes the required control for its own action using a pyramid based interpolation. Solutions to the Boundary Value Problems, once solved, does not have to be recomputed. When the centers of the eyes in the binocular system remain stationary and on the Listing's plane, we have illustrated in this paper and also in [5] that the optimal control function has a Linear structure, described in (42). The angular velocities of the two eyes increase to a peak value before reducing back to zero, as is typically observed in eye saccades [33]. The optimal control to the binocular system is entirely parameterized by its initial value at t = 0. It has not been shown rigorously, all the possible  16 Solving an initial value problem with using the external torques obtained Both sensors are on the y-axis located at (0, 0, 0) and (0, 1, 0). The gaze point shifts from (7,2,4) to (3,2,8) conditions under which the linear structure of the optimal control is maintained 18 .
Either for the purpose of optimal control or otherwise, it is important to mimic human movements and is a subject of research in social robotics in particular [61]. Specific to human gaze and eye movement, it is possible to control a pair of mechanical pan/tilt system even though the underlying configuration spaces of the human eye and the mechanical pan/tilt system are different 19 . In spite of the difference, the part of the two configuration spaces, that control the gaze/pointing direction can be identified and controlled (as depicted in Fig. 12). The main point is that the controller T FG for the mechanical system does not have to be recomputed but instead it can be generated from the control input to the human eye, by a suitable dynamic feedback (64), called theφ 3 -generator. When the axes of the eye movement are restricted to the Listing's plane, the feedback structure can be simplified to a static feedback (65) and the implementation details are shown in Fig. 13.

Conclusion
This paper has revisited the optimal binocular gaze control problem recently introduced in the Riemannian setting by the authors. It is shown that the optimal control function can be approximated by a pyramid based interpolation scheme, hence does not need to be solved in real time. Some level of discretization for the final target point will be allowed. The paper also introduces a new Biomimetic pan-tilt rotation control, where a pair of mechanical eyes are tasked to follow the human binocular gaze trajectory. The paper shows that if the optimal control to the human binocular system is known, the same can be used with a mechanical pair of eyes, with an appropriate nonlinear static feedback.