# Creating a new coordinate system

**URL:** <https://discourse.slicer.org/t/creating-a-new-coordinate-system/19446>\
**Category:** Support\
**Tags:** transforms, coordinates\
**Created:** [August 31, 2021, 6:27pm UTC](https://discourse.slicer.org/t/creating-a-new-coordinate-system/19446 "2021-08-31T18:27:37Z")\
**Posts on this page:** 1\
**Showing post:** 19

<div class="post-metadata">

**Author:** ![mau\_igna\_06](https://sea2.discourse-cdn.com/flex002/user_avatar/discourse.slicer.org/mau_igna_06/32/9056_2.png) [@mau\_igna\_06](https://discourse.slicer.org/u/mau_igna_06)\
**Post date:** [November 5, 2021, 12:58pm UTC](https://discourse.slicer.org/t/creating-a-new-coordinate-system/19446/19 "2021-11-05T12:58:08Z")

</div>

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.

 ![frameOfTibiaInPCDirections](https://us1.discourse-cdn.com/flex002/uploads/slicer/original/3X/e/a/ea5a68f1130fc2f4ec8c1344442ecd571e4a3861.png)

Here is the code:

```auto
#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

---

_[View the full topic](https://discourse.slicer.org/t/creating-a-new-coordinate-system/19446)._