Hi
Could I ask if there someone has a code to calculate the axis of rotation please.
I have two STL files, bones, which has moved from time point 1 to time point 2.
I want to calculate the displacement. So I calculate the transformation matrix by rigid registration. I would like to find out the axis of rotation, between the two.
Many thanks
I think you are looking for the Euler axis. Here is a formula to calculate it given a transformation matrix A:
Thanks @mau_igna_06
I am a medical doctor, so no experience in coding
I am after “axis of rotation” (yes I am calculating Euler angles for rotation)
Is there a code that I can get it
Thanks
(the euler angle with euler axis formula is different than yaw, pitch, row angles formula. I think you need the euler axis)
Having into account that you have no programming experience I suggest you open Transforms module and look for the transform matrix you just created. You only need its elements to use the formula. In the formula where it says A13 (for example) it means you need to use the element in first row and third column of matrix A (the matrix you just created). You can calculate each ecuation typing the values and operations in the python interactor. To calculate arccos you can write the line import math, press enter then write math.acos(X) where X is a number and press enter.
Divisions are done with / character and multiplications are done with * character.
So you can write 12/4 and the python interactor will print 3 as the result.
The axis of rotation you want is a vector which components are e1, e2 and e3
To calculate sin you can use math.sin(X) for example in the python interactor you can write math.sin(X) and press enter and it will print 0.8414709848078965 that is the result.
Here you have example code of how to make basic arithmetic operations in the python interactor (or you can just use a calculator):
https://www.tutorialspoint.com/python/arithmetic_operators_example.htm
What do you mean by this exactly? Are you interested in rotation angles? Or, would you like to get the direction of the axis that can orient one model to match the orientation of the other model by a single rotation?
What is the clinical problem you would like to solve/question you want to answer?
Thanks
Yes, I am able to calculate the rotation angles in a given coordinate system
Instantaneous centre of rotation and instantaneous axis of rotation?- axis fixed to a body undergoing planar movement that has zero velocity at a particular instant of time.
Clincal problem- I have wrist bones that I am studying, rotating as the wrist moves.I have 2 models (same bone) , one at time point 1 and second at time point 2. So I am able to calculate rotation angles using transforms module.
But I would like to calculate the "Instantaneous axis of rotation- from time point 1 to point 2) I believe it is the axis around which the bone rotates. Hope this makes sense
If you have a single relative transform then copy-paste this rotation_from_matrix function to the Python console, and then to get rotation axis from a transform node you can run this:
transformNode = getNode('Transform') # replace 'Transform' with the actual name of your transform node
import numpy as np
import math
transformMatrix = arrayFromTransformMatrix(transformNode)
angle, axisDirection, axisPosition = rotation_from_matrix(transformMatrix)
# Display rotation axis
rotationAxisMarkupsNode = slicer.mrmlScene.AddNewNodeByClass("vtkMRMLMarkupsLineNode", "RotationAxis")
rotationAxisMarkupsNode.AddControlPoint(vtk.vtkVector3d(axisPosition[:3]-50*axisDirection))
rotationAxisMarkupsNode.AddControlPoint(vtk.vtkVector3d(axisPosition[:3]+50*axisDirection))
Getting rotation axis from a single relative transform is probably not very robust or accurate. It would be better to compute an “average” axis from multiple transforms, similarly how Pivot calibration module in SlicerIGT extension does it (see a recent discussion about the algorithm that the module uses here).
Thanks @lassoan
I get an error message NameError: name ‘rotation_from_matrix’ is not defined
Thanks
Copy-paste this rotation_from_matrix function to the Python console before running the code snippet above.
Thanks, @lassoan
Really appreciate your advice, Not sure what I am doing wrong
I get this
def rotation_from_matrix(matrix):
… “”“Return rotation angle and axis from rotation matrix.
…
… >>> angle = (random.random() - 0.5) * (2*math.pi)
… >>> direc = np.random.random(3) - 0.5
… >>> point = np.random.random(3) - 0.5
… >>> R0 = rotation_matrix(angle, direc, point)
… >>> angle, direc, point = rotation_from_matrix(R0)
… >>> R1 = rotation_matrix(angle, direc, point)
… >>> is_same_transform(R0, R1)
… True
…
… “””
… R = np.array(matrix, dtype=np.float64, copy=False)
… R33 = R[:3, :3]
… # direction: unit eigenvector of R33 corresponding to eigenvalue of 1
… l, W = np.linalg.eig(R33.T)
… i = np.where(abs(np.real(l) - 1.0) < 1e-8)[0]
… if not len(i):
… raise ValueError(“no unit eigenvector corresponding to eigenvalue 1”)
… direction = np.real(W[:, i[-1]]).squeeze()
… # point: unit eigenvector of R33 corresponding to eigenvalue of 1
… l, Q = np.linalg.eig(R)
… i = np.where(abs(np.real(l) - 1.0) < 1e-8)[0]
… if not len(i):
… raise ValueError(“no unit eigenvector corresponding to eigenvalue 1”)
… point = np.real(Q[:, i[-1]]).squeeze()
… point /= point[3]
… # rotation angle depending on direction
… cosa = (np.trace(R33) - 1.0) / 2.0
… if abs(direction[2]) > 1e-8:
… sina = (R[1, 0] + (cosa - 1.0) * direction[0] * direction[1]) / direction[2]
… elif abs(direction[1]) > 1e-8:
… sina = (R[0, 2] + (cosa - 1.0) * direction[0] * direction[2]) / direction[1]
… else:
… sina = (R[2, 1] + (cosa - 1.0) * direction[1] * direction[2]) / direction[0]
… angle = math.atan2(sina, cosa)
… return angle, direction, point
…
transformNode = getNode(‘LinearTransform_3’)
import numpy as np
import math
transformMatrix = arrayFromTransformMatrix(transformNode)
angle, axisDirection, axisPosition = rotation_from_matrix(transformMatrix)
Traceback (most recent call last):
File “”, line 1, in
NameError: name ‘rotation_from_matrix’ is not defined
Copy-paste the rotation_from_matrix function implementation to the Python console first (that you included in the top of your previous post). After you copy-paste and hit Enter key once or twice, you should see the command prompt (>>>) again and no error message. That means that the function is now available.
You can then copy-paste the rest of the code to the Python console.
Thanks a lot.
It worked fine, I trialled with two models (create model module) with known rotation. Thank you very much. It draws the AoR on the 3D viewer.
I have two questions.
How to get a numerical value for this axis. (With the code, in the console, it types 0 or 1 only, I do not get a value)
I have pasted an image of my transform (resulted by model registration). _this is after registering bone position 1 to position two using ICP. When I apply the same code for this transformation matrix- Linear transform 5 in this case- All models disappear. (I have positioned my models centred on the coordinate system o, o, o)
I think, this is expected for real-life motion, where you cannot ensure that the rotation happens exactly along a single axis. In general, an orientation difference cannot be described by rotation along a single axis, but you need to rotate along two axes (or around a point). If you have several (preferably dozens) of different poses then from those you can compute an average instantaneous axis of rotation, even if there are slight additional rotations (similarly how Pivot calibration module computes the rotation axis in SlicerIGT extension), but I’m not sure if this can be applied to just two poses.
@mholden8 Do you have any recommendation?
Thanks @lassoan
I elicited the same behaviour with the created models. Whenever I introduce few mm translation to a model, this happens.
I think that is the reason with the bones. If I strictly take two STL files (of the same bone), which only is slightly displaced between time points (10-15 degrees of rotation in all three planes and 1-3 mm translation ) I believe, I still give a ‘theoretical’ instantaneous axis of rotation.
What are your thoughts? Thanks
The function that I’ve found and linked above can only compute the instantaneous axis of rotation if it exists (if the transformation can be reproduced by rotation along a single axis). If you also have a translation component then this computation fails. There are infinite number of solutions for decomposing a linear transformation to rotation along a single axis and some additional rotations and/or translations and you need to figure out that in what sense you would like an optimal solution. I would suggest to find papers in your field that compute the rotation axis and try to implement the method that they describe (or if it is not clear how exactly they do it then contact them and ask for the implementation). If you have an exact computation formula or implementation in any language then we can help with porting it to Python/Slicer. You may also contract an R&D company, such as Kitware, or a computational geometry expert to work this out for you.
Thanks for advice @lassoan
I will define the calculation a bit more from previous work
(Also will update slicer support so that if we can have this functionality, it will be useful for the users, As any bone that moves in human body is likely to have a translational component, even its few mm)
Thanks for kind support @lassoan
Found this, after some research
"Similarly we can also calculate backwards and find the intantanous axis of rotation at some distance
−
⃗
r
from any point if we know the linear velocity of the point and the angular velocity. This is given by:
[\begin{equation} \label{eq:2} -\vec{r}= \frac{\vec{\omega} \times \vec{v}_P }{ \vec{\omega}\cdot\vec{\omega} } \end{equation}]
Thus if know the position of a point
⃗
r
P
, its linear velocity
⃗
v
P
we can find closest point
⃗
r
C
on the rotation axis as:
[\vec{r}_C = \vec{r}_P - \vec{r} = \vec{r}_P + \frac{\vec{\omega}\times\vec{v}_P}{\vec{\omega}\cdot\vec{\omega} }]
Script
AnyKinRotational Rotational ={
AngVelOnOff = On;
AnyRefFrame& Ref1= .ReferenceFrame;
};
AnyKinLinear Linear ={
Ref=-1;
AnyRefFrame& Ref1= .ReferenceFrame;
};
/// Direction of the instantaneous axis of rotation
AnyVec3 e_iaor = Rotational.Vel/vnorm(Rotational.Vel)
/// The point on the rotation axis closest to ReferenceFrame origin
AnyVec3 r_iaor = Linear.Pos + cross(Rotational.Vel, Linear.Vel)/(vnorm(Rotational.Vel)^2);
Unfortunately I cant use the same software, as it only runs on Windows
Is there any possibility that slicer could enable this function, Thanks
I cannot decipher neither the formulae nor the code. Can you provide a link to the source of this code snippet?