There are two main problems in tracking and navigation. The first involves using data from targets rigidly attached to a body to determine the orientation and location of the body relative to a reference attitude and location. In many applications, the targets are optical and the data are recorded with a computer vision system. For example, in a variety of biomechanics applications, sets of optical targets are placed on two anatomical segments and the relative displacements of the targets are then used to measure the relative rotation and translation of the two segments. It has long been known that the measurement of the rotation is difficult and prone to error. Several studies in the biomechanics community have discussed and quantified these errors for a variety of situations (e.g., see [1, 2, 3]). We devote this page to discussing various methods used to estimate and from sets of experimental data. These estimates are denoted by an asterisk, i.e., and . This problem is sometimes known as the Wahba problem [4] or the Procrustes problem [5].

Contents

## Measurements and a cost function

Suppose we have a series of current position measurement vectors associated with vectors of reference position measurements, . As shown in Figure 1, the reference and current measurements typically pertain to a set of targets, or landmarks, that are rigidly attached to the bodies of interest, and the motion of these objects is detected by cameras. The targets often come in fours so that they can be used to form a triad of vectors.

Typically, the position measurements are all taken with respect to the same frame, which is often known as the laboratory frame:

We wish to determine the rotation tensor and translation that solves the equations

Because of measurement errors, no unique solution to these equations exist. The best we can do is to find the solution that minimizes the cost function

where are weights assigned to the individual measurements. Because is a rotation tensor, and thus proper-orthogonal, the cost function (3) is subject to the constraints

An alternative decomposition of the motion given by (2) expresses this transformation in terms of the parameters of a helical (or screw) motion. As denoted in Figure 1, the four parameters for this motion are the rotation axis that is parallel to the screw axis; the angle of rotation, , about this axis; the position vector of a point on the screw axis, ; and the translation along this axis. It is straightforward to show that and can be determined by solving the equation

However, the tensor is not invertible, so the solution for is not uniquely defined and indeed can be indeterminate when . One possibility is to choose to be the position vector from the origin to the screw axis that is normal to . Then, from (5), we find that

where

This solution to the parameters of the screw motion is equivalent to the one chosen by Spoor and Veldpaus in [6].

## Solution methods

We now discuss four well-known solutions to the problem of estimating the rotation tensor and translation. We start with a naive approach that assumes there is no noise, measurement error, or finite precision issues. Next, a method from the satellite dynamics community known as the TRIAD method is discussed [7]. We then present a popular algorithm that uses the singular value decomposition of a matrix [8, 9]. We close with a discussion of the -method, which employs unit quaternions (i.e., Euler-Rodrigues symmetric parameters) and features a novel eigenvalue problem [10, 11]. Comparisons of the methods we present can be found in works by Eggert et al. [12] and Metzger et al. [13]. We do not discuss error estimates for any of these methods. Instead, we refer the interested reader to Dorst’s paper [14]; his work enables one to correlate errors in the estimates and to noise and target distributions. Additional works in this area with specific application to biomechanics include studies by Panjabi [1], Spoor and Veldpaus [6], and Woltring et al. [2], among others.

### A naive method

A naive method is to ignore errors in the measurement data and assume that one has four targets. First, we define the following augmented arrays for each target’s reference and current position measurement vectors:

Assuming no errors in either the rotation or translation estimates, one can rewrite each of the four motion equations given by (2) in matrix-vector form as

where

Given four non-coplanar targets, (9) can be used to solve for and :

Because there will always be some amount of noise and measurement error, the matrix in (11) will typically not be proper-orthogonal. This fact contributes substantially to the errors observed when using this naive method.

### The TRIAD method

The TRIAD algorithm is credited to Shuster and Oh [7] and dates to the early 1980s. It is often used with magnetometer, gravity, and star-pointing data, but here we consider the case where it is used with data obtained from optical targets. Given three sets of reference and current target position data, say, and , we can use these sets of vectors to define the right-handed orthonormal bases and , respectively:

and

These unit vectors can be used to define two proper-orthogonal matrices:

With the help of the motion equations (2), we can manipulate and to solve for the matrix :

The translation vector is then calculated using the mean positions of the targets before and after the motion:

where

For a set of four targets, there are 24 possible constructions of the matrices and , and thus 24 possible estimates of and . To obtain the optimal and , we examine the error between the true and calculated positions of the targets after the motion:

The rotation tensor and translation that yield the smallest error are considered the optimal results .

### A method based on the singular value decomposition

The method featuring the singular value decomposition (SVD) is credited to Arun et al. [8] and to Hanson and Norris [9] in the 1980s. Additional discussions of this method can be found in [15, 16, 17]. The first step of the SVD method involves computing the mean positions for a set of four targets before and after the motion:

Next, a matrix is formed using position vectors relative to these mean positions:

where

We then calculate the singular value decomposition of :

where and are orthogonal matrices, and is a 3 3 diagonal matrix. Finally, the optimal estimates of the rotation and translation are given by

### The q-method

The -method dates to the latter half of the 20th century. The method uses unit quaternions (i.e., Euler-Rodrigues symmetric parameters) to parameterize the rotation and features a novel eigenvalue problem to arrive at an estimate for the optimal rotation. This solution was developed independently by Davenport [10] and later in a well-cited paper by Horn [11]. To show how the -method works, we start by defining the parameter and the “centers” of the reference and current position measurements based on assigned weights :

Expanding the cost function given in (3), we find with some straightforward manipulations that

We are interested in extremizing as a function of the rotation tensor and translation , so we can ignore the underbraced terms in (25) that are independent of and . In short, we wish to find and that maximize

We now parameterize using a set of Euler-Rodrigues symmetric parameters, :

This representation enables us to express (26) in a more compact form:

where the 4 4 symmetric matrix and the column vector are given by

Here, is the 3 3 identity matrix, the column vector

and the 3 3 matrix , for which the components

Note that the matrices and are composed entirely of measurements. Our definition of the column vector differs by a minus sign from that found in the literature for the case when (e.g., see [7, 18]). This difference is expected and can be traced to our convention for the skew-symmetric part of the matrix in (27).

We are now in a position to solve for the optimal and . Specifically, we seek and that maximize (28), subject to the Euler parameter constraint . With the help of a Lagrange multiplier , one seeks and that maximize

It is now straightforward to calculate the optimal solution by differentiating with respect to its arguments. Doing so, one finds

That is, the optimal is an eigenvector of . For to be maximized, must be the eigenvector corresponding to the largest eigenvalue of . Once this eigenvector is found, can be constructed using (27). Notice that is a simple mean value, but it can only be obtained after has been determined. You should also notice that because the 4 4 matrix is real and symmetric, it has four real eigenvalues and its eigenvectors are orthogonal. One advantage of the -method is that determining the rotation axis and the corresponding rotation angle for the screw axis representation of the motion follows easily once is known:

## An example of estimating a rigid-body motion

To illustrate the four methods for estimating rotations and translations from measurement data, we consider the example of a tossed book. Referring to Figure 2, we choose four landmark points on a uniform rectangular prism of length , width , and thickness such that their position vectors in the reference state are

For the example at hand, we choose the dimensions , , and (with units of inches) to mimic a small book. Note that in the reference state, the body’s center of mass is located at the origin and the rotation tensor can be taken as the identity transformation, i.e., in the reference configuration.

Suppose we rotate and translate the book such that at time , the book’s orientation and translation are described by, respectively,

where, using a 3-2-1 set of Euler angles to parameterize , the rotation axis . A lengthy but straightforward calculation shows that the components of at time are

In terms of Euler’s representation, the rotation is equivalent to a counterclockwise rotation of about an axis

We assume there is no noise in the measurements, in which case, using (2), the landmark points’ position vectors at time are

We now examine how the four methods described earlier estimate and at time using the position measurements (35) and (39). The forthcoming results that we present were calculated in MATLAB; the associated code is available for download here.

### Solution using the naive method

The solution procedure for the naive method reveals that the augmented matrix in (11) is

and so

Within numerical precision, the resulting predictions of and are accurate.

### Solution using the TRIAD method

Out of 24 possible permutations, the matrix that has the least error, as defined by (18), is given by

The resulting is identical to that estimated by the naive method; the same is true of the translation.

### Solution using the method based on the singular value decomposition

The matrix in (20) is calculated as

To compute the singular value decomposition of in MATLAB, we use the built-in singular value decomposition function , where , , and . In this case,

Both and have unit determinant, and so, from (23)_{1},

which agrees with the estimates from the naive and TRIAD methods. A subsequent calculation involving (23)_{2} shows that the translation estimate also agrees with the two earlier methods.

### Solution using the q-method

Using (29), (30), and (31) with unit weights , the matrices and are given by

The eigenvalues of are all real. The eigenvector that corresponds to the largest eigenvalue is

When we construct from the components of using (27), we find that the estimate of the rotation is the same as those computed by the naive, TRIAD, and singular value decomposition methods; the same is true of the translation when we calculate it from (33)_{2} after obtaining . This equivalence can be attributed to the noiseless nature of the measurements. Lastly, when the rotation has the representation , it can be shown from using (34) that the angle and axis of rotation are, respectively,

which agree with their expected values from earlier.

### The effect of measurement noise

In the absence of noise, it is well known that the four methods give equally accurate results, and this was observed in the previous calculations. However, in the presence of noise, this situation changes. To illustrate this point, we numerically added noise to the landmark points’ current position measurements and recomputed the rotation and translation estimates. We found that the naive and TRIAD methods were susceptible to noise, while the SVD method and the -method were equally accurate. Figure 3 depicts the difference between the book’s actual orientation and location at time and its configuration estimated by the TRIAD method when are subject to noise.

## Downloads

The MATLAB files used to perform calculations and generate Figure 3 for the example we discussed in this page are available here and here. The former file is the main script that executes the calculations and plots the actual and estimated configurations; the latter is a function called by the main script to draw the book in these plots.

## References

- Panjabi, M. M., Centers and angles of rotation of body joints: A study of errors and optimization,
*Journal of Biomechanics***12**(12) 911-920 (1979). - Woltring, H. J., Huiskes, R., de Lange, A., and Veldpaus, F. E., Finite centroid and helical axis estimation from noisy landmark measurements in the study of human joint kinematics,
*Journal of Biomechanics***18**(5) 379–389 (1985). - Yuan, X., Ryd, L., and Blankevoort, L., Error propagation for relative motion determined from marker positions,
*Journal of Biomechanics***30**(9) 989-992 (1997). - Wahba, G., Problem 65-1, A least squares estimate of satellite attitude,
*SIAM Review***8**(3) 384–386 (1966). - Schönemann, P. H., A generalized solution of the orthogonal Procrustes problem,
*Psychometrika***31**(1) 1–10 (1966). - Spoor, C. W., and Veldpaus, F. E., Rigid body motion calculated from spatial co-ordinates of markers,
*Journal of Biomechanics***13**(4) 391-393 (1980). - Shuster, M. D., and Oh, S. D., Three-axis attitude determination from vector observations,
*Journal of Guidance, Control, and Dynamics***4**(1) 70–77 (1981). - Arun, K. S., Huang, T. S., and Blostein, S. D., Least-squares fitting of two 3-D point sets,
*IEEE Transactions on Pattern Analysis and Machine Intelligence***9**(5) 698-700 (1987). - Hanson, R. J., and Norris, M. J., Analysis of measurements based on the singular value decomposition,
*SIAM Journal on Scientific and Statistical Computing***2**(3) 363-373 (1981). - Keat, J.,
*Analysis of Least-Squares Attitude Determination Routine DOAOP*, CSC Report CSC/TM-77/6034 (February 1977). - Horn, B. K. P., Closed-form solution of absolute orientation using unit quaternions,
*Journal of the Optical Society of America A***4**(4) 629–642 (1987). - Eggert, D. W., Lorusso, A., and Fisher, R. B., Estimating 3-D rigid body transformations: A comparison of four major algorithms,
*Machine Vision and Applications***9**(5–6) 272–290 (1997). - Metzger, M. F., Faruk Senan, N. A., O’Reilly, O. M., and Lotz, J. C., Minimizing errors associated with calculating the helical axis for spinal motions,
*Journal of Biomechanics***43**(14) 2822-2829 (2010). - Dorst, L., First order error propagation of the Procrustes method for 3D attitude estimation,
*IEEE Transactions on Pattern Analysis & Machine Intelligence***27**(2) 221–229 (2005). - Markley, F. L., Attitude determination using vector observations and the singular value decomposition,
*Journal of the Astronautical Sciences***36**(3) 245–258 (1988). - Markley, F. L., Attitude determination using vector observations: A fast optimal matrix algorithm,
*Journal of the Astronautical Sciences***41**(2) 261-280 (1993). - Söderkvist, I., and Wedin, P. A., Determining the movements of the skeleton using well-conﬁgured markers,
*Journal of Biomechanics***26**(12) 1473–1477 (1993). - Bar-Itzhack, I. Y., REQUEST: A recursive QUEST algorithm for sequential attitude determination,
*Journal of Guidance, Control, and Dynamics***19**(5) 1034-1038 (1996).