Optical Enhancement of Exoskeleton-Based Estimation of Glenohumeral Angles.

In Robot-Assisted Rehabilitation (RAR) the accurate estimation of the patient limb joint angles is critical for assessing therapy efficacy. In RAR, the use of classic motion capture systems (MOCAPs) (e.g., optical and electromagnetic) to estimate the Glenohumeral (GH) joint angles is hindered by the exoskeleton body, which causes occlusions and magnetic disturbances. Moreover, the exoskeleton posture does not accurately reflect limb posture, as their kinematic models differ. To address the said limitations in posture estimation, we propose installing the cameras of an optical marker-based MOCAP in the rehabilitation exoskeleton. Then, the GH joint angles are estimated by combining the estimated marker poses and exoskeleton Forward Kinematics. Such hybrid system prevents problems related to marker occlusions, reduced camera detection volume, and imprecise joint angle estimation due to the kinematic mismatch of the patient and exoskeleton models. This paper presents the formulation, simulation, and accuracy quantification of the proposed method with simulated human movements. In addition, a sensitivity analysis of the method accuracy to marker position estimation errors, due to system calibration errors and marker drifts, has been carried out. The results show that, even with significant errors in the marker position estimation, method accuracy is adequate for RAR.

1. Introduction

The application of robotics and Virtual Reality (VR) to motor neurorehabilitation (Figure 1) has been beneficial for patients, as they receive intensive, repetitive, task-specific, and interactive treatment [1-4].

Figure 1

Robotic and VR-based rehabilitation.

The assessment of (a) patient movement compliance with the prescribed exercises and (b) patient long-term improvement is critical when planning and evaluating the efficacy of RAR therapies. In order to obtain the patient motion data to conduct the said assessments, one has to estimate patient posture (i.e., the joint angles of the limbs). Patient posture estimation methods need to be practical and easy to set up for the physician, so that the said assessments can indeed be an integral part of the therapy.

Current methods for estimating patient posture are either cumbersome or not accurate enough in exoskeleton-based therapies. In order to overcome such limitations, we propose a method where low-cost RGB-D cameras (which render color and depth images) are directly installed in the exoskeleton and colored planar markers are attached to the patient's limb to estimate the angles of the GH joint, thereby overcoming the individual limitations of each of these systems.

2. Literature Review

Optical, electromagnetic, and inertial MOCAPs have been used in many rehabilitation scenarios for accurate posture estimation [5]. However, the use of the said MOCAPs in exoskeleton-based rehabilitation is limited by the factors discussed below:

Optical marker-based systems (e.g., Optotrak, CODA, Vicon) are considered the most accurate for human motion capture [5]. Reference [6] reports Optotrak errors of 0.1–0.15 mm. However, in the specific case of exoskeleton-based therapy, these systems require redundant sensors and markers to cope with occlusions caused by the exoskeletal body. Therefore, their specific usage for therapy is limited. Besides, the cost of these systems is high (50 K–300 K USD [7]) compared to nonoptical MOCAPs.

Electromagnetic systems do not suffer from optical occlusions. However, they are easily perturbed by surrounding metallic objects (e.g., exoskeletal body) and electric/magnetic fields [5]. An additional drawback of these systems is their limited detection volume when compared to optical systems.

Inertial and Magnetic Measurement Systems are robust, handy, and economical for full-body human motion detection (upper limb tracking in [8, 9]). With the use of advanced filtering techniques, inertial sensor drift errors are reduced and a dynamic accuracy of 3 deg. RMS [5] is achieved. However, these systems require patients to perform calibration motions/postures, which may not be suitable for those with neuromotor impairments.

In exoskeleton-based rehabilitation, the prevailing approach to estimate human limb joint angles (e.g., [10-13]) is to approximate them with the angles of the exoskeleton joints. However, misalignment between the axes of the exoskeleton and human joints may produce large estimation errors [14, 15]. Accurate estimation of GH joint angles is hard to achieve using this approach, since it requires an exoskeleton with a complex kinematic structure that considers the concurrent motion of the sternoclavicular and acromioclavicular joints.

Recognizing the differences in the kinematic structures of the limb and exoskeleton, [16] presents a computational method which considers the limb and exoskeleton parallel kinematic chains related by the cuff constraints joining them together. Then, the IK problem of the parallel kinematic chain can be solved to find the limb joint angles. A limitation of this method is that its performance has been demonstrated solely for analytic (1-DOF) movements of the elbow and wrist joints. The estimation accuracy of the GH joint angles has yet to be determined.

Reference [17] presents a computational method based on the estimation of the arm swivel angle (which parametrizes arm posture) for exoskeleton-based therapy. The arm IK is solved with a redundancy resolution criterion that chooses a swivel angle that allows the subject to retract the palm to the head efficiently. The approach in [17] extends their previous work in [18, 19] by considering the influence of the wrist orientation on the swivel angle estimation. Although the error of the swivel angle estimation (mean error ≈ 4 deg.) has been reported for compound movements [17], individual errors in the wrist, elbow, and GH joint angles are not indicated.

Reference [20] extends the method in [17] to estimate the wrist angles and assesses its performance for compound movements (mean RMSE ≈ 10 deg. in the swivel angle estimation). Reference [20] reports the individual errors of the arm joint angles solely for the movement task where the swivel angle was best estimated (mean RMSE ≈ 5 deg. in the swivel angle estimation). No errors of the arm joint angles were discussed for the other cases. A limitation of the work in [20] is that the MOCAP used to obtain the reference angles to assess their method performance is a custom-made inertial system with no reported measurement accuracy.

2.1. Conclusions of the Literature Review

We remind the reader that the general context of this paper is the estimation of the GH joint angles.

As per our literature review, no MOCAPs have been developed for the specific scenario of exoskeleton-based rehabilitation. Even if current MOCAPs and the exoskeleton could be set up for simultaneous use (e.g., [15, 16]), the setup protocol and operation are intricate and conflicting with the usual time and resources available for patient treatment.

Exoskeleton-based posture estimations present limitations in their accuracy due to kinematic mismatch of the limb and exoskeleton [15, 16].

The accuracy of the GH joint angle estimations provided by computational methods in [16, 17] is unknown. Reference [20] extends the work in [17] by estimating the wrist angles. Reference [20] solely reports the estimation accuracy of the GH angles for the best-case scenario and the precision of its ground-truth is not indicated.

2.2. Contributions of This Paper

In response to the limitations discussed in the estimation of patient joint angles in exoskeleton-based therapy (Sections 2 and 2.1), this paper introduces a hybrid approach to estimate, in real-time, the GH joint angles. This hybrid system is composed of a low-cost marker-based vision system and the rehabilitation robot, overcoming the individual limitations of its constitutive subsystems:

Occlusions are minimized, which are a major limitation of optical systems.

Accuracy of joint angle estimation is improved, which is a major limitation of exoskeleton-based systems.

This paper presents the implementation and assessment of our method using simulated human motion data. In addition, a sensitivity analysis of our method accuracy to marker position estimation errors is carried out.

We have considered the following scenarios of application for the proposed method in the RAR domain:

Precise estimation of GH joint angles during rehabilitation or evaluation sessions of GH joint analytic movements.

Acquisition of GH joint movement data enabling validation and improvement of other posture estimation methods without using expensive redundant optical MOCAPs.

3. Methods

3.1. Problem Definition

This section presents the problem of estimating the patient limb GH joint angles during the GH joint RAR using the proposed hybrid motion capture system (a detailed version of the problem definition is presented in the Appendix). This problem can be stated as follows.

Given. Consider the following:

Patient: (a) the kinematic model (e.g., the Denavit-Hartenberg parameters [21]) of the human upper limb (H) (Figure 2(a)).

Figure 2

Components of the GH joint angles estimation system: (a) human kinematic model, (b) exoskeleton kinematic model, (c) marker-based optical motion capture system, and (d) hybrid GH joint angles estimation system.

Exoskeleton: (a) the kinematic model of the exoskeleton (E) and (b) the exoskeleton joint angles at any instant of the therapy (v(t)) (Figure 2(b)).

Marker-based optical motion capture system (R): (a) color and depth information captured by the RGB-D cameras installed in the exoskeleton links and (b) geometry and color of the markers attached to the patient upper limb (Figure 2(c)).

Goal. The goal is to estimate the patient GH joint angles (v(t)) with minimum error during the GH joint rehabilitation exercises.

3.2. Kinematic Models

This section discusses the main features of the kinematic models of the human limb and exoskeleton used for the posture estimation method.

3.2.1. Kinematic Model of the Human Upper Body

The human kinematic model is denoted by H(L, J), where L and J are the sets of links and joints, respectively. We use the human upper body model presented in [16] (Figure 2(a)), which includes joints of the spine, scapuloclavicular system, and arm. The upper limb is modeled with 9 DOFs: 2 DOFs of the scapuloclavicular system, 3 DOFs of the GH joint (spherical joint), 2 DOFs of the elbow, and 2 DOFs of the wrist (see further details in the Appendix). This model presents the following advantages:

It can be easily implemented in robotic simulators and similar tools.

It is suitable for simulating human-robot interaction in real-time [16].

The spherical model of the GH joint avoids limitations of other representations of such joint, like the Gimbal lock that occurs when using the three concurrent and orthogonal 1-DOF revolute joints' model [22].

3.2.2. Kinematic Model of the Exoskeleton

The exoskeleton kinematic model is denoted by E(L, J), where L and J are the sets of links and joints, respectively. In this research, the rehabilitation exoskeleton used is the Armeo Spring (Figure 2(b)), which is a passive system that supports the weight of the patient's arm [23] with springs. The Armeo kinematic structure includes rotational joints (equipped with encoders [24, 25]) and prismatic joints (enabling exoskeleton adjustment to the size of each patient). We use the Armeo Spring kinematic model presented in [16], which includes both types of joints (see further details in the Appendix).

3.3. GH Joint Angles Estimation Method

The aim of the method is to estimate the GH joint angles with respect to (w.r.t.) a coordinate system (CS) attached to the scapuloclavicular system. Figure 2(d) shows the proposed system for the GH joint angle estimation. Our approach is based on the estimation of the upper arm orientation w.r.t. the acromion (Figure 3(a)). According to such requirements, the rationale to install the markers of the optical MOCAP R is as follows:

Figure 3

(a) Schematic diagram of the hybrid GH joint angles estimation system and (b) high-level operation of the system.

Marker m0 is rigidly installed in the acromion, so the estimated upper arm orientation can be expressed w.r.t. the m0 CS (and therefore w.r.t. the scapuloclavicular system).

Marker m1 is rigidly installed in the upper arm, so that all the rotations of the upper arm are captured by m1. The region that was chosen to attach m1 to the upper arm by using a custom-made fixation (Figure 2(d)) is the distal part of the humerus (near the elbow). Elbow rotations do not affect the orientation of m1.

Reference [26] reports a five-marker installation procedure. This reference explicitly mentions five markers as an acceptable number for clinical upper limb tracking. In this paper, we report the usage of two markers for upper arm tracking. It is not possible to compare the performance of the marker placement protocol proposed here with the one in [26] because the work in [26] addresses (a) non-RAR scenarios, (b) tracking of the entire upper limb, and (c) protocol sensitivity w.r.t. its application on the dominant/nondominant arm and w.r.t. the age of test subjects. However, the work in [26] helps to establish the number of markers compatible with the clinical application of upper limb tracking.

The cameras of the optical motion capture system R are rigidly attached (using custom-made supports) to exoskeleton links so that camera r0 detects marker m0 and camera r1 detects marker m1 during the GH joint training. Camera r0 is mounted on link l0 and camera r1 is mounted on link l8 (Figure 3(a)).

The cameras used in our system are of low cost. Commercial cameras that present similar specifications to the ones simulated here (Table 1) are Intel® SR300 (99 USD) [27, 28], DepthSense® 525 (164 USD) [29, 30], and CamBoard picoS (690 USD) [28, 31].

Table 1

Vision sensor features.

Color camera resolution (px)	128 × 128
Depth camera resolution (px)	128 × 128
Field of view (deg.)	Horizontal = 45; vertical = 45
Minimum sensing distance (meters)	0.05
Maximum sensing distance (meters)	0.3

Figure 3(b) shows an overview of the operation of the estimation method. In order to estimate the upper arm pose, the poses of the markers need to be expressed w.r.t. a common CS. A suitable CS to conduct such estimation is the exoskeleton base.

A summary of the steps to estimate the GH joint angles is as follows:

Estimate the pose of the markers w.r.t. the cameras.

Estimate the pose of the cameras w.r.t. the exoskeleton.

Estimate the pose of the markers w.r.t. the exoskeleton.

Estimate the upper arm pose w.r.t. the exoskeleton.

Refer the GH joint angles w.r.t. the acromion (marker m0 CS).

The details of the mentioned steps are presented in the following sections.

3.3.1. Estimation of the Pose of the Markers w.r.t. the Cameras

The purpose of this step is to estimate the position and orientation of the markers (Figure 4) w.r.t. the CSs of the cameras using the color and depth images provided by each camera r:

Figure 4

Schematic diagram of the iterative estimation of the pose of the markers.

The RGB image is Ic (A × B pixels). The pixel coordinates (u, v) take values 0 ≤ u ≤ A − 1 and 0 ≤ v ≤ B − 1. C (1 × 3∗A∗B) contains the RGB color associated with each pixel (u, v) ∈ Ic.

The depth image associated with the scene in Ic is Id (L × N pixels); L ≤ A and N ≤ B. The pixel coordinates (u, v) in Id take values 0 ≤ u ≤ L − 1 and 0 ≤ v ≤ N − 1. The CS of images Ic and Id is coincident. D (1 × L∗N∗3) contains the (x, y, z) coordinates of the object in each pixel (u, v) ∈ Id w.r.t. the r CS.

The pose estimation of the markers w.r.t. the cameras is based on the reconstruction of the 3D position of the colored disks on the markers. The following steps are taken to estimate the marker pose:

Estimation of disk coordinates in color image (Figure 5): the purpose of this step is to find the approximate (u, v) coordinates of the centers of the marker disks in image Ic. The following steps are carried out:

Figure 5

Estimation of disk coordinates in color image. (a) Simulated RBG image, (b) result of the color segmentation (zoomed image), and (c) result of the blob extraction (zoomed image).

Color segmentation in image Ic: image regions containing the colors of the marker disks are preserved and the other regions are colored in white. The resulting image is defined as Isc.

Blob extraction on image Isc: blob extraction consists of finding the connected regions in the image Isc sharing the same color and labeling them according to their color.

Disk center coordinates estimation: for each j (j = 0,…, n) blob extracted from Ic, the position of the center of a bounding box for the blob is obtained. This point approximates the actual center of disk p (Figure 5). The resulting set of the approximate coordinates of disk centers in Ic is . The ℤ2 center coordinates are referenced w.r.t. the internal image CS. Blobs are extracted with standard connected-component labeling algorithms.

Estimation of disk coordinates in the camera r CS: this step converts disk coordinates in the internal image CS into the ℝ3 ones w.r.t. the r sensor CS, as follows:

Convert the positions (u, v) of the disk centers in set into the image Id CS. The CSs of images Ic and Id match. Hence,

Compute the indices a of the (x, y, z) coordinates of point in array D as follows:

The point contains the (x, y, z) coordinates of point w.r.t. the r CS. The coordinates of point are obtained as follows:

The approximate marker disk centers detected by camera r form the set .

Computation of the marker m CS in the r camera CS: an SO(3) coordinate frame is attached to each marker:

Make

Use the four disk centers in the marker (Figure 5) as follows:

The submatrix is normalized to guarantee its SO(3) nature. The frame describes the estimated pose of marker m w.r.t. the CS of the camera r.

3.3.2. Estimation of the Pose of the Cameras w.r.t. the Exoskeleton

The goal of this step is to find the transformation T, which expresses the pose of the camera r w.r.t. the base of the exoskeleton (Figure 6).

Figure 6

Schematic diagram of the iterative estimation of the pose of the cameras.

The rigid transformation matrices T and T ∈ ℝ4×4, which describe the pose of the cameras r w.r.t. the CS of the link where they are installed, are estimated during system calibration (the calibration matrix can be obtained by camera detection of a 2D/3D calibration object mounted on a known location of the exoskeleton). The poses T and T of the exoskeleton links l0 and l8 w.r.t. to the exoskeleton base CS are computed using the Forward Kinematics of exoskeleton E. Then, T and T are estimated as follows:

3.3.3. Estimation of the Pose of the Markers w.r.t. the Exoskeleton

The objective of this step is to estimate the transformation (T) that describes the pose of marker m w.r.t. the exoskeleton base CS (Figure 7). Transformations T are estimated as follows:

Figure 7

Schematic diagram of the iterative estimation of the pose of the markers w.r.t. the exoskeleton CS.

3.3.4. Estimation of the Upper Arm Pose w.r.t. the Exoskeleton

The purpose of this step is to estimate the upper arm pose (Tarm) w.r.t. the exoskeleton base CS using the marker poses T (Figure 8). The upper arm direction vector is computed from the estimated position of the end-points of the upper arm (GH and elbow joint centers) as follows (CSs in Figure 9):

Figure 8

Schematic diagram of the iterative estimation of the upper arm pose.

Figure 9

Coordinate systems for the upper arm pose estimation.

Estimate the position of the GH joint center: the rigid transformation matrix T, which expresses the pose of the GH joint CS w.r.t. the m0 CS, is estimated during the calibration process of the system. Hence, the GH joint center is estimated as follows:

Estimate T, which is the pose of the GH joint CS w.r.t. the exoskeleton E base CS (see (8)).

Extract p from T. The point p is the position of the center of the GH joint seen from the E CS:

Estimate the position of the elbow joint center: the rigid transformation matrix Telw (elbow joint CS w.r.t. the m1 CS) is estimated during the calibration process of the system. Hence, the elbow joint center is computed as follows:

Estimate Telw, which is the pose of the elbow joint CS w.r.t. the exoskeleton E base CS (see (9)).

Extract pelw from Telw. The point pelw is the position of the center of the elbow joint seen from the E CS:

Estimate the upper arm position:

Estimate the arm direction vector as .

Estimate the origin of the upper arm CS as .

Estimate the upper arm orientation: the estimated orientation of the upper arm is computed using Euler angle yxz decomposition w.r.t. the base CS of exoskeleton E:

Estimate the rotation of the arm around the y-axis of the E CS using the projection of on the x-z plane of the fixed E CS.

Compute the rotation of the arm around the mobile x-axis of E CS from the inner product of with the mobile z-axis of E CS.

Estimate the rotation of the upper arm around its longitudinal axis as the rotation of the marker m1 around vector . This angle is the one between (i) the mobile x-axis of E CS and (ii) the projection of x-axis of marker m1 CS onto the x-y plane of E CS.

Express the pose of the upper arm w.r.t. the E base CS as the 4 × 4 rigid transformation Tarm.

3.3.5. Refer the Angles of the GH Joint w.r.t. the Acromion

Since m0 is rigidly attached to the acromion, the upper arm orientation can be expressed w.r.t. the acromion by using the inverse of T:

3.4. Implementation and Simulation

The arm posture estimation method was implemented by using the V-REP robotics simulator [32]. In the simulator, the scene in Figure 2(d) is created, which includes the models of (a) a human patient, (b) an Armeo Spring, (c) the RGB-D vision sensors with the couplings to attach them to the exoskeleton, and (d) the planar markers with the couplings to attach them to the human arm. The configuration of the simulated vision sensors is summarized in Table 1.

For the estimation of the coordinates of disk centers P in the image Ic, color segmentation and blob detection algorithms available in the simulator were used. Additional code was written to sort blob centers by color. All additional code was written in LUA (lightweight embeddable scripting language) scripts.

3.4.1. Generation of the Ground-Truth Poses of the Patient Upper Limb during RAR

The accuracy of the proposed method is determined by comparing its estimations of the upper arm poses with the ones of the simulated human patient (ground-truth values of Tarm). To generate movements of the simulated patient that resemble the ones of therapy, we performed the next steps:

Armeo movement generation: we recorded 4 time sequence datasets of the actual Armeo joint measurements (sampled at 66.6 Hz) while performing the following shoulder movements (Figure 10): (a) shoulder horizontal abduction-adduction (SAbAd), (b) shoulder flexion-extension (SFE), (c) shoulder internal rotation (SIR), and (d) a combination of all the mentioned movements (COMB). These movement history datasets are used to guide a simulation of the Armeo model.

Figure 10

GH joint movements: (a) shoulder flexion-extension (SFE), (b) shoulder horizontal abduction-adduction (SAbAd), and (c) shoulder internal rotation (SIR).

Patient movement generation: the movements of the patient upper limb that correspond to the recorded movements of the Armeo are computer-generated with the method in [16]. The said method provides an estimation of the patient posture given the joint angles of the exoskeleton by using an inverse kinematics approach.

In this way, four sets (one per movement dataset) of known poses of the upper arm are obtained by simulating patient movement and compared here against those estimated with our method. Our method accuracy is assessed without compensating any time offsets between the reference and estimated angles. In this way, real-time accuracy of the method is assessed. Table 2 presents the approximate amplitudes of the yxz Euler angle decomposition of the GH joint movements of the simulated patient w.r.t. its local CS.

Table 2

Movement dataset features.

Movement dataset	Amplitude (deg.)	Samples
SAbAd	(6°, 31°, 10°)	1000
SFE	(31°, 8°, 1°)	1000
SIR	(3°, 3°, 34°)	1000
COMB	(40°, 90°, 60°)	2000

3.4.2. Measurement of the Estimation Performance

Error in the estimation of the markers position: the error in the position estimation of markers m is computed as the RMS of expression , where i ∈ {0,1}.

Error in the estimation of the arm pose: the error in the arm position estimation for a GH joint movement dataset (eposarm) is computed as the RMS of for all samples in the movement dataset.

To quantify the error in the arm orientation estimation (eoriarm), the next steps are carried out:

Compute the matrix of rotation error , where Rotarm and are the rotation submatrices of transformation matrices Tarm and , respectively.

Express Roterror in exponential map notation [22] as .

Compute eoriarm as the RMS of for all samples in the movement dataset.

3.5. Sensitivity Analysis

A sensitivity analysis is carried out to study the influence of relevant parameters on the method accuracy. Formally, the sensitivity analysis determines the effect of the perturbation of the parameter Q on the objective function F(Q). The relative sensitivity of F(Q) w.r.t. Q, S, is given by (11) [33]. The value of S is the ratio (dimensionless) between the percentual changes in F and Q:

The upper arm pose accuracy (and, therefore, that of the GH joint angles) relies on the precise estimation of the position of the centers of the elbow and GH joints ( and ) (Section 3.3.4), which ultimately depend on the following transformations involving the markers:

and (markers w.r.t. exoskeleton).

T and Telw (GH and elbow joints w.r.t. markers).

The conducted sensitivity analysis focuses on errors in T and Telw, given that errors in the estimation of and (Section 4.2) are small. Possible causes of errors in T and Telw are as follows:

Inaccurate computation of T and Telw during the system calibration.

Relative displacement of the markers w.r.t. the GH and elbow joints due to skin movement.

In the sensitivity analysis, translations errors in matrices T and Telw are induced by disturbing the location of the markers m (k = [0,1]) w.r.t. the CSs of the GH and elbow joints. Since orientation information in T and Telw is not used to estimate the upper arm pose, it is excluded from the sensitivity analysis.

For the sensitivity analysis (see (11)), the vector-valued function F(q) quantifies the estimation error of the arm position and orientation (see (12)) and the parameter set q represents the marker translation errors. The parameter set q is defined as q = {q1, q2, q3, q4, q5, q6}, where each q ∈ q is a scalar representing the magnitude of a translation of a specific marker along a prescribed direction. Table 3 describes the meaning of each parameter in set q:

Table 3

Parameters of function F(q) (error in the position and orientation estimation of the upper arm (see (12))) to study in the sensitivity analysis.

Parameter	Meaning	CS of reference
q 1	Translation with magnitude ‖q1‖ of m0 along x-axis	GH joint
q 2	Translation with magnitude ‖q2‖ of m0 along y-axis	GH joint
q 3	Translation with magnitude ‖q3‖ of m0 along z-axis	GH joint
q 4	Translation with magnitude ‖q4‖ of m1 along x-axis	Elbow joint
q 5	Translation with magnitude ‖q5‖ of m1 along y-axis	Elbow joint
q 6	Translation with magnitude ‖q6‖ of m1 along z-axis	Elbow joint

The sensitivity analysis procedure (Figure 11) entails the following steps:

Figure 11

Sensitivity analysis steps.

Load the movement dataset of the GH joint to test (SFE, SAbAd, SIR, and COMB).

Select the parameter q ∈ q to perturb (selection of a marker and a direction of translation). Marker m0 translates along axes of the GH joint CS. Marker m1 translates along axes of the elbow joint CS (Figure 12).

Figure 12

Sensitivity analysis. Coordinate systems of reference for the translations of (a) marker m0 and (b) marker m1.

Apply the translation indicated by q to the corresponding marker. The marker perturbation q is applied for the complete movement dataset.

Compute the estimation errors of the upper arm position and orientation F(q) = (eposarm(q), eoriarm(q)) as the simulated patient moves according to the chosen GH joint movement dataset. The current iteration of the process is indicated by index i.

Compute the position and orientation components of S as per (11). The derivative of F(q) w.r.t. q is given by (13). The required derivatives are computed numerically [34, 35]:

Increment q by Δq and go to step (3). Repeat the process until the desired number of iterations i of the procedure is reached.

The complete sensitivity analysis was performed for each movement dataset (SFE, SAbAd, SIR, and COMB). The directions in which marker translations occur (Table 3) are chosen so that the markers do not leave the detection volume of the cameras. Table 4 summarizes the parameters of the sensitivity analysis. Translation units are in meters (mts).

Table 4

Parameters of the sensitivity analysis.

Minimum marker translation qmin (mts)	0
Maximum iterations of the sensitivity analysis imax	10
Increment of marker translation in each iteration Δq (mts)	0.002
Movement datasets evaluated	4

4. Results and Discussion

This section presents and discusses the results of (a) estimation accuracy of the marker 3D position, (b) estimation accuracy of the upper arm pose, and (c) sensitivity analysis of the estimation accuracy of the upper arm pose w.r.t. translation errors in T and Telw.

4.1. Results of Marker Position Estimation

Table 5 presents the RMS of the estimation errors of the position of the markers m per movement dataset. The mean RMS errors of the position estimation of m0 and m1 for all movement datasets are 0.00083 and 0.00208 mts, respectively.

Table 5

RMS of errors (and standard deviation in parentheses) in the position estimation of markers m in the datasets of GH joint movements.

Movement	m 0 [mts]	m 1 [mts]
SAbAd	0.00089 (0.0001)	0.00175 (0.001)
SFE	0.00060 (0.0002)	0.00197 (0.0008)
SIR	0.00088 (0.0001)	0.00135 (0.0007)
COMB	0.00097 (0.0003)	0.00324 (0.002)

Figure 13 shows the box plots of the estimation errors in the marker positions for all movement datasets. A greater variation in the position estimation accuracy of marker m1, in comparison to that of m0, is observed. We have attributed this to (a) the higher linear and rotational velocities and likewise (b) the larger translations and rotations that m1 undergoes compared to m0.

Figure 13

Box plots of estimation errors in markers position and upper arm position and orientation for all movement datasets.

4.2. Results of Upper Arm Pose Estimation

The RMS of errors in the upper arm pose estimation are presented in Table 6. By averaging the results of all movement datasets, errors of 0.00110 mts and 0.88921 deg. in the upper arm position and orientation estimation are obtained. Figure 13 shows the box plots of the estimation errors in the upper arm position and orientation for all movement datasets.

Table 6

RMS (and standard deviation in parentheses) of errors in the upper arm position and orientation estimation in the assessed movement datasets.

Movement	Position [mts]	Orientation [deg.]
SAbAd	0.00109 (0.0005)	0.92039 (0.4842)
SFE	0.00094 (0.0004)	0.83796 (0.3763)
SIR	0.00091 (0.0002)	0.73465 (0.4156)
COMB	0.00145 (0.0008)	1.0638 (0.5238)

In motor rehabilitation, angular errors in the range of 3–5 degrees are considered acceptable for mobility evaluation of patients [6, 36, 37]. Figure 13 shows that our arm orientation estimation accuracy is adequate for exoskeleton-assisted rehabilitation.

4.3. Results of the Sensitivity Analysis

The results of the sensitivity analysis per movement dataset of the shoulder are presented in Figures 14, 15, 16, and 17. In each figure, the following subfigures are presented:

Figure 14

Results of the sensitivity analysis with the SAdAd movement dataset (q: m0-x movement/m0-y movement/m0-z movement/m1-x movement/m1-y movement/m1-z movement).

Figure 15

Results of the sensitivity analysis with the SFE movement dataset (q: m0-x movement/m0-y movement/m0-z movement/m1-x movement/m1-y movement/m1-z movement).

Figure 16

Results of the sensitivity analysis with the SIR movement dataset (q: m0-x movement/m0-y movement/m0-z movement/m1-x movement/m1-y movement/m1-z movement).

Figure 17

Results of the sensitivity analysis with the COMB movement dataset (q: m0-x movement/m0-y movement/m0-z movement/m1-x movement/m1-y movement/m1-z movement).

Error in upper arm position estimation (eposarm) versus total marker translation (q): this figure shows the evolution of the absolute error in the upper arm position estimation as the error in the translation components of matrices T and Telw increases.

Error in upper arm orientation estimation (eoriarm) versus total marker translation (q): this figure shows the evolution of the absolute error in the upper arm orientation estimation as the error in the translation components of matrices T and Telw increases.

Position component of S versus total marker translation (q): this figure shows the evolution of the relative sensitivity metric corresponding to the error in the upper arm position estimation as the error in the translation components of matrices T and Telw increases.

Orientation components of S versus total marker translation (q): this figure shows the evolution of the relative sensitivity metric corresponding to the error in the upper arm orientation estimation as the error in the translation components of matrices T and Telw increases.

4.3.1. Sensitivity in Arm Position Estimation

Regarding the arm position estimation, one can observe that translations of marker m0 produce larger absolute errors than translations of marker m1. This difference is due to the fact that the translations of m0 produce a larger change in when compared to the one produced by translations of m1. Note that since is computed by using , any modification in directly affects the accuracy of .

Observing the behavior of the position component of S, one can conclude that all translations of the markers m0 and m1 contribute similarly to the error in the arm position estimation. The curves obtained for the position component of S resemble a logarithmic function with an asymptote along the value 1 of the ordinate axis. A value of 1 in the magnitude of the position component of S means that a percentage change in the magnitude of the marker translation produced the same percentage change (also matching the sign) in the magnitude of the error in the arm position estimation.

4.3.2. Sensitivity in Arm Orientation Estimation

In Figures 14, 15, 16, and 17, one can observe that the translations of marker m1 produce larger absolute errors in the upper arm orientation estimation when compared to those produced by translations of marker m0. Notice that the x-axis and z-axis of the elbow joint CS are always perpendicular to the upper arm vector () (Figure 12(b)). When the position of m1 is perturbed along the said axes, the angle between (i) the actual upper arm vector () and (ii) the estimated upper arm vector () (which is inaccurate due to the perturbation of the marker position) is maximal.

A side effect of the marker position perturbation is that the marker m suffers modifications of scale and changes in the level of perspective distortion in the images of camera r, affecting the accuracy of the system. This situation can be observed in Figures 14(b), 15(b), 16(b), and 17(b), where translations of m1 along the y-axis of the elbow joint CS should not produce variations in the orientation estimation error. However, on the contrary, slight variations in the accuracy of the orientation estimation are indeed present in the mentioned figures.

4.3.3. Robustness of the Upper Arm Pose Estimation Method

In Figures 14(c), 14(d), 15(c), 15(d), 16(c), 16(d), 17(c), and 17(d) one can observe that the position component of S increases faster than the orientation component of S. The behavior of S observed remains across the datasets used. Hence, the orientation estimation of the upper arm is more robust than the position estimation w.r.t. errors in the translational components of matrices T and Telw.

The results of the sensitivity analysis show that the assumption that transformations T and Telw are rigid is reasonable. Even with marker drifts of 0.02 mts, the GH joint angles can be estimated with an accuracy (RMSE 3.6 deg.) appropriate for the mobility evaluation of patients (in the range of 3–5 deg.).

Marker drifts must be mitigated by the marker attachments to the human body. Furthermore, marker attachments should be designed to minimize the effect of errors in T and Telw on the method accuracy. For example, notice how the attachment of marker m1 (Figure 12(b)) locates marker m1 with an offset w.r.t. the elbow joint center along the direction which least affects the upper arm orientation estimation.

The results presented suggest that the method we implemented is a feasible alternative for estimating the GH joint angles in a RAR scenario.

4.4. Comparison to Related Works

The literature review provided no references other than [16-20] for upper limb posture estimation (including the GH joint) in exoskeleton-based rehabilitation using computational methods. Among the mentioned works, only [20] reports the errors (mean RMSE 4.8 deg.) in the GH joint angles estimation. Reference [20] reports RMSE values of the GH joint angles only for the best-case scenario (swivel angle mean RMSE 5 deg.). For all the movement tasks tested, the method in [20] presents a mean RMSE of 10 deg. for the swivel angle estimations. Given that global errors of the swivel angle double those of the best-case scenario, a report of global errors of GH joint angle estimations of the method in [20] is required to reach a conclusion regarding its suitability for clinical use.

Table 7 summarizes the comparison of our contributions w.r.t. comparable works (i.e., [20]).

Table 7

Contributions of this paper w.r.t. comparable works.

Work	Method	Method evaluation	Accuracy of GH joint angles
[20]	IK-based swivel angle estimation	(1) Studied angles: swivel angle plus the shoulder, elbow, and wrist joint angles(2) Reference angles: obtained from custom-made inertial MOCAP; homologation-calibration of the readings is not reported(3) Movements: compound movements(4) Sensitivity analysis: no	Mean RMSE: 4.8 deg. (best-case scenario)
This paper	Hybrid exoskeleton-optical MOCAP	(1) Studied angles: shoulder angles(2) Reference angles: simulated(3) Movements: 1-DOF and multi-DOF shoulder movements(4) Sensitivity analysis: method accuracy w.r.t. marker position errors produced by marker drift or calibration errors	(a) Mean RMSE: 0.9 deg. (assuming no marker drift or calibration errors) (b) Mean RMSE: 3.6 deg. (with marker drift or calibration errors up to 20 mm)

5. Conclusions and Future Work

In the context of RAR, this paper presents the formulation, implementation, and assessment, in silico, of a novel accurate method to estimate the patient GH joint angles during therapy. Our method does not require redundant markers or cameras and relies on simple geometric relationships and tools of standard robotics and computer vision libraries. These characteristics make it economical and readily applicable in RAR.

The accuracy and the robustness of our method are evaluated using computer-generated human movement data corresponding to actual movement datasets of the Armeo Spring. We present a formal sensitivity analysis of the pose estimation accuracy w.r.t. marker position estimation errors produced by (a) system calibration errors and (b) marker drifts (due to skin artifacts). This analysis indicates that even in the presence of large marker position errors our method presents an accuracy that is acceptable for patient mobility appraisal.

Future work includes (a) implementation of the method using commercially available RGB-D vision sensors, (b) evaluation of the method accuracy with actual human movement data, (c) adaptation of the method using solely RGB cameras, and (d) extension of our method to address other limbs.