Function geo_nd::vector::axis_of_rotation3

pub fn axis_of_rotation3<V: Float>(rotation: &[V; 9]) -> [V; 3]

Find the axis of rotation of a Matrix3

A rotation matrix R has the property that R.v = k.v for any vector v along the axis of rotation; this is what defines it to be the axis of rotation. This is to say that v is an Eigenvector of R. For pure rotations the value of k is 1; it must be a real number in any case.

Hence to find the axis of rotation of a Matrix3 R requires finding the Eigenvector of R that has a real Eigenvalue - there will always be at least one (more than one implies a pure scaling matrix).

The algorithm used is iterative; it requires the angle of rotation to be non-tiny (i.e. > 1/1000 radians, or about 1/16 of a degree)

Note that R . axis = axis; also we will consider p being any vector perpendicular to the axis.

Consider the matrix RI = R - 99999/100000 * I (where I is the Identity)

Then RI . axis = R.axis - 99999/100000*I.axis = axis * 1/100000

And RI . p = R.p - 99999/100000*I.p, and note then that if the rotation is > 1/1000 radians then |RI.p| > 1/1000

Inverting RI yields RI_inv. Consider (RI * RI_inv) * axis:

axis = RI * RI_inv * axis = RI * (RI_inv * axis)

Since RI * axis = axis / 100000, RI_inv * axis must be 100000*axis

Consider any vector p perpendicular to the axis, and then (RI * RI_inv) * p:

p = RI * (RI_inv * p); this tells us that |RI_inv * p| < 1000

Hence the off-axis component is scaled by less than 1000, where the axis component is scaled by 100000; normalizing will provide a new vector that has the off-axis component reduced by at least a factor of 100.

This provides us with an iterative approach:

 let x = any vector p + k*axis
 let y1 = RI_inv * x / 100000
 // y1 is now k*axis + (p'/100) for some p' such that |p'|=|p|
 let y2 = RI_inv * y1 / 100000
 // y2 is now k*axis + (p''/10000) for some p'' such that |p''|=|p|
 let y3 = RI_inv * y2 / 100000
 // yn will approach the axis of rotation as n increases

etc.

The rate of convergence depends on the angle of rotation; for a rotation of > 1/10 radians the off-axis amount will reduce by around 2^12 each iteration, but for 1/1000 radians it will be only about 2^6 per iteration. For a 52-bit mantissa ten iterations suffices in all cases.

A starting guess for the axis can be (1,0,0); if the actual axis of rotation is the Y-axis or the Z-axis then the algorithm does not converge; and a second guess of (0,1,0) can be used, or even (0,0,1) if that does not converge.

If the rotation is very close to 0 then any axis will do :-)

Another option for this is to consider any vector x and R.x; the cross-product of the two should be the axis of rotation. This clearly has problems if the rotation is by 180 degrees, and is unstable for small rotations (as is any determination of the axis).

Another approach is to just consider it as a rotation matrix

Unit axis vector = (x,y,z); c=cos(angle), t=1-c, s=sin(angle)

Matrix:

t.x^2 + c t.xy - s.z t.xz + s.y t.yx + s.z t.y^2 + c t.yz - s.x t.zx - s.y t.zy + s.x t.z^2 + c

Adding the digaonal yields 3c + t = 2+2c Then subtracting across the diagonal yields 2sx, 2sy, 2sz

However, if 2+2c is close to 0 then this fails