MATLAB Central Newsreader - Find instant axis of rotation

Find instant axis of rotation

Henning Olsson — Tue, 27 Oct 2009 20:52:03 -0400

I need to estimate/calculate the location of a rotation axis for a number of points in space. I have the position of the all of points (x,y,z) at all of the different times. The relative location of all the points remain the same as they are all part of a rigid body.

I started out by creating a rotation matrix around the origin (by using 3 points in step n and step n+1 and then solving for the rotation matrix). I then used an Axis-Angle approach to calculate the axis around which they rotate.

The problem is that I am pretty sure that for the actual motion the axis doesn't go through the origin but is translated somehow in space. But I'm not sure how to find this translation. Anyone have any ideas?

Another idea I had was to to "fit" a cylinder to the points. If the points from two steps are all on the surface on the cylinder, the cylinder axis must be the instant axis of rotation. Any ideas on how to do such a calculation (preferably analytically).

Re: Find instant axis of rotation

Jan Simon — Tue, 27 Oct 2009 22:17:04 -0400

Dear Henning!

> I need to estimate/calculate the location of a rotation axis for a number of points in space. I have the position of the all of points (x,y,z) at all of the different times. The relative location of all the points remain the same as they are all part of a rigid body.
>
> I started out by creating a rotation matrix around the origin (by using 3 points in step n and step n+1 and then solving for the rotation matrix). I then used an Axis-Angle approach to calculate the axis around which they rotate.
>
> The problem is that I am pretty sure that for the actual motion the axis doesn't go through the origin but is translated somehow in space. But I'm not sure how to find this translation. Anyone have any ideas?
>
> Another idea I had was to to "fit" a cylinder to the points. If the points from two steps are all on the surface on the cylinder, the cylinder axis must be the instant axis of rotation. Any ideas on how to do such a calculation (preferably analytically).

Ask Google for "helical axis" or "screw parameters".
On the KWON3D site, the orientation of the axis is not clearly defined: If you assume a positive angle, the anti-symmetrical H-H' defines the orientation of the axes.
A simple method is shown in:
http://robotics.caltech.edu/~jwb/courses/ME115/handouts/rodriguez.pdf

Be aware that there are some numerical problems, e.g. if the coordinates are noisy or the rotational angle is tiny...

Good luck, Jan

Re: Find instant axis of rotation

Bruno Luong — Wed, 28 Oct 2009 07:18:00 -0400

"Henning Olsson" <asd@prutcom.com> wrote in message <hc7mhj$69k$1@fred.mathworks.com>...
> I need to estimate/calculate the location of a rotation axis for a number of points in space. I have the position of the all of points (x,y,z) at all of the different times. The relative location of all the points remain the same as they are all part of a rigid body.
>
> I started out by creating a rotation matrix around the origin (by using 3 points in step n and step n+1 and then solving for the rotation matrix). I then used an Axis-Angle approach to calculate the axis around which they rotate.
>
> The problem is that I am pretty sure that for the actual motion the axis doesn't go through the origin but is translated somehow in space. But I'm not sure how to find this translation. Anyone have any ideas?

1. Rotation axis (direction only) R is perpendicular to the velocity field. So estimate R is easy, for example as following:
Minimize Sum | Projection R on V_i |^2 such that |R|=1
| Projection R on V |^2 = <R,V>^2 / | V |^2

2. Project the velocity to the plane perpendicular to R, V projected is tangential to the circle centered about the axis center. The magnitude is proportional to the radius. Write this as equation than solve it for he axis center.

Bruno

Re: Find instant axis of rotation

Henning Olsson — Wed, 28 Oct 2009 22:54:02 -0400

Thank you both for your suggestions,

I did some work with based on that pdf that was posted (mainly because I have very limited experience of vector fields).

It seems like I can successfully calculate the direction of the axis (eq 8 in pdf). But I have problems getting the translation of the axis correct (eq 10).

I've tried with many different motions/data points but always get more or less the same translation (but different directions).

Here is my code with two example points (I've also have tried increased accuracy on the data points with out success).

Once again, I would be very grateful for all comments/suggestions!

Thanks!

%points before
p0=[-199.177 399.997 172.887]';
q0 = [-199.177 -400.003 172.872]';
r0 = [-2479.583 244.996 225.422]';
%points after
p1 = [-199.056 399.997 172.048]';
q1 = [-199.056 -400.003 172.034]';
r1 = [-2479.425 244.996 226.179]';
% %

w = cross(((q1-q0)-(r1-r0)),((p1-p0)-(r1-r0)))/dot(((q1-q0)-(r1-r0)),((p1+p0)-(r1+r0))); %Calculates rotation vector

rot_angle = norm(w) %this is tan(theta/2) where theta is actuall rot angle (this
angle matches really well with data output)

w=w/norm(w); %normalizes rotation axis for further calculations

p = 0.5*((cross(w,(p1-p0))/rot_angle)-(dot(w,p1+p0).*w+p0+p1)); %This is where something goes wrong

Re: Find instant axis of rotation

Jan Simon — Thu, 29 Oct 2009 10:43:03 -0400

Dear Henning Olsson!

> p0=[-199.177 399.997 172.887]';
> q0 = [-199.177 -400.003 172.872]';
> r0 = [-2479.583 244.996 225.422]';
> %points after
> p1 = [-199.056 399.997 172.048]';
> q1 = [-199.056 -400.003 172.034]';
> r1 = [-2479.425 244.996 226.179]';

What values do you expect?
I get:
  Angle: 0.000699
  Translation: -0.000261
  Helical axis: [-0.00178, 0.999998, 2.792e-5]
  Origin of axis: [-1397.32, -2.4959, -04475]
This set of values is not unique, e.g. all points on the axis can be used as origin. I've defined the point with the minimal distance to the global origin.

Please be aware, that the rotational angle of your example is tiny and rounding errors can dominate the result! There is no difference between a rotation around a point far away and a pure translation...

Jan