 Research
 Open Access
 Published:
In vivo estimation of the shoulder joint center of rotation using magnetoinertial sensors: MRIbased accuracy and repeatability assessment
BioMedical Engineering OnLinevolume 16, Article number: 34 (2017)
Abstract
Background
The human glenohumeral joint is normally represented as a spherical hinge and its center of rotation is used to construct humerus anatomical axes and as reduction point for the computation of the internal joint moments. The position of the glenohumeral joint center (GHJC) can be estimated by recording ad hoc shoulder joint movement following a functional approach. In the last years, extensive research has been conducted to improve GHJC estimate as obtained from positioning systems such as stereophotogrammetry or electromagnetic tracking. Conversely, despite the growing interest for wearable technologies in the field of human movement analysis, no studies investigated the problem of GHJC estimation using miniaturized magnetoinertial measurement units (MIMUs). The aim of this study was to evaluate both accuracy and precision of the GHJC estimation as obtained using a MIMUbased methodology and a functional approach.
Methods
Five different functional methods were implemented and comparatively assessed under different experimental conditions (two types of shoulder motions: cross and star type motion; two joint velocities: ω_{max} = 90°/s, 180°/s; two ranges of motion: Ɵ = 45°, 90°). Validation was conducted on five healthy subjects and true GHJC locations were obtained using magnetic resonance imaging.
Results
The best performing methods (NAP and SAC) showed an accuracy in the estimate of the GHJC between 20.6 and 21.9 mm and repeatability values between 9.4 and 10.4 mm. Methods performance did not show significant differences for the type of arm motion analyzed or a reduction of the arm angular velocity (180°/s and 90°/s). In addition, a reduction of the joint range of motion (90° and 45°) did not seem to influence significantly the GHJC position estimate except in a few subjectmethod combinations.
Conclusions
MIMUbased functional methods can be used to estimate the GHJC position in vivo with errors of the same order of magnitude than those obtained using traditionally stereophotogrammetric techniques. The methodology proposed seemed to be robust under different experimental conditions. The present paper was awarded as “SIAMOC Best Methodological Paper 2016”.
Background
An accurate estimation of the glenohumeral joint center (GHJC) is of primary importance in biomechanics for the upper limb motion analysis [1, 2], the design of robotic arm exoskeletons for shoulder rehabilitation [3], surgical navigation procedures and bones alignment during surgery [4]. The GHJC is commonly employed for defining the humerus anatomical coordinate system [5] and to compute internal joint moments [6]. Being an internal anatomical landmark (AL), the GHJC cannot be identified by palpation as usually done for superficial ALs [7, 8]. The GHJC position is commonly estimated in vivo by using either regressive [7, 9, 10] or functional methods [11–13]. Regressive methods estimate GHJC position from empirical geometrical relations between specific ALs. Functional methods require the subject to perform ad hoc joint movements while recording the motion of the adjacent segments [8]. The assumption underlying these methods is that the glenohumeral joint (GHJ) is well described by a ballandsocket joint [14, 15] and its center of rotation (CoR) can be made to coincide with the geometrical center of the medialsuperior portion of the humeral head [16].
Several studies have been proposed to estimate the GHJC position using stereophotogrammetry, electromagnetic or ultrasound positioning systems [7, 9, 11]. These measurement technologies are used to track the threedimensional positions of selected active or passive markers/sensors attached on the subject’s body. In particular, stereophotogrammetry markerbased methods have been extensively investigated in terms of formal description [17, 18], using mathematical simulations [19], through experiments using mechanical devices [19, 20], ex vivo [21] and in vivo [22] experiments. While few studies have estimated the precision associated to the GHJC identification in vivo [7, 9, 11], only three works evaluated the accuracy in vivo using medical imaging techniques as gold standard [10, 12, 13]. Campell et al. [10] compared the accuracy and reliability of eight regressive methods, Lampereu et al. [12] compared five different functional methods and Nikooyan et al. [13] tested two functional methods in patients with shoulder hemiarthroplasty.
However, methods based on the use of multiplecamera setup suffer from usability problems in clinical settings and in telerehabilitation scenarios. In fact, they require the adoption of expensive equipment, an adequate space for acquisitions, rarely available in the clinical facilities or homebased environment. Furthermore, the expertise of the operator is fundamental for the marker placement, data acquisition and processing. In addition, optical occlusion problems may occur, especially during passive exercises.
The above listed problems could be potentially solved by using miniaturized magnetoinertial measurement units (MIMUs). Sensors wearability, along with the absence of complex external setup and occlusion problems, make this technology a promising alternative to opticalbased approaches. Differently from the stereophotogrammetry or alternative positioning systems, MIMUs do not provide reliable position measurements, but only accelerations, angular velocities and magnetic field intensity [23]. The measured quantities can be then combined by using sensorfusion algorithms to obtain MIMU orientation [23].
Previous studies demonstrated that it is possible to estimate the CoR position of a spherical joint using inertial sensors [24–26]. However, the methods validation was carried out on a mechanical analogue of a generic spherical human joint by manually moving a segment with respect to a stationary frame. The abovementioned methods assumed that a unique stationary CoR exists during the functional calibration movement. Conversely, Seel et al. [27] presented a method for the functional estimation of joint axes direction and CoR position, in hingelike and ballandsocketlike joints, respectively, which does not require the existence of a stationary CoR. They applied the method for the estimation of the hip joint CoR from data recorded during gait. However, no quantitative information about the errors in the CoR identification was reported. To the best of our knowledge, no studies have been presented for the in vivo functional identification of the GHJC position using MIMU technology.
The aim of this study is to evaluate both accuracy and precision of the GHJC estimation using MIMUs and following a functional approach. In vivo magnetic resonance imaging (MRI) was used as a gold standard for the validation of the methodology on five healthy subjects. Different experimental conditions were considered (joint motion type, movement speed and range of motion, ROM) and five different functional GHJC estimation methods were comparatively evaluated. Selected methods differed for the number of MIMU employed and the algorithms implementation.
Methods
The following algorithms were implemented and compared.

1a.
The original algorithm proposed by Crabolu et al. [25, 26], based on the use of a single MIMU, hereafter named Null Acceleration Point (\({\text{NAP}}^{\left( 1 \right)}\)), where the superscript indicates the number of MIMUs used;

1b.
A variant of the algorithm \({\text{NAP}}^{\left( 1 \right)}\), which includes a criterion based on the analysis of the angular velocity for selecting the samples to be used for the estimation of the GHJC, hereafter named \({\text{NAP}}_{\upomega}^{\left( 1 \right)}\);

1c.
A variant of algorithm \({\text{NAP}}_{\upomega}^{\left( 1 \right)}\), which uses two MIMUs and exploits the scapulae acceleration information for selecting the samples to be used for the estimation of the GHJC, hereafter named \({\text{NAP}}_{{\upomega {\text{a}}}}^{\left( 2 \right)}\);

2.
The original algorithm proposed by Seel et al. [27], based on the use of two MIMUs, referred to as Symmetrical Specific Force Center (\({\text{SSFC}}^{\left( 2 \right)}\));

3.
A combination of the two methods 1a and 2, based on the use of two MIMUs, that we named Symmetrical Acceleration Center (\({\text{SAC}}^{\left( 2 \right)}\)).
All the algorithms exploit rigid body kinematics equations. In this regards, the acceleration of an arbitrary point P of a rigid body B during free motion is given by
where a _{ t } is the acceleration of a point O on the same body, r is the position vector from O to P, ω is the angular velocity, and \(\frac{{d{\varvec{\upomega}}}}{dt}\) the angular acceleration. When a rigid body is constrained through a ballandsocket joint, it can only experience pure rotational motion around the CoR. If O coincides with the CoR, at becomes null, thus Eq. (1) can be computed as:
After some algebraic manipulations, Eq. (2) can be rearranged as:
where K is equals to
Equation (3) is linear with respect to the unknown vector r, which represents the CoR position.
Description of the algorithms
Algorithm \(NAP^{\left( 1 \right)}\)
This algorithm requires a single sensor and assumes that the CoR point has a null acceleration. The vector r can be estimated according to Eq. (3) by measuring the observable variables a and ω, which can be obtained from the data recorded by a MIMU attached to the rotating body. At each sampling time, specific force f, angular velocity, and magnetic field vectors are measured with respect to the MIMU coordinate system (MCS). The MCS orientation with respect to a global coordinate system (GCS) can be estimated using a Kalman filter from the measures recorded by the triaxial accelerometer, gyroscope, and magnetometer [23]. In particular, MIMU orientation is required to subtract the gravitational acceleration g from the specific force f to obtain the coordinate acceleration a in Eq. (3). Applying Eq. (3) at each of the N sampling time points recorded during a pure rotational motion of the rigid body, a leastsquares solution for r can be computed.
Algorithm \(NAP_{\omega }^{\left( 1 \right)}\)
It was proved that to reduce SNR associated to angular velocity signals [26], slow calibration movements should be avoided since they lead to worse results compared to fast movements [25, 26]. To take advantage of this previous observation, only data samples characterized by a magnitude of the angular velocity higher than an empirically chosen threshold, equal to 0.5 rad/s, were included in the leastsquare solution computation.
Algorithm \(NAP_{{\omega {\text{a}}}}^{\left( 2 \right)}\)
To compensate for moderate violations of the assumption of a null acceleration of the GHJC during the shoulder movement, only data samples for which the magnitude of the acceleration vector of the MIMU placed on the “quasi stationary” body segment was lower than an adaptive threshold were included in the leastsquare solution computation. The threshold was set equal to the mean magnitude of the acceleration vector of the “quasi stationary” segment. To preserve a minimum number of samples for the computation, when the number of selected samples was less than 300, the threshold was iteratively incremented of 0.1 m/s^{2} until this criterion was met. These values were empirically chosen after some preliminary tests in order to obtain the best results while limiting the number of samples. The applied procedure can be described as follows:

1.
Set the threshold “γ” to the mean magnitude of the acceleration vectors measured by the MIMU placed on the “quasi stationary” body segment, a _{ 2i }:

2.
Retain the M samples i for which \(\left\ {{\mathbf{a}}_{{2\varvec{i}}} } \right\ < \gamma\);

3.
If M < 300, set γ = γ + 0.1 and repeat point 2, else compute the CoR using the selected M samples.
Algorithm \(SSFC^{(2)}\) [27]
Let us consider two rigid segments connected by a spherical joint and two MIMUs firmly attached to the segments. Then, the specific force of the CoR as recorded by the two MIMUs should have the same magnitude. It implies that the magnitude of the accelerations experienced by the two MIMUs must be equal each other:
where the subscripts denote the quantities measured respectively by the two MIMUs.^{Footnote 1} Calculating the gradients of Eq. (5) leads to:
Then, it is possible to find those r _{1} and r _{2} that minimize the quantity in Eq. (5) with the Gauss–Newton algorithm as follows:

1.
Define and initialize x = [r _{1}, r _{2}]^{T};

2.
Calculate the error vector e defined by:

3.
Use Eq. (6) to calculate the Jacobian (J) de/dx and its Moore–Penrosepseudoinverse pinv(J);

4.
Update x by: \({\mathbf{x}} = {\mathbf{x}}  pinv({\mathbf{J}})\) and repeat from step 1.
The vector x was initialized to [0 0 0 0 0 0]^{T}. The algorithm converged for all the trials analysed in 30 iterations. Convergence was assumed as the achievement of a stable condition in which the solution improvement is < 0.01 m/s^{2}.
Algorithm \(SAC^{\left( 2 \right)}\)
If both segments move during the functional movement, the CoR could experience a coordinate acceleration different from zero. However, by rigidly attaching a MIMU on each segment, it is possible to estimate the position of the relative CoR by considering its acceleration a _{ t }. Taking in account a _{ t }, Eq. (3) becomes:
By computing Eq. (8) for both segments, and substituting a _{ t }, it yields:
where \(\mathop {\mathbf{R}}\nolimits_{2}^{1}\) is the rotation matrix from MCS_{2} to MCS_{1}.
Equation (9) represents an indeterminate linear system with six unknowns and three equations, but applying Eq. (9) at each of the N sampling instants recorded during the joint movement, an oversized linear system is obtained:
and a leastsquares solution for r _{ 1 } and r _{ 2 } can be computed.
Experimental design
Population
Five subjects (3 M and 2 F) without upper limb disorders were enrolled in the study. The study has been performed following the principles outlined in the Helsinki Declaration of 1975, as revised in 2000, on healthy subjects. The volunteers’ age, height, and body mass index (BMI) were respectively: 36 ± 4 years, 1.7 ± 0.1 m, and 20.7 ± 1.7 kg/m^{2}. All participants signed an informed consent form prior to start the recordings, after being informed about the aims and procedures of the experiments. Each subject first underwent an MRI examination and, immediately after, followed a protocol for the evaluation of the GHJC based on a functional approach.
Measurement of the GHJC reference position using MRI
The MIMUbased GHJC position estimate is expressed with respect to the MCS [26], which is aligned to the edges of the MIMU housing. In order to obtain the gold standard position of the GHJC with respect to the MCS using MRI, a nonferromagnetic phantom of the MIMU was designed (Fig. 1), and realized by 3D printing (MakerBot Replicator X2) in acrilonitrile–butadiene–stirene (ABS), with a tolerance of 0.1 mm. The phantom included three nonaligned spherical holes filled with copper sulfate in order to be MRI visible, in positions enabling the reconstruction of a local frame coincident to the MCS of the actual MIMU (Fig. 2). One of the three holes was in the same 3D position of the origin of the triaxial accelerometer, with respect to the housing edges. The MIMU, or its phantom, could be fixed to the arm by means of a 3Dprinted custom clip and an elastic band.
The MRI data of subjects’ right scapula and humerus were acquired while the phantom was attached to the arm (Fig. 3a). MR scans of the whole right humerus and scapula were obtained by using a 1.5 T MR scanner (Philips Intera Achieva version 1.7). Spin Echo imaging sequences were used (axial T1W: TR 660 ms; TE 18 ms; flip angle 90 deg; Contiguous Slice Thickness 4 mm, FoV 280 mm). Bone contours were identified using a semiautomatic segmentation procedure. 3D reconstructions of the entire scapular and humeral bones were obtained using the AMIRA image processing software (Visualization Sciences Group, v.5.4). The gold standard GHJC position was estimated as the center of the best fitting sphere to the humeral head of the reconstructed humeral bone [16] (Fig. 2).
Immediately following the MRI acquisition, the phantom was replaced by the actual MIMU, paying attention to preserve the position of the clip (Fig. 3b). In this way, the MCS of the MIMU coincided with the coordinate system of the MRvisible phantom.
Experimental protocol
The MIMUs was attached at approximately the same anatomical location for all subjects. A first MIMU was attached laterally to the third distal portion of the arm through the clip and a Velcro strap as shown in Fig. 3. With the arm in anatomical position, the MIMU was mounted with x axis approximately directed superiorly along the long axis of the humerus, the z axis pointing laterally and the y axis posteriorly (Fig. 3). A second MIMU was attached on the scapula with the lower edge along the cranial edge of the spina scapulae using a doublesided tape [28]. A third MIMU was fixed on the thorax at the sternum level. This MIMU was not used to compute the GHJC but only to monitor the amplitude of the thorax motions during the functional movement. Each MIMU comprises triaxial accelerometer, gyroscope and magnetometer (MTw2 Awinda wireless motion tracker system, Xsens, sample frequency: 100 Hz; dynamic accuracy: Roll/pitch = 0.75° RMS; Heading = 1.5° RMS). The sensor position with respect to the housing was provided by the manufacturer. A 15min warmup period was included before starting the data collection. A preliminary spotcheck of the MIMUs orientation estimates was performed according to the guidelines proposed in [29]. Briefly, the MIMUs were aligned to each other and attached to a rigid plastic plate; then, the differences in the orientation estimates provided by each MIMU were computed while the orientation of the plastic plate was being manually varied about the three directions. This dynamic test had a duration of 30 s. In addition, the gyroscope bias and the residual acceleration after compensation of the gravity during a 30 s static recording were computed.
Subjects were instructed to perform arm movements trying to minimize the movement of the trunk and were allowed to practise few times before starting the recordings (Fig. 3c). Data were recorded under six different experimental conditions: two types of shoulder motions (cross and star type motion as in Table 1), two joint velocities (fast movements: ω_{max} = 180°/s, slow movements: ω_{max} = 90°/s) and two ranges of motion (Ɵ_{1} = 45°; Ɵ_{2} = 90°). Each trial comprised a static acquisition (5 s) followed by the selected shoulder motion. Three repetitions for each condition were recorded, for a total of 24 trials per subject.
Acquisitions were performed on each subject by the same examiner.
Data processing and analysis
For each MIMU, the bias affecting the gyroscopes during the joint motion was partially compensated by subtracting the mean angular velocity value obtained during the spotcheck. Prior to the differentiation, the angular velocity was filtered using a decimated wavelet denoising approach exploiting the Bior3.3 mother wavelet and softthresholding. Fixed leveldependent thresholds over the chosen four decomposition levels were used, whereas level 1 was completely cleared out from the signal, achieving a lowpass behavior and limiting the signal bandwidth up to a quarter of the sampling frequency. The angular acceleration was computed by differentiating the angular velocity signal with a threepoint central difference operator.
For each trial and algorithm, the vector difference e _{ i } between the gold standard GHJC location and the estimated one, and the scalar difference e _{ ri } between the true radius length (r) and the estimated one, were computed. For each subject, the overall algorithm accuracy was computed in terms of E and E _{ r } values obtained by averaging the module of e _{ i } and e _{ ri } over the 24 trials (2 type of motions × 2 two joint velocities × two ranges of motion × 3 trials repetitions), as follows:
where \(\left\ \cdot \right\\) is the Euclidean distance.
The repeatability associated to the GHJC estimation was computed for each algorithm in terms of E _{ SD } value:
where SD _{ j }, with j = x, y or z, is the standard deviation of each estimated GHJC position coordinate.
Preliminary normality test (Shapiro–Wilk test) was performed to choose the most appropriate statistical analysis. For each algorithm, we investigated whether the method’s accuracy was affected by the following independent factors: (a) joint angular velocity, (b) type of the joint motion, (c) joint angular ROM. To this purpose, a twotailed Student’s t test was performed for each one of the independent factors analysis above (individual subject comparisons).
In order to assess if there was a significant difference between the five algorithms (\({\text{NAP}}^{\left( 1 \right)}\), \({\text{NAP}}_{\upomega }^{\left( 1 \right)}\), \({\text{NAP}}_{{\upomega {\text{a}}}}^{\left( 2 \right)}\), \({\text{SSFC}}^{\left( 2 \right)}\), \({\text{SAC}}^{\left( 2 \right)}\)), oneway ANOVA test for normal samples distribution was used (individual subject comparisons). When a significant difference was detected (p < 0.05), a Student’s t test was performed between every pair of methods.
We also evaluated whether differences in the GHJC location error component along each anatomical direction were present, by performing a Student’s t test for each pair of coordinates (individual subject comparisons).
A Bonferroni’s correction was used when multiple comparisons were studied.
Results
MIMU spot check
Results relative to the spotcheck performed on the MIMUs attached over the humerus and scapula are shown in Table 2.
Effects of joint angular velocity, type of joint motion, joint angular ROM
For each algorithm and each subject, no statistically significant differences were observed between slow and fast joint angular velocities and between cross and star joint motion (p > 0.05, Student’s t test). Even for the different joint ROM, no statistically significant differences were observed except in the following cases: error E _{ r }, subject 2, methods NAP and \({\text{SSFC}}^{\left( 2 \right)}\); error E, subjects 4 and 5, method \({\text{SAC}}^{\left( 2 \right)}\) (p < 0.05, Student’s t test). For the sake of brevity, the results relative to each independent factor are shown only for the algorithm \({\text{NAP}}_{\upomega }^{\left( 1 \right)}\) (Table 3).
Accuracy and repeatability of the tested algorithms
The values of E (mean and standard deviation, STD), obtained for the different five algorithms for each subject, are reported in Table 4 and in the bar chart of Fig. 4. The smallest mean E value was obtained with the \({\text{NAP}}_{\upomega }^{\left( 1 \right)}\) algorithm (20.6 ± 10 mm). With this algorithm, E over each subject ranges between a minimum of 11.2 mm (subject 3) and 37.5 mm (subject 2). The lowest accuracy was found with the \({\text{SSFC}}^{\left( 2 \right)}\) algorithm (29.9 ± 10 mm, ranging across subjects from 22.2 to 37.9 mm).
The values of E _{ r } (mean and STD) obtained for each algorithm for each subject are reported in Table 5 and in the bar chart of Fig. 5. Also for the E _{ r }, the smallest value was obtained for \({\text{NAP}}_{\upomega}^{\left( 1 \right)}\) algorithm (8.7 ± 6 mm). The minimum value was detected for the subject 4 (4.1 ± 2 mm), whereas the highest value was again obtained on subject 2 (12.6 ± 9 mm). The largest mean E _{ r } was obtained with \({\text{SSFC}}^{\left( 2 \right)}\) algorithm (17.1 ± 11 mm, ranging from 13.9 to 23.2 mm).
The repeatability results (E_{SD}) are displayed for each subject and algorithm in Table 6. The \({\text{NAP}}_{\upomega }^{\left( 1 \right)}\) algorithm showed the smallest E_{SD} mean value compare to the other algorithms (9.4 mm, ranging from 5.3 to 16.3 mm). The \({\text{SSFC}}^{\left( 2 \right)}\) algorithm returned the largest E_{SD} values (13.8 mm, ranging from 9.2 to 17.3 mm).
Algorithms comparison
The ANOVA test showed that the resulting E and E _{ r } were significantly different among algorithms. In Tables 4 and 5, the algorithm results are marked with one star when a statistically significant difference was detected between that algorithm (for a given subject) and the others (p < 0.005, Student’s t test with Bonferroni’s correction). A significant difference was observed only between the \({\text{SSFC}}^{\left( 2 \right)}\) and the other algorithms for both E and E _{ r } with the exception of E values in the subject 2. No significant differences were detected among the other four algorithms (\({\text{NAP}}^{\left( 1 \right)}\), \({\text{NAP}}_{\upomega}^{\left( 1 \right)}\), \({\text{NAP}}_{{\upomega{\text{a}}}}^{\left( 2 \right)}\), \({\text{SAC}}^{\left( 2 \right)}\)).
GHJC error directionality
For all subjects and all methods, significant differences were observed between the errors along the x and z directions (marked with a circle in Table 7) as well as x and y except for the subject 5; between y and z directions a significant difference was detected only for the subjects 1 and 5 (p < 0.016, Student’s t test and Bonferroni’s correction). For the sake of brevity, the results along each direction in estimating GHJC are reported only for the algorithm \({\text{NAP}}_{\upomega}^{\left( 1 \right)}\) (Table 7).
Discussion
In the present study, we performed an evaluative comparison, in terms of accuracy and repeatability, of five different algorithms for the in vivo identification of the GHJC using magnetoinertial sensing technology. The methods accuracy was evaluated on five healthy subjects, using as gold standard the GHJC positions determined from bioimaging (MRI).
The analysis performed on the five participants seemed to indicate that the accuracy of the tested methods was not significantly affected by the type of arm motion analyzed (star vs. cross) or the reduction of the arm angular velocity from 180°/s to 90°/s. Similarly, a reduction of the joint ROM from 90° to 45° did not seem to influence significantly the GHJC position estimate except in a few subjectmethod combinations.
These results seemed to confirm the robustness of the methodology proposed under different experimental conditions. However, when defining the experimental protocol for the functional identification of the GHJC position, good practice guidelines should be followed (slow joint movement should be avoided, and these movements should involve at least two different axes of rotation and allowing acquisition of sufficient number of samples) [26].
The null acceleration point methods (\({\text{NAP}}^{\left( 1 \right)}\), \({\text{NAP}}_{\upomega}^{\left( 1 \right)}\) and \({\text{NAP}}_{{\upomega{\text{a}}}}^{\left( 2 \right)}\)) and the Symmetrical Acceleration Center \({\text{SAC}}^{\left( 2 \right)}\) method outperformed the Symmetrical Specific Force Center (\({\text{SSFC}}^{\left( 2 \right)}\)) method. By analysing each of the five participants, we could not find any significant differences among the different NAP method implementations. This would suggest that the sample selection procedures implemented in both \({\text{NAP}}_{\upomega}^{\left( 1 \right)}\) and \({\text{NAP}}_{{\upomega{\text{ a}}}}^{\left( 2 \right)}\) algorithms were not effective under the specific experimental conditions analyzed. However, it cannot be excluded that differences might arise for slower movements and/or larger movement of the trunk segment. Considering the highest repeatability and accuracy of the algorithm \({\text{NAP}}_{\upomega }^{\left( 1 \right)}\) and its advantage in using a single sensor, we suggest the use of this algorithm when estimating the GHJC position in similar experimental conditions.
GHJC position errors computed over all subjects slightly varied among NAP and \(SAC^{\left( 2 \right)}\) methods (20.6 to 21.9 mm), whereas repeatability varied between 9.4 and 10.4 mm.
Interestingly, the use of a second MIMU attached to the scapula (\({\text{NAP}}_{{\upomega{\text{a}}}}^{\left( 2 \right)}\) and \(SAC^{\left( 2 \right)}\)) seemed to offer no advantages with respect to single unit methods (\({\text{NAP}}^{\left( 1 \right)}\), \({\text{NAP}}_{\upomega }^{\left( 1 \right)}\)). This results can be explained by the intrinsic difficulties in tracking the scapular motion due to the soft tissue artefacts [30] and the fact that the trunk movements observed for all subjects during the functional exercises were quite limited (ROM < 30, trunk angular velocity RMS < 30°/s).
The \({\text{SSFC}}^{\left( 2 \right)}\) algorithm was the less accurate and repeatable method despite the fact it does not require to remove the gravity contribution and was conceived to take into account also for a translation of the CoR. An explanation of the poorer performance compared to the other methods can be provided by the small amplitude of the signals recorded by the second MIMU positioned on the scapula (ω_{max} < 30°/s) and the difficulty related to the scapula tracking.
For all methods, errors along the x axis (approximately longitudinal anatomical direction of the humerus) resulted slightly lower than those in the other directions. This might be either due to anisotropy of the soft tissue artefact components or to the shoulder motions performed.
To the best of authors’ knowledge, this is the first study investigating functional methods for the GHJC estimation in vivo using MIMUs, exploiting a medical imaging approach as gold standard. The errors associated to the best performing methods were of the same order of magnitude than those obtained using a functional approach from markerbased stereophotogrammetric systems [10, 12, 13] and smaller than those found by using the regressive method recommended by the International Society of Biomechanics (ISB) [7] (about 30 mm, on average) [10]. Lampereu et al. [12] obtained on four subjects errors ranging between 11 and 17 mm with five different functional methods. Nikooyan et al. tested the symmetrical CoR estimation (Score) [18] and instantaneous helical axis [15] methods in five patients with shoulder hemiarthroplasty, finding accuracy errors respectively of 20.7 and 14.7 mm [13]. GHJC errors found in the present in vivo study were about five times higher than those found on a mechanical analogue [25, 26]. These differences were expected and they could be mainly ascribed to: (1) the presence of soft tissue artefacts, which can introduce a relative movement between the body and device which is influenced by the mass of the MIMU and by the fixing technique; (2) the fact that the GHJC is not perfectly stationary during the calibration movements potentially affecting methods \({\text{NAP}}^{\left( 1 \right)}\), \({\text{NAP}}_{\upomega }^{\left( 1 \right)}\), \({\text{NAP}}_{{\upomega{\text{a}}}}^{\left( 2 \right)}\). In fact, due to the scapulahumeral rhythm, the scapula remains stationary only for the first 30° of humerus abduction. Furthermore, it is difficult for the subject to completely avoid trunk movements during the calibration exercise.
Interestingly, we found that both in the experiments using a mechanical linkage and a mathematical simulation [24] and the present in vivo study, the type and range of arm motion were not critical. Conversely, while in the study on the mechanical linkage we observed an increase of the errors when the angular velocity was decreased from 180°/s to 60°/s, in the present study we could not find any significant trend. This discrepancy could be partially explained by the larger angular velocity tested in the in vivo experiments (slow movements characterized by ω_{max} = 90°/s). However, due to the limited number of subjects analysed (five), caution is required when generalizing the study findings. Future studies including a larger number of participants, with different anthropometric characteristics, from both healthy and pathological populations, are required to further validate the proposed methodologies.
Conclusions
Based on the study evidences, it is possible to estimate in vivo the GHJC by using magnetoinertial sensors and a functional approach, obtaining accuracy and repeatability of the same order of magnitude of those achievable with more expensive and complex stereophotogrammetry techniques. Under the different experimental conditions analyzed and in presence of limited trunk movements, GHJC errors of about 20 mm were achieved using a single MIMU. The positioning of an additional MIMU did not improve the methods performance. However, further analyses are needed to investigate the advantage associated to the use of two MIMUs when large trunk movements are observed.
Notes
 1.
The accelerations in Eq. (5) (a _{ 1 }, a _{ 2 }) represent the specific forces as they are recorded by the MIMUs. For the sake of clarity, we preserved the same notation.
Abbreviations
 GHJC:

Glenohumeral joint center
 GHJ:

Glenohumeral joint
 MIMUs:

magnetoinertial measurement units
 AL:

anatomical landmark
 CoR:

center of rotation
 ROM:

range of motion
 MRI:

magnetic resonance imaging
 NAP:

null acceleration point
 SSFC:

Symmetrical Specific Force Center
 SAC:

Symmetrical Acceleration Center
 MCS:

MIMU coordinate system
 GCS:

global coordinate system
 BMI:

body mass index
 ABS:

acrilonitrile–butadiene–stirene
 e :

vector difference e between the true GHJC location and the estimated one
 e _{ r } :

scalar difference between the true radius length (r) and the estimated one
 E :

mean Euclidean distance
 E _{ r } :

mean e _{ r }
 E _{ SD } :

repeatability
 STD:

standard deviation of the error
 ISB:

International Society of Biomechanics
 Score:

symmetrical CoR estimation
References
 1.
Cutti AG, Giovanardi A, Rocchi L, Davalli A, Sacchetti R. Ambulatory measurement of shoulder and elbow kinematics through inertial and magnetic sensors. Med Biol Eng Comput. 2008;46:169–78.
 2.
Tranquilli C, Bernetti A, Picerno P. Ambulatory joint mobility and muscle strength assessment during rehabilitation using a single wearable inertial sensor. Med Sport. 2013;66:583–97.
 3.
Kumar A, Jobe C, Saha S. Location of the instantaneous center of rotation for shoulder motion. In: Proceedings of the 1995 fourteenth southern biomedical engineering conference, Shreveport, LA, 1995, p. 245–7.
 4.
Nam D, Weeks KD, Reinhardt KR, Nawabi DH, Cross MB, Mayman DJ. Accelerometerbased, portable navigation vs imageless, largeconsole computerassisted navigation in total knee arthroplasty. A comparison of radiographic results. J Arthroplasty. 2013;28:255–61.
 5.
Wu G, Van Der Helm FCT, Veeger HEJ, Makhsous M, Van Roy P, Anglin C, et al. ISB recommendation on definitions of joint coordinate systems of various joints for the reporting of human joint motion—part II: shoulder, elbow, wrist and hand. J Biomech. 2005;38:981–92.
 6.
Giroux B, Lamontagne M. Net shoulder joint moment and muscular activity during light weighthandling at different displacements and frequencies. Ergonomics. 1992;35:385–403.
 7.
Meskers CGM, Van Der Helm FCT, Rozendaal LA, Rozing PM. In vivo estimation of the glenohumeral joint rotation center from scapular bony landmarks by linear regression. J Biomech. 1997;31:93–6.
 8.
Cappozzo A. Gait analysis methodology. Hum Mov Sci. 1984;3:27–50.
 9.
Stokdijk M, Nagels J, Rozing PM. The glenohumeral joint rotation centre in vivo. J Biomech. 2000;33:1629–36.
 10.
Campbell AC, Lloyd DG, Alderson JA, Elliott BC. MRI development and validation of two new predictive methods of glenohumeral joint centre location identification and comparison with established techniques. J Biomech. 2009;42:1527–32.
 11.
Monnet T, Desailly E, Begon M, Vallée C, Lacouture P. Comparison of the SCoRE and HA methods for locating in vivo the glenohumeral joint centre. J Biomech. 2007;40:3487–92.
 12.
Lempereur M, Leboeuf F, Brochard S, Rousset J, Burdin V, RémyNéris O. In vivo estimation of the glenohumeral joint centre by functional methods: accuracy and repeatability assessment. J Biomech. 2010;43:370–4.
 13.
Nikooyan AA, van der Helm FCT, Westerhoff P, Graichen F, Bergmann G, Veeger HEJ. Comparison of two methods for in vivo estimation of the glenohumeral joint rotation center (GHJRC) of the patients with shoulder hemiarthroplasty. PLoS ONE. 2011;6:e18488.
 14.
Poppen NK, Walker PS. Normal and abnormal motion of the shoulder. J Bone Jt Surg Am. 1976;58:195–201.
 15.
Veeger HEJ. The position of the rotation center of the glenohumeral joint. J Biomech. 2000;33:1711–5.
 16.
Calderone M, Cereatti A, Conti M, Della Croce U. Comparative evaluation of scapular and humeral coordinate systems based on biomedical images of the glenohumeral joint. J Biomech. 2014;47:736–41.
 17.
Cereatti A, Camomilla V, Cappozzo A. Estimation of the centre of rotation: a methodological contribution. J Biomech. 2004;37:413–6.
 18.
Ehrig RM, Taylor WR, Duda GN, Heller MO. A survey of formal methods for determining the centre of rotation of ball joints. J Biomech. 2006;39:2798–809.
 19.
Camomilla V, Cereatti A, Vannozzi G, Cappozzo A. An optimized protocol for hip joint centre determination using the functional method. J Biomech. 2006;39:1096–106.
 20.
Piazza SJ, Okita N, Cavanagh PR. Accuracy of the functional method of hip joint center location: effects of limited motion and varied implementation. J Biomech. 2001;34:967–73.
 21.
Cereatti A, Donati M, Camomilla V, Margheritini F, Cappozzo A. Hip joint centre location: an ex vivo study. J Biomech. 2009;42:818–23.
 22.
Hicks JL, Richards JG. Clinical applicability of using spherical fitting to find hip joint centers. Gait Posture. 2005;22:138–45.
 23.
Sabatini AM. Estimating threedimensional orientation of human body parts by inertial/magnetic sensing. Sensors. 2011;11:1489–525.
 24.
McGinnis RS, Perkins NC. Inertial sensor based method for identifying spherical joint center of rotation. J Biomech. 2013;46:2546–9.
 25.
Crabolu M, Pani D, Cereatti A. Evaluation of the accuracy in the determination of the center of rotation by magnetoinertial sensors. In: 2016 IEEE sensors applications symposium (SAS), Catania, 2016, p. 1–5.
 26.
Crabolu M, Pani D, Raffo L, Cereatti A. Estimation of the center of rotation using wearable magnetoinertial sensors. J Biomech. 2016;49:3928.
 27.
Seel T, Schauer T, Raisch J. Joint axis and position estimation from inertial measurement data by exploiting kinematic constraints. In: 2012 IEEE international conference on control applications, Dubrovnik, 2012, p. 45–9.
 28.
van den Noort JC, Wiertsema SH, Hekman KMC, Schönhuth CP, Dekker J, Harlaar J. Measurement of scapular dyskinesis using wireless inertial and magnetic sensors: importance of scapula calibration. J Biomech. 2015;48:3460–8.
 29.
Picerno P, Cereatti A, Cappozzo A. A spot check for assessing static orientation consistency of inertial and magnetic sensing units. Gait Posture. 2011;33:373–8.
 30.
Cereatti A, Rosso C, Nazarian A, DeAngelis JP, Ramappa AJ, Della Croce U. Scapular motion tracking using acromion skin marker cluster: in vitro accuracy assessment. J Med Biol Eng. 2015;35(1):94–103.
Authors’ contributions
All the authors substantially contributed to the manuscript. In particular, MCr coded the algorithms and performed data analysis, AC, DP and MCr defined the scientific methodology to be followed, AC supervised the biomechanical aspects whereas DP and LR were responsible for the signal processing. MCo and PC designed and managed the MRI acquisition and preprocessing. All the authors drafted, reviewed the final version of the manuscript. All authors read and approved the final manuscript.
Acknowledgements
The authors gratefully acknowledge Dr. Manuela Calderone for the 3D reconstructions of the scapular and humeral bones by AMIRA image processing software.
Competing interests
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential competing interests.
Availability of data and materials
The signals form MIMUs analyzed in the study are available upon reasonable request along with the measurements of the GHJC from MRI images. MRI images are not directly provided.
Consent for publication
The participants acknowledged their consent to publish the acquired data and pictures.
Ethics approval and consent to participate
The study has been performed following the principles outlined in the Helsinki Declaration of 1975, as revised in 2000, on healthy subjects. All the volunteers gave their informed consent to the measurements.
Author information
Rights and permissions
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.
About this article
Received
Accepted
Published
DOI
Keywords
 Human movement
 Center of rotation
 Magnetoinertial sensing
 Accelerometers
 Gyroscope
 Wearable devices
 Shoulder
 Glenohumeral joint
 Functional method