Compute bone axis in 3D Slicer by specifying two circles.

BoneAxis.py

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