3D image alignment: How to extract rotation angles about "static" x,y,z axes from a 3x3 rotation matrix in matlab

Question

I am carrying out 3D imaging alignment using (as an overview) the following coding scheme:

ref = imref3d(size(img(:,:,:,1)),[0 32] ,[0 32] ,[0 20] );
[optimizer, metric] = imregconfig( 'monomodal');
tform=imregtform(img(:,:,:,2), ref, img(:,:,:,1), ref, 'rigid', optimizer, metric);
transform_mat=tform.T

For clarification, the variable img is a time series of 3D images (where each 3D image is a stack of 2D cross sections). At different time points, the object I am imaging can move. Thus, the purpose of the above code is to generate a geometric transformation that optimally aligns the second time point of the 3D volume ( img(:,:,:,2) ) to a reference image, which is the first time point of the 3D volume ( img(:,:,:,1)). The geometric transformation that is finally chosen by the optimization algorithm is output to the object tform. tform contains a 4x4 matrix as one of its properties (T); it is this 4x4 matrix that encodes the translation and rotation information. As you can see from my above code, I stored this 4x4 matrix in the variable transform_mat.

After reading through a lot of scattered mathworks documentation, I determined that this transform_mat variable (representing an affine rigid body transformation matrix) is in a "post-multiplied" form, which, as I understand it, just means that it is the transposed version of what one would traditionally see in a linear algebra text book.

For my purposes, I am interested in extracting specific rotation information from transform_mat. Here is what I have done so far:

rot_postmultiply=transform_mat(1:3,1:3); %extracting the elements that encode rotation-related info
rot_premultiply=rot_postmultply'; %transposing 

Just to quickly interject, I created the premultiplied version of the rotation matrix because I believe that many of the functions that act on rotation matrices assume that it is in its premultiplied form.

From this 3x3 rotation matrix, I want to extract the rotations (in radians) about the STATIC x-axis, y-axis, and z-axis of the reference image. My initial attempt at carrying this out is as follows:

eul = rotm2eul(rot_premultiply);

The rotm2eul function provides me the 3 Euler Angles associated with this 3x3 rotation matrix. eul is a 1x3 vector and, according to documentation, the "default order for Euler angle rotations is 'ZYX' ". However, I am not certain that Euler Angles are actually describing the information I want to extract (i.e. rotations about the static x, y, z axes of the reference image).

I do not have a strong linear algebra/geometric transformation background, but my understanding of Euler Angles is that each rotation that takes place changes the coordinate system. For example, after the rotation about the Z axis (the first value in the eul vector), we now have new X and Y axes (call them X' and Y'). Then the "Y-axis rotation" (the second value in eul) is actually the rotation about Y'...not Y. Repeat this argument for the final "X-rotation" (which would really be about an X'' axis).

If anyone could offer some insight about how to proceed (and if my concerns about Euler Angles are correct), I would greatly appreciate it!

Also, sorry if the word "static" is the incorrect terminology. Hopefully, I have provided sufficient context so that no confusion arises.

Answer 1

If I understand you correctly, all you need is rotm2axang that converts a rotation given as an orthonormal rotation matrix, rotm, to the corresponding axis-angle representation, axang. The input rotation matrix must be in the premultiply form for rotations.

The output of rot2axang is an n-by-4 matrix of n axis-angle rotations. The first three elements of every row specify the rotation axis, and the last element defines the rotation angle (in radians).

If you dont have access to the robotics system toolbox that has that function, consider this alternative instead.

EDIT:

Rotation of a point in 3D space by an angle alpha about an arbitrary axis defined by a line between two points (x1,y1,z1) and (x2,y2,z2) can be achieved by the following steps:

1. translation so that the rotation axis passes through the origin
2. rotation around x axis so that the rotation axis lies in the xz plane (or projection on the xz plane)
3. rotation around y axis so that the rotation axis lies along the z axis (or projection on the z axis)
4. rotation by alpha around the z axis
5. the inverse of step 3
6. the inverse of step 2
7. the inverse of step 1

Are you asking about the angles in steps 2,3?

Otherwise, are you already in the origin for the axis around which you rotate?

If indeed you need the angles in step 2,3 then:

step #1 is done via a the translation matrix T:

T =  [ 1    0   0   -x1
      0     1   0   -y1
      0     0   1   -z1
      0     0   0   1  ];

    

Step #2 is done via the following: Given that v = (a,b,c) is the unit vector along the rotation axis, then d = sqrt(b^2 + c^2) as the length of the projection onto the yz plane. Then rotate by the angle between the projection of rotation axis in the yz plane and the z axis given by the dot product of the z component of the unit vector v and its yz projection. This angle (beta) is determine by: cos(beta) = dot([0,0,c],[0,b,c]) /(c*d) = c/d , sin(beta) = cross( [0,0,c],[0,b,c]) /(c*d) = b/d

The rotation matrix Rx is therefor:

Rx = [ 1    0     0     0
       0    c/d  -b/d   0
       0    b/d   c/d   0
       0    0     0     1 ] ;

Step 3 we rotate around the y axis so that the rotation axis lies along the positive z axis. Using dot and cross product relationships the cosine of the angle is d, the sine of the angle is a. The rotation matrix about the y axis Ry is

Ry = [ d    0   -a  0
       0    1   0   0
       a    0   d   0
       0    0   0   1];

Last, we rotate around z simply by the angle alpha

Rz = [ cos(alpha)   sin(alpha)  0   0
      -sin(alpha)   cos(alpha)  0   0
       0            0           1   0
       0            0           0   1 ] ;

Answer 2

I think the answers to this question should help you out: to decompose a 3D rotation into Cartesian components.