Frankfurt Horizontal Plane alignment of one model, angle measurements

[esomjai (Eszter Maja Somjai)]

Hi there,
I am currently starting a masters project in facial approximation with a landmark-heavy plan.
I apologise if my questions are answered elsewhere, I truly am new to 3D Slicer.

1. I would like to establish the Frankfort horizontal plane (anatomical plane based on 4 landmarks) as my midline for the models (see example pictures ) - they do not need to be aligned, their positions just need to be adjusted according to the plane. I established the plane before, having a hard time adjusting the view.
   

   
   
   

2. This plane is the basis of my angle measurements - I have lines in between landmarks crossing the plane and I will need the coordinates of these to be precise when measuring angles in relation to the plane. I am using this study: Guyomarc’h, P. and Stephan, C.N. (2012), The Validity of Ear Prediction Guidelines Used in Facial Approximation. Journal of Forensic Sciences, 57: 1427-1441. https://doi.org/10.1111/j.1556-4029.2012.02181.x

Thank you in advance for your time and assistance!

[muratmaga]

We use this tutorial as an example of scripting in Slicer. It uses three points to define the Frankfort, but it should be trivial to adjust it to use your four points:

W_2020/Lab09_Slicer_Scripting at master · SlicerMorph/W_2020 - Example 1: Aligning volumes using landmarks (Frankfort alignment plane).

[esomjai (Eszter Maja Somjai)]

Thank you so much for this! I have adjusted the code to suit version 5.2.2 and it works perfectly!
I decided to do the 3-landmark approach anyway.

in case anyone reading these posts is interested in this version:

import numpy
scene = slicer.mrmlScene
F = getNodesByClass('vtkMRMLMarkupsFiducialNode')
F = F[0]

# Get the fiducial IDs of porions and zygomatico - orbitale  suture from the fiducial list by name
po1_id = -1; po2_id = -1; zyo_id = -1;

for i in range(0, F.GetNumberOfControlPoints()):
	if F.GetNthControlPointLabel(i) == 'poR' :
		po1_id = i
	if F.GetNthControlPointLabel(i) == 'poL' :
		po2_id = i
	if F.GetNthControlPointLabel(i) == 'zyoL' :
		zyo_id = i


# Get the coordinates of landmarks
po1 = [0, 0, 0]
po2 = [0, 0, 0]
zyo = [0, 0, 0]

F.GetNthControlPointPosition(po1_id,po1)
F.GetNthControlPointPosition(po2_id,po2)
F.GetNthControlPointPosition(zyo_id,zyo)

# The vector between two porions that we will align to LR axis by calculating the yaw angle
po =[po1[0] - po2[0], po1[1] -po2[1], po1[2]-po2[2]]
vTransform = vtk.vtkTransform()
vTransform.RotateZ(-numpy.arctan2(po[1], po[0])*180/numpy.pi)
transform = slicer.mrmlScene.AddNewNodeByClass('vtkMRMLLinearTransformNode')
transform.SetMatrixTransformToParent(vTransform.GetMatrix())

# Apply the transform to the fiducials and the volume
transform = slicer.vtkMRMLLinearTransformNode()

scene.AddNode(transform) 
transform.SetMatrixTransformToParent(vTransform.GetMatrix())
V = getNodesByClass('vtkMRMLScalarVolumeNode')
V = V[0]
V.SetAndObserveTransformNodeID(transform.GetID())
F.SetAndObserveTransformNodeID(transform.GetID())

# Get the updated (transformed) coordinates from the list
po1 = [0, 0, 0]
po2 = [0, 0, 0]
zyo = [0, 0, 0]

F.GetNthControlPointPosition(po1_id,po1)
F.GetNthControlPointPosition(po2_id,po2)
F.GetNthControlPointPosition(zyo_id,zyo)

# Apply the transform to the coordinates
po1 = vTransform.GetMatrix().MultiplyPoint([po1[0], po1[1], po1[2], 0])
po2 = vTransform.GetMatrix().MultiplyPoint([po2[0], po2[1], po2[2], 0])
zyo = vTransform.GetMatrix().MultiplyPoint([zyo[0], zyo[1], zyo[2], 0])
po =[po1[0]-po2[0], po1[1]-po2[1], po1[2]-po2[2]]

# Calculate the rotation for the roll 
vTransform2 = vtk.vtkTransform()

vTransform2.RotateY(numpy.arctan2(po[2], po[0])*180/numpy.pi)

# Apply the transform to previous transform hierarchy
transform2 = slicer.vtkMRMLLinearTransformNode()
scene.AddNode(transform2) 
transform2.SetMatrixTransformToParent(vTransform2.GetMatrix())
transform.SetAndObserveTransformNodeID(transform2.GetID())

# Get the coordinates again
po1 = [0, 0, 0]
po2 = [0, 0, 0]
zyo = [0, 0, 0]

F.GetNthControlPointPosition(po1_id,po1)
F.GetNthControlPointPosition(po2_id,po2)
F.GetNthControlPointPosition(zyo_id,zyo)

# Apply transforms to points to get updated coordinates
po1 = vTransform.GetMatrix().MultiplyPoint([po1[0], po1[1], po1[2], 0])
po2 = vTransform.GetMatrix().MultiplyPoint([po2[0], po2[1], po2[2], 0])
zyo = vTransform.GetMatrix().MultiplyPoint([zyo[0], zyo[1], zyo[2], 0])
po1 = vTransform2.GetMatrix().MultiplyPoint([po1[0], po1[1], po1[2], 0])
po2 = vTransform2.GetMatrix().MultiplyPoint([po2[0], po2[1], po2[2], 0])
zyo = vTransform2.GetMatrix().MultiplyPoint([zyo[0], zyo[1], zyo[2], 0])

# The vector for pitch angle
po_zyo = [zyo[0]-(po1[0]+po2[0])/2, zyo[1]-(po1[1]+po2[1])/2, zyo[2]-(po1[2]+po2[2])/2]

# Calculate the transform for aligning pitch
vTransform3 = vtk.vtkTransform()
vTransform3.RotateX(-numpy.arctan2(po_zyo[2], po_zyo[1])*180/numpy.pi)

# Apply the transform to both fiducial list and the volume
transform3 = slicer.vtkMRMLLinearTransformNode()
scene.AddNode(transform3) 
transform3.SetMatrixTransformToParent(vTransform3.GetMatrix())
transform2.SetAndObserveTransformNodeID(transform3.GetID())

slicer.vtkSlicerTransformLogic().hardenTransform(V)

Now that part 1 of my loaded question is solved, the second part remains - How do I get the coordinates of the intersection point where my measurement line crosses the plane? I could do this with approximation but that would be imprecise…

Thank you again for your time and assistance!

[lili-yu22]

can you help me?
I don’t know pathon please give me script to produce a sagittal plane through N point and vertical to FHplane.

[esomjai (Eszter Maja Somjai)]

Dear lili-yu22,

I think you’d need a different approach for that…

I would not do this by python scripting if you already have your FHPlane in the environment - I’d do the following:
1.

1. Duplicate the FHP ONCE (right click on Markups list > Clone)

• Rename it to Sag

Time to rotate the new planes

• Open Transforms

• Click Sag, then the >○ arrow

1. Adjust Rotation LR to 90 degrees

2. Click again on Sag and on ‘Harden Transform’

3. Adjust the planes to meet the dots described

4. Toggle the visibility of the point you need the plane to meet (N)

5. Then enable under Display¬ Interaction Visibility enable translation for all and manually toggle until they cross the landmarks described above

6. If they don’t move, unlock their control points!1. I hope this helps

Best wishes,
Eszter

[lili-yu22]

I copy the script ，but I also want the segment of the skull adjust as the volume 。so I modify the script as the picture ， but after I usethe modified scipt， the volume and the segment of the skullnot overlap，please help me​:rose:

[esomjai (Eszter Maja Somjai)]

Hi,
I believe the issue is in the order in which you execute the steps.
I re-orientate my scan (script) BEFORE any segmentation takes place - this way, they will be oriented the same.
I hope this helps.

so

1. Open the DICOM file,
2. Put the orientating landmarks for FHP on the image
3. Run the script
4. Create an segmentations needed

[lili-yu22]

thank you very much