How to Calculate the Axis of Rotation in 3D Space?

[yasith]

hey guys,

I have a problem I need to solve for my masters thesis. I have an object that rotates in 3D space. There is no translation. Given that I have coordinates of one landmark at 2 instances of time I need to calculate axis of rotation and angle of rotation.

is this sufficient information?
if so can someone please help me with the working.

Regards
Yasith (Medical Doctor)

 

[M Quack]

If you also know the center of rotation that is enough.

(1) calculate the vector from the center of rotation to the initial and final points.

(2) take the cross product of these two vectors. That is the axis of rotation. (you may need to normalize)

(3) normalize the two vectors and take the dot product. That is the cosine of the angle of rotation.

(4) The cos does not really tell you if the angle is positive or negative (i.e. it does not change when you exchange the two vectors). However, the cross product changes sign. So the angle and direction do provide the complete information. (Sorry, but I don't remember off the top of my head if you have to choose the positive or negative value of the arccos.)

(5) If you don't know the pivot point, then you need at least 3 points.

Hope this helps.

 

[yasith]

thank you for the reply,

yes the pivot point is not known

i can get 3 landmark points.

if this is the case how can I approach this problem?

 

[D H]

No, it doesn't, not even if you know the center of rotation.

Suppose you use M Quack's procedure and determine that the object has rotated through 140 degrees between your two time points. How do you know it hasn't rotated through 500 degrees (140+360), or 860 degrees, or even more? There's aliasing going on here thanks to the wrapping function. Another aliasing problem arises due to the axis of rotation having a sign ambiguity by this method. Flip the sign and that 140 degrees becomes 220 degrees (or 580, or ...).

Yet another problem is that it is rarely, if ever, a good idea to use the minimal number of measurements.

A final problem is that the object may not have a fixed axis of rotation if this is an object rotating in free space. You will at best get a mean axis of rotation for the time interval of question. The Euler torque (aka "the polhode rolls without slipping on the herpolhode") means that such objects will have a different axis of rotation at some later time.

 

[yasith]

thanks D H for your input.

The problem I have at hand is to map the motion of a bone.
There is no translation component to this problem. It is pure rotation.

I am mapping landmarks on the bone at each 1 second interval.

therefore between each time point I want to map the instant axis of rotation and the angle of rotation.

Any ideas on how I can approach this problem?

 

[D H]

Any ideas on how I can approach this problem?

Search for hand+pose and hand+motion at Google scholar or CiteSeer. Eun-Jung Holden's 1997 PhD thesis, "Visual Recognition of Hand Motion" appears to be directly applicable.

Addendum
Another good one is Erol et al, "Vision-based hand pose estimation: A review", Computer Vision and Image Understanding 108 (2007).

You have a pose estimation and motion detection problem here. There is lots and lots of technical literature on these subjects.

 

[Vargo]

One point rotating around an axis traverses a circle. If you observe it at point A, B, C then those are three points on the circle. Let's say we want to know the measure of the arc from A to C. According to Euclid, it is 360 minus double the measure of the angle ABC. The axis of rotation is simple to calculate because it is perpendicular to both AB and BC, so you can take the cross product to find the right direction. Let N be the direction of the axis of rotation.

The axis of rotation is the intersection of the perpendicular bisectors of the segments AB and BC. The equations for the first plane is:
$(A-B)\cdot X = (1/2)(A\cdot A - B\cdot B)$.

The second plane has a similar equation. If you find just one point P that satisfies both equations, then you can easily find the direction of the axis by taking the cross product of A-B and C-B (since it has to be perpendicular to both of those). A point and direction characterize a line, so you're done.

 

[M Quack]

No, it doesn't, not even if you know the center of rotation.

Suppose you use M Quack's procedure and determine that the object has rotated through 140 degrees between your two time points. How do you know it hasn't rotated through 500 degrees (140+360), or 860 degrees, or even more? There's aliasing going on here thanks to the wrapping function. Another aliasing problem arises due to the axis of rotation having a sign ambiguity by this method. Flip the sign and that 140 degrees becomes 220 degrees (or 580, or ...).

True in the general case, but probably not a problem here. But you are right, I should have pointed out this limitation.

Yet another problem is that it is rarely, if ever, a good idea to use the minimal number of measurements.

I fully agree with this. No redundancy = no idea about the errors of your measurement.

A final problem is that the object may not have a fixed axis of rotation if this is an object rotating in free space. You will at best get a mean axis of rotation for the time interval of question. The Euler torque (aka "the polhode rolls without slipping on the herpolhode") means that such objects will have a different axis of rotation at some later time.

I was indeed thinking about 1 point (mapping landmark) at 3 different times to define an arc. For a bone that can rotate in several directions this is not a good idea, you will get false results.

It would be much better to identify as many landmark points as possible and then calculate the center of rotation, axis and angle from these (including some estimate of the error). I am sure that there more sophisticated and elaborate methods available in the literature than I can come up with here :-)

BTW, yasith, what kind of input data are you using?

 

[yasith]

this is a diagram to the problem
I have a perfect 3D volume rendered model of the wrist in 3D motion which I take measurements from

I have managed to map all the points of interest to a 3D coordinate system.
Therefore my input is p = [x,y,z].

The problem is to calculate instant axis of rotation. Therefore Vargo I can't assume the center of rotation will stay constant across point A -> B -> C.

 

[Tinyboss]

Mathematically speaking (in other words, being anal about not making assumptions, to the point that my answer might not be relevant to your problem!), you do not have enough information to determine the rotation axis and angle. For a given pair of points, there are infinitely many (axis, angle) pairs that move point A to point B.

Think of a globe, and imagine that you have a rotation that moves the north pole (point A) to some specific point B on the equator. The axis itself could go through the equator (90 degrees east and west from point B), in which case point A moves along a line of longitude. Or the axis could go through that line of longitude (exactly halfway between A and B), so that point A moves along a smaller circle. And there are infinitely many more axes that would work, each with a different angle.

If you first drew an arrow on the globe, pointing in some direction from the north pole, then observing the direction of the arrow after the rotation would allow you to deduce the axis and angle. And maybe the fact that your object isn't a sphere implicitly gives you this "arrow"? Anyway, you either need that, or a second landmark point (which is not antipodal to the first) to get the information you want.

 

[yasith]

Mathematically speaking (in other words, being anal about not making assumptions, to the point that my answer might not be relevant to your problem!), you do not have enough information to determine the rotation axis and angle. For a given pair of points, there are infinitely many (axis, angle) pairs that move point A to point B.

Think of a globe, and imagine that you have a rotation that moves the north pole (point A) to some specific point B on the equator. The axis itself could go through the equator (90 degrees east and west from point B), in which case point A moves along a line of longitude. Or the axis could go through that line of longitude (exactly halfway between A and B), so that point A moves along a smaller circle. And there are infinitely many more axes that would work, each with a different angle.

If you first drew an arrow on the globe, pointing in some direction from the north pole, then observing the direction of the arrow after the rotation would allow you to deduce the axis and angle. And maybe the fact that your object isn't a sphere implicitly gives you this "arrow"? Anyway, you either need that, or a second landmark point (which is not antipodal to the first) to get the information you want.

thank you tinyboss for your response.

I can get track as many points as I need to solve this problem.

So the solution to this problem can be achieved by tracking 2 points at 2 time instances?

if so can someone be kind enough to show me the working for this solution.

also is there any software/matlab coding to approach this problem?

Thank you.

 

[Tinyboss]

Assuming you also don't know the center of rotation, I'd proceed this way:

1. If point A moves to point A', then you know that the center of rotation is somewhere on the equidistant plane (the plane of points in 3D space which are at equal distance from A and A').
2. Now use another pair of points B and B' to get another plane. Intersecting these planes gives you a line, and you know the center of rotation is somewhere on that line.
3. Use a third pair C and C' to get a third plane, and the intersection of it with your line will be a single point, which is the center of rotation.
4. Now that you know the center, you can find the axis and angle using at least two before/after pairs. I would have to work out the recipe for that; I don't know it off the top of my head, sorry!

Edit: Of course, I'm assuming there is no measurement error. Three nonparallel planes will always intersect in a single point, so that part will be okay. But when you start to work out the axis and angle, you might find there is no solution that moves all point pairs exactly as required. I don't have any experience in dealing with that kind of issue.

Good luck!

 

[Vargo]

Earlier you said you could take 1 measurement per second... Can you speed up those measurements so that you are taking them every .1 s? If so, then we can do some numerical calculus here and solve these problems.

If not then we could expand on Tinyboss' suggestion. With 4 measurements, you could find the location O of the joint. All axes of rotation must pass through O, so given two measurements A and B, the axis of rotation would lie in the plane that is the perpendicular bisector of AB and it would pass through the point O. Unfortunately that leaves us with an infinite number of axes of rotation... So even knowing O is not enough information to pin down the axis of rotation from just two measurements. There is also the problem that the wrist joint is hardly a point, so whether there is a stable point O for all motions of the wrist is doubtful in the first place.

I think what you need Yasith is measurements over short timescales. Then it would make sense to discuss "instantaneous" quantities and it would be practically feasible to calculate them using calculus.

 

[D H]

yasith, you mentioned this is a project for a masters thesis. What literature searches have you done? I quickly found several and posted a couple, one from over a decade ago, the other from about four years ago.

 

[M Quack]

I thought I had typed this before, but apparently it got lost...

You have a set of 3D coordinates of points before and after the rotation,
p --> p'. You want to find the center of rotation COR, c, the axis of rotation (e.g. defined by Euler angles), and the angle of rotation.

That makes 3 variables for the x,y,z coordinates of the COR, one for the angle, and 3 for the Euler angles (but I think you can get away with 2).

The equation describing the transformation for each points is

$\vec{p}' - \vec{c} = U \cdot R(\theta) \cdot U^T \cdot (\vec{p} - \vec{c})$

$R(\theta) = \left( \begin{array}{ccc} \cos(\theta) & -\sin(\theta) & 0 \\ \sin(\theta) & \cos(\theta) & 0 \\ 0 & 0 & 1 \end{array}\right)$
describes a rotation about the z-axis by the angle $\theta$

$U = (\hat{u}_1, \hat{u}_2, \hat{u}_3)$ is a coordinate transformation from the z-axis to the acutal axis of rotation (u_1,2,3 are the basis vectors, the rotation will be around u_3). I have to look up how to best express U in terms of Euler angles.

What you need to do is fit the free variables (coordinates of the COR, θ, and the orientation of the axis of rotation) of this expression to the data. It might be helpful to use a simpler procedure to produce a good initial guess that assures the rapid convergence of the fit.

 

[meldraft]

An easier(?) and dirtier solution would be this:

Your (instantaneous) center of rotation is defined as the point (or axis) that has the following property:

A line connecting the center of rotation to any point on the body is perpendicular to the velocity vector at that point. If the center of rotation is an axis, this line must be perpendicular to the axis as well.

Therefore, you can get 3 velocity vectors on your body and define 3 disks in 3d space with the velocities' normal vectors. These will intersect at least at 1 point (or define the axis directly if you're lucky). Repeat this one more time for 3 more vectors, and you end up with 2 points in 3D space. Connect them and you will have the (instantaneous) axis of rotation. You can repeat this process for more sets of points if you want to verify your result.

One pitfall here is that the second point you define can be the same as the first. This means that all of your new velocity vectors are parallel to the first ones, meaning you have to calculate the velocity at different points

edit: if you want to do it without tracking too many points, 4 in total should be enough. You can then do the process above by using only 3 of them at a time.

 

[M Quack]

meldraft, it sounds as if this could work. However, I don't see an easy systematic way to use the redundancy and estimate errors using this method.

Least-squares fits are standard in most data analysis packages these days, so my method should be easy to implement.


$U = \left(\begin{array}{ccc} \cos(\phi) & -\sin(\phi) & 0 \\ \sin(\phi) & \cos(\phi) & 0 \\ 0 & 0 & 1 \end{array}\right) \cdot \left(\begin{array}{ccc} \cos(\chi) & 0 & \sin(\chi) \\ 0 & 1 & 0 \\ -\sin(\chi) & 0 & \cos(\chi) \end{array}\right) = \left(\begin{array}{ccc} \cos(\phi) \cos(\chi) & -\sin(\phi) & \cos(\phi)\sin(\chi) \\ \sin(\phi) \cos(\chi) & \cos(\phi) & \sin(\phi)\sin(\chi) \\ -\sin(\chi) & 0 & \cos(\chi) \end{array}\right)$,

where $0 \leq \phi < 2\pi$ and $0 \leq \chi \leq \pi$

The axis of rotation is the $\left(\begin{array}{c} \cos(\phi) \sin(\chi) \\ \sin(\phi) \sin(\chi) \\ \cos(\chi) \end{array}\right)$