2.1.8. Tetrahedrons B and C Apexes

In this subsection, we show that, by knowing Ap, the point B<sup>p</sup> needs two sensors to be found. To determine the result of the tetrahedron TB, we consider the base of a triangle △ - B1, B3, A<sup>p</sup> in Figure 11a .

We compute the angle *α* from the XY plane to a normal vector *n* ˆ *ApB*:

$$\mathfrak{H}\_{ApB} = \frac{(B\_3 - Ap) \times (B\_3 - Ap)}{|| (B\_3 - Ap) \times (B\_3 - Ap) ||} \, \text{} \tag{35}$$

and, the angle *α* is:

$$\mathfrak{a} = \operatorname{acos}(\mathfrak{h}\_{ApB} \cdot \mathfrak{h}\_{\mathbb{Z}}) . \tag{36}$$

where *n* ˆ *<sup>z</sup>* is the unitary vector normal to the XY plane.

The tetrahedron sides are the lengths *lB*1 = k*Bp* − *B*1k, *lB*3 = k*Bp* − *B*3k, and *dApBp* = k*Bp* − *Ap*k. The rotation axis is in the direction *B*1 − *B*3. The *Bpr* is *Bps* rotated *α* in angle about this axis. In Figure 11b, we show how to find the *Bpr* apex, similarly to

that of a tetrahedron TA. Finally, when *Bpr* is found, the contrary rotation about the axis *B*1 − *B*3 gives the *Bps*.

**Figure 11.** Rotation of *α* angle about the axis *B*<sup>1</sup> − *B*3: (**a**) original tetrahedron, *T<sup>B</sup>* (**b**) rotated tetrahedron.

There are two possible apex values: *Bps*1 over, and *Bps*2 below of the XY plane. We show the *Bpr* apex below the XY plane in Figure 12.

**Figure 12.** Finding the apex Bpr.

We use the same method to solve the T<sup>C</sup> apex. For the correct apex selection, the condition when the side of the platform distance *dCpBp* is:

$$d\_{\mathbb{C}pBp} = \|Bps - \mathbb{C}ps\|.\tag{37}$$

2.1.9. Procedure for Found Platform Positions

We must fix the shank and the foot to the base and platform. Then we mark the MMP and the MLP. To do so, we design a detachable reference point from module A. Initially, we attach the foot and the shank to the device, and then we mark and record the MMP and MLP; Figure 13 shows the detailed view.

**Figure 13.** Adjusting the foot, the shank and the Most Medial Point reference.

We compute the platform position from the seven sensor lengths. The main steps for capture data series are:


First, we capture the sensor lengths by activating a button in the computer software. Every time, we compute the absolute difference from IMU2 to the IMU1 readings. If the differences are constants, then there is no platform and base relative movement. We compute the jerk by relative acceleration differentiation. The data capturing process ends when the acceleration change crosses zero. Jerk changes activate the capture of IMU data.

The symbolic equations to find Ap, Bp, and C<sup>p</sup> from the captured data, were found by the SageMath CAS. By using the prototype dimensions and the sensor lengths, we compute the platform's position and orientation. Here, the origin is from the initial DWS lengths *lMi*<sup>0</sup> , where M is the module A, B, or C; and *i* is the sensor number *i* = 1, 2, 3.

After MLP and MMP registering, we attach the apex of module A to the platform, define the sagittal plane perpendicular to the ABC base plane, and intersect point A. By implementing the trilateration method mentioned before, we compute the points *A*0, *B*0, and *C*0.

Figure 14a illustrate the point positions with the device in the initial portable configuration. The apexes' computation are in Figure 14b–d.

**Figure 14.** Computed positions from sensor lengths at portable configuration: (**a**) the rest position, (**b**) apex A, (**c**) apex B, and (**d**) apex C.

#### 2.1.10. Computing the Axis Position and Direction

From the anthropometric values [40], we put a mean model in the Turmell-Meter. The TC axis will be defined by *M*1<sup>0</sup> and *M*2<sup>0</sup> . The sagittal plane intersection with the *M*1*M*<sup>2</sup> segment is *r*1. For example, the TC axis approximation is computed from most lateral point (MLP), the most medial point (MMP) and L, K, P, and O:

$$M\_{1\_0} = MLP - [L, 0, K],\tag{38}$$

and:

$$M\_{\text{2}\_0} = MMP - [-P, 0, O]\_{\prime} \tag{39}$$

from these values, we solve for *r*<sup>1</sup> from the plane *y* = 0 intersection with the line *LTC*:

$$L\_{T\mathbb{C}} = V - M1 = \rho \frac{(M\_2 - M\_1)}{\|M\_2 - M\_1\|} \, \tag{40}$$

where V is a point pertaining to *LTC*.

The ST axis sagittal intersection *r*<sup>2</sup> initial point is:

$$r\_2 = r\_1 + \frac{\sqrt{2}}{2} [\mathbb{Q}, 0, \mathbb{Q}].\tag{41}$$

Here, *r*<sup>1</sup> and *r*<sup>2</sup> are reference values computed from the previously mentioned anthropometric mean values. Such initial points are for reference, comparison and validation of the trilateration and regression method. The tracked trajectory data set is processed offline. We use the least squares normal vector to the plane, this direction is similar to the circle approximation. From here, we compute the TC axis first, and then the ST axis. To do so, we compute the TC axis position by employing dorsiflexion and plantarflexion, because the TC axis is the most dominant in such movements. The method used is circle fitting in a plane containing the trajectory points. A further model refinement can be made with optimization, and machine learning methods, such as gradient descent and the symbolic product of exponential formula.

First, we found the TC axis orientation *ω*<sup>1</sup> by registering several trajectories. For each trajectory, we have a list of data points *P* = [*x*, *y*, *z*], which pertain to a plane:

$$a\mathbf{x} + by + c\mathbf{z} + d = \mathbf{0},\tag{42}$$

where *a*, *b*, and *c* are the components of a direction vector perpendicular to a plane containing the points. Solving out for *z*, we have the system:

$$
\begin{bmatrix} x\_0 & y\_0 & 1 \\ x\_1 & y\_1 & 1 \\ \vdots & \vdots & \vdots \\ x\_{n-1} & y\_{n-1} & 1 \end{bmatrix} \begin{bmatrix} a \\ b \\ d \end{bmatrix} = -c \begin{bmatrix} z\_0 \\ z\_1 \\ \vdots \\ z\_{n-1} \end{bmatrix} \tag{43}
$$

which has the form:

$$A\mathfrak{x} = B \tag{44}$$

there are more equations than unknowns. From linear algebra and least squares we knew that the pseudo inverse is *A* <sup>+</sup> = (*A <sup>T</sup>A*) <sup>−</sup>1*A T* , then a normal vector is:

$$
\begin{bmatrix} a \\ b \\ d \end{bmatrix} = (A^T A)^{-1} A^T B \tag{45}
$$

Now, we compute *c* by replacing *a*, *b*, *d* in the plane equation, and finally we get *n*ˆ = [*a*, *b*, *c*] *T* . We found the angle between the normal plane and the X-Y plane, after knowing the normal vector by applying the Rodrigues' formula, *<sup>v</sup>*<sup>ˆ</sup> <sup>=</sup> *<sup>n</sup>*<sup>ˆ</sup> <sup>×</sup> <sup>ˆ</sup>*k*, with <sup>ˆ</sup>*<sup>k</sup>* = [0, 0, 1] *T*

$$P\_r = P\cos(\theta) + (\vartheta \times P)\sin(\theta) + \vartheta(\vartheta \cdot P)(1 - \cos\theta). \tag{46}$$

where *<sup>θ</sup>* <sup>=</sup> arccos *n*ˆ· ˆ*k* k*n*ˆk .

After this, we estimate the plane, and rotate all the data points onto the X-Y plane. We search for a circle in the X-Y plane, and rearrange the equation for least squares estimation by using a variable substitution.

$$\begin{array}{rcl}(\mathbf{x} - \mathbf{x}\_c)^2 + (\mathbf{y} - \mathbf{y}\_c)^2 &=& r^2\\(2\mathbf{x}\_c)\mathbf{x} + (2\mathbf{y}\_c)\mathbf{y} + (r^2 - \mathbf{x}\_c^2 - \mathbf{y}\_c^2) &=& \mathbf{x}^2 + \mathbf{y}^2\\c\_0\mathbf{x} + c\_1\mathbf{y} + c\_2 &=& \mathbf{x}^2 + \mathbf{y}^2\end{array} \tag{47}$$

where *c* = [*c*0, *c*1, *c*2] *<sup>T</sup>* with *c*<sup>0</sup> = 2*xc*, *c*<sup>1</sup> = 2*yc*, and *c*<sup>2</sup> = *r* <sup>2</sup> <sup>−</sup> *<sup>x</sup>* 2 *<sup>c</sup>* − *y* 2 *c* By taking the rotated points, *P<sup>r</sup>* we have a linear system:

$$
\begin{bmatrix} x\_0 & y\_0 & 1 \\ x\_1 & y\_1 & 1 \\ \vdots & \vdots & \vdots \\ x\_{n-1} & y\_{n-1} & 1 \end{bmatrix} \begin{bmatrix} c\_0 \\ c\_1 \\ c\_2 \end{bmatrix} = \begin{bmatrix} x\_0^2 + y\_0^2 \\ x\_1^2 + y\_1^2 \\ \vdots \\ x\_{n-1}^2 + y\_{n-1}^2 \end{bmatrix}. \tag{48}
$$

that has the form:

$$A\mathcal{c} = b \tag{49}$$

.

In this system, we have more equations than unknowns, then, we search for the c values that minimize the squared difference k*b* − *Ac*k 2 .

$$\underset{c \in \mathbb{R}^3}{\arg\min} \|b - Ac\|^2. \tag{50}$$

We found the center point *C<sup>p</sup>* = [*xc*, *yc*] and radius *r* by solving:

$$\begin{array}{ccccc} 2\mathfrak{x}\_{\mathfrak{c}} & = & c\_0 \\ 2\mathfrak{y}\_{\mathfrak{c}} & = & c\_1 \\ r^2 - \mathfrak{x}\_{\mathfrak{c}}^2 - \mathfrak{y}\_{\mathfrak{c}}^2 & = & c\_2 \end{array} \tag{51}$$

Finally, we apply a rotation to the center in respect to the original plane. This point pertains to the TC axis. For each trajectory A, B, C, we get three planes, and three centers, the TC line direction is parallel to the planes' normal vectors. The information is complete by determining the plane orientation.

The ST axis estimation is similar, but employs trajectories from inversion and eversion movements.

This is a basic estimation; by conducting optimization on the product of the exponential formula, we enhance the accuracy of the axis position estimation.

#### 2.1.11. Ankle Joint Movements as a Manifold

In this subsection, we explain how the centers *r*1, *r*<sup>2</sup> and directions *ω*1, *ω*<sup>2</sup> define a manifold representing the HAJ movements. The circle center points calculated pertain to the TC and ST axes; they are the initial data to fit the product of the exponential formula. In Figure 15a, we show the complete platform's center point manifold. It is topologically similar to a torus.

**Figure 15.** Simulation of the platform central point with variations in the mean statistical values: (**a**) platform's center point manifold, (**b**) manifold chart and a geodesic.

A manifold chart represents the range of motion limits, we show an example of the geodesic as a trajectory on the manifold in Figure 15b; this explains how to map ankle coordinates, and a straight trajectory with initial velocity and no external force action. We have the data necessary for the line intersection with the sagittal plane, the center points, and the direction gives a line:

$$
\mathfrak{p}\_l = \mathfrak{l}\_0 + \mathfrak{l}d\mathfrak{l} \tag{52}
$$

where *p*ˆ*<sup>l</sup>* is the parametric line, ˆ *l* is a parallel vector to it, ˆ *l*<sup>0</sup> is a known vector in such line, and *d* ∈ R, replacing the parametric equation in the plane equation:

$$(\mathfrak{p}\_l - \mathfrak{p}\_0) \cdot (\mathfrak{n}\_p) = 0,\tag{53}$$

where *p* ˆ <sup>0</sup> is a known vector in the plane, and *n* ˆ *<sup>p</sup>* is the plane's normal vector, solving for *d*, gives:

> *d* = (*p* ˆ <sup>0</sup> <sup>−</sup> <sup>ˆ</sup> *l*0) · *n* ˆ *p* ˆ *l* · *n* ˆ *p* , (54)

and replacing in the TC axis line equation:

$$
\sigma\_1 = \sigma\_1 + \omega\_1 d\_\prime \tag{55}
$$

where *r*<sup>1</sup> is the TC axis intersection with the sagittal plane. The point *c*<sup>1</sup> is the center, and the axis direction *ω*1, both were found by circle fitting. Furthermore, packing in six dimensional Plücker line coordinates, we have:

$$
\mathfrak{M}\_1 = r\_1 \times \omega\_1. \tag{56}
$$

and the *l*<sup>1</sup> six dimensional vector is:

$$l\_1 = [\omega\_{1\_x} : \omega\_{1\_y} : \omega\_{1\_z} : m\_{1\_x} : m\_{1\_y} : m\_{1\_z}].\tag{57}$$

We include those data for the PoE formula simulation and the manifold representation.

#### *2.2. Mechatronic System Design*

In this section, we design DWS to measure the lengths of the tetrahedron sides; they are arranged as structural parts. Their maximal length estimation is from the forward kinematics simulation. We design the shank attachment from the dimensions, proportions, and statistical data.

#### 2.2.1. Draw-Wire Sensor

We use flat springs. They are not exposed to a high load against gravity, and are in two or three concurrent groups. In Figure 16, we depict the design, composed of three 3D printed parts, potentiometer, flat spring, bolts, and nuts.

**Figure 16.** Draw-wire sensor design.

A two-coil winch drives the potentiometer; a flat spring retracts a wire attached to the winch. When we pull the wire, the spring retracts it. The value of each turn is from the nominal value of the potentiometer, R<sup>n</sup> = 2.2 kΩ, divided into ten turns, that is 220 Ω per turn. The diameter is D = 3.8 cm, the spring could be compressed in four turns. The maximal length is as follows:

$$l\_{\max} = 4 \cdot \mathbf{D} \cdot \boldsymbol{\pi} \tag{58}$$

Which is 47.75 cm approximately, this value is greater than *lmax* for all groups of movements.

#### 2.2.2. Mechanical Parts

The attachment on the calf has a size according to the simulation. We use the mesh model of a leg to guide the shape of the calf support, as in Figure 17a. We also scale and divide this structure into seven parts for 3D printing. An aluminum tube is the support structure, as in Figure 17b, and a neoprene band attaches the shank to the support with Velcro fabric.

**Figure 17.** Mechanical attachment: (**a**) calf support and (**b**) aluminum tube structure.

All the DWS modules are in a plate, the A module has three DWS, B and C modules has two DWS, as in Figure 18a. The design of the foot attachment is from standard measurements to adjust the foot's length and width, as in Figure 18b.

**Figure 18.** Base and platform: (**a**) DWS modules support (**b**) platform with foot's size adjustment.
