Generation of realistic saddle trajectories from captured horseback motion

Hippotherapy, riding a horse in the context of rehabilitation, is a medical treatment that successfully has been employed in various fields, e.g. for improving locomotion performance of patients with movement disorders. Robotic systems enable the application of hippotherapy in clinical environments with additional benefits, like adjustable speed and high repeatability. Fundamental for a therapy outcome equivalent to classical hippotherapy is that the trajectory of the robotic system is as realistic as possible. This paper introduces a method for the synthesis of horseback motions using motion capture data of various horse gaits. Based on complete gait cycles, an analytical, time-continuous and periodic motion description is constructed. Measured 3D marker positions are reconstructed with a mean error not exceeding 8.6 mm. If the motion capture data of several gait cycles are considered, a more generalized trajectory is generated. An adjustable time dilatation parameter enables the adaptation of the generated motion according to the physical abilities of the patient or the capabilities of the robotic system. This method allows for motion synthesis with arbitrary time span and time resolution, generating realistic trajectories effectively applicable to robotic systems for riding simulation in general and robotic hippotherapy in particular.


Introduction
Riding a horse is a complex, cognitive and physically demanding movement with significant similarities to human walking [15,36]. The rhythmic, dynamic interaction between horse and rider generates sensory stimuli with positive rehabilitational effects on the human body. Therapeutic horse riding, also known as hippotherapy and a type of equine assisted therapy (EAT), is therefore employed as a medical treatment, e.g. for persons with movement disorders. It has successfully been applied to improve gross motor function [35], postural control [41], balance [12], muscle symmetry [5], trunk/head stability and functional reach B J. Ziegler jakob.ziegler@jku.at 1 Institute of Robotics, Johannes Kepler University Linz, Altenberger Strasse 69, 4040 Linz, Austria [33], as well as for post-stroke gait training [4]. To support traditional therapy methods, attention is given more and more to robotic rehabilitation approaches. Concerning EAT, the development of robotic hippotherapy systems is motivated by the thereby gained possibility to facilitate hippotherapy in clinical environments or generally avoiding the necessity of a real horse. It has been shown that the therapy outcome of robotic systems is comparable to classical EAT, where real horses are involved [13,22,26]. In contrast to classical hippotherapy, robotic hippotherapy allows for arbitrary adjustable cadence or speed (within the system limits) with high repeatability and precision of the executed movements. In [21] the authors mention the association between greater motion variability and decreased therapeutic effect. Robotic systems are not dependent on the temperament of a real horse, and do therefore not contain the risk of unpredictable motions. Additionally, they enable the reproduction of various motion patterns related to different horses. Depending on their body type, horses show differing suitability for diverse applications, as reported in [29]. Furthermore, robotic hippotherapy has the potential to increase the accessibility of this treatment for the patients. A review on the benefits of classical hippotherapy and horse riding simulation can be found in [16].
To the best knowledge of the authors there is no study investigating the reproduction of horseback motion using a comprehensive set of optical markers for the transfer to a robotic system. Nevertheless, there are some studies focusing on the description and assessment of horse motions and the development of riding simulators. The authors of [24] developed a robotic riding simulator where the riding motion is modeled as periodic sinusoidal signals that have dominant up-and-down motions with two, three and four beats. In [14] the dynamics of the horseback during motion is captured by an inertial measurement unit rigidly attached to the saddle. The measured signals are then transformed to input signals for a hydraulic Gough-Stewart platform by digital filtering. A motion capture system was used in [7] to measure the position of reflective markers attached to horse, saddle and rider while the horse was walking on a treadmill. In [27] the authors introduce a robotic horseback-riding simulator consisting of a base motion platform with 4 degrees of freedom and an additional revolving system with 2 degrees of freedom mounted on top. The saddle motion trajectories were based on motion capture data and expressed as the sum of two sine functions. Mapping real movements to robotic systems is also investigated in other fields of research, especially regarding human gait, exoskeletons and gait rehabilitation [9,39] or biped locomotion [28].
We assume that, in order to improve the sensation of real horseback riding and to achieve a therapeutic success equivalent to classical hippotherapy, the robotic horse trajectories should be as close to the movements of real horses as possible. To this end, this paper introduces a motion synthesis method for the horseback, modeled as a rigid body, during typical horse gaits based on motion capture data. Improvements to a previously proposed approach in [42], already in use in robotic hippotherapy [20], are presented. An analytical, time-continuous description of the motion of a virtual saddle, or a hypothetical therapy sit, is constructed so to best reproduce the captured horseback movements. Horseback motion can be synthesized with arbitrary time span and time resolution, resulting in smooth, cyclic and realistic trajectories, immediately applicable to robotic systems. With an adjustable time dilatation parameter the speed of the generated movement can easily be adjusted to the current physical abilities of the patient or the capabilities of the robotic system. Note that this work does not investigate the transitions between the individual horse gaits or possible speed-dependent changes of the motion patterns, which might be part of future research. The presented approach consists of the three main steps: 1. The positional variance of the individual markers with respect to each other during motion is investigated. A subset of markers with a relative movement below a certain thresh- old, and therefore more appropriate to describe the recorded movement as rigid body motion, is selected for further calculations. 2. Based on the measured marker positions of the selected markers of the horse in steady standing position, a kinematic model of the horseback is generated. This rigid body model of the horseback represents the saddle of a riding simulator or a robotic hippotherapy system imitating the motion of a real horse. 3. The parameters of a cyclic trajectory, describing the rigid body motion of the horseback model, are identified in an optimization process using recorded marker data of complete gait cycles. This method allows for generating general high-fidelity saddle motions since the algorithm allows for considering data of an arbitrary number of complete gait cycles simultaneously.

Motion data
Basis for the analysis of movements is the application of a (kinematic) model to motion capture data, which enables a description of the underlying motion, i.e. the time evolution of position, orientation and, if joints are defined, joint angles of the model. Motion capture thereby refers to the recording of the motion of general objects, including humans and animals, using primarily vision-based or inertial systems. It is applied in various fields such as computer animation [25], robotics [31], human activity recognition [1] or biomechanics [11]. Particularly, measuring the 3D position of passive (purely reflective) markers is a frequently used method [10]. Horse motions were already analyzed on the basis of motion capture data [7], acceleration data [14,36] or photographs [19].
Within the framework of this paper, the authors had access to motion capture data recorded with passive markers. Two horses, Lambada and Doolittle, were equipped with sets of reflective markers whose 3D positions were measured during the three typical horse gaits walk, trot and gallop (Fig. 1). Walk and trot measurements were recorded for both horses, while gallop motion was only available for one horse. In total, the raw marker data of five movement trials together with snapshots of each horse in standing position were provided to the authors.
Pre-processing As none of the data sets contains more than one complete gait cycle, all motion data sets are cropped to the samples of one single movement cycle per gait. The marker data of the standing horses are used to define the reference pose, i.e. the shape of the respective horseback, or a virtual saddle. To this end, the x-axis of the horse coordinate system is defined to point anteriorly (forwards), the y-axes to the left side of the horse body and the z-axis superiorly (upwards). Its origin is positioned at the centroid of the x-and zcomponents of the marker positions of the reference pose. The markers nearest to the horse spine are aligned with the x-axis, as they approximately represent the center line of the horse body. Additionally, as the horses were not necessarily moving parallel to the axes of the inertial coordinate system and the generated trajectories should be applicable on a stationary robotic system, a trend correction of the recorded marker positions is applied. To this end, the position of the horseback during the motion cycle is calculated according to the least squares rigid motion approach described in the appendix. For every position component the trend is now calculated as the linear progress from the first to the last position and is subtracted from the according measurement data. This results in the horseback position deviating almost around zero, or another fixed point by adding some constant offset. The latter is important as the executing robotic system is assumed to be stationary and its workspace to be limited.

Marker selection
Modeling the horseback as a rigid body contradicts the fact that the markers attached to the horse are subject to relative movement with respect to each other. Nevertheless, as a relatively large part of the horses is covered by the markers it is assumed that some markers show more relative movement than others. The markers on the back and the middle of the flank of the horse are suspected to have smaller positional variance than the markers near the fore and hind leg, which would be in accordance with the traditional placement of a horse saddle. Therefore, a marker selection based on the positional variance of the markers is performed.
Points on a rigid body are not moving with respect to each other over time, which means the marker to marker distances are constant and the standard deviation of the distances is zero. By calculating the standard deviation of the marker to marker distances for each marker over the samples of a recorded horse gait, information about the 'rigidness' of marker regions can be gained. For the n = 1, . . . , N samples of a gait measurement the distance from marker i to marker j is calculated as with m i,n and m j,n being the 3D positions of the corresponding markers. The standard deviation of the distance of a marker pair concerning all N measurement samples is denoted as whered i,j is the mean distance between markers i and j during the considered measurement. All s i,j are summarized in a symmetric matrix S, where one row or column represents the relative movement from a specific marker to all markers. In order to visualize potential marker clusters, the rows and columns of S are ordered based on the median standard deviation in marker distances of the individual markers. This results in a new matrixS, which now has rows and columns with increasing median values. With the method of spectral clustering [30], the marker set can be segmented into groups. This method was already used to identify human body segments and positions of the connecting joints from motion capture data [23]. Nevertheless, in our case all markers with an average standard deviation of the marker to marker distances below a certain threshold, i.e. the median value with respect to all markers, were selected. Figure 2 shows the original matrix S of standard deviations and its version with sorted rows and columnsS for one of the data sets. As can be seen in Fig. 2(a), there is a specific subset of markers (towards the lower right corner of the matrixS) with a small standard deviation of marker-to-marker distances within this group but with a high standard deviation with respect to the rest. This subset corresponds to the small marker group on or near the hind leg of the horse. The selected markers correspond to the upper left quarter of this matrix. Figure 3 shows the selected markers for each of the data sets. Some markers are missing, since there was no corresponding positional information during the specific measurement trial.

Marker kinematics
The position of the coordinate frame of the reference pose, in the following referred to as horse frame, is described by the vector I r H ∈ R 3 represented in the inertial frame, while its orientation relative to the inertial frame is described by the rotation matrix R IH ∈ SO(3). Notice that the left subscript of a vector denotes the coordinate frame in which the vector is represented, I for inertial and H for horse. In summary, the pose of the horse is represented by the homogeneous transformation matrix Denote with H p j the (constant) position of the j th horse marker measured in the horse frame as determined from the reference pose, i.e. the standing horse. Its position measured in the inertial frame, denoted I m j , is then determined by which can be written in terms of homogeneous coordinates The pose of the horseback model is described by the generalized coordinates, summarized in the vector q = [q 1 , q 2 , . . . , q 6 ] T ∈ V 6 , which comprises the position vector, as well as the Cardan angles α, β and γ (rotation sequence X-Y'-Z") used to parameterize the rotation matrix R IH . I r H (q) depends on the variables describing the position of the horseback [q 1 , q 2 , q 3 ] T while the rotation matrix R IH (q) depends on the rotational degrees of freedom [q 4 , q 5 , q 6 ] T with q usually being time-dependent, i.e. q = q(t). Thus, the homogeneous transformation matrix is T IH = T IH (q(t)).

Time-continuous description of trajectories
Due to the cyclic nature of locomotion movements, the repetitive motion patterns of the individual horse gaits are described as periodic signals expressed by truncated Fourier series.
The Fourier series of order l for the generalized coordinate q i (t), with cycle duration T , is where a i = [a i,0 , . . . , a i,l ] T and b i = [b i,1 , . . . , b i,l ] T are the 2l + 1 Fourier coefficients. If a parameterization as described is chosen to represent the motion of the presented model, it follows that the generalized coordinates q = q(t, π ) are now dependent on trajectory parameters denoted with π = [a T 1 , . . . , a T 6 , b T 1 , . . . , b T 6 ] T . Typically, locomotion movements are quasi-periodic, which means that the individual gait cycles are similar but not identical. Shape and period duration are subject to variations and therefore, e.g. considering gait, single steps or strides do not necessarily take the same amount of time. The proposed method enables the possibility to find an optimal trajectory based on an arbitrary number of complete motion cycles. As each cycle consists of n = 1, . . . , N c samples depending on the sample rate of the motion capture system, a time normalization, based on a unified time-scale, is applied. A defined normalized period duration ofT = 1 results, for every cycle c, in a sample vectort = [t 0 = 0, . . . ,t Nc =T ]. Thus, the generalized coordinates are nowq =q(t, π ).
Time-scaling Extending Eq. (6) by a dilatation parameter ξ , with 0 < ξ ≤ 1, easily allows for an adaption of the original trajectory concerning the movement speed according to the current physical abilities of the patient or to the capabilities of the robotic system. The timescaled version of the approximated trajectory and its first two derivatives can be calculated as As can be noted in Eq. (8) and Eq. (9), the velocity decreases linearly, while the acceleration decreases quadratically with ξ , which might be important considering that the robotic system has to be capable of producing the required actuator forces, e.g. due to the varying body mass of different patients.

Optimization
Problem statement Goal of the optimization is to identify the motion of the model from marker data of a specific gait measurement. This amounts to find an optimal set of trajectory parameters π * reproducing the measured 3D Initialization In order to improve the performance of the optimization in terms of speed and robustness, it is advisable to provide a reasonable initial guess of the optimization variables. The method used in our previous work on this topic, as described in [42], can be applied to this end. Thereby, the position and orientation of the horseback, also modeled as a rigid body, is determined for every sample of a measured motion cycle by finding the rigid body transformation (R IH,n , I t H,n ) ∈ SE ( The problem in (11), also called weighted Procrustes problem [32], has a closed form solution [3,18], which is briefly summarized in the appendix. The extracted trajectory of the horseback model in terms of position and orientation can now be approximated with a Fourier series and used as initial guess for the optimization.

Results
The CasADi framework [2] was used to implement the nonlinear optimization problem (10) which was solved with IPOPT [40]. All algorithms and tests were executed on a computer equipped with an i7 CPU (3.60 GHz) and 8 GB RAM. The functionality of the method was validated with artificial marker data without the effects of noise, soft tissue artifacts or other inaccuracies. To this end, the reference pose was used together with a Fourier series parameterization of a prescribed motion to generate a synthetic 3D marker set. For validation the trajectory parameters were recovered with the proposed optimization method based on the artificial data, which could be achieved with arbitrary precision. Provided that the number of parameters describing the motion is sufficiently large, the accuracy of the optimization result depends only on the used time resolution of the simulation.
Concerning measurement data, an order of l = 6 was used for the motion parameterization with Fourier series, which results in a total of 78 optimization variables. Figure 4 shows the optimized trajectories of the horseback position and orientation for the three gaits walk, trot and gallop of one horse, the according velocities and angular velocities are computed as v = ẋ,ẏ,ż T and and are shown in Fig. 5. The maximum errors between estimated and measured marker positions were 21.9 mm for the walking, 22.8 mm for the trotting and 30.0 mm for the galloping motion. The respective mean error values were 5.5 mm (walk), 6.8 mm (trot) and 8.6 mm (gallop), respectively. As formulated in Eq. (10), an arbitrary number of complete gait cycles can be considered at once in order to calculate a more general version of the respective movement. To demonstrate this functionality, the walking gait cycles of the two horses were used to generate a common horseback motion that best fits both of the individual marker data sets (Fig. 6). This resulted in an increased marker position error with a maximum value of 85.5 mm and a mean value of 33.5 mm. For a typical motion sequence, consisting of about 180 samples acquired with a frequency of 200 Hz, an optimization of the trajectory parameters took 0.78 s on average. Considering multiple motion cycles simultaneously in order to obtain a more general solution of the motion pattern, the computational time increased with a rate slightly greater than linear with respect to the number of motion sequences.

Discussion
In this study a method for the synthesis of time-continuous rigid body motions of a horseback is presented. An analytical, time-continuous and periodic description of the movement allows for the generation of realistic trajectories with arbitrary time span and time resolution, effectively executable by robotic systems for riding simulation and robotic hippotherapy. Motion capture data in the form of 3D marker positions measured during three typical horse gaits, walk, trot and gallop, of two horses without rider are used to calculate the trajectories of position and orientation of a rigid body model of the horseback.
The reference pose of the model, defining the shape of the horseback, is obtained from marker data of the horse in steady standing position. As the horseback is modeled as a rigid body, e.g. representing a virtual saddle or the seat of a therapy system, a marker selection is performed in order to identify the markers most significantly describing the measured motion of the horseback.
Fourier series are used to parameterize the positional and rotational degrees of freedom of the model, which inherently results in time-continuous and periodic trajectories. The Fourier coefficients are optimized on the basis of the measured marker positions of complete gait cycles, where a higher number of cycles of the same movement results in a more general solution. As there are notable differences in the motion patterns of different horses, this can be used to determine a typical, characteristic (horseback) motion pattern of a specific horse gait. Additionally, the implementation of a time dilatation parameter allows for a convenient option to adjust the generated motion according to the current physical abilities of the patient or the capabilities of the robotic system.
The method was validated using artificially generated marker data, which could be reproduced with arbitrary precision only depending on the used time resolution of the simulation.   Prerequisite to do so is a motion parameterization with a sufficiently large number of parameters that can be adapted. After the validation the proposed approach was tested with measurement data. For all used motion data sets the marker positions could be recovered within a mean error not exceeding 8.6 mm. Although no publications on soft tissue artifacts (STA), i.e. artifacts caused by the relative movement between bones and soft tissues like muscles and skin, concerning the back of a horse could be found, this value is within the range of typical skin displacements in the limbs [37,38] and in accordance with values reported in the literature [6].
It is well known that horses without a rider have measurably different motion patterns compared to the riding situation [8]. Nevertheless, we do not assume that this significantly affects the therapy outcome in the context of robotic hippotherapy, although it might be worthwhile to do additional investigations in this direction. We are not aware of an existing publication investigating the therapeutic effect of trajectories from horses with versus without a rider. Note that transitions between the different horse gaits or speed-dependent alterations of the motion patterns are not addressed within this work. Moreover, motion capture data of only one complete motion cycle per horse and gait were available, additional investigations on this topic would benefit from larger data sets. Future research will focus on a novel kinematic design of a robotic hippotherapy system and the generation of generalized horseback motions based on more comprehensive motion capture data.
Authors' contributions All authors were fully involved in and have made substantial contributions to this work, including conception and implementation of the methods, as well as preparation and revision of the manuscript.
Funding Open access funding provided by Johannes Kepler University Linz. This work has been supported by the COMET-K2 Center for Symbiotic Mechatronics of the Linz Center of Mechatronics (LCM) funded by the Austrian federal government and the federal state of Upper Austria.

Conflict of interest
The authors declare that they do not have any financial or personal relationship with other people or organizations that could inappropriately influence this work.
Publisher's Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

Appendix: least squares rigid motion
As described in Sect. 5 the problem stated in Eq. (11) amounts to find a rotation matrix R IH,n and a translation vector I t H,n for data sample n, mapping the m = 1, . . . , M marker coordinates { H p 1 , . . . , H p M } of the reference pose to the marker positions { I m 1,n , . . . , I m M,n } measured in the inertial coordinate system. Note that the subscripts I and H are omitted hereinafter for better readability. To briefly summarize an approach for a closed form solution using singular value decomposition [17,34], let us assume the rotation to be fixed, and denote The optimal translation can be found by searching for the roots of the derivative of L with respect to t n : wherep is the weighted centroid of the reference marker positions andm n is the weighted centroid of the measured marker positions of frame n, and rearranging Eq. (14), the optimal translation is given as w m x T m R T n R n x m − y T m,n R n x m − x T m R T n y m,n + y T m,n y m,n .
As rotation matrices are orthogonal and therefore R T n R n = I, where I is the identity matrix, and y T m,n R n x m is a scalar, which implies y T m,n R n x m = (y T m,n R n x m ) T = x T m R T n y m,n , the optimization problem becomes R n = arg min w m y T m,n R n x m = Tr WY T n R n X = Tr R n XWY T n = Tr (R n S n ) , (19) where W = diag[w 1 , . . . , w M ] contains the weights and S n = XWY T n is the weighted covariance matrix of the centered vectors x m and y m,n . By applying the singular value decomposition on S n , Eq. (19) can be reformulated as Tr (R n S n ) = Tr R n U n n V T n = Tr n V T n R n U n .
As n is a diagonal matrix with non-negative singular values and V T n R n U n is orthogonal, Eq. (20) is maximized when V T n R n U n is the identity matrix, that is, The optimal orthogonal matrix is calculated with Eq. (22), which covers rotations as well as reflections, where det V n U T n = +1 for rotations and det V n U T n = −1 for reflections. Therefore, if R n is found to have a determinant of −1, this means that the global maximum is not achievable with a rotation (but only with a reflection). Since any reflection matrix can be converted to a rotation matrix by changing the sign of a row and since the singular values in n are arranged in decreasing order, the next best thing to the global maximum is n diag (1, 1, . . . , −1). Thus, if we restrict R n to be a rotation matrix, it is required that V T n R n U n = diag(1, 1, . . . , 1) if det V n U T n = 1 and V T n R n U n = diag(1, 1, . . . , −1) if det V n U T n = −1. Therefore, the optimal rotation matrix can finally be calculated as