Unphysical properties of the rotation tensor estimated by least squares optimization with specific application to biomechanics

Analysis of the transformation of one data set into another is a ubiquitous problem in many fields of science. Many works approximate the transformation of a reference cluster of n vectors X i (i=1,2,..,n) into another cluster of n vectors x i by a translation and a rotation using a least squares optimization to obtain the rotation tensor Q . The objective of this work is to prove that this rotation tensor Q exhibits unphysical dependence on the shape and orientation of the reference cluster. In contrast, when the transformation is approximated by a translation and a general non-singular tensor F , which includes deformations, then the associated rotation tensor R does not exhibit these unphysical properties. An example in biomechanics quantifies the errors of these unphysical properties.


Introduction
Analysis of the transformation of one data set into another is a ubiquitous problem in many fields of science. For examples: behavioral science analysis (Hurley and Catell, 1962;Schonemann,1966); satellite attitude estimation (Wahba, 1966); registration and motion of 3-D shapes (Laub and Shiflett, 1982;Arun et al., 1987;Besl and McKay, 1992); anthropometric scaling (Lew and Lewis,1977;Sommer et al., 1982); and biomechanical motion analysis (Spoor and Veldpaus, 1980;Veldpaus et al., 1988;Soderkvist and Wedin,1993;Challis, 1995;Cappozzo et al., 1996;Ball and Pierrynowski, 1998 In the applications discussed above the vectors x i include inhomogeneous deformations relative to X i due to a number of sources associated with measurement error and actual inhomogeneous deformations. Specifically, in biomechanical motion analysis the forces and moments on body joints can be estimated by knowing the rigid motion of bones in the body. However, piercing the skin by placing pins in the bone to determine actual bone position cannot be done for general patient diagnosis. Therefore, estimates of the bone pose using markers on the deformable skin are essential. From a continuum mechanics point of view, it is obvious that the rotation of a material line element in a deformable body depends on the deformation field and on the specific orientation of the line element in the reference configuration. If the deformations are not too large then it is reasonable to use a rigid body approximation. Often (e.g., Schonemann, 1966;Wahba, 1966;Spoor and Veldpaus, 1980;Arun et al., 1987;Veldpaus et al., 1988;Besl and McKay, 1992;Soderkvist and Wedin, 1993;Challis, 1995) the transformation of X i into x i is approximated as a translation and rotation using least squares optimization to determine the rotation tensor Q. The main objective of this work is to prove that this rotation tensor Q exhibits an unphysical dependence on the orientation and shape of the reference cluster X i . In contrast, when the transformation between these data sets is approximated by a translation and a general non-singular tensor F, which includes deformations, then the associated rotation tensor R is uninfluenced by shape and orientation changes of the reference cluster. For biomechanical motion analysis this means that the estimates of the underlying bone pose using Q will include errors due the STA as well as additional unphysical errors which depend on the placement of the markers. These additional unphysical errors can be removed using the analysis based on F.

The affine approximation
Within a general context, the objective is to determine a simple approximate relationship between the reference cluster of vectors X i and another cluster of vectors x i .
To this end, it is convenient to define the centroids {X, x} of {X i , x i } by the expressions and define the difference vectors, X i and x i , such that Then, the estimates x i * of x i based on an affine transformation of X i are defined by where t is the approximate translation vector of X and F is a second order non-singular transformation tensor.
Using a least squares procedure with an affine approximation (2.3), define the function of the sum of squared errors (e.g. Plackett, 1960) (2.4) where (  ) denotes the inner product between the vectors. If {X i , x i } are vectors of dimension m, then t has dimension m and F has dimension mm. In continuum mechanics the vectors are in 3-space with m=3, X i represent material line elements in the reference configuration, x i represent material line elements in the present configuration and F is called the deformation gradient. In the following discussion use will be made of the terms translation, deformation, rotation and stretch from continuum mechanics even though other names for the same mathematical quantities are used in other fields.
Next, taking the variation f of f(t,F) with respect to the variations {t, F} of {t, F} yields where ab denotes the tensor (outer) product of the vectors {a,b}. Since t and F are independent, critical values of f are determined by the condition that the coefficients of {t, F} vanish, which yields where ( T ) denotes the transpose operator and it has been assumed that the tensor H is nonsingular. It follows that t is the translation of the centroids and it is noted that F includes both rotation and stretching of X i since the orientation and length of x i can be different from those of X i .
Using a generalization of the polar decomposition theorem (Malvern, 1969) and assuming that F has a positive determinant, F can be represented in the form where R is a unique proper orthogonal rotation tensor, M is a unique positive definite, symmetric stretch tensor and I is the identity tensor.

The rotation tensor based on the rigid body approximation
For the rotation tensor based on the rigid body approximation, F in (2.3) is restricted to be a proper orthogonal rotation tensor Q, which satisfies the conditions where  is a skew-symmetric tensor.
Since  is skew-symmetric and H is symmetric, critical values of f require Now, assuming that -F has a positive determinant, use can be made of the polar decomposition theorem to express -F in the form where -R is a unique proper orthogonal rotation tensor and -M is a unique positive definite stretch tensor. It then follows that the unique solution of (3.3b) requires

The rotation tensor based on the affine approximation
In Section 3 it was shown that when the least squares approximation is used together with the constraint that F be a rotation tensor Q, the optimal value of Q is given by the rotation tensor -R associated with the polar decomposition of the auxiliary tensor -F defined in (2.7). An alternative specification of Q can be obtained using the least squares optimization based on the affine approximation to obtain the tensor F defined in (2.7).
Then, the value of Q can be specified by the rotation tensor R in the polar decomposition Next, using (2.7), (2.8) and (3.4) it follows that which proves that the rotation tensor R is different from -R whenever -MH -1 is not symmetric. This means that the solutions (3.5) and (4.1) are different, with the estimated rotation being influenced by deformations due to F.
In order to quantify the influence of deformations on the estimate of the rotation tensor, it is noted that a general line element (i.e. vector) in the reference configuration with unit direction S has stretch  and is rotated to the unit direction s in the present configuration, where A  B = tr(AB T ) is the inner product of two second order tensors {A, B}. It then follows that the angle  between MS and S is given by Since S is a unit vector its variation S can be expressed in the form where  is an arbitrary skew-symmetric tensor. Now, for fixed M the variation of (4.4) which vanishes for all  provided that the skew-symmetric tensor A vanishes In order to solve this equation for S it is convenient to write M in its spectral form where  i are the positive ordered eigenvalues, P i are the associated orthonormal eigenvectors of M and  ij is the generalized Kronecker delta. If S is equal to one of the principle directions P i then A vanishes and  = 0, which means that R is the exact rotation of the principle directions of M. Furthermore, it can be shown that the maximum value so that (4.7) requires Thus, the values of  and  max are given by ] . (4.11) It can easily be seen that this expression for  max gives the maximum value of  since it is based on the minimum value of the ratio  i / 1 . Moreover, this value of  max determines the maximum error in the rotation of material line elements relative to that estimated by the rotation tensor R. Figure 1

Unphysical properties of the rotation tensor
The objective of this section is to discuss some unphysical properties of the rotation tensor Q = -R determined in Section 3 using the least squares approximation with the constraint that F be a rotation tensor. To this end, consider another set of vectors X i which are defined by applying a pure rotation  to X i  which is not valid for arbitrary rotation tensors  so R does not equal -R in general. This means that the estimation of the rotation tensor (3.5) exhibits an unphysical dependence on the orientation of the vectors X i through the rotation tensor . In the next section an example will be considered which shows that this rotation tensor also exhibits an unphysical dependence on the shape of the vectors X i through the tensor H.
In contrast, it follows from (2.7) and (5.3) that the tensors -FH -1 and FH -1 yield the correct affine transformation tensor F F = -FH -1 = FH -1 , (5.9) so the alternative rotation tensor Q proposed in (4.1), which equals the rotation tensor R of F, is unaffected by shape and orientation changes of the vectors X i .

A quantitative example
Consider the simple planar case with a triangular cluster of three points defined relative to the fixed rectangular Cartesian base vectors {e 1 , e 2 } by where L is a length measure and  controls the shape of the triangle (see Fig. 2). Next, consider the two-dimensional rotation tensor  defined by the angle    = (cos e 1 + sin e 2 )e 1 + (-sin e 1 + cos e 2 )e 2 , (6.2) with the rotated triangle in Fig. 2 defined by the points Moreover, the rotated triangle is attached to a body which is distorted (deformed with no area change) by the homogeneous two-dimensional deformation gradient F  FM F = M =  (e 1 e 1 ) + 1  (e 2 e 2 ) , R = I , (6.4) where  is the stretch in the e 1 direction and I is the two-dimensional identity tensor. Using

Discussion
The analysis in the previous sections shows that when the transformation of n vectors X i into another set of n vectors x i is approximated by a pure rotation, the least squares estimate of the rotation tensor Q exhibits unphysical dependence on the shape and orientation of the reference data X i . In contrast, when the transformation is approximated by a general non-singular tensor F, which includes deformation, the rotation tensor remains unaffected by shape and orientation changes of X i . Specifically, the average transformation tensor F is defined by where it is assumed that H is non-singular. Further, assuming that the determinant of F is positive, the polar decomposition theorem (2.8) is used to determine the unique proper orthogonal rotation tensor R and the unique symmetric positive definite stretch tensor M. can also be used to analyze the inhomogeneity in a deformation field which can reveal relative motion between different segments of the body. More specifically, the TCPE approach adds a vector normal to the triangle in order define a non-singular F for full three-  However, they did not recognize the stronger implications of using an average deformation gradient F to obtain an average rotation tensor R which is unaffected by shape and orientation changes of points in the reference configuration.
The main new aspect of the work in this paper is to prove that when the transformation of the data set X i into x i is approximated by a translation and rotation, then the rotation tensor Q obtained by least squares optimization exhibits an unphysical dependence on the orientation and shape of the reference cluster X i , as shown in the example in Section 6.
Also, the value  max in (4.11) is a new measure which gives an upper bound for the error in estimating the orientation of the underlying bone pose in a bone-tissue structure experiencing homogeneous deformations.