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.
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.
2 Likes
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
I think this it is a common need to define a frame that follows some principal directions of the model. I’ve created a script for this. Also, it allows interaction to transform the model and give it different position/orientation.
Here is the code:
#Create a frame for the model in the principal components directions,
#set up a plane with interaction handles to move it.
#The script can be reused to allow interactions with n models
import numpy as np
#Change this number for more accuracy (more points, take more time to compute)
numberOfSampledPointsOfModel = 2000
model = getNode('mymodel')
maskPointsFilter = vtk.vtkMaskPoints()
maskPointsFilter.SetInputData(model.GetPolyData())
if model.GetPolyData().GetNumberOfPoints() > numberOfSampledPointsOfModel:
ratio = round(model.GetPolyData().GetNumberOfPoints()/numberOfSampledPointsOfModel)
else:
ratio = 1
#This works but the sampling could be biased spatially I think
#If you have vtk9 you should use UNIFORM_SPATIAL_SURFACE option
maskPointsFilter.SetOnRatio(ratio)
maskPointsFilter.RandomModeOn()
maskPointsFilter.Update()
polydata = vtk.vtkPolyData()
polydata.ShallowCopy(maskPointsFilter.GetOutput())
#Calculate the SVD and mean
from vtk.util.numpy_support import vtk_to_numpy
modelPoints = vtk_to_numpy(polydata.GetPoints().GetData())
# Calculate the mean of the points, i.e. the 'center' of the cloud
modelPointsMean = modelPoints.mean(axis=0)
# Do an SVD on the mean-centered data.
uu0, dd0, eigenvectors0 = np.linalg.svd(modelPoints - modelPointsMean)
# Create a frame for the model
modelZ = np.zeros(3)
modelX = eigenvectors0[0]
modelY = eigenvectors0[1]
vtk.vtkMath.Cross(modelX, modelY, modelZ)
modelZ = modelZ/np.linalg.norm(modelZ)
modelOrigin = modelPointsMean
#Create the plane to move the model
planeNode = slicer.mrmlScene.AddNewNodeByClass("vtkMRMLMarkupsPlaneNode", slicer.mrmlScene.GetUniqueNameByString("modelToWorldPlane"))
slicer.modules.markups.logic().AddNewDisplayNodeForMarkupsNode(planeNode)
displayNode = planeNode.GetDisplayNode()
displayNode.SetGlyphScale(2.5)
displayNode.SetOpacity(0)
displayNode.HandlesInteractiveOn()
planeNode.SetAxes(modelX,modelY,modelZ)
planeNode.SetOrigin(modelOrigin)
#Create the transformNode that the model will observe
planeToWorldTransformNode = slicer.vtkMRMLLinearTransformNode()
planeToWorldTransformNode.SetName(slicer.mrmlScene.GetUniqueNameByString("planeToWorld"))
slicer.mrmlScene.AddNode(planeToWorldTransformNode)
#Save initial transform of the plane
worldToInitialModelTransformNode = slicer.vtkMRMLLinearTransformNode()
worldToInitialModelTransformNode.SetName(slicer.mrmlScene.GetUniqueNameByString("worldToInitialModel"))
slicer.mrmlScene.AddNode(worldToInitialModelTransformNode)
SLICER_CHANGE_OF_API_REVISION = '29927'
planeToWorldMatrix = vtk.vtkMatrix4x4()
if int(slicer.app.revision) > int(SLICER_CHANGE_OF_API_REVISION):
planeNode.GetObjectToWorldMatrix(planeToWorldMatrix)
else:
planeNode.GetPlaneToWorldMatrix(planeToWorldMatrix)
planeToWorldTransformNode.SetMatrixTransformToParent(planeToWorldMatrix)
planeNode.SetAttribute('planeNodeToWorldTransformNodeID',planeToWorldTransformNode.GetID())
worldToInitialModelMatrix = vtk.vtkMatrix4x4()
worldToInitialModelMatrix.DeepCopy(planeToWorldMatrix)
worldToInitialModelMatrix.Invert()
worldToInitialModelTransformNode.SetMatrixTransformToParent(worldToInitialModelMatrix)
planeNode.SetAttribute('worldToInitialModelTransformNodeID',worldToInitialModelTransformNode.GetID())
model.SetAndObserveTransformNodeID(worldToInitialModelTransformNode.GetID())
worldToInitialModelTransformNode.SetAndObserveTransformNodeID(planeToWorldTransformNode.GetID())
###THIS PART SHOULD ONLY BE COPYPASTED ONCE TILL MARKED
#update model's position/orientation interactively
def onPlaneModified(sourceNode,event):
SLICER_CHANGE_OF_API_REVISION = '29927'
planeToWorldMatrix = vtk.vtkMatrix4x4()
if int(slicer.app.revision) > int(SLICER_CHANGE_OF_API_REVISION):
sourceNode.GetObjectToWorldMatrix(planeToWorldMatrix)
else:
sourceNode.GetPlaneToWorldMatrix(planeToWorldMatrix)
planeNodeToWorldTransformNode = slicer.mrmlScene.GetNodeByID(sourceNode.GetAttribute('planeNodeToWorldTransformNodeID'))
planeNodeToWorldTransformNode.SetMatrixTransformToParent(planeToWorldMatrix)
planeNodesAndObservers = []
###TILL HERE SHOULD ONLY BE COPYPASTED ONCE
planeNodesAndObservers.append(
[planeNode,planeNode.AddObserver(slicer.vtkMRMLMarkupsNode.PointModifiedEvent,onPlaneModified)]
)
#If you want to stop interactions for plane i, do this:
i = 0#change this number
planeNodesAndObservers[i][0].RemoveObserver(planeNodesAndObservers[i][1])
#If you want the model to world transform AFTER YOU MOVE IT for plane i, do this:
i = 0#change this number
worldToInitialModelTransformNode = slicer.mrmlScene.GetNodeByID(planeNodesAndObservers[i][0].GetAttribute('worldToInitialModelTransformNodeID'))
modelToWorldMatrix = vtk.vtkMatrix4x4()
worldToInitialModelTransformNode.GetMatrixTransformToWorld(modelToWorldMatrix)
print(modelToWorldMatrix)
#If you just want to define a frame for the model (assigned to plane i), do this:
i = 0#change this number
worldToInitialModelTransformNode = slicer.mrmlScene.GetNodeByID(planeNodesAndObservers[i][0].GetAttribute('worldToInitialModelTransformNodeID'))
initialModelToWorldMatrix = vtk.vtkMatrix4x4()
worldToInitialModelTransformNode.GetMatrixTransformToParent(initialModelToWorldMatrix)
initialModelToWorldMatrix.Invert()
print(initialModelToWorldMatrix)
@lassoan maybe this can be added to the script repository
Thank you @mau_igna_06 this is example works very nicely. It would be probably too long for including it entirely in the script repository but you can upload it into a gist and we can link it from the script repository (as we do it for other longer examples). Gists are version controlled, so you can update the gist and the same link will automatically post to the latest version (and user has access to previous versions of the gist, too).
There are also a few things to fix/improve, such as:
- Node IDs are not persistent (when you save the scene and load it again a different ID may be assigned to the node), therefore they must never stored in node attributes. Instead, node references can be used, which maintain link between nodes, even if the node IDs change. Since parent transform is already stored using node references, you don’t even need to add another node reference.
- Simplify the sampling. If ratio=1 then just use the model output and skip the filter. It would be also more relevant to use the VTK9 filter, as VTK8 will not be around for long.
- Always leave a space after the
#character - Instead of hardcoding any specific revision, put the method call in a try/except block (put the more recent API in the main path and the old API in the except path), as it is more pythonic and may be faster (not that speed would matter here, just as a general note)
- Instead of updating the transform in the main path and then implementing it again in onPlaneModified function, just call onPlaneModified once.
- I would remove the multiple planes support and the
ivariable to keep things simple. - Remove the
THIS PART SHOULD ONLY BE COPYPASTED ONCE TILL MARKEDcomments. It is understood that the entire code snippet is to be executed once. It would just confuse users that they are not allowed to copy-paste some parts. - Comment out the …RemoveObserver… line so that if the user copy-pastes the entire code snippet the axes are kept interactive
Two notes for using principal axes of surface mesh points as object coordinate system:
-
For solid objects Segment Statistics module gives more relevant axis directions. Principal X, Y, Z axes computed by Segment Statistics provides the object’s true principal axes of inertia, based on its volume. PCA of surface mesh points can be quite similar for many shapes, but they are not the same. Even if you sample the surface uniformly, a small feature with a large surface area can throw off the axis directions computed by surface mesh points PCA.
-
For models that have flat sides (such as CAD models) Segment Statistics module’s oriented bounding box axis directions are more relevant than principal axes. The oriented bounding box axes will not be slightly off if you have small attachments or engravings in the object, but the axes will be generally aligned with flat sides.
All changes done except this one:
It would be also more relevant to use the VTK9 filter, as VTK8 will not be around for long.
I’ll update the gist when Slicer changes VTK version.
Thanks for the notes.
Here is the gist: Create a frame for the model in the principal components directions, set up a plane with interaction handles to move it. · GitHub
Do you know if the option UNIFORM_SPATIAL_VOLUME for the maskPointsFilter will make it return points inside the volume of the mesh? We could try that for this script later, maybe it would give more robust results…
I tried to use vtkPoissonDiskSampler to update the script that creates a frame for a model using the principal directions of the surface but I couldn’t find the filter on latest preview release of Slicer.
maskFilter = vtk.vtkPoissonDiskSampler() gives this exception: AttributeError: module 'vtkmodules.all' has no attribute 'vtkPoissonDiskSampler'
The vtkPoissonDiskSampler class was added earlier this year and is not included in SLicer’s VTK yet. This issue tracks the next VTK update.
However, uniform sampling will probably not make much difference, unless your surface has huge variations in surface area of cells and you don’t use many samples.
For significant change you could use volumetric sampling (to get true mechanical moments) or oriented bounding box (so that small indentations or extrusions do not affect the axis directions in objects with flat sides).
Hi,
I have been trying to “define a coordinate system using a markups plane node” as suggested by @Iassoan. I have unfortunately no knowledge in coding and have not used the python interactor before.
I was wondering if there is any way to generate the plane linked to the 3 defining points without using the python interactor?
Is there a module that would do this?
I mainly use 3D Slicer to plan screw positioning in vertebrae and the first step of the workflow is to manually align the bone model to the RAS coordinate system using a transform. It would be much easier to place 3 points defining the sagittal plane assuming I can then define my screw positions in that coordinate system.
Many thanks,
Guillaume
Try this:
PlaneNode = getNode('P')
PlaneNode.SetPlaneType(PlaneNode.PlaneTypePlaneFit)
And restrict yourself to using 3 points while creating the plane.
1 Like
Thanks @mau_igna_06
I am afraid I am not sure I know how to.
Is there a basic tutorial on the steps to run a script like this?
I don’t want to waste your time as I have no experience at all in any scripting/coding so maybe I need to start somewhere else?
I tried to copy and paste these 2 lines in the python interactor using an empty Slicer scene, press enter and obtained a few error lines:
PlaneNode = getNode(‘P’)
Traceback (most recent call last):
File “”, line 1, in
File “/Applications/Slicer.app/Contents/bin/Python/slicer/util.py”, line 1312, in getNode
raise MRMLNodeNotFoundException(“could not find nodes in the scene by name or id ‘%s’” % (pattern if (isinstance(pattern, str)) else “”))
slicer.util.MRMLNodeNotFoundException: could not find nodes in the scene by name or id ‘P’
PlaneNode.SetPlaneType(PlaneNode.PlaneTypePlaneFit)
Traceback (most recent call last):
File “”, line 1, in
NameError: name ‘PlaneNode’ is not defined
I suspect I need a plane object in my scene maybe, which should be named P? Or is the P node the Fiducial list which is F by default in my 3D Slicer version? Do I need 3 points in the scene before I run the code or after? Will the plane appear as I place the 3 points? Is there any other prerequisite in the the scene before the script is run? Like a 3D object or a CT study?
Very sorry if these questions are extremely basic. I completely understand if this scripting module is meant for people with basic programming skills and maybe I should stick to the modules already created.
Guillaume
I suspect I need a plane object in my scene maybe, which should be named P?
You have to do that, yes.
Do I need 3 points in the scene before I run the code or after?
You should have drawn a plane before running the code I think
Is there any other prerequisite in the the scene before the script is run? Like a 3D object or a CT study?
No
planeToWorld = vtk.vtkMatrix4x4()
getNode('P').GetObjectToWorldMatrix(planeToWorld)
This part does not appear to work for me, Plane node is named “P”, however i cant get the matrix.
Has something changed or am i doing something wrong?
thank you
I’ve just tested this in the current stable release (5.0.3) and it works well. You can print matrix values:
>>> planeToWorld = vtk.vtkMatrix4x4()
>>> getNode('P').GetObjectToWorldMatrix(planeToWorld)
>>> print(planeToWorld)
vtkMatrix4x4 (000001D1F157F7B0)
Debug: Off
Modified Time: 4320154
Reference Count: 1
Registered Events: (none)
Elements:
1 0 0 87.172
0 1 0 -74.5663
0 0 1 -75.25
0 0 0 1
1 Like
Thank you so much, i still make a lot of small misstakes so i truly appreciate you helping me. I dont know how i did it but i messed up my print command.
Last question, what would be the best way to measure the model registration from point 1? could i transform the model to the world coordinate system using the inverse of the matrix i now have to allign my model with the world coordinate system?
again, thank you.