VTK Tutorial 13: Fiducial Localization Error with Bias - Introduction to Computer Assisted Surgeries

In this ( $13^{th}$ ) tutorial, we are going to investigate how Fiducial Localization Error (FLE) with bias influence the rigid-body registration accuracy.

More examples and tutorials can be found here

Fiducial Localization Errors with bias¶

[We previously discussed how bias may be introduced into the Fiducial Localization Error (FLE)](#ref_FLE_with_bias). Here, we are going visualize this process.

Python Setup¶

Most of the methods for manipulating this transformation, e.g. Translate, Rotate, and Concatenate, can operate in either PreMultiply (the default) or PostMultiply mode. In PreMultiply mode, the translation, concatenation, etc. will occur before any transformations which are represented by the current matrix. In PostMultiply mode, the additional transformation will occur after any transformations represented by the current matrix.

The presentation we use adheres to the PostMultiply mode.

To begin with, we need to import all the necessary Python functions and external libraries:

import math
import numpy as np

from vtkmodules.vtkCommonTransforms import vtkTransform
from vtkmodules.vtkCommonTransforms import vtkLandmarkTransform
from vtkmodules.vtkCommonCore import vtkPoints

Geometry¶

Using this previous example of optical tracking as the context, let’s define 4 registration fiducial points (i.e. the location of the retroreflective spheres) and a single target fiducial point (i.e. tip of the stylus):

Marker	x (mm)	y (mm)	z (mm)
A	0.00	0.00	0.00
B	0.00	0.00	50.00
C	0.00	25.00	100.00
D	0.00	-25.00	135.00
tip	10.00	0.00	-150.00

RegistrationFiducial = vtkPoints()
RegistrationFiducial.SetDataTypeToDouble()

RegistrationFiducial.InsertNextPoint( 0, 0, 0 )
RegistrationFiducial.InsertNextPoint( 0, 0, 50 )
RegistrationFiducial.InsertNextPoint( 0, 25, 100 )
RegistrationFiducial.InsertNextPoint( 0, -25, 135 )

TargetFiducial = vtkPoints()
TargetFiducial.InsertNextPoint( 10, 0, -150 )

Define a transformation. For the purpose of this tutorial, the exact transformation does not matter as long as it is not identity:

alpha = 60 # 60 degree
beta = 45
transform = vtkTransform()
transform.PostMultiply()
transform.Identity()
transform.RotateX( alpha )
transform.RotateY( beta )
transform.Translate( 1, 2, 3)
transform.Update()
print( "\nThe ground truth transformation is:\n", transform.GetMatrix() )

Perform the transformation to both the registration and target fiducial points:

transformedRegistrationFiducial = vtkPoints()
transformedTargetFiducial = vtkPoints()

transform.TransformPoints( RegistrationFiducial, transformedRegistrationFiducial )
transform.TransformPoints( TargetFiducial, transformedTargetFiducial )

Ground Truth¶

Now we have established a ground truth, i.e. we explicitly specified a known transformation, and applied the model points with the transformation without any errors.

Let’s print these ground truth on the screen.

print( "The model registration fiducial points are: \n")
for i in range(RegistrationFiducial.GetNumberOfPoints() ):
    print( RegistrationFiducial.GetPoint(i) )
print( "\nThe error-free transformed fiducial points are: \n" )
for i in range( transformedRegistrationFiducial.GetNumberOfPoints() ):
    print(transformedRegistrationFiducial.GetPoint(i))

print( "\nThe model target fiducial point is: \n" )
print( TargetFiducial.GetPoint(0) )
print( "\nThe transformed target fiducial point is: \n" )
print(transformedTargetFiducial.GetPoint(0))

print( "These are the ground truth!")

Transform via OPA¶

Given the model and measured fiducial point sets, can we recover the underlying transformation?

OPA = vtkLandmarkTransform()
OPA.SetSourceLandmarks( RegistrationFiducial )
OPA.SetTargetLandmarks( transformedRegistrationFiducial )
OPA.SetModeToRigidBody()
OPA.Update()
print( OPA.GetMatrix() )

Indeed we can recover the underlying transformation using vtkLandmarkTransform!

Introducing Fiducial Localization Errors¶

Let’s contaminate the measured point fiducial with some FLE. Again, for the purpose of this tutorial, the exact magnitude does not matter.

We will be manually introduce a small bias: the exact values does not matter.

measuredFiducial = vtkPoints()

FLE = np.array([1,2,3])

for i in range( transformedRegistrationFiducial.GetNumberOfPoints() ):
    point = transformedRegistrationFiducial.GetPoint( i )
    measuredFiducial.InsertNextPoint( point+FLE )

    
for i in range( measuredFiducial.GetNumberOfPoints() ):
    print( measuredFiducial.GetPoint(i) )
    print( transformedRegistrationFiducial.GetPoint(i), "\n" )

Transformation with FLE?¶

Now our measurements are contaminated with FLE, can we still recover the underlying transformation? If so, how do we evaluate its fitness?

OPA_FLE = vtkLandmarkTransform()
OPA_FLE.SetSourceLandmarks( RegistrationFiducial )
OPA_FLE.SetTargetLandmarks( measuredFiducial )
OPA_FLE.SetModeToRigidBody()
OPA_FLE.Update()
print( OPA_FLE.GetMatrix() )

Calculate FRE¶

transformedRegistrationFiducial_with_FLE = vtkPoints()

OPA_FLE.TransformPoints( RegistrationFiducial, transformedRegistrationFiducial_with_FLE )

FRE = 0
for i in range( transformedRegistrationFiducial_with_FLE.GetNumberOfPoints() ):
    point = measuredFiducial.GetPoint(i)
    point_FLE = transformedRegistrationFiducial_with_FLE.GetPoint( i )
    FRE = FRE + np.linalg.norm( np.array(point) - np.array(point_FLE) )

print( "FRE with bias: ", FRE )

Target Registration Error¶

We now see that the transformation we derived using the contaminated measurements is different from the ground truth, which is expected, but how good (or how bad) is this registration?

One way to quantify is via Target Registration Error (TRE).

transformedTargetFiducial_with_FLE = vtkPoints()

OPA_FLE.TransformPoints( TargetFiducial, transformedTargetFiducial_with_FLE )

print( "The TRE is: \n", 
      np.linalg.norm( np.array(transformedTargetFiducial.GetPoint(0)) - 
                     np.array(transformedTargetFiducial_with_FLE.GetPoint(0))))

Final Thoughts¶

You should perhaps wrap this calculation of target registration error in a loop and run it through, perhaps, thousands of times and calculate the mean TRE instead. Based on the Law of large numbers, what do you expect the mean TRE to be?

Notice the magnitude of the FRE, it is non-existant! What happens here?