Soft Tissue Artifact Compensation Using Triangular Cosserat Point Elements (TCPEs)

. Existing methods which compensate for the Soft Tissue Artifact (STA) in optoelectronic motion measurements estimate the rigid motion of a nearly rigid underlying body segment based on analysis of the motion of all fiducial markers. The objective of the proposed Triangular Cosserat Point Elements (TCPE) method is to estimate the motion of the underlying body segment even when the STA in the entire cluster of markers can be large. This is accomplished by characterizing the cluster of markers with TCPEs defined by triangles based on all combinations of three markers. Then, scalar deformation measures characterizing the magnitudes of strain and relative rotation of pairs of TCPEs are defined for each TCPE. These deformation measures are used to define a filtered group of TCPEs which best represents the motion of the underlying body segment. The method was tested using an experimental setup that consists of a rigid pendulum with a deformable 300ml silicone breast implant attached to it as a simulation of the soft tissue around a bony segment. The rotation angles extracted from markers on the deformable implant were compared with simultaneous measurements of the rigid pendulum using an optoelectronic system. Analysis of the experimental data shows that this filtering process substantially reduces the error due to the STA even though the data set includes large deformations. In particular, the analysis shows that the error reduction using the TCPE approach is larger than the reductions obtained using standard least-squares minimization methods.

The most widespread technique used for the reconstruction of human or animal skeletal system 3D kinematics is optoelectronic streophotogrammetry.This method involves placing markers on the skin of the analysed segment and capturing their locations in space as a function of time.Three non-collinear markers are adequate to determine the 3D pose of a rigid bony segment.If n markers   3 n  on a segment move as a rigid body, then the position vector i x of each marker i at any time t can be mapped from its position i X at an arbitrary reference time by the transformation ii  x t QX (1) where Q is a proper orthogonal second order rotation tensor and t is a translation vector, both function of time t only.However, since the markers are attached to the soft tissues their motion relative to the underlying bony segment is that of a deformable body and equation (1) cannot be satisfied exactly.This deformation causes a limitation on this technique that is referred to as the Soft Tissue Artifact (STA).Since body segments are not rigid bodies, and the STA at joint landmarks where markers are usually placed is even more prominent than at other locations on the body segments (Cappozzo et al., 1996), many methods employ redundant skin markers and propose different techniques to reduce the effect of STA on the prediction of the underlying motion of the bony segment.One method to compensate for the STA is based on a kinematic assumption and the least-squares minimization technique (Ball and Pierrynowski, 1998;Cheze et al., 1995;Soderkkvist and Wedin, 1993;Spoor and Veldpaus, 1980;Veldpaus et al., 1988).Another method assigns weights to each marker in the cluster of markers to define an effective inertia tensor and then adjusts these weights to minimize differences between the eigenvalues of the inertia tensor at each time step to estimate rigid body motion (Andriacchi et al., 1998).The relationship between this method and another SVD-based least-squares method has been discussed by (Cereatti et al., 2006).A more general approximation i x of the vectors i x at time t assumes homogeneous deformation, such that ii  x t FX (2) where F is a non-singular second order tensor with a positive determinant which is the deformation gradient of the homogeneous deformation and t is a translation vector, both function of time t only.Since the measured vectors i x include measurement errors and the STA the transformation (2) is also only an approximation of the actual motion.If the number of markers in the cluster is 4 n  , then the best-fit values of F and t can be found by solving the unweighted least-squares minimization problem characterized by (3) as was shown previously by (Sommer III et al., 1982) for the purpose of osteometric scaling of skeletal segments.Applying this method to STA compensation, the rigid rotation tensor can be extracted from F using the polar decomposition (Dumas and Cheze, 2009) or by defining alternative rotation tensors based on analysis of deformations associated with F (Ball and Pierrynowski, 1998).In numerous other works, F was required to be an orthogonal rotation tensor, with the mapping (2) corresponding to rigid body motion.In this case, the same minimization problem (3) was solved by requiring F to satisfy the constraint
(Wolf et al., 2010) described an alternative approach in which the cluster of markers was characterized by a group of tetrahedrons based on all combinations of four markers.The kinematics were analysed using the nonlinear continuum theory of a Cosserat point (Rubin, 1985) by treating each tetrahedron as a Cosserat Point Element (CPE).The magnitude of strain in each CPE at each time frame was then used to define a filtered group of CPEs which was used to estimate the underlying rigid-body motion.The method was validated by experimental data.
In this paper, several modifications and improvements are presented.Here, the tetrahedral CPEs are replaced by triangular CPEs (TCPEs) which reduce potential problems with poor aspect ratios of relatively flat tetrahedrons.Furthermore, in addition to the magnitude of strain in each TCPE, the relative rotation between pairs of TCPEs is used as a measure of deformation of the body segment.
In both the homogeneous deformation and the rigid body methods the entire cluster of markers on a body segment is analysed, including markers attached to potentially highly deformed parts of the segment or markers with large measurement errors.This tends to corrupt the resulting estimate of rigid body motion.In contrast, in the proposed TCPE method the deformation measures are used to define a filtered group of TCPEs which best represents the motion of the underlying body segment.Moreover, the objective of the proposed TCPE method is to estimate the motion of the underlying body segment even when the STA in the entire cluster of markers is large.
Recently (Průša et al., 2013) have analyzed errors in the estimation of the local value of the deformation gradient from the motion of markers on a deformable body.They developed an upper bound for this error.In particular, their results indicate that this local error can be surprisingly large for a general inhomogeneous deformation field.The work here focuses on estimating the rigid body motion of a body part by using an average homogeneous deformation based on a quality group of TCPEs.Since the spatially averaged value of the deformation gradient is relatively robust to spatial variations, it is expected that the quality group of TCPEs can be used to obtain an accurate estimation of rigid body motion.
An experimental setup using a silicone breast implant as a non-rigid body attached to a rigid pendulum has been constructed for validation of the new method.The effectiveness of the TCPE method is tested by comparing the rotation of the pendulum predicted by the kinematic measurements from markers placed on the implant with direct measurements of the pendulum motion.Analysis of the experimental data shows that the filtering process in the TCPE method substantially reduces the error due to the STA even though the data set includes large magnitudes of strain and relative rotations.The experimental results are also compared with predictions of the homogeneous deformation and the rigid body least-squares methods.

Materials and Methods
The proposed method is based on a cluster of markers randomly distributed on the surface of a deformable segment.The cluster of n markers is characterized by a group of TCPEs based on all combinations of three markers.For markers placed on a rigid body the kinematics of the TCPEs should have zero strains and zero relative rotations (i.e. they should all rotate with the same rotation tensor), apart from inaccuracies in the measurement system.For markers placed on a deformable body the TCPEs can experience large strains as well as large relative rotations.In the following sub-sections the analysis of the non-rigid kinematics of a TCPE will be used to determine two scalar deformation measures; the magnitude of the strain and the relative rotations between pairs of TCPEs.Most importantly, the method uses these deformation measures to filter the data to define groups of those TCPEs that are most likely to represent the underlying rigid body motion.

Kinematics of a TCPE
The present configuration of a TCPE at time t is characterized by the position vectors 0 1 2 } , {, x x x of the three nodes (vertices) of the TCPE relative to a fixed origin as shown in Figure 1.The values of the positions 0 1 2 } , {, x x x in the reference configuration are denoted by the vectors where d is a scalar that denotes the square root of the determinant of the metric defined by ij dd.Note that 3 d is a unit vector that is perpendicular to the plane of the triangle defined by 1 d and 2 d .Moreover, the reference reciprocal director vectors are defined by where i j  is the Kronecker delta symbol.Furthermore, the deformation gradient F of the TCPE is defined as a second order tensor using the tensor product (outer product) operator  in the The deformation gradient F is a nonsingular second order tensor so that the polar decomposition theorem (Malvern, 1969) can be used to determine the unique proper orthogonal second order rotation tensor R and the unique positive definite symmetric second order stretch tensor M, such that where I is the second order unity tensor.Moreover, the associated symmetric right Cauchy-Green deformation C and the Lagrangian strain E are second order tensors defined by: The rotation tensor R can then be used to obtain the rotation angle  and the associated rotation axis (a unit vector) e by (11) (Rubin, 2007) where ij R are the components of the rotation tensor R relative to the fixed rectangular Cartesian base vectors i e .
Often in biomechanics it is of interest to determine the angle of rotation about a user specified joint axis parallel to a unit vector 3 a .If 3 a is not an eigenvector of R then this angle cannot be defined uniquely.However, an approximate average angle of rotation 3 ()  a about the axis 3 a can be defined by the following procedure.Let i a be an orthonormal set of base vectors, such that and define the auxiliary vectors i b and i b by  a of the relative rotation about the axis 3 a by the expression

Strain and rotational measures of deformation
For ideal rigid body motion all TCPEs have the same rotation tensor R and their strains E vanish.However, when markers are placed on a deformable body, it is convenient at each time step to define two scalar deformation measures I E and that characterize the magnitude of the strain in a TCPE and a measure of the rotation of the TCPE relative to all other TCPEs.Specifically, the magnitude of the strain of the I th TCPE is defined by ( ) and the relative rotation angle between the I th and J th TCPEs is defined by  An average rotation tensor avg R is determined for each of these groups using the method described in (Sharf et al., 2010) based on the minimization of the Riemannian bi-invariant metric for rotation matrices in terms of the Frobenius norm such that the corresponding mean 12 ( , ,..., ) (3), arg min log( ) Once an average rotation tensor avg R is obtained, the rotation angle 3 ()  a about a specified axis 3 a can be determined using equations.( 12)-( 13), with R replaced by avg R .Also, it is noted that other simpler definitions for the average rotation tensor can be used, as described in (Sharf et al., 2010).

Least-Squares Methods for Comparison
Two least-squares minimization techniques for the determination of the average rotation from the position vectors of a cluster of markers on a deformable body segment are presented in this sub-section.The first technique is obtained by forcing the transformation to conform to a rigid body motion and the second one is obtained by allowing homogeneous deformation of the cluster of markers.
From the definitions of ,, xxX in the introduction, it is convenient to define the following vectors It then follows that equation ( 2) can be re-written in the form Using the assumption of homogeneous deformation, it can be shown that the best-fit F that satisfies (3) can be written as (Sommer III et al., 1982).The corresponding unique rotation tensor R and Lagrangian strain tensor E associated with F can be obtained using equations ( 8)-( 10).Moreover, the magnitude of the strain in the entire cluster of markers can be obtained from E using equation ( 14).
Furthermore, using the assumption of rigid body, it can be shown that the best-fit that satisfies (3) with constraint (4) can be calculated from the polar decomposition of F as follows   (23) (Spoor and Veldpaus, 1980;Veldpaus et al., 1988).It is important to emphasize that the expressions ( 22) and ( 23) yield different rotation tensors for general motions.Then, the rotation angles 3 ()  a about a specified axis 3 a can be determined using equations ( 12)-( 13), with R replaced by R or R , respectively.

Experimental setup
An experimental setup was constructed in order to test the proposed method.The setup was composed of a rigid pendulum with a 300cc silicon breast implant attached to its end (see Figure 3).Seven markers (markers 5-11) were randomly placed on the deformable silicone implant, which create 35 different triangles.Four additional markers were placed on the rigid pendulum to obtain an accurate direct measurement of the rigid motion of the pendulum; two of them (markers 1-2) were used to determine the pendulum's axis of rotation and two (markers 3-4) were used to determine the pendulum's motion about its axis of rotation.All the markers were 14mm diameter reflective markers and their positions as functions of time were extracted using the 1.8.1 version of Vicon Nexus software (Vicon, Oxford, UK) connected to a Vicon MX13 8-camera motion capture system with a frequency of 120 Hz.A series of experiments were conducted with different initial amplitudes of rotation of the pendulum's rod.The pendulum was allowed a free passive movement on one side of the motion, and on the other side it was stopped and pushed manually before reaching maximum height in order to induce larger accelerations and deformations of the implant.Although the system does not replicate the mechanical properties of the human body, it can be used to examine the effectiveness and accuracy of the proposed method to compensate for the STA.
The rotations angles i  about the pendulum axis were determined using the rotation tensors of each TCPE.Then, following the filtering procedure described in sub-section 2.3, four different average rotation tensors were calculated at each time step using equation ( 19) about the axis of the pendulum were calculated using equations ( 12)-( 13).In addition, the rotation angle All  was calculated using 11 the average rotation tensor from the group of all possible TCPEs (35 in our case) without any filtering.Furthermore, the rotation angles {,} HD RB   were obtained using the least-squares minimization procedures described in sub-section 2.4 for the assumptions of homogeneous deformation and rigid body motion, respectively.All the angles above were then compared with the "true" pendulum angle true  which was calculated using the markers attached to the rigid pendulum.
The objective of this analysis was to examine the effect of filtering based on the two different deformation measures on predictions of the rotation angles of the pendulum from data for a highly deformable body using the average rotation tensors.Specifically, the predictions of the filtered group of TCPEs are compared with those based on the least-squares minimization techniques that analyse the entire cluster of markers.

Results
The results for the experiment described in Section 3 were recorded for a total of 25 cycles of the pendulum in five different trials.The predictions of the rotation angles for three cycles of one trial of are plotted in Figure 4a.The lower part of the graph (negative angles) corresponds to the side of the motion where the pendulum was allowed a free passive movement, and the upper part (positive angles) corresponds to the side where the pendulum was stopped and pushed manually before reaching maximum height.
Each TCPE represents a slightly different rotation, and the difference between them is especially large at the top of the graph where manual perturbations were induced to increase the accelerations and deformations of the implant.The same data is plotted in Figure 4b where t N is the number of time steps.The results are recorded in Table 1.The left column in Table 1 presents the range of errors recorded in the entire simulation; the middle column presents the RMSE for the entire motion, while the right column presents results for only the highly deformable regions of motion where manual perturbations were introduced to increase the magnitude of the deformation.Table 1 also   the implant; all 35 individual TCPEs (solid grey lines) , the group of all TCPEs (dashed green line), the group of filtered TCPEs (dashed black line), the group of all markers using the homogeneous deformation least-squares method (dashed cyan line) and the group of all markers using the rigid body least-squares method (dashed magenta line).These predictions are compared with the true angle (solid red line).(a) For three cycles of the pendulum and (b) for a region of motion with large deformations.

Discussion
This paper presents a new method for compensation of the Soft Tissue Artifact (STA) using Triangular Cosserat Point Elements (TCPEs) based on all combinations of three markers from a cluster of markers randomly placed on the surface of a body segment.Two scalar deformation measures of the TCPE have been defined to characterize the strain of each TCPE and the relative rotation of pairs of TCPEs.Standard least-squares minimization methods use the entire cluster of markers to estimate the rigid motions of an underlying body segment based on assumptions of either homogeneous deformation or rigid body motion.In contrast, the TCPE approach filters the data to determine quality groups of TCPEs that best represent the underlying rigid body motion of the body segment.The experimental results discussed in the last section demonstrate that this filtering process substantially reduces the error due to the STA despite the large magnitudes of strain and relative rotations associated with the motion of the entire cluster of markers.Furthermore, the results were analysed to reveal that filtering based on both measures of deformations produces more accurate results than those predicted by filtering based on only one of the deformation measures, with filtering based on only the strain measure providing higher quality TCPEs than filtering based on only the relative rotation measure.All three filtered groups of TCPEs provided more accurate predictions than obtained by either of the two least-squares methods, which showed similar predictions.
The accuracy of the proposed TCPE method most likely will be compromised if the number of TCPEs is too small due to a small number of visible markers.In this case, the filtering process may require the use of poor quality TCPEs for predicting the underlying rigid body motion.However, it is likely that the accuracy of techniques using the entire cluster of markers will also be compromised with a small number of visible markers.
The TCPE method provides local information about strains and rotations that can be used for clusters of markers with large deformations and relative rotations.Consequently, it can be used to quantify the motion of skin and muscles relative to the underlying body segment and it can be applied to clusters of markers placed on obese patients.Moreover, since the deformation measures of strain and relative angles are physically different, additional work is needed to optimize the filtering process.

Figure 1 :
Figure 1: Positions of a TCPE in its reference (t0) and present (t) configurations.
are unit vectors that are parallel to the projections of   12 , bb in the 12  aa plane which is the plane normal to 3 a (see Figure2).In general, the angles between   are not equal.However, it is possible to define the average angle 3 () strain tensor and rotation tensor of the I th TCPE, respectively.Note that / Δ IJ  is limited to the range   0,(Rubin, 2007).By definition / TCPEs have the same rotation tensor.Then, the scalar measure of relative rotation total number of TCPEs defined by the cluster of markers.It is assumed that the rotation tensor of a TCPE having small values of the rigid body motion of the underlying body part than that for a TCPE with large values of these measures.Therefore, these deformation measures can be used to define filtered groups of quality TCPEs that exclude TCPEs having large values of these deformation measures.

Figure 2 :
Figure 2: The definition of the unit vectors   12 , bb and the angles   12 ,  from the unit the groups of TCPEs associated with the parameters{ ,

Figure 3 :
Figure 3: The experimental setup.A rigid pendulum with a 300cc silicon breast implant attached to its end.
for a shorter time period which corresponds to a region of the motion with high deformations.The Root Mean Square Errors (RMSE) { , , simulation containing 25 cycles of the pendulum are calculated using formulas of the type includes the percentage (in parentheses) of the error reduction relative to the values for All  determined by each of the methods.The maximum values of the deformation measures in the different groups provide an estimate of the non-rigidity in the motion of the cluster.Figures 5 and 6 focus on the same three cycles as in Figure 4a. Figure 5 plots the maximum value of I E as a function of time for the group of all TCPEs and the maximum value of I E for the filtered group of TCPEs associated with Combined N .For comparison, the value of HD E for the entire cluster of markers obtained using the assumption of homogeneous deformation by equation (22) is also plotted.Similarly, the maximum value of I   for the group of all TCPEs and for the filtered group of TCPEs associated with Combined N are plotted in Figure 6.The maximum value of I E of all TCPEs during the entire simulation of 25 cycles was 1.27, while the maximum value of I E for all TCPEs in the filtered group associated with Combined N was 0.09 and the maximum value of HD E 13 associated with the homogeneous deformation gradient was 0.43.The maximum value of I   of all TCPEs during the simulation was 99.75 degrees, while the maximum value of I   for all TCPEs in the filtered group associated with Combined N was 20.87.These values indicate that the implant experienced large deformations during the simulation which simulates a large STA.

Figure 4 :
Figure 4: Rotation angles of the pendulum as functions of time predicted from markers on Table 1: The Error range and Root Mean Square Errors (RMSE) { , of the pendulum.The percentages of the error reduction relative to the values for All  are presented in parentheses.

Figure 5 :
Figure 5: The maximum value of strain norm

Figure 6 :
Figure 6: The maximum value of relative rotation measure