VTK Tutorial 8: Intersection between a line and an ellipsoid - Introduction to Computer Assisted Surgeries

In this ( $8^{th}$ ) tutorial, we are going to compute the intersection between a line and an ellipsoid, and visualize using VTK.

More examples and tutorials can be found here

Line-Ellipsoid Intersection¶

We previously discussed how to compute the intersection between a line and an ellipsoid. Here, we are going visualize this process.

Python Setup¶

To utilize VTK built-in function of an ellipsoid, add the following

from vtkmodules.vtkCommonComputationalGeometry import vtkParametricEllipsoid
from vtkmodules.vtkFiltersSources import vtkParametricFunctionSource

import math
import numpy as np

# First access the VTK module (and any other needed modules) by importing them.
# noinspection PyUnresolvedReferences
import vtkmodules.vtkRenderingOpenGL2
from vtkmodules.vtkCommonColor import vtkNamedColors
from vtkmodules.vtkCommonTransforms import vtkTransform
from vtkmodules.vtkFiltersSources import vtkSphereSource
from vtkmodules.vtkFiltersSources import vtkLineSource
from vtkmodules.vtkFiltersCore import vtkTubeFilter
from vtkmodules.vtkCommonComputationalGeometry import vtkParametricEllipsoid
from vtkmodules.vtkFiltersSources import vtkParametricFunctionSource
from vtkmodules.vtkRenderingAnnotation import vtkAxesActor
from vtkmodules.vtkInteractionStyle import vtkInteractorStyleTrackballCamera
from vtkmodules.vtkInteractionWidgets import vtkBoxWidget
from vtkmodules.vtkRenderingCore import (
    vtkActor,
    vtkPolyDataMapper,
    vtkRenderWindow,
    vtkRenderWindowInteractor,
    vtkRenderer
)

Geometry¶

An ellipsoid is specified by its center, and the magnitude of the principal axes.

colors = vtkNamedColors()

# Geometry of an ellipsoid, specified by its centre and the length of the principal axes
centre = np.array( [0, 4, 0] )
radius = np.array( [3, 4, 5] )
ellipsoid = vtkParametricEllipsoid()
ellipsoid.SetXRadius( radius[0] )
ellipsoid.SetYRadius( radius[1] )
ellipsoid.SetZRadius( radius[2] )

# in VTK, an ellipsoid is assumed to be located at the origin. We'll need to find a way to move it later

# generate polydata
ellipsoidSource = vtkParametricFunctionSource()
ellipsoidSource.SetParametricFunction( ellipsoid )
ellipsoidSource.Update()

# mappers and actors
ellipsoidMapper = vtkPolyDataMapper()
ellipsoidMapper.SetInputConnection( ellipsoidSource.GetOutputPort() )

ellipsoidActor = vtkActor()
ellipsoidActor.SetMapper( ellipsoidMapper )
ellipsoidActor.GetProperty().SetColor( colors.GetColor3d( 'Red') )
ellipsoidActor.SetPosition( centre.tolist() )   # this is how we are going to set the centre of an ellipsoid

A line segment is specified by its 2 end points. Feel free to change these locations.

Note that for our numerical solution, we need to compute the direction (orientation) of the line, which is represented as a unit vector.

line = vtkLineSource()
point1 = np.array([ 0, 4, -10])
point2 = np.array([ 0, 4, 10])
v = point2 - point1
v = v/np.linalg.norm( v )

line.SetPoint1( point1.tolist() )
line.SetPoint2( point2.tolist() )

tube = vtkTubeFilter()
tube.SetInputConnection( line.GetOutputPort() )
tube.SetRadius( 0.1 )
tube.SetNumberOfSides(10)

lineMapper = vtkPolyDataMapper()
lineMapper.SetInputConnection( tube.GetOutputPort() )
 
lineActor = vtkActor()
lineActor.SetMapper( lineMapper )

The quadratic equation that allows us to solve for the intersection. Since this is a quadratic equation, if the term inside the square-root (i.e. determinant) is:

positive, there are 2 solutions (i.e. 2 intersections),
zero, there is exactly 1 solution (i.e. the line is tangent to the sphere), or
negative, there are 0 solution (i.e. there is no intersection)

# compute the coefficient of the quadratic equation
#
# a t^2 + b t + c = 0

a = ( ( v[0]/radius[0] )**2 + 
      ( v[1]/radius[1] )**2 + 
      ( v[2]/radius[2] )**2 )
b = 2 * ( (point1[0] - centre[0]) * v[0] / (radius[0]**2) +
          (point1[1] - centre[1]) * v[1] / (radius[1]**2) +
          (point1[2] - centre[2]) * v[2] / (radius[2]**2) )
c = ( ( (point1[0] - centre[0])**2/(radius[0])**2 ) +
      ( (point1[1] - centre[1])**2/(radius[1])**2 ) + 
      ( (point1[2] - centre[2])**2/(radius[2])**2 ) - 1 )
determinant = b*b - 4*a*c

## calculate the intersection

if determinant > 0:                     # 2 intersections
    t1 = (-1 * b + math.sqrt( determinant ) ) / (2*a)
    L1 = point1 + t1 * v
    t2 = (-1 * b - math.sqrt( determinant ) ) / (2*a)
    L2 = point1 + t2 * v
elif determinant == 0:                  # 1 intersection
    L1 = point1 + (-1*b/a) * v

In terms of visualization, we need to either 1 or 2 (or 0) actors to the scene, depending on what this determinant is:

if determinant > 0:                         # 2 intersections
    s1 = vtkSphereSource() 
    s1.SetCenter( L1.tolist() )
    s1.SetRadius( 0.4 )
    
    s1Mapper = vtkPolyDataMapper()
    s1Mapper.SetInputConnection( s1.GetOutputPort())
    
    s1Actor = vtkActor()
    s1Actor.SetMapper( s1Mapper )
    s1Actor.GetProperty().SetColor( colors.GetColor3d('Yellow'))
    
    s2 = vtkSphereSource()
    s2.SetCenter( L2.tolist() )
    s2.SetRadius( 0.4 )
    
    s2Mapper = vtkPolyDataMapper()
    s2Mapper.SetInputConnection( s2.GetOutputPort())
    
    s2Actor = vtkActor()
    s2Actor.SetMapper( s2Mapper )
    s2Actor.GetProperty().SetColor( colors.GetColor3d('Yellow'))
elif determinant == 0:                      # 1 intersection
    s1 = vtkSphereSource() 
    s1.SetCenter( L1.tolist() )
    s1.SetRadius( 0.4 )
    
    s1Mapper = vtkPolyDataMapper()
    s1Mapper.SetInputConnection( s1.GetOutputPort())
    
    s1Actor = vtkActor()
    s1Actor.SetMapper( s1Mapper )
    
# if the determinant is 0, there is no intersection
#
# do NOTHING

Add an Axes to help us orient.

# add an 3D Axes
axes = vtkAxesActor()
axes.SetTotalLength( 10, 10, 10 )
axes.AxisLabelsOff()

# properties of the axes labels can be set as follows
# this sets the x axis label to red
axes.GetXAxisCaptionActor2D().GetCaptionTextProperty().SetColor(colors.GetColor3d('Red'));
axes.GetYAxisCaptionActor2D().GetCaptionTextProperty().SetColor(colors.GetColor3d('Green'));
axes.GetZAxisCaptionActor2D().GetCaptionTextProperty().SetColor(colors.GetColor3d('Blue'));

Rendering

ren = vtkRenderer()
ren.AddActor( ellipsoidActor )
ren.AddActor( lineActor )
ren.AddActor( axes )

if determinant > 0:
    ren.AddActor( s1Actor )
    ren.AddActor( s2Actor )
elif determinant == 0:
    ren.AddActor( s1Actor )

ren.SetBackground( colors.GetColor3d('MidnightBlue'))

renWin = vtkRenderWindow()
renWin.AddRenderer(ren)
renWin.SetSize(640, 480)
renWin.SetWindowName('AISE 4025, Intersection between a line and an ellipsoid')

iren = vtkRenderWindowInteractor()
iren.SetRenderWindow(renWin)

Rendering and user interaction:

iren.Initialize()
iren.Start()

Figure 1:Intersection between a line and an ellipsoid.