Creating a new coordinate system

Hello all,

I am trying to create a local coordinate system within the global coordinate system in order to compute the displacement of a bone block between two different scans. The reason why the global coordinate system does not suffice is because I would like the rotational displacement of the bone block to be in terms of the orientation of the reference bone block rather than the global space.

Could someone please provide some guidance on how to create a local coordinate system?

Thank you very much in advance,

Michele

@matsuba8,

I am not sure about creating a ‘local’ coordinate system, but here are two things that have worked for me in the recent past;

  1. Center all of your models to the origin of the scene itself. You can do this through the Python interactor by calculating the centroid of your model, then the distance between the centroid and the origin, and then translating the model to the origin. If you do this, then your rotation will be about the centroid of the model, as you would hope. If you have a secondary model that can be placed at an arbitrary distance from the origin, now you know that arc-length of its rotation. This may take a bit of programming, but it has helped me in a few workflows I am developing :sweat_smile:

  2. You can rotate nodes in slicer about ‘points’ or ‘lines’ using the Python interactor also. Maybe you only need one of your objects to rotate about a long bone, where the length of the bone is defined as axis of rotation…?

Perhaps you can expand a little on your intended goal and current problem. What are these ‘bone blocks’? Are you trying to recreate the rotation of a joint from bones segmented from different scans?

Let me know!

Hi Fluvio,

Thank you very much for your reply!

  1. I’ve used the same method as you for translating the centroid to the origin, but I didn’t know there was a way to use the arc-length to calculate the rotation! This is good to know. I am quite an inexperienced coder (and by that I mean 0 experience, haha), but this seems very useful!
  2. I have not tried this method, but it may not be ideal, as it is crucial in my project that the three axes are orthogonal to one another.

To explain my project further, I have attached images below.


The model seen in the 3D view is the bone block of interest, and the blue model is the neutral condition while the red model is the condition (models are overlaid on top of each other). While they appear to be registered, they are not, as the difference in rotation we are looking for is expected to be ≤ 1˚. In the images attached, the centroids of both the neutral and the condition bone block have been moved to the origin.
Our objective is to find the amount of translation and rotation of the condition bone block relative to the neutral bone block.
In method #1 which you have explained in your reply, did you move the centroids of both models to the origin? or just the reference model?
Further, it would be awesome if you could provide me with some guidance as to where I could learn the coding for using the arc length to calculate the rotation.

Thank you!

Michele

P.S. Please note that the red axes in both images were added to a screenshot rather than created in 3d slicer. The top image represents the global coordinate system, and the second image is what I was hoping for the local coordinate system to look like.

1 Like

@matsuba8,

Glad I could help!
I am also happy you shared more details, as I think there are better strategies that can help you.

Our objective is to find the amount of translation and rotation of the condition bone block relative to the neutral bone block.

Two approaches;

  1. Model-to-Model Distance. Since the “bone blocks” are similar in shape, you could use this module to calculate the point-to-point displacement between the two. I am pretty sure you can then use this displacement map to approximate an overall transformation matrix that describes the overall translation and rotation matrix that describes the displacement of the part.

  2. Fiducial Registration Wizard. You can use fiducials to determine the rigid transformation that describes the displacement of the part. This module returns the transformation matrix itself! You have to download the SlicerIGT module to use this tool.

In method #1 which you have explained in your reply, did you move the centroids of both models to the origin? or just the reference model?
Further, it would be awesome if you could provide me with some guidance as to where I could learn the coding for using the arc length to calculate the rotation.

If you use the two new solutions I propose here, you wont need to program at all :rofl:. Assuming the result is what you want of course.

Now, I guess the other question is… are the frames of reference aligned? In other words, are the two imaging datasets (DICOM) aligned? If they are not, you may need to align them to remove that additional transformation from your final displacement. Does that make sense?

Hi again,

I have not tried using

yet, but I have tried using the IGT model registration module. The only problem with this is that when I use the matrix to extract the rotational values, it’s in terms of the global coordinate system. I have attached an image with planes running through the model, to demonstrate how I would like the coordinate system to be oriented.


In terms of

I am unsure whether the bone blocks have enough reliable landmarks for it to be accurate. Further, I believe the same problem regarding the coordinate system still exists for this method.

I have aligned the two datasets already:)

1 Like

This bone looks like a thick elliptical disk. If I were you I would make PCA to the points of the bone mesh. With that you would get two vectors that you could use to create a frame (coordinate system) I think. Hope it helps

@matsuba8,

What about;

  1. Moving both models to the center of the scene (origin)
  2. Creating three orthogonal planes that delineate the local coordinate system (or more like the basis) on both models (e.g. xy1, xz1, yz1, xy2, xz2, yz2)… similar to what you did on your screenshot…
  3. Calculate the angles between corresponding planes using the Python interactor and this example… this should be equivalent to a dot product between vectors, so the angle will not be tied to the overall coordinate system as in “angle planes” module you are using.

The challenge here is making sure those planes are made with the same features/landmarks on the model. My guess is that you already thought of this though.

Hello Mauro,

Thank you for your response. I am unfamiliar with the term PCA, but is that “Principal Components Analysis”?
I will look into this, thanks again.

Michele

1 Like

Hi again Fluvio,

Yes, I have thought about this but I haven’t found a way to make sure that the planes are orthogonal to each other… Do you know how to do this?
Perhaps if I can find a way to ensure this, calculating the angles between the corresponding planes is my best bet!

Michele

1 Like

Yes. I’m referring to Principal Component Analysis

2 Likes

If you have a linear transformation matrix which takes the position of one bone block to the other (as can be produced as an output from any Slicer registration technique), you can easily extract the rotation component of that transform using a module I developed for Slicer called CharacterizeTransformMatrix (available at GitHub - mikebind/SlicerCharacterizeTransformMatrix). I’m not sure what you want to use the new coordinate system for, but if you just want to know how rotated one block is relative to the other, this module will answer that.

@matsuba8,

Give PCA a try. I am not familiar with the process myself but I think this will ensure your analysis process carries less user error. Also, PCA seems like the more scalable solution.

If you still want to do something with planes, let me know!

1 Like

Your module looks good and I think it will be very useful for lots of Slicer users.

Although I didn’t catch how you calculate the scaleMatrices from the eigenvectors and eigenValues here: https://github.com/mikebind/SlicerCharacterizeTransformMatrix/blob/bde60887b7291c5ffd337e382f469ed5aadaf569/CharacterizeTransformMatrix/CharacterizeTransformMatrix.py#L323-L325

I suppose you have some formula to create a symmetric matrix but it’s just a guess.

1 Like

Yes, each scale matrix is symmetric, and is constructed by computing the outer product of the eigenvector with itself (which is symmetric by construction), multiplying that by the (eigenvalue - 1) and adding the result to the identity matrix. Each scale matrix has the effect of stretching vectors by a factor of the eigenvalue in the direction of the eigenvector (for the eigenvalue and eigenvector used to construct it).

1 Like

Thanks. That was a good explanation. I didn’t know that formula

1 Like

You can define a coordinate system using a markups plane node:

The advantage of using this instead of a transform that you can conveniently place the node using 3 points and then further adjust its position and orientation using interaction handles.

Once you have defined the plane, you can get the the transform between this coordinate system and the RAS coordinate system like this:

planeToWorld = vtk.vtkMatrix4x4()
getNode('P').GetObjectToWorldMatrix(planeToWorld)

If you concatenate this transform with the result of the model registration then you’ll get the rotation and translation between the bones in the plane coordinate system.

1 Like

Hello Andras,

Thank you very much for your response!
I just had a question while I was attempting this method: Is there a code to place a fiducial at a specific coordinate?
For example, as seen in the image below, the centroid of the red segment is at (-110.676, -105.632, 31.2008), and I would like to set my ‘P’ fiducial at this point.

Thank you for your help in advance.
Best,

Michele

You can use markupsNode.SetNthControlPointPositionWorld(i, r, a, s) method.