Skip to content

Instantly share code, notes, and snippets.

@lassoan
Last active October 31, 2023 17:40
Show Gist options
  • Save lassoan/99514fb5edf98e1e3dda926253dc9742 to your computer and use it in GitHub Desktop.
Save lassoan/99514fb5edf98e1e3dda926253dc9742 to your computer and use it in GitHub Desktop.
Compute bone axis in 3D Slicer by specifying two circles.
# Create a markups fiducial node and place 6 points then copy paste this code into 3D Slicer's Python console.
# Two circles will be created (one for each 3 points) and the axis will be computed by connecting their centers.
# The circles and the axis is updated as the fiducial point positions are adjusted.
def sphereFrom3Points(markupsNode, startPointIndex):
"""Compute center and radius of 3-point sphere from 3 fiducial points
source: https://stackoverflow.com/questions/20314306/find-arc-circle-equation-given-three-points-in-space-3d
"""
import numpy as np
A = np.zeros(3)
B = np.zeros(3)
C = np.zeros(3)
markupsNode.GetNthFiducialPosition(startPointIndex,A)
markupsNode.GetNthFiducialPosition(startPointIndex+1,B)
markupsNode.GetNthFiducialPosition(startPointIndex+2,C)
a = np.linalg.norm(C - B)
b = np.linalg.norm(C - A)
c = np.linalg.norm(B - A)
s = (a + b + c) / 2
R = a*b*c / 4 / np.sqrt(s * (s - a) * (s - b) * (s - c))
b1 = a*a * (b*b + c*c - a*a)
b2 = b*b * (a*a + c*c - b*b)
b3 = c*c * (a*a + b*b - c*c)
P = np.column_stack((A, B, C)).dot(np.hstack((b1, b2, b3)))
P /= b1 + b2 + b3
return P, R
def UpdateModels(param1=None, param2=None):
"""Update the sphere and line models from the fiducial points"""
import math
center0, radius0 = sphereFrom3Points(markups, 0)
center1, radius1 = sphereFrom3Points(markups, 3)
spheres[0].SetCenter(center0)
spheres[0].SetRadius(radius0)
spheres[1].SetCenter(center1)
spheres[1].SetRadius(radius1)
line.SetPoint1(center0+3*(center0-center1))
line.SetPoint2(center1+3*(center1-center0))
tube.Modified()
slicer.util.forceRenderAllViews()
# Get markup node from scene
markups = slicer.mrmlScene.GetFirstNodeByClass('vtkMRMLMarkupsFiducialNode')
# Create model nodes add to scene, and enable slice intersection display
modelsLogic = slicer.modules.models.logic()
# Spheres
spheres = [vtk.vtkSphereSource(), vtk.vtkSphereSource()]
for itemIndex, sphere in enumerate(spheres):
sphere.SetPhiResolution(30)
sphere.SetThetaResolution(30)
model = modelsLogic.AddModel(sphere.GetOutputPort())
model.GetDisplayNode().SetSliceIntersectionVisibility(True)
model.GetDisplayNode().SetSliceIntersectionThickness(3)
model.GetDisplayNode().SetColor(1,1-itemIndex*0.3,itemIndex*0.3)
# Line
line = vtk.vtkLineSource()
tube = vtk.vtkTubeFilter()
tube.SetInputConnection(line.GetOutputPort())
model = modelsLogic.AddModel(tube.GetOutputPort())
model.GetDisplayNode().SetSliceIntersectionVisibility(True)
model.GetDisplayNode().SetSliceIntersectionThickness(2)
model.GetDisplayNode().SetColor(1, 0.4, 0.6)
UpdateModels()
# Call UpdateSphere whenever the fiducials are changed
markups.AddObserver(slicer.vtkMRMLMarkupsNode.PointModifiedEvent, UpdateModels, 2)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment