Creating a new coordinate system

[mau_igna_06 (Mauro I. Dominguez)]

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…

[mau_igna_06 (Mauro I. Dominguez)]

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'

[lassoan (Andras Lasso)]

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).

[GL1984 (Guillaume Leblond)]

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

[mau_igna_06 (Mauro I. Dominguez)]

Try this:

PlaneNode = getNode('P')
PlaneNode.SetPlaneType(PlaneNode.PlaneTypePlaneFit)

And restrict yourself to using 3 points while creating the plane.

[GL1984 (Guillaume Leblond)]

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

[mau_igna_06 (Mauro I. Dominguez)]

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

[jegberink (jeroen egberink)]

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

[lassoan (Andras Lasso)]

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 

[jegberink (jeroen egberink)]

@lassoan

print(planeToWorld)

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.

[jegberink (jeroen egberink)]

I figured i might aswell try, took the inverse matrix and it worked. no more help needed for now, thank you!

[GL1984 (Guillaume Leblond)]

Apologies for not following up on this thread after guidance provided. I only have limited opportunities to work on this project and it has been several months since I have tried again.
The recent discussions motivated me to try again and I have made progress. The reason I could not figure out how to use mark-up plane nodes is that I was using an older version of Slicer (4.11) and I believe that function was not fully available then.
I am now able to to generate the plane with 3 points and could obtain the transformation matrix as suggested.
My residual problem is that transformation seem to cause a 90 degree rotation around the A axis when I try to align my model to the world coordinate. I am not sure I understand why? Maybe it has to do with the way the plane coordinates are defined?
I am attaching a screen shot so you can visualise what I am talking about (green model before transformation, red model post transformation).

Anyone knows what I am doing wrong?

[jegberink (jeroen egberink)]

hey!
You can see if this helps:

Defining a new coordinate system using markups plane node

The axes of a MarkupsPlane node are always orthogonal, even if the triangle of points used to define the plane is not a right triangle. If you have a MarkupsFiducial node with your three defining points, this function will produce a MarkupsPlane node def makePlaneMarkupFromFiducial(self, FNode, planeName): """Create MarkupsPlane using first three control points of the input fiducial node. """ if FNode.GetNumberOfControlPoints()<3: logging.warning('Not enough control points to…

Your point 1 to point 2 is your positive X-axis, which should correspond to the LR axis in Slicer.

[GL1984 (Guillaume Leblond)]

Thanks @jegberink
This explains well my problem. I am not sure how to resolve it though.
The way I want to use this coordinate system is by defining 3 points in the Sagittal plane so unfortunately I don’t have the option to use point 1 and 2 in the left to right direction. I want to use this to plan dog’s neurosurgeries so in the spine the easiest plane to define is sagittal (plane of symmetry). I could however use point 1 and 2 in either the inferior to superior or posterior to anterior direction.
Would there be a way to modify the code so I get a transform matrix in the proper orientation?
Otherwise I guess I would have to add a 90 degree rotation manually so I get the proper anatomical alignment.

[jegberink (jeroen egberink)]

I had the same problem, a 90 degree transform is the easiest solution. Create a linear transform and rotate it 90 degrees towards the the axis you prefer.
Do keep in mind that point 1 will remain in 0.

Good luck!

[GL1984 (Guillaume Leblond)]

Thanks,
The method is working now.
I would like to find a way to simplify it in the future but need to explore other steps in my workflow to decide the best way forward. I suspect setting up a module will be the end goal.

[lili-yu22]

I have doubts. I registered two objects in the world coordinate system and got a matrix 1. Then I can know how an object moves and rotates in the world coordinate system to register with another object. If I create a local coordinate system using the makeupplane node and want to know the movement in the local coordinate system, How do I deal with these two matrices next, as shown in the figure, is it right?

[lili-yu22]

can you tell me how to make the plane parallel to the world’s x or y or z axis using the script

[lassoan (Andras Lasso)]

You can create a plane with patient RAS axes if you set plane to world to identity matrix. Change order of columns to change what plane axis is aligned with which patient axis.

[lili-yu22]

@jegberink

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

first,thanks for your reply.
according to the above script:
planeToWorld = vtk.vtkMatrix4x4()
getNode(‘P’).GetObjectToWorldMatrix(planeToWorld)
I can get the matix，then I use the matrix on the plane.Do I understand correctly that ？