I am working on a project to quantify how closely a surface model generated from a CT scan matches the true bone geometry. This has been done quite extensively in literature but we have a need to quantify accuracy specifically with our tools and scan protocols etc.
The plan at the moment is to get a piece of animal bone and mount it on a small stand. We will then put the bone in a coordinate measuring machine and touch some points on the bone to get a low density point cloud of XYZ points which are exactly on the edge of the bone surface.
We will then put the bone in a container and fill it with gelatin to simulate soft tissue and put the container through the CT scanner to get an image set. This image set will then be used to generate an surface model of the bone with commercial segmentation software.
The next question is how to assess the accuracy of the surface model to the true surface of the bone. Thus far my thinking is to put all the XYZ points into the computer as fiducials, then load the STL model in as well like so. (pictures are just a practice example)
The fiducial to model registration module gives a readout of the average distance which is great but it would also be very handy to know the distance between each fiducial point and the closest point on the model surface say in a table or something.
Is there an easy way to do this or does it involve writing python scripts? I have found some scripts which look like they do similar things but not quite the same in the script repository.
you should first register the models. You have number of tools in 3D slicer,
IGT Fiducial registration wizard
CMFRegistration - Surface registration, ROI registration and Landmark registration.
For me the best results i have obtained was with CMFReg- ROI registration. But first try surface registration with Max iterations and max landmarks as it is the easiest method.
Next you need to apply the transformation to the moving model and then use Model to model distance module to calculate the distance. and you can export the values as a file too if you want to.
or you can use ShapePopulation viewer to look at the areas of difference the aid of a color map which would be best in this i case for visualization purposes.
The question is how do you know the closest point registration found on the model derived from CT is indeed the point you recorded with your coordinate reading machine.
If you go with sparse manual sampling, perhaps stick garnet sands (you can find them about 50-100 micron diameter), then record their coordinates with your digitizer, and then scan them. They will show up like bright spots in the CT scan.
Or you can use a well-calibrated microCT to obtain a dense mesh and use that as your reference and use model-to-model distance module.
Manjula - thanks for the suggestions, Unfortunatly model-model distance is not useful in this case. I would need a fiducial to model distance module but I havent seen that one exists.
Muratmaga - Yes good suggestion about the markers. I suppose I could put another set of fiducials on the markers as seen on the scan and register the two sets of fiducials. It is also simpler for me to measure the distance between two fiducials than between fiducials and a surface. Other studies use surface scans of the bone as the gold standard to compare the CT to which would also be easier to analyse using model-model distance, however the others I am working with are keen to use the method described above.
Wow thanks Andras, that is really useful. I was trying to figure it out on python last night but I am more or less a complete noob at programming. I am trying to learn though as more and more I am realising how useful python is in slicer. I can try and make a module it just might take me a while…
Thanks again. Do you have any pointers on where a complete novice could learn to write a script like the one you just made?
You need to be comfortable with Python language, which you can learn in a few days from tutorials on the web. Then the rest is about getting to know libraries. Slicer mainly relies on VTK (that you can learn from VTK textbook and examples), ITK (for image processing, that you can learn from the ITK software guide), and Qt (for GUI; which you can learn from Qt documentation and examples). Slicer adds a little extra layer and logic, which you can learn from studying Slicer programming tutorials, script repository, and source code of Slicer core and extensions.
The tutorial that I linked above is a good starting point and you can always ask questions here. You can also come to one of the Slicer project weeks where you bring your own problem that you work on with help from Slicer experts.
I have been trying to learn some python 3 using codeacademy.com which seems good.
I added some extra lines to your code to print out the min, max, mean and RMS values from the distances in the table.
I am not sure about the RMS calculation. I was looking at the code for the fiducial to model registration module in Slicer IGT here and trying to copy what the ComputeMeanDistance function does. However their closest point function must work differently to your one as your code seems to give a negative distance when the fiducial is inside the surface of the model (which is good), however it is not possible to take the square root of a negative number.
To get around this I have taken the absolute value of the numbers in the table but I cannot get the same result as is printed in the fiducial to model registration. I have tried applying the transfrom to the model, I have tried inverting the transform and applying it to the model and I have tried with no transform applied. I cannot get my RMS printout to match the value in the fiducials-model registration output box.
Do you know what I might be doing wrong? Here is my code. I don’t understand how to use the vtk arrays so I have made a parallel list to do operations on.
import statistics
import math
markupsNode = getNode('F')
modelNode = getNode('mymodel')
# Transform model polydata to world coordinate system
if modelNode.GetParentTransformNode():
transformModelToWorld = vtk.vtkGeneralTransform()
slicer.vtkMRMLTransformNode.GetTransformBetweenNodes(modelNode.GetParentTransformNode(), None, transformModelToWorld)
polyTransformToWorld = vtk.vtkTransformPolyDataFilter()
polyTransformToWorld.SetTransform(transformModelToWorld)
polyTransformToWorld.SetInputData(modelNode.GetPolyData())
polyTransformToWorld.Update()
surface_World = polyTransformToWorld.GetOutput()
else:
surface_World = modelNode.GetPolyData()
# Create arrays to store results
indexCol = vtk.vtkIntArray()
indexCol.SetName("Index")
labelCol = vtk.vtkStringArray()
labelCol.SetName("Name")
distanceCol = vtk.vtkDoubleArray()
distanceCol.SetName("Distance")
# Create blank list becuase I dont know how to use vtk arrays yet
dist_list = []
distanceFilter = vtk.vtkImplicitPolyDataDistance()
distanceFilter.SetInput(surface_World);
nOfFiduciallPoints = markupsNode.GetNumberOfFiducials()
for i in range(0, nOfFiduciallPoints):
point_World = [0,0,0]
markupsNode.GetNthControlPointPositionWorld(i, point_World)
closestPointOnSurface_World = [0,0,0]
closestPointDistance = distanceFilter.EvaluateFunctionAndGetClosestPoint(point_World, closestPointOnSurface_World)
indexCol.InsertNextValue(i)
labelCol.InsertNextValue(markupsNode.GetNthControlPointLabel(i))
distanceCol.InsertNextValue(closestPointDistance)
#fill the list with closest point distances
dist_list.append(closestPointDistance)
# Create a table from result arrays
resultTableNode = slicer.mrmlScene.AddNewNodeByClass("vtkMRMLTableNode", "Points from surface distance")
resultTableNode.AddColumn(indexCol)
resultTableNode.AddColumn(labelCol)
resultTableNode.AddColumn(distanceCol)
# Show table in view layout
slicer.app.layoutManager().setLayout(slicer.vtkMRMLLayoutNode.SlicerLayoutFourUpTableView)
slicer.app.applicationLogic().GetSelectionNode().SetReferenceActiveTableID(resultTableNode.GetID())
slicer.app.applicationLogic().PropagateTableSelection()
# math for min, max, mean and RMS values
max_value = max(dist_list)
min_value = min(dist_list)
mean_value = statistics.mean(dist_list)
sqrt_list = [math.sqrt(abs(i)) for i in dist_list]
rms_value = statistics.mean(sqrt_list)
# print min, max, mean and RMS values
print("The min value is " + "%.3f" % min_value + ".")
print("The max value is " + "%.3f" % max_value + ".")
print("The mean value is " + "%.3f" % mean_value + ".")
print("The Root Mean Square Value is " + "%.3f" % rms_value + ".")
Ok, this new code seems to work and get the same value. I had to take the absolute value of all the entries on the table and take the average of them.
I got confused by the sqrt in the fiducials-model code and thought I had to take the sqrt of all the values in the list but apparently not.
import statistics
import math
markupsNode = getNode('F')
modelNode = getNode('mymodel')
# Transform model polydata to world coordinate system
if modelNode.GetParentTransformNode():
transformModelToWorld = vtk.vtkGeneralTransform()
slicer.vtkMRMLTransformNode.GetTransformBetweenNodes(modelNode.GetParentTransformNode(), None, transformModelToWorld)
polyTransformToWorld = vtk.vtkTransformPolyDataFilter()
polyTransformToWorld.SetTransform(transformModelToWorld)
polyTransformToWorld.SetInputData(modelNode.GetPolyData())
polyTransformToWorld.Update()
surface_World = polyTransformToWorld.GetOutput()
else:
surface_World = modelNode.GetPolyData()
# Create arrays to store results
indexCol = vtk.vtkIntArray()
indexCol.SetName("Index")
labelCol = vtk.vtkStringArray()
labelCol.SetName("Name")
distanceCol = vtk.vtkDoubleArray()
distanceCol.SetName("Distance")
# Create blank list becuase I dont know how to use vtk arrays yet
dist_list = []
distanceFilter = vtk.vtkImplicitPolyDataDistance()
distanceFilter.SetInput(surface_World);
nOfFiduciallPoints = markupsNode.GetNumberOfFiducials()
for i in range(0, nOfFiduciallPoints):
point_World = [0,0,0]
markupsNode.GetNthControlPointPositionWorld(i, point_World)
closestPointOnSurface_World = [0,0,0]
closestPointDistance = distanceFilter.EvaluateFunctionAndGetClosestPoint(point_World, closestPointOnSurface_World)
indexCol.InsertNextValue(i)
labelCol.InsertNextValue(markupsNode.GetNthControlPointLabel(i))
distanceCol.InsertNextValue(closestPointDistance)
#fill the list with closest point distances
dist_list.append(closestPointDistance)
# Create a table from result arrays
resultTableNode = slicer.mrmlScene.AddNewNodeByClass("vtkMRMLTableNode", "Points from surface distance")
resultTableNode.AddColumn(indexCol)
resultTableNode.AddColumn(labelCol)
resultTableNode.AddColumn(distanceCol)
# Show table in view layout
slicer.app.layoutManager().setLayout(slicer.vtkMRMLLayoutNode.SlicerLayoutFourUpTableView)
slicer.app.applicationLogic().GetSelectionNode().SetReferenceActiveTableID(resultTableNode.GetID())
slicer.app.applicationLogic().PropagateTableSelection()
# math for min, max, mean
abs_list = [abs(i) for i in dist_list]
absolute_max_value = max(abs_list)
absolute_min_value = min(abs_list)
absolute_mean_value = statistics.mean(abs_list)
# print min, max, mean
print("The absolute min value is " + "%.3f" % absolute_min_value + ".")
print("The absolute max value is " + "%.3f" % absolute_max_value + ".")
print("The absolute mean value is " + "%.3f" % absolute_mean_value + ".")
I made this code to compute the RMS how it is defined in wikipedia vs how it is computed in the fiducial-model distance code. As you can see there is a difference.
import math
dist_list = [-0.0778, 0.0521, -0.0069, -0.0620, 0.0526, 0.0078]
# How it is defined in wikipedia
squared_list = [i ** 2 for i in dist_list]
sum_of_squared_list = sum(squared_list)
squared_mean = sum_of_squared_list / len(squared_list)
rms_wiki = math.sqrt(squared_mean)
# How it is computed in fiducial-model module
squared_list2 = [i ** 2 for i in dist_list]
root_list = [math.sqrt(j) for j in squared_list2]
rms_fid_mod = sum(root_list) / len(root_list)
print (rms_wiki)
print (rms_fid_mod)
Definition of RMSE above is correct (root of the mean of squared differences). It is up to you if you use RMS or MAD (mean of absolute differences) as error metric. RMSE is more sensitive to outliers.
I have made your code above into a scripted module. I have also added some output boxes in the UI to show some other data such as mean of absolute differences, root mean square, min, max etc.
I have got it working but some of my code feels a bit dubious. I am still very newb at programming. I am also not really sure what to put in the test area.
I have made an icon graphic and put it in the the icon folder but I am not sure where to put it for displaying in the extension store/manager thing.
do you mind taking a look at the files in github here?
Hi Juicy,
Andras already reply a efficient answer to you, and I almost miss the boat, but it was also possible for you to use this function of GeodesicSlicer to have the measurement of your fiducial points by geodesic distance.
Best,
Frederic
For me, it seems that @Juicy is interested in length of the straight line connecting a point and the closest surface point, while GeodesicSlicer computes length of the curve that connects two points on the surface. Both of them are important and interesting metrics, but they are different.
@Juicy it is very useful to have this module that computes point-to-surface error metrics. As a test, you can create a small test function that creates a model node (e.g., using vtkCylinderSource), adds a few points (just hardcode their positions), run the logic, and check if the computed values match what you expect (values that you computed and verified before).
It would be nice to extend it to be able to compute point-to-point error metrics, too (for people who use landmark point based registration).
You can submit this module as a separate extension, but it would be probably simpler to just add it to SlicerIGT extension.
Hi @lassoan, I have done some more work on the Fiducial to Model Distance scripted module.
I have added the ability to compute the same set of error metrics between two point clouds made of fiducials.
I have added a simple test which adds a 100 x 100 x 100mm vtk cube and then adds a set of six “fixed” reference fiducial points on the surface of the cube. Then another set of 6 “moving” fiducials are added around the cube before running both the Fiducial to Model and then Fiducial to Fiducial algorithms. 3 of the “moving” fiducial points are 1mm from the cube surface and 3 more are 2mm from the cube surface so the average distance will be 1.5mm.
I was reading through your powerpoint and one thing I don’t really understand is how to make a parameter node to save the settings in the GUI?
Do you mind having a quick look over the code here? Let me know if you have any suggestions. Again, sorry about any dubious coding.
