MATLAB: Extracting a rotation axis from a rotation matrix

Question

I am working with 3 x 3 rotation matrices, R(ij), of cubic bicrystal misorientations, so there are up to 24 symmetry related descriptions (i.e. different permutations of the rotation matrix) for any starting misorientation. Usually I extract the axis/angle description from each one by standard methods – acos (theta) = acos(((traceR)-1)/2), and the axis as [R32-R23, R13-R31, R21-R12]. This works fine, except when theta = 180 degrees, for which this formula gives the axis as [0 0 0], as we know. So to get the axes for the 180 case, it seems to me I can do the following:

1) Use [V,D] = eig® to get a matrix D whose columns are the eigenvalues of R. One of these will be equal to +1, the other two will be -1.

2) Identify which column of D contains the +1, and then take the equivalent column of V as the eigenvector, i.e. the rotation axis. So far so good.

I am trying to use a for-(else)-end loop to do this, by using a counter i to examine each column of D in turn, and then summing the three elements of that column. If the sum comes to +1, I want to use that corresponding i to obtain the eigenvector as:

[V(1,i), V(2,1), V(3,1)], and then come out of the loop. I am relatively new to MATLAB, and I cannot get a loop to work. Can anyone help point me to a working code to do this loop?

Thanks.

Answer 1

tol = 100*eps;
i = find(abs(diag(D)-1)<tol);
rotAxis = V(:,i);

(Edited to take rounding error into account. Increase tol if necessary.)