Using Screw Theory for the Kinematics and Statics Design of the Turmell-Bot: a Cable-Driven, Reconfigurable and Compliant Ankle Rehabilitation Parallel Robot

1. Introduction

The ankle is a complex joint with high injury incidence. There are growing research interest in robotics to assist ankle rehabilitation. Our purpose is to create a device that can aid in the movement of the human ankle joint by restoring its range of motion. We focused on the reconfigurable kinematics, statics, and adjustable design.

In this work, we have taken a human-centered approach [1,2,3,4,5,6,7,8,9,10,11]. We propose an adjustable device for different foot sizes and be able to switch between limbs. The design is for the patients in lying down positions. Our approach is to mimicry the tendons involved in ankle movements while avoiding any forces that could pull the foot away from the ankle. We propose reciprocal tendons that apply pressure against the plantar surface of the foot, which will induce ankle joint integration forces to generate motion.

1.1. Ankle Rehabilitation Robotics

There are many state-of-the-art articles, literature reviews, and design considerations related to ankle rehabilitation robotics. For example, we found designs based on platforms in [12], state-of-the-art techniques in [13], and robot-assisted therapy in [14]. We also found the state-of-the-art robot-assisted ankle neurorehabilitation and challenges reviewed in [15], a mechanical design review in [16], robot-assisted techniques reviewed in [17], and a home-based review in [18]. Furthermore, we found anatomy-based designs regarding ankle-foot robotics considered in [19], and research progress in [20]. In general, most of the efforts are focused on parallel mechanisms, which have advantages over serial counterparts[21] and can be analyzed by using screw theory [22,23], as is explained in modern robotics [24,25].

1.2. Cable-Driven Robotics

Cable-driven parallel mechanisms have lightweight designs acting in tension against gravity or mutual tension configurations. They are reviewed in [26] and [27] and studied in [28]. A state-of-the-art is discussed in [29]. Cable-driven robots were applied to human arms [30], with Bowden guides [31] in lower-limbs [32]. The idea of attaching the shank and the foot was discussed in [33]. Additionally, wires [34] and a cable-driven ankle rehabilitation robot, were used in [35]. Such mechanisms are related to tendon-driven mechanisms [36], musculoskeletal robots, and artificial muscles.

2. Materials and Methods

First, we study the ankle anatomy in [37] and the International Society of Biomechanics (ISB)[38]. Then, we use the open-source digital model (z-anatomy) [39] to identify the bones, insertions, tendons, and muscles related to the ankle movement. We introduced the dimensions and geometrical model in a Jupyter Notebook using a SageMath 10 kernel [40] in [41]. We provided the CAD model in [42] and the sources in GitHub [43]. We identified the tendons involved in ankle movement acting as an over-actuated bio-mechanism (two degrees-of-freedom and four tendon groups). Then, we choose antagonistic tendon groups. We identified the pressure forces from the platform against the plantar surface of the foot. Then, we drew a schematic with the tension forces between the base attached to the shank and a platform attached to the foot. We adapted the robot dimensions from [44], the proportions from [45], and the statistical data from [46].

We start the design by analyzing a simple cable-driven two-axis serial chain, which helps us to understand the reciprocal products between the cable and axis screw representations. We create artificial data from the ankle model using the forward kinematics computed from the “product of exponentials” (PoE). The data serve us to validate how to get the ankle model in practice from trajectory measurements received from a modification of the Turmell-Meter system [47]. When synthesizing the robot, we enlarged the device size to prevent cable contact with the body from the platform and base. We analyzed the cable-driven antagonistic actuation in the perpendicular axes model and use it for the initial calibration. Afterward, we propose a method for the robot reconfiguration, the workspace, and the static simulation using MuJoCo 2.3.1 [48] to validate the stability by changing the axis and tendon positions. Finally, we proposed a CAD model using SolidWorks 2022.

2.1. Robot Inspired on the Analysis of the Ankle Joint

The model has two rotational joints and four tendons involved in ankle movement. First, we analyzed the tension forces and then the compression forces. Finally, we propose a simplified schematic for the robot.

2.1.1. Ankle-Foot Tendons

The two-axis representation and the tendon insertions in the bones at the mid-foot are in Figure 1. The Figure 1 illustrates the main tendons from the lateral-anterior view of the right foot. As we show in Figure 1, the fifth metatarsal has two insertions, the peroneus tertius, with components for dorsiflexion and eversion, and the peroneus brevis, for eversion and plantarflexion. The calcaneus has the insertion of the Achilles tendon for inversion and plantarflexion. In Figure 1, we show the medial-bottom view of the foot, and the insertion points are between the first metatarsal and the medial cuneiform. The tibialis anterior contributes to inversion and dorsiflexion, the peroneus longus contributes to plantar flexion and eversion, and the tibialis posterior contributes to inversion and plantar flexion.

2.1.2. Foot Compression Forces

Some requirements for ankle rehabilitation are security and comfort. The foot supports high compression forces in some regions of the plantar surface. In contrast, the dorsal skin is thin, soft, and less tolerant of compression forces. We propose a device that mainly applies pressure at the plantar surface to integrate the ankle joint and generate motion. From the bottom view of the foot, we note that the pressure points are involved in the subtalar and talocrural axes motion. The Figure 2 shows the axes projection and the triangle vertices of the transversal plane. The three main points are the tubercle of the calcaneus and the first and fifth metatarsal heads. The pressure points are the projection of the triangular-shaped dome, limited by the lateral, medial-longitudinal arches, and the anterior transverse arch. These contact points transmit a compression force to the platform.

2.1.3. Robot Based on the Ankle-Foot Model

We designed a reconfigurable cable-driven robot to actuate on synergy with the ankle tendons. Cables from the base pull the platform attached to the foot. The cable anchor points have a relocation mechanism to equilibrate the rotation forces on the ankle joint. We align the centers to the ankle axes’ intersection on the transverse plane. We show the schematic design in the Figure 3. In Figure 3, we sketch the approximated schematic and show the tendon directions regarding the subtalar (ST) and the talocrural (TC) axis. We represent the foot platform with a circle with a radius $r_p$, centered on $P{M}_0$, and on the same plane of the anchor points $a{p}_1$, $a{p}_2$, $a{p}_3$, and $a{p}_4$. The base is a circle with a radius $r_b$ in the same plane as $P_O$, and the anchor points $a{b}_1$, $a{b}_2$, $a{b}_3$, and $a{b}_4$. We simplified the plantar surface of the foot with three contact points $f{s}_1$, $f{s}_2$, and $f{s}_3$. Finally, $T_1$, $T_2$, $T_3$ and $T_4$ are tendons. The serial kinematic chain is RR, starting at the origin $P_O$, followed by a rotational joint at $r_1$ on the TC axis, followed by the $r_1-r_2$ link to a rotational joint at $r_2$ on the ST axis. A link from $r_2$ to $P{M}_0$ connects the platform with the ankle. The tendons from $a{p}_1$, $a{p}_2$, $a{p}_3$, and $a{p}_4$ on the platform to the corresponding $a{b}_1$, $a{b}_2$, $a{b}_3$, and $a{b}_4$ on the platform complete the parallel closed-loop structure. Each human ankle has different axis positions; thus, we propose reconfigurable cable attachments. In Figure 3, the green arrows illustrate the concept of reconfiguration for the right foot by displacing the anchor endpoints $a{p}_1$, $a{p}_2$, $a{p}_3$, and $a{p}_4$ on the platform and the anchors $a{b}_1$, $a{b}_2$, $a{b}_3$, and $a{b}_4$ on the base. We designed a mechanism for centering pivot points $P_{pb}$ and $P_{pp}$.

We computed the center points from the intersection of the talocrural and subtalar axes projected on the transverse plane parallel to the base and the platform. Such pivot points are not the same as $P_0$ and $P{M}_0$, which pertain to the position sensors reference system.

2.1.4. Dimensions and Initial Configuration

In this section, we estimate the initial robot size and configuration. We use the ankle model and measurements from [44], human proportions from [45] and statistical data from [46]. We show the proportions in Figure 4. The TC axis is dominant, and we use its statistical value for the cable-body collision analysis. We simplified the body as a sphere centered in half of the ankle’s most medial point (MMP) and most lateral point (MLP); between the two malleoli. We compare this dimension with a leg model in Figure 4b, illustrating the distance between the ankle and PM on the platform. Additionally, we show that the cables must not be in contact with the foot or the malleolus. The sphere represents the radius between the ankle and the plantar surface of the foot.

2.2. Ankle Kinematic Model

For the ankle kinematics representation, we used data from a two-axis model representation of the ankle joint explained in [47]. There, we described a system for platform pose capture and a method for model approximation using circle fitting, the Figure 5 shows the concept.

The points $M_1$, $M_2$, and $r_1$ are references on the talocrural axis. Similarly, $N_1$, $N_2$, and $r_2$ are in the subtalar axis. We used it to compute the following vectors and matrices:

$$\begin{array}{cc}\hfill \widehat{r_1}=r_1-P_O,& \widehat{r_2}=r_2-P_O\end{array}$$

(1)

$$\begin{array}{cc}\hfill {\widehat{\omega }}_1=\frac{M_2-M_1}{\parallel M_2-M_1\parallel },& {\widehat{\omega }}_2=\frac{N_2-N_1}{\parallel N_2-N_1\parallel }\end{array}$$

(2)

$$\begin{array}{cc}\hfill {\widehat{\nu }}_1=-{\widehat{\omega }}_1\times {\widehat{r}}_1,& {\widehat{\nu }}_2=-{\widehat{\omega }}_2\times {\widehat{r}}_2\end{array}$$

(3)

$$\begin{array}{cc}\hfill {\widehat{\omega }}_1=({\omega }_{1x},{\omega }_{1y},{\omega }_{1z}),& {\mathsf{\Omega }}_1=\left[\begin{array}{ccc}0& -{\omega }_{1z}& {\omega }_{1y}\\ {\omega }_{1z}& 0& -{\omega }_{1x}\\ -{\omega }_{1y}& {\omega }_{1x}& 0\end{array}\right]\end{array}$$

(4)

$$\begin{array}{cc}\hfill {\widehat{\omega }}_2=({\omega }_{2x},{\omega }_{2y},{\omega }_{2z}),& {\mathsf{\Omega }}_2=\left[\begin{array}{ccc}0& -{\omega }_{2z}& {\omega }_{2y}\\ {\omega }_{2z}& 0& -{\omega }_{2x}\\ -{\omega }_{2y}& {\omega }_{2x}& 0\end{array}\right]\end{array}$$

(5)

$$\begin{array}{cc}\hfill {\xi }_1=\left(\begin{array}{c}{\widehat{\nu }}_1\\ {\widehat{\omega }}_1\end{array}\right),& {\xi }_2=\left(\begin{array}{c}{\widehat{\nu }}_2\\ {\widehat{\omega }}_2\end{array}\right)\end{array}$$

(6)

Then, we compute the Rodrigues’ formulas:

$$\begin{array}{c}\hfill e^{{\widehat{\omega }}_1{\theta }_1}=I_{3\times 3}+{\mathsf{\Omega }}_1{\theta }_2+{\mathsf{\Omega }}_1^2\left(1-cos{\theta }_1\right)\end{array}$$

(7)

$$\begin{array}{c}\hfill e^{{\widehat{\omega }}_2{\theta }_2}=I_{3\times 3}+{\mathsf{\Omega }}_2{\theta }_2+{\mathsf{\Omega }}_2^2\left(1-cos{\theta }_2\right)\end{array}$$

(8)

representing the ${\theta }_1$ and ${\theta }_2$ rotations about ${\widehat{\omega }}_1$ and ${\widehat{\omega }}_2$, respectively. Also, we compute the matrices:

$$\begin{array}{cc}\hfill e^{{\widehat{\xi }}_1{\theta }_1}=\left[\begin{array}{cc}{e}^{{\widehat{\omega }}_1{\theta }_1}& {\widehat{\tau }}_1\\ 0_{1\times 3}& 1\end{array}\right],& e^{{\widehat{\xi }}_2{\theta }_2}=\left[\begin{array}{cc}{e}^{{\widehat{\omega }}_2{\theta }_2}& {\widehat{\tau }}_2\\ 0_{1\times 3}& 1\end{array}\right]\end{array}$$

(9)

where:

$$\begin{array}{c}\hfill {\widehat{\tau }}_1=\left(I_{3\times 3}-e^{{\widehat{\omega }}_1{\theta }_1}\right){\widehat{\omega }}_1\times {\widehat{\nu }}_1+{\widehat{\omega }}_1·{{\widehat{\omega }}_1}^T{\widehat{\nu }}_1{\theta }_1\end{array}$$

(10)

$$\begin{array}{c}\hfill {\widehat{\tau }}_2=\left(I_{3\times 3}-e^{{\widehat{\omega }}_2{\theta }_2}\right){\widehat{\omega }}_2\times {\widehat{\nu }}_2+{\widehat{\omega }}_2·{{\widehat{\omega }}_2}^T{\widehat{\nu }}_2{\theta }_2\end{array}$$

(11)

If we define the initial pose representation as:

$$g_P\left(0\right)=\left[\begin{array}{cc}{R}_0& P_0\\ 0_{1\times 3}& 1\end{array}\right]$$

(12)

where:

$$R_0=\left[{\widehat{s}}_0{\widehat{n}}_0{\widehat{a}}_0\right]$$

(13)

where each column vector is:

$$\begin{array}{c}\hfill {\widehat{s}}_0=\frac{B_0+C_0-2A_0}{\parallel B_0+C_0-2A_0\parallel }\end{array}$$

(14)

$$\begin{array}{c}\hfill {\widehat{n}}_0=\frac{(B_0-A_0)\times (C_0-A_0)}{\parallel (B_0-A_0)\times (C_0-A_0)\parallel }\end{array}$$

(15)

$$\begin{array}{c}\hfill {\widehat{a}}_0={\widehat{s}}_0\times {\widehat{n}}_0\end{array}$$

(16)

Then, we finally get the product of exponentials (PoE) representation for the serial chain with two hinge joints for all the P points on the platform:

$$g_P=e^{{\widehat{\xi }}_1{\theta }_1}{e}^{{\widehat{\xi }}_2{\theta }_2}{g}_P\left(0\right)=\left[\begin{array}{cc}{R}_T& {\tau }_T\\ 0_{1\times 3}& 1\end{array}\right]$$

(17)

where $R_T$ is the total rotation matrix and ${\tau }_T$ is the total translation vector.

The model formulas are in the Appendix A.1.

2.3. Synthesis of the Parallel Tendon-Driven Robot

The cables must not be in contact with the body and be smaller than the base and platform used by the draw-wire sensors. To simplify the collision study, we analyzed a coronal section of the foot and ankle, and the contact body is a circle containing the ankle, as illustrated in Figure 6. We use this coronal section to approximate the platform and base sizes. In Figure 6, two tendons drive a hinge joint, and the platform anchor points trace two concentric circular trajectories. The segments $T_1$ and $T_2$ represent two antagonistic tendons, $r_s$ is the radius of a solid body containing the axis of rotation, and the interior circle represents a solid body. The collision contact of the cable with the body depends on the base and platform proportions and the radius $r_s$. By observing, we note that if the base is larger than the platform, we can enhance the range of motion of the hinge joint. The points $M_1$, $M_2$, and $r_1$ are references on the talocrural axis. Similarly, $N_1$, $N_2$, and $r_2$ are in the subtalar axis. We used it to compute the following vectors and matrices:

We use this coronal section to approximate the platform and base sizes. From the initial position, we define a maximum and minimum reached angle. In Figure 7, we illustrate the full range of motion of the platform.

The limits are due to collision between the cables and the base body. When $T_1$ is extending and $T_2$ is contracting, Tendon $T_1$ touches the foot in $C{P}_1$. The minimum angle is limited for such a collision. The maximum angle occurs when $T_1$ is contracting and $T_2$ is extending. We select a radius $r_s$ greater than the foot width, and then we evaluate a platform radius $r_p$ greater than the radius. By selecting a base radius $r_b$, we can evaluate the range of movement. The positive arc and its derivate are:

$$\begin{array}{ccc}\hfill y& =& \hfill \sqrt{r_s^2-x^2}\hfill \end{array}$$

(18)

$$\begin{array}{ccc}\hfill \frac{\partial y}{\partial x}& =& \hfill \frac{-x}{\sqrt{r_s^2-x^2}}\hfill \end{array}$$

(19)

The derivative is the slope of the line $T_1$ in the tangent point $P_T$:

$$\frac{-x}{\sqrt{r_s^2-x^2}}=\frac{y-y_0}{x-x_0},$$

(20)

solving for y yields:

$$y=-\frac{x^2-x·x_0-\sqrt{r_s^2-x^2}·y_0}{\sqrt{r_s^2-x^2}}$$

(21)

We found the tangential point by substituting y from Equation (18) in Equation (21) (22), and simplifying yields (23):

$$\begin{array}{c}\hfill \sqrt{r_s^2-x_T^2}=-\frac{x_T^2-x_T·x_0-\sqrt{r_s^2-x_T^2}·y_0}{\sqrt{r_s^2-x_T^2}}\end{array}$$

(22)

$$\begin{array}{c}\hfill x_T·x_0-(r^2-\sqrt{r_s^2-x_T^2}·y_0)=0\end{array}$$

(23)

To solve for $x_T$, we use the software SageMath 10.0. Then, we replace the positive value in Equation (18) to find $y_T$.

Finally, we found the intersection $P_2$ with the circular trajectory that has a radius $r_2$. To solve for $x_2$, we use the following equation:

$$\frac{y_T-y_0}{x_T-x_0}(x_2-x_0)+y_0=\sqrt{r_2^2-x_2^2}$$

(24)

We found $y_2$ by replacing $x_2$ in Equation (18).

To find the arch length, we consider the initial position of anchor point $P_1$ on the platform. We use the following Equation:

$$d_{af}=\sqrt{r_2^2-x_1^2}$$

(25)

where $d_{af}$ is the length from the circle center to the platform central point and defines the trajectory radius regarding the positive semicircle. We solve for $x_1$ yields two values. By selecting the positive value, we obtain $y_1=d_{af}$. The arc length is given by the absolute difference of the two corresponding angles. We computed the angles by using the Equation:

$$\gamma =|{arctan}_2\left(y_2,x_2\right)-{arctan}_2\left(y_1,x_1\right)|$$

(26)

2.3.1. Cable-driven Two-Rotational Serial Chain.

In this subsection, we study a model in three dimensions that is analogous to the ankle joint. The model uses the screws in an antagonistic configuration to achieve tension forces at the cables. In the Figure 8, the view is normal to the proximal axis, and the Figure 8 shows the normal view of the distal axis.

We use this simplified model to compute the cable locations. Also, we can use the model for testing and calibrating the mechanism before use in humans. With the purpose of designing a reconfigurable physical model, we represent the rotational joints as axes located coincident in opposite edges on a tetrahedral structure, as shown in Figure 9. We choose two different configurations by rotating the anchor points by ${45}^{\circ }$ (half of ${90}^{\circ }$) from their initial position on the base and the platform. The first coincides with the rotation axis as in Figure 9, the second configuration is similar to that observed in the ankle, and we show this configuration in Figure 9.

The first representation allows us to visually identify the intersection points between the tendons’ lines of action.

$$\begin{array}{cc}\hfill L_{{\$}_1}\cap L_{{\$}_{12}}=P_1,& L_{{\$}_2}\cap L_{{\$}_{12}}=P_2,\end{array}$$

(27)

$$\begin{array}{cc}\hfill L_{{\$}_3}\cap L_{{\$}_{34}}=P_3,& L_{{\$}_4}\cap L_{{\$}_{34}}=P_4\end{array}$$

(28)

Such a condition results in null reciprocal twists because they are coplanar [49]. So the reciprocal products are:

$$\begin{array}{cc}\hfill {\$}_1\circ {\$}_{12}=0,& {\$}_2\circ {\$}_{12}=0,\end{array}$$

(29)

$$\begin{array}{cc}\hfill {\$}_3\circ {\$}_{34}=0,& {\$}_4\circ {\$}_{34}=0\end{array}$$

(30)

This configuration must be avoided because it leads to a singularity.

The relations for all the tendons on the rotary joints are as follows:

$$\begin{array}{cc}\hfill {\$}_{12}\circ \left({\$}_1+{\$}_2+{\$}_3+{\$}_4\right),& {\$}_{34}\circ \left({\$}_1+{\$}_2+{\$}_3+{\$}_4\right)\end{array}$$

(31)

In the Figure 9, we show the following:

$$\begin{array}{cc}\hfill {\$}_1\circ {\$}_{34}=-{\$}_1\circ {\$}_{34},& {\$}_3\circ {\$}_{12}=-{\$}_4\circ {\$}_{12}\end{array}$$

(32)

Then, the sum of all reciprocal products is zero, which means that the platform is in a static position. However, it is unstable because little variation in the anchor position suddenly changes the product’s sign. In the second configuration, the twist pair ${\$}_A$ and ${\$}_C$ is antagonistic with regards to ${\$}_B$ and ${\$}_D$ about the proximal twist ${\$}_x$. Also, ${\$}_C$ and ${\$}_B$ are antagonistic with respect to ${\$}_A$ and ${\$}_D$ about the twist ${\$}_t$. We will use this in our design.

In summary, two conditions are important to avoid: singular configurations and cable collisions. We address the first condition by changing the angle position of the anchor points and the second condition by selecting the base and platform radiuses bigger than the foot standard size.

2.4. Robot Configuration

For simulation and validation, we generate trajectories and capture the axes pose using the formulas in the Appendix A.2 and Appendix A.3. The method uses the formulas in the Appendix A.4 for selecting the trajectories. The axis approximation formulas are in the Appendix A.5. Finally, the range of motion and the common perpendicular line are in the Appendix A.6 and Appendix A.7 complete the captured bi-axial model characterization. This data are important for the minimum and maximum cable lenghts computation.

We assume that the base plane is parallel to the transversal plane of the shank. With the reference systems aligned and $P_O$ and $P{M}_0$ centered on the base and the platform. The shank line is perpendicular to the foot plane, and the neutral position occurs when the four tendons have similar tension. Therefore, the actuators used in antagonistic operation will have seamless forces. When the platform of the foot is parallel to the transversal plane of the shank, the talocrural and subtalar axes projected to the plane in symmetric form are:

$$\begin{array}{cc}\hfill \frac{x-r_{1x}}{{\omega }_{1x}}=\frac{y-r_{1y}}{{\omega }_{1y}},& \frac{x-r_{2x}}{{\omega }_{2x}}=\frac{y-r_{2y}}{{\omega }_{2y}}\end{array}$$

(33)

Solving for y yields:

$$\begin{array}{c}\hfill y=\frac{{\omega }_{1y}}{{\omega }_{1x}}x-\frac{{\omega }_{1y}}{{\omega }_{1x}}{r}_{1x}+\frac{r_{1y}}{{\omega }_{1y}}\end{array}$$

(34)

$$\begin{array}{c}\hfill y=\frac{{\omega }_{2y}}{{\omega }_{2x}}x-\frac{{\omega }_{2y}}{{\omega }_{2x}}{r}_{2x}+\frac{r_{2y}}{{\omega }_{2y}}\end{array}$$

(35)

Subtracting (35) from (34), we have:

$$x\left(\frac{{\omega }_{1y}}{{\omega }_{1x}}-\frac{{\omega }_{2y}}{{\omega }_{2x}}\right)-\left(\frac{{\omega }_{1y}}{{\omega }_{1x}}{r}_{1x}-\frac{{\omega }_{2y}}{{\omega }_{2x}}{r}_{2x}\right)+\frac{r_{1y}}{{\omega }_{1y}}-\frac{r_{2y}}{{\omega }_{2y}}=0$$

(36)

Solving for x yields:

$$x_p=-\frac{({\omega }_{2y}{r}_{2x}-{\omega }_{2x}{r}_{2y}){\omega }_{1x}+({\omega }_{1x}{\omega }_{2x}-{\omega }_{1y}{\omega }_{2x})r_{1x}}{{\omega }_{1y}{\omega }_{2x}-{\omega }_{1x}{\omega }_{2y}}$$

(37)

Replacing x in (33) and solving for y yields:

$$y_p=-\frac{({\omega }_{2y}{r}_{2x}-{\omega }_{2x}{r}_{2y}){\omega }_{1y}+({\omega }_{1x}{\omega }_{2y}-{\omega }_{1y}{\omega }_{2y})r_{1x}}{{\omega }_{1y}{\omega }_{2x}-{\omega }_{1x}{\omega }_{2y}}$$

(38)

We create two planes parallel to the z-axis, coincident with the intersecting point $P_{ip}=(x_p,y_p)$:

$$\begin{array}{cc}\hfill {\mathsf{\Pi }}_1:(P-P_{ip})·{\widehat{n}}_{1\parallel z}=0,& {\mathsf{\Pi }}_2:(P-P_{ip})·{\widehat{n}}_{2\parallel z}=0\end{array}$$

(39)

where the normal vectors are:

$$\begin{array}{cc}\hfill {\widehat{n}}_{1\parallel z}={\widehat{\omega }}_1\times [0,0,1],& {\widehat{n}}_{2\parallel z}={\widehat{\omega }}_2\times [0,0,1]\end{array}$$

(40)

We compute the angle between the planes. This equation can be used to calculate the minimum and maximum angles, and the angles can be obtained as follows:

$$\begin{array}{c}\hfill {\gamma }_{12}=\frac{1}{2}arccos\left(\frac{{\widehat{n}}_{1\parallel z}·{\widehat{n}}_{2\parallel z}}{|n_1\left|\right|n_2|}\right),\end{array}$$

(41)

$$\begin{array}{c}\hfill {\gamma }_{21}=\frac{1}{2}\left[\pi -arccos\left(\frac{{\widehat{n}}_{1\parallel z}·{\widehat{n}}_{2\parallel z}}{|n_1\left|\right|n_2|}\right)\right]\end{array}$$

(42)

We rotate the unitary vector normal to the planes ${\gamma }_{12}$ and ${\gamma }_{21}$ about an axis parallel to the z-axis and pass it through $P_{ip}$. The resulting vectors are:

$$\begin{array}{cc}\hfill {\omega }_{12}=\mathrm{rot}(\widehat{k},-{\gamma }_{12}),& {\omega }_{21}=\mathrm{rot}(\widehat{k},-{\gamma }_{21})\end{array}$$

(43)

where $\widehat{k}$ is the unitary vector in the direction of the positive z-axis.

The resulting lines in symmetric form are:

$$\begin{array}{cc}\hfill \frac{x-x_p}{{\omega }_{12x}}=\frac{y-y_p}{{\omega }_{12y}},& \frac{x-x_p}{{\omega }_{21x}}=\frac{y-y_p}{{\omega }_{21y}}\end{array}$$

(44)

Then, we obtain the pivot point $P_{pp}=\left[x_{pp},y_{pp},z_{pp}\right]$ on the robot platform, and the pivot point $P_{pb}=\left[x_{pb},y_{pb},z_{pb}\right]$ on the base, corresponding to $\left[x_p,y_p\right]$ projected on the platform and the base. The platform and base attachment points are at the intersection between the circle centered on the pivot points $P_{pp}$ and $P_{pb}$ with the lines of Equations (44), the circle equations are:

$$\begin{array}{cc}\hfill y^2=r_p^2-{(x-x_{pp})}^2+y_{pp},& y^2=r_b^2-{(x-x_{pb})}^2+y_{pb}\end{array}$$

(45)

The platform radius is $r_p$, and the base radius is $r_b$. By solving for y in (44), substituting it in (45), and solving for x, we can obtain the following two values:

$$\frac{\sigma \mp \sqrt{\left({\omega }_{12x}^2+{\omega }_{12y}^2\right)\left(y_{pp}+r_p^2\right)-{\omega }_{12x}^2{y}_{pp}^2}}{{\omega }_{12x}^2+{\omega }_{12y}^2}$$

(46)

where:

$$\sigma =\left({\omega }_{12x}^2+{\omega }_{12y}^2\right)x_{pp}-{\omega }_{12x}{\omega }_{12y}{y}_{pp}$$

(47)

Substituting these values back into Equation (44) yields:

$$y_{a{p}_1,a{p}_2}=y_{pp}-\frac{{\omega }_{12y}}{{\omega }_{12x}}\left(x_{pp}-x_{a{p}_1,a{p}_2}\right)$$

(48)

• Then, by substituting ${\omega }_{12}$ by ${\omega }_{21}$, we can obtain the other two points. Finally, we obtained four points for the platform. Obtaining the base anchor points is similar to computing the platform anchor points, replacing $P_{pp}$ by $P_{bp}$, and $r_p$ by $r_b$ we compute four base corresponding points. The reconfigurable structure is easy to set up by changing the angle position given by:

  $$\begin{array}{cc}\hfill {\theta }_{pi}=arctan(y_{a{p}_i},x_{a{p}_i}),& {\theta }_{bi}=arctan(y_{a{b}_i},x_{a{b}_i})\end{array}$$

  (49)

  where ${\theta }_{pi}$ and ${\theta }_{bi}$ are the corresponding angles of each anchor point on the platform and the base related to the x direction around the intersection projected points $P_{pp}$ and $P_{bp}$.

3. Workspace From the Product of Exponentials

By knowing the platform’s initial configuration, we can plot the anchor point group of movements from the product of exponential matrices by applying the Equation 17 to each platform anchor point. The product of exponential matrices for each anchor point on the platform is:

$$g_{a{p}_i}=e^{{\widehat{\xi }}_1{\theta }_1}{e}^{{\widehat{\xi }}_2{\theta }_2}{g}_{a{p}_i}\left(0\right)=\left[\begin{array}{cc}{R}_{a{p}_i}& {\tau }_{a{p}_i}\\ 0_{1\times 3}& 1\end{array}\right]$$

(50)

Moreover, the cable lengths can be calculated as follows:

$$l_{ci}=\parallel {\tau }_{a{p}_i}-a{b}_i\parallel$$

(51)

The range of the talocrural angle is ${\theta }_1\in \left[-{20}^{\circ },{20}^{\circ }\right]$, and the subtalar angle range is ${\theta }_2\in \left[-{15}^{\circ },{15}^{\circ }\right]$ The surfaces represent the group of movements for each anchor point on the platform.

3.1. Reconfiguration and Statics Simulation

The robot configuration depends on the ankle axis location. We compared the analogous two-axis system and the ankle biaxial model on MuJoCo 2.3.1. We used this software because it is open-source and free to use. Additionally, it is used for model-based control. We provided the XML source in [43]. The simulation process is easy, we edited the tendon lengths in the text editor and reloaded the simulation in MuJoCo 2.3.1. We stopped the simulation, we changed the hinge joint angles, then we ran the simulation with the computed tendon lengths. The simulation stops in a static position at previously computed angles.

4. Results

4.1. Base and Platform Dimensions.

First, we show the resulting base radius in Table 1 with a relation $r_b=1.3r_p$. Every module has radius $r_m$.

The draw-wire sensor modules are in the base, we compute the position of each module by translating the origin centered triangle with vertexes $P_1$, $P_2$, and $P_3$ in Table 2.

The base reference vertexes are in Table 3.

We can generate sensor lengths for the synthetic data generation by subtracting points pertaining to the groups from the corresponding module vertexes. For example, we get the distance between the initial position from $A_0$ to the module $P_{1A}$.

$$l_{1A}=\parallel A_0-P_{1A}\parallel ,$$

(52)

We do the same for all points, then we have the lengths from each initial point tho the corresponding sensors in Table 4.

From the sensor lengths, we compute the initial position for A, B and C from the Equations (A15), (A16) and (A17). Then we compare the original with the computed sensor lengths in Table 5.

4.2. Trajectories From the Group of Movements.

In this subsection, we use synthetic trajectories generated from the product of exponential matrices. We take 50 samples for each generated trajectory and obtain the plane containing the trajectories in Table 6.

The resulting circle-fitting estimation is in Table 7.

We compute the mean rotation axis and the mean axis and center. The rotation axis and trajectory center are in Table 8.

We apply a similar method for the subtalar axis estimation, then we summarize the resulting normal vector, angle, and axis regarding the plane parallel to x-y in Table 9.

The resulting circle-fitting estimation is in Table 10.

The mean rotation axis and mean trajectory center for the talocrural axis are in Table 11.

We summarize the resulting characterization compared to the provided by the synthetic data and the minimal distance error in Table 12.

4.3. Talocrural Range of Motion.

In this subsection we use decimated data for illustrative purposes. First, we divided the trajectories in two groups. From the initial position reference we separate the negative and the positive angles. The Table 13 shows the negative angle group.

The points differences for negative angles are in Table 14.

The rotation matrix vectors components are in Table 15.

The rotation matrices and the angles are in Table 16.

The greatest negative value is 0.34906 radians, corresponding to $-{20}^{\circ }$. The positive values start from the initial position tho the trajectory end, are in Table 17.

The points differences for the positive angles are in Table 18.

The rotation matrix components are in the Table 19.

The rotation matrices for the positive angles are in Table 20.

The greatest positive value is 0.34906, corresponding to $+{20}^{\circ }$. Then the total range of motion is ${40}^{\circ }$.

4.4. Robot’s Platform and Base Sizes.

In this section, we select the platform, the base sizes and the location of the anchor points for mean values. We use Figure 7 and the Equation (26). The body segment is proportional to the human height H=1752 mm. The robot platform and base radius are smaller than the platform and base of the draw-wire sensor system.

We multiply the robot platform radius $r_{pr}=P_f·r_p$ and robot base radius $r_{br}=P_f·r_b$, where Pf is a factor from 50% to 100% of the sensor base and sensor platform.

The range of motion for different proportional values are in Table 21.

For H maximum, the angles are in Table 22.

For H minimum, the resulting angles are in Table 23.

For all three cases, we need at least ${30}^{\circ }$ for the angle, also a proportion near to the sensors or the foot is not desirable, we select a proportion of 80% from the base (0.8rb) and platform (0.8rp) of the draw-wire sensor system, then and , which are enough for the mean height H.

The Figure 10 shows the range for the mean, maximum and minimum heights.

4.5. Screws and Reciprocity Initial Analysis.

In this subsection, we analyze the platform screws system. We use Figure 9 as a reference, and then we build a model based on the cable anchor with the platform and base. Then, we reconfigure the anchor point places to ${45}^{\circ }$ rotation about the z unitary axis $\widehat{k}$. The base and platform corresponding points are provided in the Table 24.

The resulting reciprocals are in Table 25.

We observe that in the 45° rotated configuration, each pair of cable-driven actuators works in the tensegrity mode. Additionally, it is possible to apply work from the cables in the tension to the rotation axes in the serial kinematic chain.

4.6. Projection of the Intersection Pivot Point.

The projected axis on the plane x-y, has a intersection point. We use the different axis attitudes for finding the intersection point. The results are in the Table 26.

4.7. Robot Configuration

We show how to align the robot at the initial projected pivot point intersection. It is imperative to align the references before placing the cable anchors on the platform and the base. For simulation, we use the ankle forward kinematics, common perpendicular, and range of motion in the Appendix B.1, Appendix B.2 Appendix B.3 appendix sections.

In the Figure 20, the horizontal and vertical guides have markers when the robot is centered with respect to the draw-wire sensor system. On the top-left side, we designed scales for fine alignment after patient-specific ankle-axis estimation. The pivot point depends on the anthropometric measurements; in this case, we use the standard deviation and the mean values for the common perpendicular distance Q=mean. The resultin projection points are in Table 27.

For Q maximum, we show the results in Table 28.

For Q minimum, the results are in Table 29.

We note that by changing the different parameters, the robot has sufficient range to reconfigure the anchor points. So the aligning system work with different axis attitudes.

4.8. Angle Computation of the Anchor Points

After we centered the robot around the intersection of the talocrural and subtalar axes projected on the base and initial platform planes, we configured the four cable anchor points. To do so, we use Equations (46), (47), and (48) to compute the anchor points. In this subsection, we only change the axis attitudes and then compute the anchor points. The Table 30 shows the resulting computed data for the axes’ mean, maximum, and minimum attitude values. We refer to the axes’ mean attitude when TC and ST have mean statistical values. Also, we call the maximal attitude when we add the standard deviation to the TC and ST means. The same goes for the minimum.

From this table, we compute the angle configuration of the attaching points in Table 31.

The angles for the base and platform anchor points for maximum attitude are in Table 32.

The angles for the base and platform anchor points for maximum attitude are in Table 33.

The Figure 11 shows the anchor points the three axes attitude groups.

4.9. Workspace

The workspace is limited by the range of motion of the ankle joint. We use Equation (50) to compute the different groups of movements shown in Figure 12.

4.10. Cable Lenghts at Initial Position

We use Equation (51) to compute the cable lengths at the initial position.

The resulting lengths for the mean values are shown in the Table 34.

For the initial lenghts, we draw the forward kinematics is the intersection between four spheres and the circle representing the platform, with the pose represented on $g_{PM\left(0\right)}$ , as in the Figure 12.

4.11. Cable Lengths at Extreme Positions

We create three tables with different axis attitudes and compute the positions. The Table 35 shows $r_{2_i}$ and $P_{M_i}$.

The Table 36 shows the anchor points $a{p}_1$, $a{p}_2$, $a{p}_3$, and $a{p}_4$.

The Table 37 lengths for the cables.

The Figure 13 shows the four configurations corresponding to axes’ mean attitude values.

The resulting platform position anchor points corresponding to the axes’ maximum attitude values are in Table 38.

The cable lengths corresponding to the axes’ maximum attitude values are in Table 39.

The Figure 14 shows the four configurations corresponding to the axes’ maximum attitude values.

The resulting platform position anchor points corresponding to the axes’ minimum attitude values are in Table 40.

The cable lengths corresponding to the axes’ minimum attitude values are in Table 41.

The Figure 15 shows the four configurations corresponding to the axes’ minimum attitude values.

4.12. MuJoCo Simulations

In this section, we replace the anchor points, the kinematic chain and the tendon lengths in the XML script included in the repository [43].

In the simulation, we rotate ${90}^{\circ }$ about the y-axis and scale the positions by a factor of 0.1. The units are in centimeters. The input values are in Table 42.

We put this model the XML document as in the Listing 4.12.

XML MuJoCo file.

<?xml version="1.0" encoding="UTF-8"?>

<mujoco model="turmellbot">

<compiler coordinate="local"/>

<default>

<geom rgba=".8 .6 .4 1"/>

</default>

<asset>

<texture type="skybox" builtin="gradient" rgb1="1 1 1" rgb2=".6 .8 1" width="256" height="256"/>

<mesh file="a-foot.stl" refpos="-10.84427 0.1037934 -17.62868" refquat="0.99619 0.034862 0.00000 0.00000"/>

</asset>

<worldbody>

<light pos="10 -10 50" dir="-10 10 -50" diffuse="1 1 1"/>

<site name="reference" pos="5.938405e-19 -0.09110710 17.62735"/>

<body>

<geom name="base" type="cylinder" fromto="-0.4000000 -0.09110710 17.62735 5.938405e-19 -0.09110710 17.62735" size="12.032592"/>

<site name="b1" pos="-6.093131e-16 6.659522 7.666812" size="0.2"/>

<site name="b2" pos="6.105008e-16 -6.841736 27.58788" size="0.2"/>

<site name="b3" pos="4.139506e-16 9.869429 24.37798" size="0.2"/>

<site name="b4" pos="-4.127630e-16 -10.05164 10.87672" size="0.2"/>

<geom name="leg" type="capsule" fromto="5.938405e-19 -0.09110710 17.62735 10.84427 -0.1037934 17.62868" size="0.3"/>

<site name="shnk" pos="10.84427 -0.1037934 17.62868"/>

<body>

<geom name="comp" type="capsule" fromto="10.84427 -0.1037934 17.62868 11.19642 -0.01170672 17.90907" size="0.3"/>

<joint name="TC" type="hinge" pos="10.84427 -0.1037934 17.62868" axis="-0.1740000 0.9790000 -0.1030000" limited="true" range="-20 20"/>

<body>

<geom name="b-foot" type="capsule" fromto="11.19642 -0.01170672 17.90907 17.61765 -0.09110710 17.62735" size="0.3"/>

<joint name="ST" type="hinge" pos="11.19642 -0.01170672 17.90907" axis="-0.6420000 0.2080000 0.7380000" limited="true" range="-15 15"/>

<geom name="ptfm" type="cylinder" fromto="17.61765 -0.09110710 17.627 18.01765 -0.09110710 17.62735" size="9.25584"/>

<site name="ptfmc" pos="18.01765 -0.09110710 17.62735" size="0.3"/>

<geom name="s-foot" type="mesh" mesh="a-foot" rgba=".85 .79 .72 0.7" density="0.01"/>

<site name="p1" pos="17.61765 5.101684 9.965397" size="0.2"/>

<site name="p2" pos="17.61765 -5.283899 25.28930" size="0.2"/>

<site name="p3" pos="17.61765 7.570844 22.82014" size="0.2"/>

<site name="p4" pos="17.61765 -7.753058 12.43456" size="0.2"/>

</body>

</body>

</body>

</worldbody>

<tendon>

<spatial name="t1" limited="true" range="17.83503 17.83513" width="0.05" stiffness="100" damping="100000">

<site site="b1"/>

<site site="p1"/>

</spatial>

<spatial name="t2" limited="true" range="17.83503 17.83513" width="0.05" stiffness="100" damping="100000">

<site site="b2"/>

<site site="p2"/>

</spatial>

<spatial name="t3" limited="true" range="17.83503 17.83513" width="0.05" stiffness="100" damping="100000">

<site site="b3"/>

<site site="p3"/>

</spatial>

<spatial name="t4" limited="true" range="17.83503 17.83513" width="0.05" stiffness="100" damping="100000">

<site site="b4"/>

<site site="p4"/>

</spatial>

</tendon>

<actuator>

<muscle name="t1" tendon="t1"/>

<muscle name="t2" tendon="t2"/>

<muscle name="t3" tendon="t3"/>

<muscle name="t4" tendon="t4"/>

</actuator>

<sensor>

<accelerometer name="ptfacc" site="ptfmc"/>

<gyro name="ptfgyro" site="ptfmc"/>

<tendonpos name="t1pos" tendon="t1"/>

<tendonpos name="t2pos" tendon="t2"/>

<tendonpos name="t3pos" tendon="t3"/>

<tendonpos name="t4pos" tendon="t4"/>

<jointpos name="TCpos" joint="ST"/>

<jointpos name="STpos" joint="TC"/>

</sensor>

</mujoco>

There are bodies like the base and the platform, three links, two hinge joints, four tendons with sensors and actuators. The data configuration is from the simulation at the initial position. We save the file with XML format, then drag and drop in Simulate, and the resulting screenshot is in Figure 16.

The sensors in the print data button from MuJoCo generate the MJDATA.txt file, which contains the sensor measurements. We put the sensors ouptut in Table 43.

Regarding the minimum and maximum values, we stop the simulation, changing the range of the tendons and moving the joints sliders to the minimum angles.

Snapshots of the different configurations are in the Figure 17.

4.13. Resulting Robot Design

In this subsection, we used the data captured from the ankle model and integrated the draw-wire sensors into the robot. The device dimensions are based on human proportions with a height of H=175 cm. The objective is to adjust the measurement device and the ankle approximated model in a configurable structure. The design is intended for the laying position. We divide the design into two main subassemblies: the platform for the foot and the base for the shank.

The resulting platform design is based on the foot anatomy by observing the Figure 1, Figure 1, and Figure 2. The platform is adaptable to various foot sizes based on proportions [45]. The length is also adaptable, and we show the assembly in Figure 18. We added a length ruler, heel support, and three force sensing resistors (FSRs). Two FSRs are used for the forefoot, and one FSR is used for the hindfoot.

To align the center of the platform to $P_{pp}$ at the initial position. We designed perpendicular sliders, as shown in Figure 19.

In Figure 20, we show the base assembly, with sensors for ankle axis estimation. We design guides for centering the shank position $P_{bp}$. The anchor points are manually adjusted in different positions, depending on the ankle model.

A main tube structure is attached to a baseplate. We use spacers and 8 mm steel bolts to fasten the two plates supporting the sensors. Finally, we place Bowden guides for the cable endpoints.

The final design includes the platform, the sensors, the base, and a possible configuration for the actuators and electronics. We also show a resulting assembly with the approximated adjustable axis mechanical model in Figure 21.

5. Discussion

In this work, we show that, by using the human-centered design of an ankle rehabilitation robot, we can use the robot in a broader group of patients following anthropometry, human proportions, and statistics.

We focus on the ankle model specific to each patient to equilibrate tension and pressure forces in an initial position. For us, ankle kinematic model identification is imperative to configure the platform and base anchor points. We show that such a design can enhance the range of motion, adapting to several sizes on both feet. Using this human-centered approach, we can limit the pressure forces acting on the plantar surface of the foot to avoid unnatural positions. We estimate the platform and base sizes. We also ensured that the cables did not touch the foot. The cable lengths can reach the ankle joint range of motion. We reconfigured the robot around the initial equilibrium position. This position is like when the human body is standing. We used projected axes on the transversal planes perpendicular to gravity. Screw theory is a powerful tool, and the results show we can effectively use it in robot geometry, kinematics, and static analysis. Using MuJoCo, we simulated the statics by editing an XML file with the computed tendon lengths and anchor points, then using sensors, the ankle joint angles measures are like the initially assumed.

6. Conclusions

Ankle sprains are a common injury, and there is abundant research on ankle rehabilitation robots. We found that Screw Theory is applied in many works. By observing the anatomy, we realized that the tendons associated with the ankle movement are attached to the bones at the base of the plantar dome, transmitting pressure to the foot surface in contact with the platform. We designed the Turmell-Bot to be used by patients lying in bed or sitting. The robot configuration depends on the patient-specific ankle model. The model is an approximation that we can refine by using piece-wise function approximations and machine learning. We designed a lightweight, low-cost, low-energy, portable, configurable, and comfortable device. We enhance the device with sensors to measure the foot pressure forces involved in ankle movements. We will use compliant actuators and ratchets to hold the desired position without energy consumption. We will search for antagonistic actuation and tension control. We also plan to integrate electromyography (EMG) and functional electrical stimulation (FES) through the shank attachment to register the activation signals when the device is making rehabilitation movements.

Appendix A.1. Ankle Model.

For obtaining the ankle joint model, we use data from anthropometric studies. The values of K, L, O, P, Q with W, and w ratio with standard deviation are in Table A1.

Table A1. Mean anthropometric values.

Variable	K (mm)	L (mm)	O (mm)	P (mm)	Q (mm)	R = W/w
Mean value	12	16	11	-1	5	0.54

Table A1. Mean anthropometric values.

Variable	K (mm)	L (mm)	O (mm)	P (mm)	Q (mm)	R = W/w
Mean value	12	16	11	-1	5	0.54

We show the graphical representation in Figure A1. The figure shows the malleolus most lateral point (MLP), and malleolus mos medial point (MMP). The other values are Q, W and w. $M_1$, $M_2$ pertain to the talocalcaneal axis; $N_1$ and $N_2$ to the subtalar axis.

The vectors representing the talocrural and talocalcaneal axis can be reflected about the plane y=0 by using the following equation.

$${\mathsf{Ref}}_{\widehat{n}}\left(\widehat{\nu }\right)=\widehat{\nu }-2\frac{\widehat{\nu }·{\widehat{n}}_j}{{\widehat{n}}_j·{\widehat{n}}_j}{\widehat{n}}_j.$$

(A1)

Where ${\widehat{n}}_j$ is the unitary column vector:

$${\widehat{n}}_j={\left[0,1,0\right]}^T,$$

(A2)

yielding

$${\mathsf{Ref}}_{\widehat{n}}\left(\widehat{\nu }\right)=\widehat{\nu }-2\left(\widehat{\nu }·{\widehat{n}}_j\right){\widehat{n}}_j={\left[{\nu }_x,-{\nu }_y,{\nu }_z\right]}^T,$$

(A3)

which have the matrix representation:

$$R_{xz}=I-2\frac{{\widehat{n}}_j·{{\widehat{n}}_j}^T}{{{\widehat{n}}_j}^T{\widehat{n}}_j}$$

(A4)

Replacing the unitary vector $\widehat{j}$, the resulting reflection matrix is:

$$R_{xz}=\left[\begin{array}{cccc}& 1& 0& 0\\ & 0& -1& 0\\ & 0& 0& 1\end{array}\right]$$

(A5)

Then, by multiplying all the right ankle axis vectors by this matrix, we get the left foot kinematic model.

Figure A1. Ankle anthropometric data.

Figure A1. Ankle anthropometric data.

Appendix A.2. Synthetic Trajectories Generation

In case of the talocrural axis estimation, the trajectories are by plantarflexion and dorsiflexion movements. This means that by blocking the subtalar mobility, the exponential factor will be constant or zero, regarding the initial position election. Assuming the product of exponential formula yields:

$$g_A{\mid }_{{\theta }_2=0}=e^{{\widehat{\xi }}_1{\theta }_1}{g}_A\left(0\right)$$

(A6)

Which is a single parameter equation. It corresponds to a circle in the space. We generate the trajectory in the range ${\theta }_1\in \left[-{20}^{\circ },{20}^{\circ }\right]$. We also generate artificial trajectories for the subtalar axis, this is done by replacing the angle , and evaluating for ${\theta }_2\in \left[-{10}^{\circ },{10}^{\circ }\right]$. Then we get:

$$g_A{\mid }_{{\theta }_1=0}=e^{{\widehat{\xi }}_2{\theta }_2}{g}_A\left(0\right)$$

(A7)

Appendix A.3. Platform Position

The schematic model is in the Figure A2. It has three base triangles and three attaching points to the foot’s platform. For each point, A, B, and C, we compute the localization in the same way.

Figure A2. Capturing the platform position.

Figure A2. Capturing the platform position.

In the figure, we represented three spheres with centers A, B, and C on the plane $z=0$. Also, we represented the distances $r_1$, $r_2$, and $r_3$ measured by each sensor. We translated each equilateral triangle to the origin, then, the vertexes are $P_1=\left[r_t,0,0\right]$, $P_2=\left[-\frac{1}{2},r_t,\frac{1}{2},\sqrt{3}{r}_t,0\right]$, and $P_3=\left[-\frac{1}{2},r_t,-\frac{1}{2},\sqrt{3}{r}_t,0\right]$. Where $r_t$ is the radius from each module center to his vertexes. The sphere’s radius in each module are $r_1$, $r_2$, and $r_3$, then the equations are:

$$\begin{array}{ccc}\hfill S_1:{\left(r_t-x\right)}^2+y^2+z^2-r_1^2& =& \hfill 0\hfill \end{array}$$

(A8)

$$\begin{array}{ccc}\hfill S_2:\frac{1}{4}{\left(r_t+2x\right)}^2+\frac{1}{4}{\left(\sqrt{3}{r}_t-2y\right)}^2+z^2-r_2^2& =& \hfill 0\hfill \end{array}$$

(A9)

$$\begin{array}{ccc}\hfill S_3:\frac{1}{4}{\left(r_t+2x\right)}^2+\frac{1}{4}{\left(\sqrt{3}{r}_t+2y\right)}^2+z^2-r_3^2& =& \hfill 0\hfill \end{array}$$

(A10)

By subtracting the sphere equations $S_1$ from $S_2$, $S_3$ from $S_1$, and $S_3$ from $S_2$, we get the planes:

$$\begin{array}{ccc}\hfill P_{12}:-3r_{t}x+\sqrt{3}{r}_{t}y-r_1^2+r_2^2& =& \hfill 0\hfill \end{array}$$

(A11)

$$\begin{array}{ccc}\hfill P_{13}:-3r_{t}x-\sqrt{3}{r}_{t}y-r_1^2+r_3^2& =& \hfill 0\hfill \end{array}$$

(A12)

$$\begin{array}{ccc}\hfill P_{23}:-2\sqrt{3}{r}_{t}y-r_2^2+r_3^2& =& \hfill 0\hfill \end{array}$$

(A13)

Which are all parallel to z-axis. Adding $P_{12}$ and $P_{13}$ we have the plane:

$$P_{123}:-2r_1^2+r_2^2+r_3^2-6r_{t}x=0$$

(A14)

Solving y in the plane $P_{23}$, we get:

$$y=-\frac{\sqrt{3}\left(r_2^2-r_3^2\right)}{6r_t}$$

(A15)

And solving for x in plane $P_{123}$ , we get:

$$x=-\frac{2r_1^2-r_2^2-r_3^2}{6r_t}$$

(A16)

Finally, replacing x and y in $S_1$, and solving for z, we have:

And solving for x in plane $P_{123}$ , we get:

$$z=\mp \frac{\sqrt{-r_1^4+r_1^2{r}_2^2-r_2^4-r_3^4-9r_t^4+\left(r_1^2+r_2^2\right)r_3^2+3\left(r_1^2+r_2^2+r_3^2\right)r_t^2}}{3r_t}$$

(A17)

In the Figure A3, we show an example when $r_t=2$, $r_1=5$, $r_2=6$, and $r_3=7$

Figure A3. Sphere intersection.

Figure A3. Sphere intersection.

We choose the z-negative value because the origin is over the sensors positions. The Figure A4 shows the geometric representation of the analytic equations.

Figure A4. Equations’ geometric representation

Figure A4. Equations’ geometric representation

Finally, we translate the origin to A, B, an C center points at the sensor modules.

Appendix A.4. Selecting the Best Circle-Fitting Trajectories

In this subsection, we select the best measured trajectory. We use draw-wire sensors (DWS). The foot position and orientation can be found from the draw-wire sensors by computing vectors from the measured points $A_M$, $B_M$, and $C_M$.

$$\begin{array}{c}\hfill {\widehat{v}}_{AB}=B_M-A_M\end{array}$$

(A18)

$$\begin{array}{c}\hfill {\widehat{v}}_{AC}=C_M-A_M\end{array}$$

(A19)

The normal unitary vector is:

$${\widehat{n}}_v=\frac{{\widehat{V}}_{AB}\times {\widehat{v}}_{AC}}{\parallel {\widehat{V}}_{AB}\times {\widehat{v}}_{AC}\parallel }$$

(A20)

A unitary vector pointing forward the foot is:

$${\widehat{s}}_v=\frac{{hatV}_{AB}+{\widehat{v}}_{AC}}{\parallel {\widehat{V}}_{AB}+{\widehat{v}}_{AC}\parallel }$$

(A21)

And a perpendicular vector to ${\widehat{n}}_v$ and ${\widehat{s}}_v$ is:

$${\widehat{a}}_v={\widehat{n}}_v\times {\widehat{s}}_v$$

(A22)

Then, for each instantaneous position we have a measured transformation matrix:

$$g_M=\left[\begin{array}{cccc}{s}_x& a_x& n_x& p_x\\ s_y& a_y& n_y& p_y\\ s_z& a_z& n_z& p_z\\ 0& 0& 0& 1\end{array}\right]$$

(A23)

We add two inertial measurement units (IMU), one in the platform and the other in the base. In our design, we use the module GY-6500 with the Invensense’s MPU-6500 chip, which integrates a 3-axis gyroscope, a 3-axis accelerometer and a digital motion processor (DMP). These units have to be first calibrated in a flat surface with the platform coordinates system aligned to the base coordinate’s system. The condition for start and finish a measured trajectory is done by comparing the absolute changes in the accelerometer and the gyroscope.

$$\begin{array}{c}\hfill \parallel a_i-a_{i-1}\parallel =\parallel \left(a_{ix},a_{iy},a_{iz}\right)-\left(a_{(i-1)x},a_{(i-1)y},a_{(i-1)z}\right)\parallel >0\end{array}$$

(A24)

$$\begin{array}{c}\hfill \parallel {\alpha }_i-{\alpha }_{i-1}\parallel =\parallel \left({\alpha }_{ix},{\alpha }_{iy},{\alpha }_{iz}\right)-\left({\alpha }_{(i-1)x},{\alpha }_{(i-1)y},{\alpha }_{(i-1)z}\right)\parallel >0\end{array}$$

(A25)

In the real world, we have noise, and the change must be greater than the noise signal. The noise can be estimated during calibration and measuring absolute maximum and minimum values with the accelerometer and the gyroscope in static positions. We subtract the foot and shank measured quaternions:

$${\widehat{q}}_{SF}=\left[q_w,q_x,q_y,q_z\right]=\left[q_{Fw},q_{Fx},q_{Fy},q_{Fz}\right]-\left[q_{Sw},q_{Sx},q_{Sy},q_{Sz}\right]$$

(A26)

And converted the data to angle-axis representation with:

$$\begin{array}{c}\hfill \theta =2\ast arccos\left(q_w\right)\end{array}$$

(A27)

$$\begin{array}{c}\hfill {\omega }_x=\frac{q_x}{\sqrt{1-q_w·q_w}}\end{array}$$

(A28)

$$\begin{array}{c}\hfill {\omega }_y=\frac{q_y}{\sqrt{1-q_w·q_w}}\end{array}$$

(A29)

$$\begin{array}{c}\hfill {\omega }_z=\frac{q_z}{\sqrt{1-q_w·q_w}}\end{array}$$

(A30)

For example, we use measurements in A to get a set of points.

$$P_{MA}=\left[P_{A0},P_{A1}\dots P_{An}\right]$$

(A31)

For each point, there are a corresponding IMU attitude:

$${\widehat{\omega }}_M=\left[{\widehat{\omega }}_0,{\widehat{\omega }}_1\dots ,{\widehat{\omega }}_n\right]$$

(A32)

Then, the mean $\overline{\omega }$ of the vector list is perpendicular to the plane containing the circular trajectory.

$${\mathsf{\Pi }}_C=({\overline{\omega }}_x,{\overline{\omega }}_y,{\overline{\omega }}_z)·(P-P_{A0})=0$$

(A33)

Appendix A.5. Axes Approximation

In this subsection, we used artificial data generated from the anthropometric studies, as an approximation based on real measurements. But the real data will be obtained by using the last subsection method. We briefly describe the method for axis estimation from a selected trajectory containing positions in space. The point positions $P_A$ pertain to a plane ${\mathsf{\Pi }}_A$ with the equation:

$$Ax+By+Cz+D=0$$

(A34)

Multiplying this equation by $-\frac{1}{C}$ , $C\ne 0$:

$$-\frac{A}{C}x-\frac{B}{C}y-z-\frac{D}{C}=0$$

(A35)

If we define:

$$\begin{array}{ccc}\hfill a& =& \hfill -\frac{A}{C}\hfill \end{array}$$

(A36)

$$\begin{array}{ccc}\hfill b& =& \hfill -\frac{B}{C}\hfill \end{array}$$

(A37)

$$\begin{array}{ccc}\hfill d& =& \hfill -\frac{D}{C}\hfill \end{array}$$

(A38)

Then, we have:

$$ax+by-z+d=0,$$

(A39)

solving for z on the right side:

$$ax+by+d=z,$$

(A40)

The list of the n sample points has the form:

$$\left[\begin{array}{ccc}{x}_0& y_0& 1\\ x_1& y_1& 1\\ ⋮& ⋮& ⋮\\ x_{n-1}& y_{n-1}& 1\end{array}\right]\left[\begin{array}{c}a\\ b\\ d\end{array}\right]=\left[\begin{array}{c}{z}_0\\ z_1\\ ⋮\\ z_{n-1}\end{array}\right]$$

(A41)

We define:

$$\begin{array}{ccc}\hfill A_{n\times 3}& =& \hfill \left[\begin{array}{ccc}{x}_0& y_0& 1\\ x_1& y_1& 1\\ ⋮& ⋮& ⋮\\ x_{n-1}& y_{n-1}& 1\end{array}\right]\hfill \end{array}$$

(A42)

$$\begin{array}{ccc}\hfill x_{3\times 1}& =& \hfill \left[\begin{array}{c}a\\ b\\ d\end{array}\right]\hfill \end{array}$$

(A43)

$$\begin{array}{ccc}\hfill B_{n\times 1}& =& \hfill \left[\begin{array}{c}{z}_0\\ z_1\\ ⋮\\ z_{n-1}\end{array}\right]\hfill \end{array}$$

(A44)

Then, the system:

$$A_{n\times 3}{x}_{3\times 1}=B_{n\times 1}$$

(A45)

has the solution:

Then, the system:

$$x_{3\times 1}=A^{+}B$$

(A46)

where:

$$A^{+}={\left(A^{T}A\right)}^{-1}{A}^T$$

(A47)

Is the pseudo-inverse of $A_{n\times 3}$.

Then we have solved a, b, and d. We replace ,and solve for A, B and D in Equations (A36), (A37), and (A38). In the Equation (A34), the vector ${\widehat{N}}_{\mathsf{\Pi }A}=\left[A,B,C\right]$ is perpendicular to the plane ${\mathsf{\Pi }}_A$, normalizing the vector, we have:

$${\widehat{n}}_{\mathsf{\Pi }A}=\frac{{\widehat{N}}_{\mathsf{\Pi }A}}{\parallel {\widehat{N}}_{\mathsf{\Pi }A}\parallel }$$

(A48)

The angle with the plane $z=0$, normal to the vector $\widehat{k}=\left[0,0,1\right]$ is given by the dot product equation:

$${\theta }_{\widehat{n}\to \widehat{k}}=arccos({\widehat{n}}_{\mathsf{\Pi }A}·\widehat{k}),$$

(A49)

the cross product:

$${\widehat{v}}_{\widehat{n}\perp \widehat{k}}=\widehat{n}\times \widehat{k}$$

(A50)

is perpendicular to ${\widehat{n}}_{\mathsf{\Pi }A}$ and $\widehat{k}$.

Rotating all the points an angle ${\theta }_r$ about the direction vector $\widehat{\nu }=\left[{\nu }_x,{\nu }_y,{\nu }_z\right]$ can be achieved by applying the rotation transformation to each point in the arc:

$$P_R=R_{\widehat{\nu },\theta }P=e^{N\Theta }P$$

(A51)

Where:

$$e^{N\Theta }=I_{3\times 3}+Nsin\Theta +N^2\left(1-cos\Theta \right)$$

(A52)

And:

$$N=\left[\begin{array}{cccc}0& & -{\nu }_z& {\nu }_y\\ {\nu }_z& 0& & -{\nu }_x\\ -{\nu }_y& {\nu }_x& 0\end{array}\right]$$

(A53)

Now the points are in a plane parallel to the plane $z=0$, then we use the two-dimensional circle equation:

$${\left(x-x_c\right)}^2+{\left(y-y_c\right)}^2=r^2$$

(A54)

Expanding the equation yields:

$$x^2-2x{x}_c+x_c^2+y^2-2y{y}_c+y_c^2=r^2$$

(A55)

Rearranging:

Expanding the equation yields:

$$\left(2x_c\right)x+\left(2y_c\right)y+(r^2-x_c^2-y_c^2)=x^2+y^2$$

(A56)

Which has the form:

$${\chi }_{0}x+{\chi }_{1}y+{\chi }_2=x^2+y^2$$

(A57)

Where:

$$\begin{array}{ccc}\hfill {\chi }_0& =& \hfill 2x_c\hfill \end{array}$$

(A58)

$$\begin{array}{ccc}\hfill {\chi }_1& =& \hfill 2y_c\hfill \end{array}$$

(A59)

$$\begin{array}{ccc}\hfill {\chi }_2& =& \hfill r^2-x_c^2-y_c^2\hfill \end{array}$$

(A60)

The $P_R$ set of points can be arranged in the system:

$$\left[\begin{array}{ccc}{x}_0& y_0& 1\\ x_1& y_1& 1\\ ⋮& ⋮& ⋮\\ x_{n-1}& y_{n-1}& 1\end{array}\right]\left[\begin{array}{c}{\chi }_0\\ {\chi }_1\\ {\chi }_2\end{array}\right]=\left[\begin{array}{c}{x}_0^2+y_0^2\\ x_1^2+y_1^2\\ ⋮\\ x_{n-1}^2+y_{n-1}^2\end{array}\right]$$

(A61)

Has the form:

$$Ax=B$$

(A62)

Which is similar to the Equation (A45). We found the vector, we found the center $c_r=\left[x_c,y_c\right]$ and the radius $r_f$ replacing the components in Equations (A58), (A59) and (A60). Then, we rotate back the center $C_{rt}$ using the negative angle $-{\Theta }_{\widehat{n}\to \widehat{k}}$, and the same axis ${\widehat{\nu }}_{\widehat{n}\perp \widehat{k}}$. We have three generated trajectories by tracking A, B, and C. Then, we averaged the resulting vectors and centers.

$$\begin{array}{ccc}\hfill {\overline{\nu }}_{rt}& =& \hfill \frac{1}{3}({\nu }_{rtA}+{\nu }_{rtB}+{\nu }_{rtC})\hfill \end{array}$$

(A63)

$$\begin{array}{ccc}\hfill {\overline{c}}_{rt}& =& \hfill \frac{1}{3}(c_{rtA}+c_{rtB}+c_{rtC})\hfill \end{array}$$

(A64)

We apply a similar method for the subtalar axis joint. Finally, the initial axis representation is in Plücker line coordinates.

Appendix A.6. Range of Motion

With the axis estimation from a selected trajectory, we compute the corresponding angles for the different positions. The measured homogeneous transformation matrix in the Equation (A.4) measured rotation matrix.

$$R_M=\left[\begin{array}{ccc}{s}_x& a_x& n_x\\ s_y& a_y& n_y\\ s_z& a_z& n_z\end{array}\right],$$

(A65)

by using the definition of logarithm matrix, we compute the rotation angles as follows. We take the measured trajectories $A_M$, $B_M$, and $C_M$ to compute the vectors $\widehat{s}$, $\widehat{a}$, and $\widehat{n}$ to get the rotation matrix for each set of measurements. For each trajectory set $A_{M\left(i\right)}$, $B_{M\left(i\right)}$, $C_{M\left(i\right)}$ with index i=1 to n measurements. The estimate the angle from the initial position is:

$${\theta }_m=arccos\left(\frac{1}{2}\left(\mathbf{tr}\left(R_M\right)-1\right)\right)\in \left[0,\pi \right)$$

(A66)

Where $\mathbf{tr}\left(R_M\right)$ is the trace of the measured matrix $R_M$. In this equation, the angles must be positive, then we choose two groups of measurements. In the set of computed angles, the zero value divides the table in negative and positive values, we defined ${\theta }_{mn}$ and ${\theta }_{mp}$ , as the negative and positive group of angles. The first group is from the zero position to the maximum value in the negative rotation regarding the rotation axis. The negative maximum angle from the initial position is:

$${\theta }_{min}=-\mathbf{max}\left({\theta }_{mn}\right)$$

(A67)

The positive maximum angle from the initial position is:

$${\theta }_{min}=\mathbf{max}\left({\theta }_{mp}\right)$$

(A68)

The range of movement (Rom) is the difference between the positive and the negative value.

$$\mathbf{R}.\mathbf{O}.\mathbf{M}.={\theta }_{mp}-{\theta }_{mn}$$

(A69)

This procedure has the same results when we add the two greatest values. We validate the data from the measured rotation matrix, with the skew-symmetric matrix:

$$\left[{\widehat{\omega }}_m\right]=\frac{1}{2sin{\theta }_m}\left(R_M-R_M^T\right)$$

(A70)

With the selected trajectories, we prove that the angle identification is correct.

Appendix A.7. Common Perpendicular and its Feet

The lines representing the axes contain two points $r_{1C}$ and $r_{2C}$, which are the feet of the common perpendicular line between the axes. We use the Plücker line coordinates to find the intersection points with the axis. Then, we represent the lines with six-dimensional vectors:

$$\begin{array}{ccc}\hfill {\$}_1& =& \hfill \left[{\widehat{\omega }}_1,{\widehat{\nu }}_1\right]\hfill \end{array}$$

(A71)

$$\begin{array}{ccc}\hfill {\$}_2& =& \hfill \left[{\widehat{\omega }}_2,{\widehat{\nu }}_2\right]\hfill \end{array}$$

(A72)

The common perpendicular components are:

$$\begin{array}{ccc}\hfill {\widehat{\omega }}_C={\widehat{\omega }}_{1\perp 2}& =& \hfill {\widehat{\omega }}_1\times {\widehat{\omega }}_2\hfill \end{array}$$

(A73)

$$\begin{array}{ccc}\hfill {\widehat{\nu }}_C={\widehat{\nu }}_{1\perp 2}& =& \hfill {\widehat{\nu }}_1\times {\omega }_2-{\widehat{\nu }}_2\times {\widehat{\omega }}_1+\frac{({\$}_1·{\$}_2)({\widehat{\omega }}_1·{\widehat{\omega }}_2)}{\parallel {\widehat{\omega }}_1\times {\widehat{\omega }}_2{\parallel }^2}({\widehat{\omega }}_1\times {\widehat{\omega }}_2)\hfill \end{array}$$

(A74)

We compute the intersection points between the common perpendicular and the axes with:

$$\begin{array}{ccc}\hfill r_{1C}=r_{1\perp 2}& =& \hfill \frac{-{\widehat{\nu }}_1\times \left[{\widehat{\omega }}_2\times \left({\widehat{\omega }}_1\times {\widehat{\omega }}_2\right)\right]+\left[{\widehat{\nu }}_2·\left({\widehat{\omega }}_1\times {\widehat{\omega }}_2\right)\right]{\widehat{\omega }}_1}{\parallel {\widehat{\omega }}_1\times {\widehat{\omega }}_2{\parallel }^2}\hfill \end{array}$$

(A75)

$$\begin{array}{ccc}\hfill r_{2C}=r_{2\perp 1}& =& \hfill \frac{{\widehat{\nu }}_2\times \left[{\widehat{\omega }}_1\times \left({\widehat{\omega }}_1\times {\widehat{\omega }}_2\right)\right]-\left[{\widehat{\nu }}_1·\left({\widehat{\omega }}_1\times {\widehat{\omega }}_2\right)\right]{\widehat{\omega }}_2}{\parallel {\widehat{\omega }}_1\times {\widehat{\omega }}_2{\parallel }^2}\hfill \end{array}$$

(A76)

Finally, we have the initial position of the two axes. We use the transverse planes perpendicular to the total tension force in the shank-ankle-foot kinematic chain. For select a central point for the attaching points configuration, we found the projection of each axis on the base plane. It is the result of multiplying the orthogonal projection matrix on the base $P_b$ by ${\widehat{\omega }}_1$:

$$P_b=\left[\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ 0& 0& 0\end{array}\right]·{\widehat{\omega }}_1=\left[{\omega }_{1}x,{\omega }_{1}y,0\right]$$

(A77)

Also, we project the vectors, ${\widehat{\omega }}_2$ , $r_{1C}$, and $r_{2C}$. The resulting projections define two lines, we define two points pertaining to the talocrural axis in homogeneous representation:

$$\begin{array}{ccc}\hfill P_{1b1}& =& \hfill P_b·{\widehat{r}}_{1C}+\widehat{k}\hfill \end{array}$$

(A78)

$$\begin{array}{ccc}\hfill P_{1b2}& =& \hfill P_b·\left({\widehat{r}}_1+{\widehat{\omega }}_1\right)+\widehat{k}\hfill \end{array}$$

(A79)

And for the subtalar axis:

$$\begin{array}{ccc}\hfill P_{2b1}& =& \hfill P_b·{\widehat{r}}_{2C}+\widehat{k}\hfill \end{array}$$

(A80)

$$\begin{array}{ccc}\hfill P_{2b2}& =& \hfill P_b·\left({\widehat{r}}_2+{\widehat{\omega }}_2\right)+\widehat{k}\hfill \end{array}$$

(A81)

With $\widehat{k}=\left[0,0,1\right]$. Then by using determinants the intersection point $P_{ixy}=\left[P_{i}x,P_{i}y,0\right]$ components are:

$$\begin{array}{ccc}\hfill P_{ix}& =& \hfill \frac{\left|\begin{array}{cc}{D}_{1xy}& D_{1x1}\\ D_{2xy}& D_{2x1}\end{array}\right|}{\left|\begin{array}{cc}{D}_{1x1}& D_{1y1}\\ D_{2x1}& D_{2y1}\end{array}\right|}\hfill \end{array}$$

(A82)

$$\begin{array}{ccc}\hfill P_{iy}& =& \hfill \frac{\left|\begin{array}{cc}{D}_{1xy}& D_{1y1}\\ D_{2xy}& D_{2y1}\end{array}\right|}{\left|\begin{array}{cc}{D}_{1x1}& D_{1y1}\\ D_{2x1}& D_{2y1}\end{array}\right|}\hfill \end{array}$$

(A83)

Where:

$$\begin{array}{ccc}\hfill D_{1xy}& =& \hfill \left|\begin{array}{cc}{P}_{1b1x}& P_{1b1y}\\ P_{1b2x}& P_{1b2y}\end{array}\right|\hfill \end{array}$$

(A84)

$$\begin{array}{ccc}\hfill D_{2xy}& =& \hfill \left|\begin{array}{cc}{P}_{2b1x}& P_{2b1y}\\ P_{2b2x}& P_{2b2y}\end{array}\right|\hfill \end{array}$$

(A85)

And:

$$\begin{array}{ccc}\hfill D_{1x1}& =& \hfill \left|\begin{array}{cc}{P}_{1b1x}& 1\\ P_{1b2x}& 1\end{array}\right|\hfill \end{array}$$

(A86)

$$\begin{array}{ccc}\hfill D_{1y1}& =& \hfill \left|\begin{array}{cc}{P}_{1b1y}& 1\\ P_{1b2y}& 1\end{array}\right|\hfill \end{array}$$

(A87)

$$\begin{array}{ccc}\hfill D_{2x1}& =& \hfill \left|\begin{array}{cc}{P}_{2b1x}& 1\\ P_{2b2x}& 1\end{array}\right|\hfill \end{array}$$

(A88)

$$\begin{array}{ccc}\hfill D_{2y1}& =& \hfill \left|\begin{array}{cc}{P}_{2b1y}& 1\\ P_{2b2y}& 1\end{array}\right|\hfill \end{array}$$

(A89)

The reference frames on the base $P_O$ and the platform $P_{M0}$ are are related to the ankle kinematic chain reference frames $P_{Oi}$ and $P_{M0i}$ by:

$$\begin{array}{ccc}\hfill P_{Oi}& =& \hfill P_O+P_{ixy}\hfill \end{array}$$

(A90)

$$\begin{array}{ccc}\hfill P_{M0i}& =& \hfill P_{M0}+P_{ixy}\hfill \end{array}$$

(A91)

We get values for $P_{Oi}$, $P_{M0i}$, $r_{1C}$, $r_{2C}$, ${\widehat{\omega }}_1$, and ${\widehat{\omega }}_2$. The kinematic chain, start from $P_{Oi}$, followed by a hinge joint at $r_{1C}$, then a perpendicular segment to a hinge joint at $r_{2C}$, ending in the platform center $PM0i$. In summary, we obtained a simplified model for the shank-ankle-foot kinematic chain.

Appendix B.1. Ankle Forward Kinematics.

There are two reference systems, the base with origin $P_O=(0,0,0)$ and the platform, with origin $P_M=(0,0,zp)$ at the platform center. The initial points are in Table A2.

Table A2. Platform’s initial position.

Platform radius ${\mathit{r}}_{\mathit{p}}$	Origin distance ${\mathit{z}}_{\mathit{p}}$	Initial Point ${\mathit{A}}_0$	Initial Point ${\mathit{B}}_0$	Initial Point ${\mathit{C}}_0$
11.57	17.62	(-11.57, 0.0000, -17.62)	(5.785, 10.02, -17.62)	(5.785, -10.02, -17.62)

* units in cm.

Table A2. Platform’s initial position.

Platform radius ${\mathit{r}}_{\mathit{p}}$	Origin distance ${\mathit{z}}_{\mathit{p}}$	Initial Point ${\mathit{A}}_0$	Initial Point ${\mathit{B}}_0$	Initial Point ${\mathit{C}}_0$
11.57	17.62	(-11.57, 0.0000, -17.62)	(5.785, 10.02, -17.62)	(5.785, -10.02, -17.62)

* units in cm.

The point $P{M}_0$ is the mean value from $A_0$, $B_0$, and $C_0$. The mean values from the sample anthropometric measurements the resulting values for $r_1$, $M_1$, $M_2$ are in Table A3.

Table A3. Talocrural axis approximation.

Foot width	${\mathit{r}}_1$	${\mathit{M}}_1$	${\mathit{M}}_2$
96.42	(0,0,-107.81)	(5.359, -50.99, -116.9)	(-4.565, 43.44, -100.1)

¹ units in mm.

Table A3. Talocrural axis approximation.

Foot width	${\mathit{r}}_1$	${\mathit{M}}_1$	${\mathit{M}}_2$
96.42	(0,0,-107.81)	(5.359, -50.99, -116.9)	(-4.565, 43.44, -100.1)

¹ units in mm.

And, for $r_2$, $N_1$, $N_2$ are in Table A4.

Table A4. Subtalar axis approximation.

Foot length	${\mathit{r}}_2$	${\mathit{N}}_1$	${\mathit{N}}_2$
266.5	(3.2803, 0, -111.58)	(-95.07, -27.70, -197.1)	(101.6, 27.70, -26.09)

¹ units in mm.

Table A4. Subtalar axis approximation.

Foot length	${\mathit{r}}_2$	${\mathit{N}}_1$	${\mathit{N}}_2$
266.5	(3.2803, 0, -111.58)	(-95.07, -27.70, -197.1)	(101.6, 27.70, -26.09)

¹ units in mm.

In this context, we have the points $r_1$, $M_1$, and $M_2$ pertaining to the talocrural axis. Also, $r_2$, $N_1$, and $N_2$ pertain to the subtalar axis. We used the mean values for ${\alpha }_1$, ${\alpha }_2$, ${\beta }_1$, ${\beta }_1$ from Table A5.

Table A5. Mean values for the axis attitude.

Angle Identifier	Mean (sexagesimal)
${\alpha }_1$	$6^{\circ }$
${\alpha }_2$	${80}^{\circ }$
${\beta }_1$	$-{23}^{\circ }$
${\beta }_2$	${41}^{\circ }$

Table A5. Mean values for the axis attitude.

Angle Identifier	Mean (sexagesimal)
${\alpha }_1$	$6^{\circ }$
${\alpha }_2$	${80}^{\circ }$
${\beta }_1$	$-{23}^{\circ }$
${\beta }_2$	${41}^{\circ }$

The resulting vectors for the product of exponentials matrix representation are in Table A6.

Table A6. Vectors for the product of exponentials matrices.

Axis and point	${\widehat{\mathit{\omega }}}_{\mathit{i}}$	${\widehat{\mathit{\nu }}}_{\mathit{i}}$
Talocrural on $r_1$	(-0.103, 0.979, 0.174)	(106., 11.1, 0)
Subtalar on $r_2$	0.738, 0.208, 0.642)	(23.2, -84.5, 0.682)

Table A6. Vectors for the product of exponentials matrices.

Axis and point	${\widehat{\mathit{\omega }}}_{\mathit{i}}$	${\widehat{\mathit{\nu }}}_{\mathit{i}}$
Talocrural on $r_1$	(-0.103, 0.979, 0.174)	(106., 11.1, 0)
Subtalar on $r_2$	0.738, 0.208, 0.642)	(23.2, -84.5, 0.682)

The components $R_{ij}$ of the rotation matrix are:

$$\begin{array}{cc}\hfill R_{11}& =\left(0.426cos\left({\theta }_1\right)-0.437sin\left({\theta }_1\right)+0.028\right)cos\left({\theta }_2\right)+\left(0.06cos\left({\theta }_1\right)-0.315sin\left({\theta }_1\right)-0.06\right)sin\left({\theta }_2\right)\hfill \\ \hfill & +0.563cos\left({\theta }_1\right)+0.437sin\left({\theta }_1\right)-0.018\hfill \\ \hfill R_{12}& =\left(-0.057cos\left({\theta }_1\right)-0.296sin\left({\theta }_1\right)-0.095\right)cos\left({\theta }_2\right)+\left(-0.621cos\left({\theta }_1\right)+0.723sin\left({\theta }_1\right)-0.019\right)sin\left({\theta }_2\right)\hfill \\ \hfill & +0.158cos\left({\theta }_1\right)+0.123sin\left({\theta }_1\right)-0.005\hfill \\ \hfill R_{13}& =\left(-0.471cos\left({\theta }_1\right)+0.599sin\left({\theta }_1\right)-0\right)cos\left({\theta }_2\right)+\left(0.131cos\left({\theta }_1\right)+0.128sin\left({\theta }_1\right)+0.076\right)sin\left({\theta }_2\right)\hfill \\ \hfill & +0.489cos\left({\theta }_1\right)+0.38sin\left({\theta }_1\right)-0.015\hfill \\ \hfill R_{21}& =\left(0.12cos\left({\theta }_1\right)+0.03sin\left({\theta }_1\right)-0.273\right)cos\left({\theta }_2\right)+\left(0.061cos\left({\theta }_1\right)-0.021sin\left({\theta }_1\right)+0.58\right)sin\left({\theta }_2\right)\hfill \\ \hfill & -0.019cos\left({\theta }_1\right)+0.143sin\left({\theta }_1\right)+0.172\hfill \\ \hfill R_{22}& =\left(0.046cos\left({\theta }_1\right)-0.04sin\left({\theta }_1\right)+0.91\right)cos\left({\theta }_2\right)+\left(-0.19cos\left({\theta }_1\right)-0.035sin\left({\theta }_1\right)+0.19\right)sin\left({\theta }_2\right)\hfill \\ \hfill & -0.005cos\left({\theta }_1\right)+0.04sin\left({\theta }_1\right)+0.048\hfill \\ \hfill R_{23}& =\left(-0.153cos\left({\theta }_1\right)-0.021sin\left({\theta }_1\right)+0.019\right)cos\left({\theta }_2\right)+\left(-0.009cos\left({\theta }_1\right)+0.036sin\left({\theta }_1\right)-0.729\right)sin\left({\theta }_2\right)\hfill \\ \hfill & -0.016cos\left({\theta }_1\right)+0.124sin\left({\theta }_1\right)+0.15\hfill \\ \hfill R_{31}& =\left(-0.425cos\left({\theta }_1\right)-0.429sin\left({\theta }_1\right)-0.048\right)cos\left({\theta }_2\right)+\left(-0.31cos\left({\theta }_1\right)-0.066sin\left({\theta }_1\right)+0.102\right)sin\left({\theta }_2\right)\hfill \\ \hfill & +0.443cos\left({\theta }_1\right)-0.549sin\left({\theta }_1\right)+0.03\hfill \\ \hfill R_{32}& =\left(-0.2948cos\left({\theta }_1\right)+0.051sin\left({\theta }_1\right)+0.161\right)cos\left({\theta }_2\right)+\left(0.704cos\left({\theta }_1\right)+0.628sin\left({\theta }_1\right)+0.033\right)sin\left({\theta }_2\right)\hfill \\ \hfill & +0.124cos\left({\theta }_1\right)-0.154sin\left({\theta }_1\right)+0.008\hfill \\ \hfill R_{33}& =\left(0.584cos\left({\theta }_1\right)+0.477sin\left({\theta }_1\right)+0.003\right)cos\left({\theta }_2\right)+\left(0.129cos\left({\theta }_1\right)-0.127sin\left({\theta }_1\right)-0.129\right)sin\left({\theta }_2\right)\hfill \\ \hfill & +0.385cos\left({\theta }_1\right)-0.477sin\left({\theta }_1\right)+0.026\hfill \end{array}$$

(A92)

The components of the translation vector ${\tau }_A$ for A reference point are:

$$\begin{array}{cc}\hfill {\tau }_{1A}& =-95.442cos\left({\theta }_1\right)-80.284sin\left({\theta }_1\right)+3.283+\left(-20.251cos\left({\theta }_1\right)+13.324sin\left({\theta }_1\right)-3.099\right)cos\left({\theta }_2\right)\hfill \\ \hfill & +\left(-15.735cos\left({\theta }_1\right)+29.205sin\left({\theta }_1\right)+2.305\right)sin\left({\theta }_2\right)\hfill \\ \hfill {\tau }_{2A}& =4.364cos\left({\theta }_1\right)-24.931sin\left({\theta }_1\right)-31.244+\left(-4.402cos\left({\theta }_1\right)-2.196sin\left({\theta }_1\right)+31.281\right)cos\left({\theta }_2\right)\hfill \\ \hfill & +\left(-6.729cos\left({\theta }_1\right)+0.214sin\left({\theta }_1\right)-21.937\right)sin\left({\theta }_2\right)\hfill \\ \hfill {\tau }_{3A}& =-81.197cos\left({\theta }_1\right)+93.028sin\left({\theta }_1\right)-113.349+\left(12.824cos\left({\theta }_1\right)+20.287sin\left({\theta }_1\right)+5.546\right)cos\left({\theta }_2\right)\hfill \\ \hfill & +\left(28.626cos\left({\theta }_1\right)+l16.104sin\left({\theta }_1\right)-3.889\right)sin\left({\theta }_2\right)\hfill \end{array}$$

(A93)

The group of movements for the point A are given by.

$$g_A\left({\theta }_1,{\theta }_2\right)=\left[\begin{array}{cc}{R}_T& {\tau }_A\\ 0_{1\times 3}& 1\end{array}\right]$$

(A94)

Where the rotation matrix is:

$$R_T=\left[\begin{array}{ccc}{R}_{11}& R_{12}& R_{13}\\ R_{21}& R_{22}& R_{23}\\ R_{31}& R_{32}& R_{33}\end{array}\right]$$

(A95)

And the translation vector is:

$${\widehat{\tau }}_A=\left[\begin{array}{c}{\tau }_{1A}\\ {\tau }_{2A}\\ {\tau }_{3A}\end{array}\right]$$

(A96)

With the range of movement for ${\theta }_1\in \left[-20,20\right]$, and ${\theta }_2\in \left[-10,10\right]$, we get the group of movements for each vertex on the platform. We also apply the same transformation to the group of movements for B, C, and $P_M$.

We change the talocrural and subtalar axes attitude from the Table A7.

Table A7. Ankle axes attitude modifications.

Angle Identifier	Maximum	Minimum
${\alpha }_1$	${13}^{\circ }$	$-1^{\circ }$
${\alpha }_2$	${76}^{\circ }$	${84}^{\circ }$
${\beta }_1$	$-{12}^{\circ }$	$-{34}^{\circ }$
${\beta }_2$	${50}^{\circ }$	${32}^{\circ }$

Table A7. Ankle axes attitude modifications.

Angle Identifier	Maximum	Minimum
${\alpha }_1$	${13}^{\circ }$	$-1^{\circ }$
${\alpha }_2$	${76}^{\circ }$	${84}^{\circ }$
${\beta }_1$	$-{12}^{\circ }$	$-{34}^{\circ }$
${\beta }_2$	${50}^{\circ }$	${32}^{\circ }$

We show the resulting ${\widehat{\omega }}_1$ and ${\widehat{\omega }}_2$ values in the Table A8.

Table A8. Direction vectors of the ankle axes.

Group of values	${\widehat{\mathit{\omega }}}_1$	${\widehat{\mathit{\omega }}}_2$
Maximum	(0.0174, 0.994, 0.105)	(0.629, 0.208, 0.749)
Minimum	(-0.218, 0.945, 0.242)	(0.703, 0.559, 0.439)

Table A8. Direction vectors of the ankle axes.

Group of values	${\widehat{\mathit{\omega }}}_1$	${\widehat{\mathit{\omega }}}_2$
Maximum	(0.0174, 0.994, 0.105)	(0.629, 0.208, 0.749)
Minimum	(-0.218, 0.945, 0.242)	(0.703, 0.559, 0.439)

The resulting ${\widehat{\nu }}_1$ and ${\widehat{\nu }}_1$ are in Table A9.

Table A9. Momentum for maximum and minimum axis attitude.

Group of values	${\widehat{\mathit{\omega }}}_1$	${\widehat{\mathit{\omega }}}_2$
Maximum	(107., -1.87, 0)	(23.1, -72.7, 0.796)
Minimum	(102., 23.5, 0)	(62.7, -79.9, 1.48)

Table A9. Momentum for maximum and minimum axis attitude.

Group of values	${\widehat{\mathit{\omega }}}_1$	${\widehat{\mathit{\omega }}}_2$
Maximum	(107., -1.87, 0)	(23.1, -72.7, 0.796)
Minimum	(102., 23.5, 0)	(62.7, -79.9, 1.48)

We also show the data generated from reflection matrix, resulting in artificial data for the left foot. The resulting reflection matrix about the $x-z$ sagittal plane is:

$$R_{xz}=\left[\begin{array}{ccc}1& 0& 0\\ 0& -1& 0\\ 0& 0& 1\end{array}\right],$$

(A97)

The resulting reflected platform initial points are in Table A10.

Table A10. Left foot platform initial points.

${\mathbf{A}}_0$ left	${\mathbf{B}}_0$ left	${\mathbf{C}}_0$ left
(-115.7, 0, -176.18)	(57.849, 100.20, -176.18)	(57.849, -100.20, -176.18)

Table A10. Left foot platform initial points.

${\mathbf{A}}_0$ left	${\mathbf{B}}_0$ left	${\mathbf{C}}_0$ left
(-115.7, 0, -176.18)	(57.849, 100.20, -176.18)	(57.849, -100.20, -176.18)

The talocrural reflected points are in Table A11.

Table A11. Talocrural axis reflected points.

${\mathbf{M}}_1$ reflected	${\mathbf{M}}_2$ reflected
(5.3595, 50.992, -116.85)	(-4.5655, -43.438, -100.11)

Table A11. Talocrural axis reflected points.

${\mathbf{M}}_1$ reflected	${\mathbf{M}}_2$ reflected
(5.3595, 50.992, -116.85)	(-4.5655, -43.438, -100.11)

And for the subtalar axis the resulting reflected points are in Table A12.

Table A12. Subtalar axis reflected points.

${\mathbf{N}}_1$ reflected	${\mathbf{N}}_2$ reflected
(-95.071, 27.700, -197.08)	(101.63, -27.700, -26.088)

Table A12. Subtalar axis reflected points.

${\mathbf{N}}_1$ reflected	${\mathbf{N}}_2$ reflected
(-95.071, 27.700, -197.08)	(101.63, -27.700, -26.088)

Appendix B.2. Axes’ Common Perpendicular and Feet Points.

After characterizing the axes and its range of motion, we compute the common perpendicular line and the intersection points (feet) with the estimated axes. The feet points will be used for represent the shank-foot-ankle kinematic chain. In this subsection, we use the results from Table A6, Table A8, and Table A9. We compute the common perpendicular from the minimum, maximum and mean axis attitude from the Table A13.

Table A13. Axes’ vectors for minimum, mean and maximum values.

Axis attitude	${\widehat{\mathit{\omega }}}_1$	${\widehat{\mathit{\nu }}}_1$	${\widehat{\mathit{\omega }}}_2$	${\widehat{\mathit{\nu }}}_2$
Minimum	(-0.218, 0.945, 0.242)	(102, 23.5, 0)	(0.703, 0.559, 0.439)	(62.7, -79.9, 1.48)
Mean	(-0.103, 0.979, 0.174)	(106., 11.1, 0)	(0.738, 0.208, 0.642)	(23.2, -84.5, 0.682)
Maximum	(0.0174, 0.994, 0.105)	(107, -1.87, 0)	(0.629, 0.208, 0.749)	(23.1, -72.7, 0.796)

Table A13. Axes’ vectors for minimum, mean and maximum values.

Axis attitude	${\widehat{\mathit{\omega }}}_1$	${\widehat{\mathit{\nu }}}_1$	${\widehat{\mathit{\omega }}}_2$	${\widehat{\mathit{\nu }}}_2$
Minimum	(-0.218, 0.945, 0.242)	(102, 23.5, 0)	(0.703, 0.559, 0.439)	(62.7, -79.9, 1.48)
Mean	(-0.103, 0.979, 0.174)	(106., 11.1, 0)	(0.738, 0.208, 0.642)	(23.2, -84.5, 0.682)
Maximum	(0.0174, 0.994, 0.105)	(107, -1.87, 0)	(0.629, 0.208, 0.749)	(23.1, -72.7, 0.796)

For each group of values, the resulting common perpendicular and feet are in Table A14.

Table A14. Common perpendicular and feet points.

Axis attitude	${\widehat{\mathit{\omega }}}_{\mathit{c}}$	${\widehat{\mathit{\nu }}}_{\mathit{c}}$	${\widehat{\mathit{r}}}_{1\mathit{c}}$	${\widehat{\mathit{r}}}_{2\mathit{c}}$
Minimum	(0.2796, 0.2658, -0.7862)	(26, -30, -2.4)	(-0.94, 4.1, -110)	(-32, -26, -19)
Mean	(0.5923, 0.1945, -0.7439)	(20, -64, -2.2)	(-0.24, 2.3, -110)	(-51, -14, -44)
Maximum	(0.7227, 0.05301, -0.6216)	(15, -79, 9.2)	(-0.24, -14, -110)	(-54, -18, -63)

Table A14. Common perpendicular and feet points.

Axis attitude	${\widehat{\mathit{\omega }}}_{\mathit{c}}$	${\widehat{\mathit{\nu }}}_{\mathit{c}}$	${\widehat{\mathit{r}}}_{1\mathit{c}}$	${\widehat{\mathit{r}}}_{2\mathit{c}}$
Minimum	(0.2796, 0.2658, -0.7862)	(26, -30, -2.4)	(-0.94, 4.1, -110)	(-32, -26, -19)
Mean	(0.5923, 0.1945, -0.7439)	(20, -64, -2.2)	(-0.24, 2.3, -110)	(-51, -14, -44)
Maximum	(0.7227, 0.05301, -0.6216)	(15, -79, 9.2)	(-0.24, -14, -110)	(-54, -18, -63)

Now, we have almost three different configurations, only by changing the axes’ attitude statistical values.

Appendix B.3. Subtalar Axis Range of Motion.

The range of motion for the subtalar axis is similar to the talocrural, the points for the negative angles are in the Table A15.

Table A15. Subtalar axis negative angles points.

Sample set i	Trajectory A	Trajectory B	Trajectory C
1	(-111.42 6.5037 -183.21)	(42.270 -118.34 -152.38)	(76.602 75.518 -189.75)
2	(-112.41 5.3489 -181.69)	(45.794 -115.24 -157.44)	(73.062 81.010 -187.46)
3	(-113.34 4.1020 -180.21)	(49.150 -111.81 -162.41)	(69.365 86.283 -184.92)
4	(-114.22 2.7668 -178.78)	(52.329 -108.09 -167.27)	(65.520 91.322 -182.13)
5	(-115.03 1.3469 -177.39)	(55.322 -104.06 -172.02)	(661.539 96.111 -179.10)

Table A15. Subtalar axis negative angles points.

Sample set i	Trajectory A	Trajectory B	Trajectory C
1	(-111.42 6.5037 -183.21)	(42.270 -118.34 -152.38)	(76.602 75.518 -189.75)
2	(-112.41 5.3489 -181.69)	(45.794 -115.24 -157.44)	(73.062 81.010 -187.46)
3	(-113.34 4.1020 -180.21)	(49.150 -111.81 -162.41)	(69.365 86.283 -184.92)
4	(-114.22 2.7668 -178.78)	(52.329 -108.09 -167.27)	(65.520 91.322 -182.13)
5	(-115.03 1.3469 -177.39)	(55.322 -104.06 -172.02)	(661.539 96.111 -179.10)

The points difference in this case is on the Table A16.

Table A16. Points differences for the subtalar negative angles.

Sample set i	$\mathit{B}-\mathit{A}$	$\mathit{C}-\mathit{A}$
1	(153.69 -124.85 30.829)	(188.02 69.014 -6.5481)
2	(158.20 -120.59 24.255)	(185.47 75.661 -5.7695)
3	(162.49 -115.92 17.808)	(182.71 82.181 -4.7012)
4	(166.54 -110.85 11.507)	(179.74 88.555 -3.3472)
5	(170.35 -105.41 5.3694)	(176.56 94.764 -1.7107)

Table A16. Points differences for the subtalar negative angles.

Sample set i	$\mathit{B}-\mathit{A}$	$\mathit{C}-\mathit{A}$
1	(153.69 -124.85 30.829)	(188.02 69.014 -6.5481)
2	(158.20 -120.59 24.255)	(185.47 75.661 -5.7695)
3	(162.49 -115.92 17.808)	(182.71 82.181 -4.7012)
4	(166.54 -110.85 11.507)	(179.74 88.555 -3.3472)
5	(170.35 -105.41 5.3694)	(176.56 94.764 -1.7107)

The components for the rotation matrices in this case are in the Table A17.

Table A17. Rotation matrices components for the subtalar axis negative angles.

Sample set i	Vector $\widehat{\mathit{n}}$	Vector $\widehat{\mathit{s}}$	Vector $\widehat{\mathit{a}}$
1	(-0.0377 0.196 0.980)	(0.984 -0.161 0.0700)	(0.171 0.967 -0.187)
2	(-0.0328 0.156 0.987)	(0.990 -0.129 0.0533)	(0.136 0.979 -0.150)
3	(-0.0264 0.116 0.993)	(0.995 -0.0972 0.0378)	(0.101 0.989 -0.112)
4	(-0.0186 0.0755 0.997)	(0.998 -0.0642 0.0235)	(0.0658 0.995 -0.0741)
5	(-0.00945 0.0356 0.999)	(0.999 -0.0307 0.0105)	(0.0310 0.999 -0.03534)

Table A17. Rotation matrices components for the subtalar axis negative angles.

Sample set i	Vector $\widehat{\mathit{n}}$	Vector $\widehat{\mathit{s}}$	Vector $\widehat{\mathit{a}}$
1	(-0.0377 0.196 0.980)	(0.984 -0.161 0.0700)	(0.171 0.967 -0.187)
2	(-0.0328 0.156 0.987)	(0.990 -0.129 0.0533)	(0.136 0.979 -0.150)
3	(-0.0264 0.116 0.993)	(0.995 -0.0972 0.0378)	(0.101 0.989 -0.112)
4	(-0.0186 0.0755 0.997)	(0.998 -0.0642 0.0235)	(0.0658 0.995 -0.0741)
5	(-0.00945 0.0356 0.999)	(0.999 -0.0307 0.0105)	(0.0310 0.999 -0.03534)

The rotation matrices for the subtalar axis negative angles are in the Table A18.

Table A18. Subtalar axis negative angles rotation matrices.

Sample set i	Rotation matrix	Angle (Rad)	Angle (Deg)
1	$\left(\begin{array}{ccc}\hfill 0.98449& \hfill 0.17132& \hfill -0.037672\\ \hfill -0.16086& \hfill 0.96740& \hfill 0.19561\\ \hfill 0.069956& \hfill -0.18652& \hfill 0.97996\end{array}\right)$	0.26181	15.001
2	$\left(\begin{array}{ccc}\hfill 0.99016& \hfill 0.13607& \hfill -0.032763\\ \hfill -0.12943& \hfill 0.97930& \hfill 0.15560\\ \hfill 0.053258& \hfill -0.14983& \hfill 0.98728\end{array}\right)$	0.20838	11.939
3	$\left(\begin{array}{ccc}\hfill 0.99455& \hfill 0.10087& \hfill -0.026412\\ \hfill -0.097195& \hfill 0.98854& \hfill 0.11552\\ \hfill 0.037762& \hfill -0.11232& \hfill 0.99295\end{array}\right)$	0.15495	8.8777
4	$\left(\begin{array}{ccc}\hfill 0.99766& \hfill 0.065822& \hfill -0.018631\\ \hfill -0.064242& \hfill 0.99507& \hfill 0.075498\\ \hfill 0.023509& \hfill -0.074124& \hfill 0.99697\end{array}\right)$	0.10153	5.8173
5	$\left(\begin{array}{ccc}\hfill 0.99947& \hfill 0.031024& \hfill -0.0094457\\ \hfill -0.030669& \hfill 0.99889& \hfill 0.035639\\ \hfill 0.010541& \hfill -0.035331& \hfill 0.99932\end{array}\right)$	0.048085	2.7551

Table A18. Subtalar axis negative angles rotation matrices.

Sample set i	Rotation matrix	Angle (Rad)	Angle (Deg)
1	$\left(\begin{array}{ccc}\hfill 0.98449& \hfill 0.17132& \hfill -0.037672\\ \hfill -0.16086& \hfill 0.96740& \hfill 0.19561\\ \hfill 0.069956& \hfill -0.18652& \hfill 0.97996\end{array}\right)$	0.26181	15.001
2	$\left(\begin{array}{ccc}\hfill 0.99016& \hfill 0.13607& \hfill -0.032763\\ \hfill -0.12943& \hfill 0.97930& \hfill 0.15560\\ \hfill 0.053258& \hfill -0.14983& \hfill 0.98728\end{array}\right)$	0.20838	11.939
3	$\left(\begin{array}{ccc}\hfill 0.99455& \hfill 0.10087& \hfill -0.026412\\ \hfill -0.097195& \hfill 0.98854& \hfill 0.11552\\ \hfill 0.037762& \hfill -0.11232& \hfill 0.99295\end{array}\right)$	0.15495	8.8777
4	$\left(\begin{array}{ccc}\hfill 0.99766& \hfill 0.065822& \hfill -0.018631\\ \hfill -0.064242& \hfill 0.99507& \hfill 0.075498\\ \hfill 0.023509& \hfill -0.074124& \hfill 0.99697\end{array}\right)$	0.10153	5.8173
5	$\left(\begin{array}{ccc}\hfill 0.99947& \hfill 0.031024& \hfill -0.0094457\\ \hfill -0.030669& \hfill 0.99889& \hfill 0.035639\\ \hfill 0.010541& \hfill -0.035331& \hfill 0.99932\end{array}\right)$	0.048085	2.7551

The greatest of the negative angles is $-{15}^{\circ }$. The Table A19 shows the points for the subtalar positive angles.

Table A19. Subtalar axis positive angles trajectories.

Sample set i	Trajectory A	Trajectory B	Trajectory C
1	(-116.32 -1.4090 -175.01)	(60.212 -96.112 -180.22)	(54.066 104.06 -173.08)
2	(-116.94 -3.0434 -173.76)	(62.640 -91.324 -184.56)	(49.766 108.08 -169.43)
3	(-117.49 -4.7460 -172.58)	(64.851 -86.289 -188.74)	(45.374 111.81 -165.59)
4	(-117.97 -6.5117 -171.46)	(66.840 -81.021 -192.73)	(40.903 115.23 -161.55)
5	(-118.37 -8.3355 -170.40)	(68.602 -75.534 -196.54)	(36.366 118.33 -157.34)

Table A19. Subtalar axis positive angles trajectories.

Sample set i	Trajectory A	Trajectory B	Trajectory C
1	(-116.32 -1.4090 -175.01)	(60.212 -96.112 -180.22)	(54.066 104.06 -173.08)
2	(-116.94 -3.0434 -173.76)	(62.640 -91.324 -184.56)	(49.766 108.08 -169.43)
3	(-117.49 -4.7460 -172.58)	(64.851 -86.289 -188.74)	(45.374 111.81 -165.59)
4	(-117.97 -6.5117 -171.46)	(66.840 -81.021 -192.73)	(40.903 115.23 -161.55)
5	(-118.37 -8.3355 -170.40)	(68.602 -75.534 -196.54)	(36.366 118.33 -157.34)

The points differences are in the Table A20.

Table A20. Points difference for the subtalar axis positive angles.

Sample set i	$\mathit{B}-\mathit{A}$	$\mathit{C}-\mathit{A}$
1	(176.53 -94.703 -5.2102)	(170.38 105.47 1.9314)
2	(179.58 -88.281 -10.798)	(166.70 111.13 4.3313)
3	(182.34 -81.543 -16.158)	(162.86 116.55 6.9919)
4	(184.81 -74.509 -21.276)	(158.87 121.74 9.9048)
5	(186.97 -67.199 -26.137)	(154.74 126.66 13.062)

Table A20. Points difference for the subtalar axis positive angles.

Sample set i	$\mathit{B}-\mathit{A}$	$\mathit{C}-\mathit{A}$
1	(176.53 -94.703 -5.2102)	(170.38 105.47 1.9314)
2	(179.58 -88.281 -10.798)	(166.70 111.13 4.3313)
3	(182.34 -81.543 -16.158)	(162.86 116.55 6.9919)
4	(184.81 -74.509 -21.276)	(158.87 121.74 9.9048)
5	(186.97 -67.199 -26.137)	(154.74 126.66 13.062)

The rotation matrices components are in the Table A21.

Table A21. Components of the rotation matrices for the subtalar axis positive angles.

Sample set i	Vector $\widehat{\mathit{n}}$	Vector $\widehat{\mathit{s}}$	Vector $\widehat{\mathit{a}}$
1	(0.0105 -0.0353 0.999)	(0.999 0.0310 -0.00945)	(-0.0307 0.999 0.0356)
2	(0.0235 -0.0741 0.997)	(0.998 0.0658 -0.0186)	(-0.0642 0.995 0.07553)
3	(0.0378 -0.112 0.993)	(0.995 0.101 -0.0264)	(-0.0972 0.989 0.116)
4	(0.0533 -0.150 0.987)	(0.990 0.136 -0.0328)	(-0.129 0.979 0.156)
5	(0.0700 -0.187 0.980)	(0.984 0.171 -0.0377)	(-0.161 0.967 0.196)

Table A21. Components of the rotation matrices for the subtalar axis positive angles.

Sample set i	Vector $\widehat{\mathit{n}}$	Vector $\widehat{\mathit{s}}$	Vector $\widehat{\mathit{a}}$
1	(0.0105 -0.0353 0.999)	(0.999 0.0310 -0.00945)	(-0.0307 0.999 0.0356)
2	(0.0235 -0.0741 0.997)	(0.998 0.0658 -0.0186)	(-0.0642 0.995 0.07553)
3	(0.0378 -0.112 0.993)	(0.995 0.101 -0.0264)	(-0.0972 0.989 0.116)
4	(0.0533 -0.150 0.987)	(0.990 0.136 -0.0328)	(-0.129 0.979 0.156)
5	(0.0700 -0.187 0.980)	(0.984 0.171 -0.0377)	(-0.161 0.967 0.196)

The rotation matrices for the subtalar axis positive angles are in Table A22.

Table A22. Rotation matrices for the subtalar axis positive angles.

Sample set i	Rotation matrix	Angle (Rad)	Angle (Deg)
1	$\left(\begin{array}{ccc}\hfill 0.99947& \hfill -0.030669& \hfill 0.010542\\ \hfill 0.031024& \hfill 0.99889& \hfill -0.035329\\ \hfill -0.0094465& \hfill 0.035638& \hfill 0.99932\end{array}\right)$	0.048085	2.7551
2	$\left(\begin{array}{ccc}\hfill 0.99766& \hfill -0.064241& \hfill 0.023509\\ \hfill 0.065821& \hfill 0.99507& \hfill -0.074124\\ \hfill -0.018631& \hfill 0.075498& \hfill 0.99697\end{array}\right)$	0.10153	5.8173
3	$\left(\begin{array}{ccc}\hfill 0.99455& \hfill -0.097194& \hfill 0.037759\\ \hfill 0.10087& \hfill 0.98854& \hfill -0.11233\\ \hfill -0.026409& \hfill 0.11552& \hfill 0.99295\end{array}\right)$	0.15495	8.8777
4	$\left(\begin{array}{ccc}\hfill 0.99016& \hfill -0.12943& \hfill 0.053257\\ \hfill 0.13607& \hfill 0.97930& \hfill -0.14983\\ \hfill -0.032762& \hfill 0.15560& \hfill 0.98728\end{array}\right)$	0.20838	11.939
5	$\left(\begin{array}{ccc}\hfill 0.98449& \hfill -0.16086& \hfill 0.069953\\ \hfill 0.17132& \hfill 0.96740& \hfill -0.18652\\ \hfill -0.037669& \hfill 0.19561& \hfill 0.97996\end{array}\right)$	0.26180	15.000

Table A22. Rotation matrices for the subtalar axis positive angles.

Sample set i	Rotation matrix	Angle (Rad)	Angle (Deg)
1	$\left(\begin{array}{ccc}\hfill 0.99947& \hfill -0.030669& \hfill 0.010542\\ \hfill 0.031024& \hfill 0.99889& \hfill -0.035329\\ \hfill -0.0094465& \hfill 0.035638& \hfill 0.99932\end{array}\right)$	0.048085	2.7551
2	$\left(\begin{array}{ccc}\hfill 0.99766& \hfill -0.064241& \hfill 0.023509\\ \hfill 0.065821& \hfill 0.99507& \hfill -0.074124\\ \hfill -0.018631& \hfill 0.075498& \hfill 0.99697\end{array}\right)$	0.10153	5.8173
3	$\left(\begin{array}{ccc}\hfill 0.99455& \hfill -0.097194& \hfill 0.037759\\ \hfill 0.10087& \hfill 0.98854& \hfill -0.11233\\ \hfill -0.026409& \hfill 0.11552& \hfill 0.99295\end{array}\right)$	0.15495	8.8777
4	$\left(\begin{array}{ccc}\hfill 0.99016& \hfill -0.12943& \hfill 0.053257\\ \hfill 0.13607& \hfill 0.97930& \hfill -0.14983\\ \hfill -0.032762& \hfill 0.15560& \hfill 0.98728\end{array}\right)$	0.20838	11.939
5	$\left(\begin{array}{ccc}\hfill 0.98449& \hfill -0.16086& \hfill 0.069953\\ \hfill 0.17132& \hfill 0.96740& \hfill -0.18652\\ \hfill -0.037669& \hfill 0.19561& \hfill 0.97996\end{array}\right)$	0.26180	15.000

The greatest positive value is $+{15}^{\circ }$. Then the range of motion is ${30}^{\circ }$.