Upper limb joint kinematics using wearable magnetic and inertial measurement units: an anatomical calibration procedure based on bony landmark identification

The estimate of a consistent and clinically meaningful joint kinematics using wearable inertial and magnetic sensors requires a sensor-to-segment coordinate system calibration. State-of-the-art calibration procedures for the upper limb are based on functional movements and/or pre-determined postures, which are difficult to implement in subjects that have impaired mobility or are bedridden in acute units. The aim of this study was to develop and validate an alternative calibration procedure based on the direct identification of palpable anatomical landmarks (ALs) for an inertial and magnetic sensor-based upper limb movement analysis protocol. The proposed calibration procedure provides an estimate of three-dimensional shoulder/elbow angular kinematics and the linear trajectory of the wrist according to the standards proposed by the International Society of Biomechanics. The validity of the method was assessed against a camera-based optoelectronic system during uniaxial joint rotations and a reach-to-grasp task. Joint angular kinematics was found as characterised by a low-biased range of motion (<−2.6°), a low root mean square deviation (RMSD) (<4.4°) and a high waveform similarity coefficient (R2 > 0.995) with respect to the gold standard. Except for the cranio–caudal direction, the linear trajectory of the wrist was characterised by a low-biased range of motion (<11 mm) together with a low RMSD (8 mm) and high waveform similarity (R2 > 0.968). The proposed method enabled the estimation of reliable joint kinematics without requiring any active involvement of the patient during the calibration procedure, complying with the metrological standards and requirements of clinical movement analysis.

sound joint kinematics requires knowledge of the orientation of the MIMU-embedded coordinate system with respect to the anatomical coordinate system of the body segment on which the MIMU is strapped. This information is retrieved during a sensor-to-segment axes calibration procedure (simply called "anatomical calibration"), which can be implemented in different ways. With special reference to upper limb kinematics, thorax, upper arm and forearm anatomical coordinate systems are generally determined according to one of the following approaches proposed in the literature: (a) matching the MIMU coordinate system with the anatomical coordinate system through manual alignment of the MIMU case 7 ; (b) asking the subject to assume a given body segment configuration ("N-pose" or "T-pose") 8 ; (c) exploiting the direction of the angular velocity vector as measured during specific monoaxial rotations of a body segment together with the direction of the gravity vector as measured while keeping a segment's axis aligned with the vertical line 9 . While simple manual alignment of the MIMU case could be critical from a repeatability perspective, in cases of severe impaired mobility, applying the second and third approaches could be difficult as the subjects may be not able to perform functional calibration tasks or assume specific postures. In addition, all three approaches may be particularly critical with bedridden acute patients as they are usually equipped with medical devices and apparatus the positioning of which might jeopardise the alignment of the MIMU with the underlying bone and limit the joint range of motion.
An alternative to overcome the aforementioned limitations is to determine the direction of the anatomical axes on the base of the bone morphology, particularly, from the position of a few selected palpable anatomical landmarks (ALs). To the best of the authors' knowledge, a similar approach has been previously used for the estimation of lower limb joint angular kinematics 10 and scapula orientation tracking 11 .
The aim of this study was to develop a MIMU-based anatomical calibration procedure to estimate linear and angular kinematics of an upper limb kinematic model and assess the procedure's validity against a camera-based optoelectronic system.

Material and Methods
Anatomical calibration. Let g R s be a 3 × 3 orientation matrix describing the 3D orientation of the MIMU coordinate system s with respect to a common global coordinate system g: where the unit vectors u are the column elements of the matrix denoting the orientation of each single axis of s with respect to g. With reference to the example shown in Fig. 1, a MIMU is strapped on the subject's upper arm (sU) for collecting movement data. An extra MIMU is mounted on a light-weight calliper-like pointing device (sC), made of aluminium and having adjustable length, such that its x-axis is aligned with the axis joining the tips of the calliper's arms ( u g x sC ). The orientation of the axis defined by two pointed ALs (in Fig. 1, the line joining the medial epicondyle, ME, to the lateral epicondyle, LE) with respect to g can be directly retrieved from the first column of the 3 × 3 direction cosine matrix output by the sensor: and subsequently expressed in the MIMU coordinate system (sU) used to track the motion of the segment through the rigid, time-invariant transformation 10 .
Note that R (0) g sU and  ME LE g are simultaneously collected and the rigid transformation is allowed as the two MIMUs share the same global coordinate system. If the same procedure is applied for a second axis, an upper arm-embedded anatomical coordinate system (aU) can be defined and rigidly associated to the MIMU strapped to the segment using a time-invariant rotation matrix Upper limb kinematic model. The adopted upper limb kinematic model comprises three rigid body segments: the thorax, humerus (upper arm), and radio-ulna (forearm) connected by spherical joints. Each body segment is associated to an anatomical coordinate system, the axes of which are determined according to the guidelines presented by the International Society of Biomechanics 12 (see Eqs 1, 3 and 7 of the Appendix for details) and are based on the identification of the following palpable ALs: processus xiphoideus, incisura jugulars, right and left acromion, lateral and medial epicondyle, and the radius and ulnar styloid (Fig. 2). The pointing device also provides the inter-anatomical landmark distance along the axis defined by a pair of selected ALs. This information, in addition to the measured anatomical axes orientation and some anthropometric assumptions, is used to estimate the orientation of the upper arm and forearm longitudinal axes together with the distances between the shoulder and elbow joints centres (i.e. upper arm length) and between the elbow joint centre and ulnar styloid (i.e. forearm length; see Eqs. 2 and 4 of the Appendix for details). Figure 2 shows the MIMU setup for collecting the right-side upper limb kinematics together with the measured (pointed) and estimated (internal) anatomical axes used for the anatomical coordinate system determination.
Angular and linear joint kinematics. The orientation of the humerus was computed with respect to the thorax, whereas the orientation of the forearm was computed with respect to the humerus. The 3D shoulder and elbow joint angular kinematics were obtained by decomposing the relevant joint orientation matrices following the recommendations of the International Society of Biomechanics 12 . Furthermore, as the distances between shoulder and elbow joint centres and between the elbow joint centre and ulnar styloid are available in addition to the time-variant orientation of the upper arm and forearm anatomical coordinate systems with respect to the thorax, a homogeneous transformation matrix can be used to solve the forward kinematics of a two-link open kinematic chain ( Fig. 3) at every sampled instant of time and obtain the instantaneous 3D position of the ulnar styloid (i.e. the end-effector of the chain) with respect to the thorax (see Eqs 8-10 of the Appendix) 13 . Table 1 summarises the model's output. experimental approach. The proposed anatomical calibration approach was validated using a six-camera optoelectronic (OPTO) stereophotogrammetric system (Smart-DX, BTS Bioengineering, Italy), considered the gold standard 14 , during two experimental sessions. In the first experimental session, we evaluated the accuracy of the proposed procedure in estimating shoulder and elbow 3D angular kinematics during monoaxial joint rotations. Fourteen young adults (four females and ten males, aged 34 ± 7 years) with no previous shoulder or elbow injury were enrolled in the first experiment to evaluate the performance of the method when wide angular rotations are considered. A second experiment was performed to assess the accuracy of the 3D trajectory of the EF during a reach-to-grasp task that was chosen for its clinical relevance in the assessment of the residual motor function, particularly, in subjects with acquired brain injury 15 . The accuracy of the estimated 3D position of the EF was tested on six young-old subjects (four females and two males, aged 62 ± 13 years) with no previous shoulder or elbow injury and with a range of age similar to that of a typical acquired brain injury population. For both experiments, sample size was determined through an a priori power analysis on pilot data as detailed in the statistical analysis paragraph.
In both experiments, each participant was equipped with three MIMUs (MTw, Xsens Technologies BV) mounted on the thorax, upper arm and forearm, as depicted in Fig. 2. Three retro-reflective markers were attached on each MIMU for validation. Figure 4 shows an example of the measurement setup utilised in this study with both measurement systems used for collecting reference (OPTO) and MIMU-based data. Prior to the execution of the tasks, the full anatomical calibration procedure was performed by a single operator who, first, identified and marked the selected palpable ALs depicted in Fig. 2, and then used the calliper-like device to point the previously marked ALs and compute the relevant quantities in Eqs 1-4 and Eqs 6-7 of the Appendix. Following the latter procedure, additional retro-reflective markers were placed on the previously marked ALs to ensure that each anatomical coordinate system was consistently defined for both the MIMU-based and OPTO-based methods used to collect the 3D markers' position. Markers located on ALs were removed after their position was recorded in a technical coordinate system defined by the three markers placed on the MIMUs during a static trial according to Figure 1. Anatomical calibration of a magnetic and inertial measurement unit (MIMU) fixed on a body segment. The direction of an anatomical axis as defined by two palpable ALs is measured in the MIMU's global coordinate system using an extra MIMU hosted on a calliper-like device to point to the two considered anatomical landmarks. As the two MIMUs share the same global reference frame, the direction of the anatomical axis can be expressed with respect to the sensor-embedded coordinate system of the MIMU fixed on the segment and used to collect motion data.
www.nature.com/scientificreports www.nature.com/scientificreports/ the CAST protocol 16 . The two measurement systems were electronically synchronised using an external trigger signal, and data were collected at a rate of 100 samples/s. MIMU orientation data were retrieved using the manufacturer's proprietary software (MT Manager, v.4.6). A spot-check of MIMU performance was conducted at the experiment's location prior to performing the tests to verify inter-MIMU global coordinate system consistency 17 . The study was approved by the Ethic Committee of the "Agostino Gemelli" University Polyclinic Foundation (Rome, Italy), and informed consent was signed by the participants. All experiments were performed in accordance with relevant guidelines and regulations.
Experiment A: Joint angular kinematics assessment. Starting from the anatomical position, subjects performed the test in the following order: (1) arm elevation in the sagittal plane, (2) arm elevation in the scapular plane, (3) arm elevation in the frontal plane, (4) elbow flexion-extension (6). Finally, with the elbow flexed at 90°, subjects performed: (5) shoulder external-internal rotation, and (6) elbow prono-supination. Three non-consecutive repetitions of each of the abovementioned motor tasks, performed at a self-selected speed, were recorded. Experiment B: Wrist trajectory assessment during reach-to-grasp movements. Fifteen non-consecutive reach-to-grasp movements were performed from a seated position with the dominant arm. Seat and table heights were adjusted such that the seated position was comfortable with the subjects' thighs parallel to the ground, feet flat on the ground and knees flexed at about 90°. Subjects started in a position with the hand on the table in neutral pronation, elbow flexed at approximately 90° and upper arm adducted along the trunk. They were instructed to reach forward at a self-selected speed, grab the object and return to their initial position. The object, centred with respect to the subject's longitudinal axis, was placed on the table at a subject-specific distance corresponding to a quasi-complete extension of the elbow avoiding scapular protraction and elevation.

Data reduction and statistical analysis.
For the first experiment, the β h range of motion (RoM) was computed during arm elevation in the sagittal plane, arm elevation in the scapular plane and arm elevation in the frontal plane, whereas γ h RoM was computed during shoulder external-internal rotation, α e RoM was computed during elbow flexion-extension and, finally, β e RoM was computed during elbow prono-supination. Shoulder www.nature.com/scientificreports www.nature.com/scientificreports/ plane of elevation values (α h ) were also computed during arm elevation in the sagittal, scapular and frontal planes. For the second experiment, the magnitude of the wrist position vector at the grasp point ( → p grasp ) and the linear RoM of the wrist along the x, y and z directions were computed at every reach-to-grasp cycle and considered for statistical analysis.
For both experiments, absolute agreement between the two systems in the considered angular and linear variables was assessed through Bland-Altman analysis corrected for the effects of repeated measurements 18 . The presence of heteroscedasticity was checked by calculating the Kendall rank correlation coefficient τ 19 : when τ > 0.1 a logarithmic transformation of the data was assessed before calculating bias and 95% upper and lower limits of agreement 20,21 . Finally, the estimated (MIMU-based) and reference (OPTO-based) angular and linear trajectories were compared at each movement cycle in terms of root mean square deviation (RMSD) and via the "linear fit method" using typical parameters: R 2 (waveform similarity), a 1 (amplitude) and a 0 (offset) 22 . In brief, R 2 and a 1 are coefficients that range from 0 to 1 (where 1 implies same shape and amplitude), whereas a 0 quantifies the offset between the two curves when a 1 tends to 1. When R 2 > 0.5, the assumption of linearity is considered valid and a 1 can be interpreted as meaningful 23 .
The normal distribution of data was assessed using the Shapiro-Wilk test. The alpha level of significance was set at 0.05. Sample size analysis (β = 0.2), based on Bland-Altman plots 24 , was performed on pilot data that included three subjects performing experiment A and one subject performing experiment B. For experiment A, power analysis was performed using β h RoM relative to arm elevation in the sagittal plane and α e RoM relative to elbow flexion-extension as these were considered the most involved degrees of freedom during the reach-to-grasp task. For experiment B, power analysis was performed using → p grasp . For experiment A, the expected mean difference between methods was 0.616° for β h RoM and 0.548° for α e RoM, whereas the maximum allowed difference was set to 5° resulting in a minimum required number of pairs equal to 7 and 27 for β h RoM and α e RoM, respectively; however, the number of pairs eventually chosen for experiment A equalled 42 (14 subjects × 3 repetitions). For experiment B, the expected mean difference between methods was 7 mm and the maximum allowed difference was set to 20 mm resulting in a minimum required number of pairs equal to 62; however, the number of pairs eventually chosen for experiment B equalled 90 (6 subjects × 15 repetitions).

Results
RMSD, R 2 , a 1 and a 0 values relative to both the angular and linear trajectories comparison were not normally distributed. Hence, these variables are presented in terms of median and inter-quartile range (IQR). Regarding the first experiment, Table 2 reports the reference and estimated angular RoM in addition to the results of the comparison between the reference and estimated angular trajectories. Reference frontal, scapular and sagittal shoulder plane of elevation values (α h ) were 8.4° ± 10.8°, 31.8° ± 10.4° and 84° ± 11.8°, respectively, against the corresponding 7.9° ± 12.1°, 33.2° ± 13.4° and 87.2° ± 12.7° estimated by our MIMU-based model. Figure 5 shows the Bland-Altman plots depicting the agreement between the reference and estimated joint angular RoM (bias and 95% lower and upper limits of agreement are reported in the graphs in Fig. 5).
Regarding the reach-to-grasp experiment, the target was reached at an average velocity of 0.181 ± 0.06 m/s. Reference and estimated → p grasp values were, on average, 516 ± 13 and 522 ± 14 mm, respectively, and their agreement is reported in Fig. 6a. Heteroscedasticity was observed via visual inspection of the Bland-Altman plot   Table 2. Mean ± sd of the reference (OPTO) and estimated (MIMU) angular RoM variables considered in the analysis. Median ± IQR are reported for RMSD (also expressed as percentage of the reference RoM, RoM OPTO ) and for the three coefficients obtained from the Linear Fit Method R 2 (coefficient of determination, i.e. waveform similarity), a 1 (slope, i.e. amplitude) and a 0 (intercept, i.e. offset) used for comparing pairs of reference and estimated angular trajectories. (2019) 9:14449 | https://doi.org/10.1038/s41598-019-50759-z www.nature.com/scientificreports www.nature.com/scientificreports/ depicted in Fig. 6a and subsequently verified using Kendall's τ = 0.309. Bias and limits of agreement depicted in Fig. 6a were adjusted accordingly. Finally, Table 3 reports the reference and estimated wrist linear RoM together with the results of the comparison between the reference and estimated linear trajectories; Fig. 6b-d show Bland-Altman plots depicting the agreement between the reference and estimated wrist linear RoM (bias and 95% lower and upper limits of agreement are reported in the graphs in Fig. 6).

Discussion
This study proposes an innovative upper limb movement analysis protocol based on the use of wearable MIMUs, which provide anatomically and clinically consistent joint kinematics estimates. The validity of the proposed approach was verified against a camera-based optoelectronic system during two experiments aimed at assessing the accuracy of a two-link open kinematic chain in estimating the 3D shoulder and elbow angular kinematics together with the 3D trajectory of the wrist during standard mobility tests and reach-to-grasp movements.
In the first experiment, the angular RoM were slightly underestimated in a systematic manner, (bias <−2.6°) and they showed a good level of agreement. Largest errors were observed for arm elevation in the frontal plane which was characterised by the highest bias and limits of agreement range (Fig. 5). During the elevation of the arm in the frontal, scapular and sagittal planes, α h was correctly estimated with <3° of error. When comparing the angular trajectories, the RMSD was <4% of the reference RoM for all considered angular RoMs (Table 2). A high waveform similarity (0.995 < R 2 < 0.999), high amplitude similarity (0.951 < a 1 < 1.02) and small offset (a 0 < 3°) were found for all considered angular variables except for the arm elevation in the frontal plane, which was characterised by the highest offset (Table 2). Unfortunately, none of the previous MIMU-based studies reported the level of accuracy associated with the estimate of upper limb joint kinematics. Assessing the accuracy of the proposed method against a gold standard joint kinematics is, indeed, a highlight of this study as the same anatomical coordinate system determination and same ALs were used both for MIMUs and the reference camera-based motion analysis system.
Regarding the second experiment, the agreement analysis revealed a systematic overestimation of both wrist position and linear RoM (Fig. 6). A low bias (<10 mm) and good level of agreement were found between the reference and estimated → p grasp (Fig. 6a). With regard to the EF trajectory, the best performance was observed in the antero-posterior and medio-lateral directions with a low bias (7 to 11 mm) and small limits of agreement ranges (Fig. 5b,c) in addition to a low RMSD (8 mm, <5.6% of the reference RoM) and offset values, high amplitude and waveform similarity coefficients (Table 3). Conversely, the cranio-caudal direction was characterised by the largest limits of agreement range (Fig. 6c), the largest RMSD and offset values and the lowest waveform and amplitude similarity coefficients (Table 3). An overall picture of the accuracy of the proposed method in tracking the 3D position of the wrist is shown in Fig. 7.
In this study, the magnitude of the errors of the wrist trajectory was comparable with those reported in previous studies based either on double numerical integration with optimisation during reach-to-grasp movements 25 or on the use of forward kinematics during standard shoulder/elbow mobility tests 26 .
The worst performance in estimating the wrist linear position was observed along the cranio-caudal direction, which was also characterised by the smallest RoM. The source of this error might be found in the task itself, which requires the elevation of the shoulder girdle to allow hand-table clearance. To reduce the number of MIMUs involved, we decided to use a simplified open chain kinematic model of the upper limb that assumes the humerus to be hinged to the thorax and no translational degrees of freedom are considered. Thus, as the scapula is not included in the open chain kinematic model, any translation of the glenohumeral joint (e.g. elevation of the shoulder girdle) would be neglected. The main peculiarity of the proposed method is that it can be applied when dealing with patients presenting severe impaired mobility and/or those equipped with obtrusive medical apparatuses. In such cases, the anatomical landmark identification approach is the most feasible method among those proposed in the literature.
Finally, an important point is that the proposed methodology relies on the assumption that all MIMUs utilised in the protocol share the same global coordinate system, and, therefore, the overall method performance is subject to the quality of the orientation estimates obtained using the sensor fusion algorithm. As this assumption declines in the presence of non-homogenous ferromagnetic disturbances and physical calibration issues of sensors, a system spot-check is necessary prior to data collection 17 .
In conclusion, the proposed anatomical calibration estimated reliable joint angular kinematics in full compliance with the metrological standards of clinical movement analysis. Compared to functional approaches, the   Table 3. Mean ± sd of reference (OPTO) and estimated (MIMU) linear RoM of the ulnar styloid ([p x , p y , p z ]). Median ± IQR are reported for RMSD (also expressed as percentage of the reference RoM, RoM OPTO ) and for the three coefficients obtained from the linear fit method R 2 (coefficient of determination, i.e. waveform similarity), a 1 (slope, i.e. amplitude) and a 0 (intercept, i.e. offset) used for comparing pairs of reference and estimated linear trajectories.
Scientific RepoRtS | (2019) 9:14449 | https://doi.org/10.1038/s41598-019-50759-z www.nature.com/scientificreports www.nature.com/scientificreports/ proposed anatomical calibration can also be used on bedridden patients who are incapable of performing segment rotations or assuming fixed postures. Additionally, the approach allows one to extend the analysis to linear joint kinematics, providing an adequate spatial resolution for clinical assessment purposes. Consequently, the current study lays the methodological foundation to delve further into the reach-to-grasp kinematics in measurement settings wherein the use of traditional motion capture technologies is not viable. From this perspective, the proposed approach could enable the measurement of the rehabilitative progress of stroke patients from the very acute phase (when patients are bedridden in stroke units) to the chronic stage of the disease.

Data Availability
Upon reasonable request, the datasets used and analysed during the current study will be made available by the corresponding author.