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

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.

[Matt J]

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**

No, the terms used to distinguish decompositions about static versus nonstatic axes are "extrinsic" (static) rotations versus instrinsic (non-static) rotations.

Technically, the term "Euler angles" is supposed to refer to a decomposition in which two of the rotation axes are the same, e.g. Z-X-Z. When all the axes are different, e.g., Z-Y-X, the angles from the decomposition are called Tait-Bryan angles. That said, the terminology is often abused, and sometimes "Euler angles" is used when really the person means Tait-Bryan.

[Matt J]

Note also that conversion between extrinsic and intrinsic angles is extremely simple. Given an intrinsic rotation about the nonstatic axis sequence Z-Y'-X'' by angles gamma then theta then phi, the same rotation can be represented as an extrinsic rotation about static axes X-Y-Z by angle sequence phi then theta then gamma. See also this thread,

https://www.mathworks.com/matlabcentral/answers/2074066-convert-quaternion-to-euler-angle-extrinsically?s_tid=srchtitle

Answer 1

[Sai Pavan]

Hello Cramer,

I can conclude from your description that your understanding of the Euler angles is right. Euler angles indeed represent intrinsic rotations about the axes of a coordinate system, and after each rotation, the coordinate system is altered, which means subsequent rotations are about the new, rotated axes. However, I understand that you are looking for is extrinsic rotation, which is a rotation about the static axes of the reference frame.

The "rotm2eul" function can be used to extract these extrinsic angles, but we need to specify the correct rotation sequence that represents extrinsic rotations. The 'ZYX' order in "rotm2eul" by default is for intrinsic rotations. For extrinsic rotations about static axes, we need to reverse this order to 'XYZ', which means a rotation first about the x-axis, then the y-axis, and finally the z-axis.

Please refer to the below code snippet that illustrates the workflow to perform extrinsic rotation:

rot_postmultiply = transform_mat(1:3, 1:3); % Extract the 3x3 rotation matrix from the 4x4 transform matrix
rot_premultiply = rot_postmultiply'; % Since the matrix is in post-multiplied form, transpose it to get the pre-multiplied form
eul = rotm2eul(rot_premultiply, 'XYZ'); % Extract the extrinsic rotation angles about static x, y, and z axes

Please refer to the below documentation to learn more about axis-rotation sequence in "rotm2eul" function: https://www.mathworks.com/help/nav/ref/rotm2eul.html#:~:text=sequence%20%E2%80%94%20Axis%2Drotation%20sequence

Hope it helps!

[S. Cramer]

Thank you for the comments! It has been some time since I wrote this question and am glad that someone could finally offer input. Quickly, though, based off of your explanation, wouldn't the XYZ specification only provide me with the extrinsic rotation about the X (static) axis? (just like how the ZYX default only proides me with extrinsic rotation information about the Z axis) Doesn't any 3 letter permutation of X, Y, and Z only provide me extrinsic rotation information about axis that is encoded as the first letter?

Answer 2

[Matt J]

Attached is the utility that I use. In alignment with what you are pursuing, the code does assume the rotation matrix is a column vector transformer, applied via pre-multiplication. In your case, if the rotation is to be first about the X axis, then about Y, and finally about Z, the syntax you would use is,

 extrinsicAngles = rot2taitbryan(transform_mat,'zyx')