Revisiting Rigid Mesh Registration

manjula · December 7, 2019, 1:43pm

I first posted the problem few months back and Prof: Lassoan help me with it. Now i have had time to read about this topic and play around I would like to have a better understanding of the problem as i am not 100% happy with the results that i am getting with CMF surface registraion module and IGT fiducial registration wizard. i am no mathematician or programmer and many parts of this will be copy pasted from other published work whici i read.

Original Problem : We have a STL files from surface scans of bone from preop and post op with changes. We have initially put landmarks on the bone.

I have attached the sample data from our original scan. Blue is the pre op and the yellow is the post op.
We need to align the two meshes perfectly so we can see the difference of bone volume. (model to model distance and shape population viewer.)

As i understand mathematical problem

Given the coordinates of a set of points measured in two Cartesian coordinate systems (left, right) find the rigid transformation, T, between the two systems so that for corresponding points Pr , Pl we have:
Pr = T(Pl )

From what i read on this,

Pairing between points is known → closed form (analytic) solutions.
Pairing between points is unknown, → iterative solutions (require initialization and only guarantee convergence to local optimum).

And the advantages and disadvantages of the each method are,

Analytic solution

requires pairing
sensitive to outliers
guarantees optimality
assumes noise is isotropic and IID (~N(0,))

Iterative solution

does not require a known pairing
sensitive to outliers
does not guarantee optimality
assumes noise is isotropic and IID (~N(0,))
requires initialization
requires use of spatial data structures to achieve reasonable running times on a reasonably large data set.

My Question 01

CMF surface registration falls under ICP right ?

My Question 02

IGT Fiducial registration wizard → Is it ICP or Analytic solution ??? or a hybird ???

My Experiment

I made 3 stl models with Blender.
Red and Green models are exact duplicates but with two locations. (2 origins)

Yellow models is also a exact duplicate but i distorted one part of the mesh and with as you can see different origin.

Then I use CMF registriaon of Red and Green models. And the results are excellent.

ExactModelsSuperimposed

However when i use it to align the distorted model (yellow) it gives very very bad results. ( with max iterations 10000 and max landmarks 4000 )

Also i used IGT Fiducial registration with 4 landmarks and the results were not good but better than CMF registration.

Then i used the cloud compare on the same and even with the default settings it gives very good results.

So my question 3

Why is 3D slicer ICP method in CMF registration module gives inferior results to cloudcompare ?

For my original study i did registration with IGT module with manual landmarks and the results were good as evident by the colour map. we can see exact superimposition on the landmarks we fixed.

But i could get the same results with the cloud compare too in which alignment was as above completely automated. ( I think with ICP registration - cloud compare fine alignment)

So Question 04

Where is the problem with CMF registration module in comparison to cloud compare ? Why is it failing ? I edited the source file of CMF slicer module and tried it with 10000 landmarks and the results were the same albeit it was lot of computer power for my laptop - had to wait a long time !.

Question 05

In IGT Fiducial registration wizard is there any part which is automatic or done by a algorithm ? or is it s just transform or move the mesh landmarks we select on to each other without any kind of matching ?

Question 06
Is there a remedy to improve this ? As i prefer a fully automated system (without selecting landmarks) to align the models with better precision.

lassoan · December 7, 2019, 3:08pm

Typically, you align surfaces approximately first (if they are not already aligned, either manually translating/rotating or using landmark registration), then perform surface-based registration (masking out or removing regions that should be ignored).

Landmark registration result can be solved without iterations and it gives the global optimum.

I haven’t used this module but it seems that it can do both landmark registration and surface registration (using ICP).

Fiducial registration wizard performs landmark registration (non-iterative, global optimum solution). You can enable automatic point matching up to 7 landmark points.

You can get accurate alignment reliably:

remove those regions from the moving surface that you want to be ignored (regions affected by the medical procedure)
perform initial alignment if surfaces are not aligned already; if orientation and region of interest approximately matches then automatic centroid-based alignment may be enough
register using ICP; it is important to choose the clipped surface as moving and the full surface as fixed surface

manjula · December 7, 2019, 3:19pm

@lassoan Many thanks again for taking the time.

I have a question that came out of this ? I don’t know it is valid or just a stupid question.

Is ICP better than land mark registration ? or which one is superior or are they the same. ?

2nd. How do i remove the region that has changed before ICP registration and then add it back to the new location to calculate model to model distance ? i can think of how to do this in dicom data with masking but these are stl files. How to do I do that ???

Thank you.

lassoan · December 8, 2019, 2:42pm

Are there any more error details shown in the application log?

manjula · December 8, 2019, 3:46pm

@lassoan

These are the results I get with the cloud compare fine registration with RMS. As you can see when i look at the model to model distance and view it in shape population viewer with the same setting… our land marks show almost exact matches ( Green ) in the viewer.

I get almost same resutls with the cloud compare iteration method also. Again very good match on the landmarks we put in surgically.

However when i used the CMF surface registration results are like this, with not very good registration.

I managed to get good results with IGT fiducial registration wizard but it was after lot of playing around with fiducial markers and changing ti manually until i could get very good alignment. This took lot of time to chaging the location of fiducials until it gave this result.

**So my major Problem is ** why is cloud compare fine alignment gives better results in comparison to CMF surface registration ? Because as i read on both the module they use the ICP algorithms and why should it then give different results ?

Today i installed nightly version because i wanted to try it with CMFreg ROI method since i did not knew how to do only to select the region of interest as you suggested in a stl file and as for some reason in the stable version in my computer ROI registration was not marking the regions of interest for some reason. So i tried with nightly version and it worked and gave excellent results with ease compared to IGT registration. But still cloud compares registration is fully automated as i do not need to select any points.

So i am out of interest and also because almost all the research that we have planned or are currently been conducted involves registration i am very much interested in learning about this and for publication purposes i think it is much convenient to stick to one software (3D slicer) to do the registration and the model to model distance and colour mapping.

lassoan · December 9, 2019, 5:15am

Maybe there is some heuristics in CloudCompare’s registration algorithms that make the registration ignore the modified parts of the surface.

Instead of relying on heuristics, it would be safer to explicitly restrict registration to only take into account neighborhood of landmarks. You can write a short Python script for this:

compute distance from landmark points on the “moving” mesh: create a polydata from the landmark points and use that in vtkDistancePolyDataFilter
clip “moving” mesh to only contain regions that are near landmark points (maximum distance should be approximate size of the landmark object, or maybe a bit larger) using vtkThreshold
register the two meshes approximately using landmark transform
refine registration by registering the clipped “moving” mesh to the “fixed” mesh (the other mesh, without any clipping) using surface registration

Maybe CMFreg already supports something like this, but if not, then it should be quite easy to add there. Or, such a registration module could be added to SlicerMorph extension, too.

@bpaniagua @muratmaga

manjula · December 9, 2019, 8:33am

Dear Prof @lassoan
I started learning python but at the moment not capable of writing that script. But i am willing to try it if someone can help me with it when your time permits it.

we are happy with the results we get with the CMFreg ROI registration and we are going to analyse the data with it for now. But if we can get the surface registration to work with the python code i am willing to learn it.

mau_igna_06 · October 14, 2021, 2:53pm

Hi @lassoan. I am thinking about developing this improved surface registration you were talking about. Although I would need a little bit more info.

A question, would the improved version work as well as this one?

Would some kind of paint selection of cells of the mesh be possible using the segment editor or using this script as base? Would that improve usability?

This blender’s ICP is currently used in blender to plan deformity correction surgeries using the normal bone (contralateral) to guide the aligment of the to-be-aligned deformed bone pieces after the closing/opening wedge osteotomy. But they lack access to the ICP transform that results from the algorithm in blender. So they can’t plan kirchner wires on the corrected bone segments and transform them to the deformed original bone to create surgical guides for cut and aligment.

Related relevant post: Robust ICP alignment tool

mau_igna_06 · October 17, 2021, 8:00pm

Meanwhile I got the faces selection working:
(copy it to a code editor for better display)

fiducialList = getNode('MarkupsFiducial')
appendFilter = vtk.vtkAppendPolyData()
for i in range(fiducialList.GetNumberOfControlPoints()):
    position = [0,0,0]
    fiducialList.GetNthControlPointPosition(i,position)
    sphere = vtk.vtkSphereSource()
    sphere.SetRadius(0.1)
    sphere.SetCenter(position)
    sphere.SetThetaResolution(5)
    sphere.SetPhiResolution(5)
    sphere.Update()
    appendFilter.AddInputData(sphere.GetOutput())

appendFilter.Update()

myModel = getNode('deformed_Right_tibia')

distanceFilter = vtk.vtkDistancePolyDataFilter()
distanceFilter.SetInputData(0, myModel.GetPolyData())
distanceFilter.SetInputData(1, appendFilter.GetOutput())
distanceFilter.SignedDistanceOff()
distanceFilter.Update()

myModel.SetAndObservePolyData(distanceFilter.GetOutput())
#Turn on distance scalars on Models module's properties

"""
Creating a points-only polydata doesn't work
points = vtk.vtkPoints()
vertsCellArray = vtk.vtkCellArray()
appendFilter = vtk.vtkAppendPolyData()
for i in range(fiducialList.GetNumberOfControlPoints()):
    position = [0,0,0]
    fiducialList.GetNthControlPointPosition(i,position)
    pointId = points.InsertNextPoint(position)
    vertsCellArray.InsertNextCell(1)
    vertsCellArray.InsertCellPoint(pointId)

polydata = vtk.vtkPolyData()
polydata.SetPoints(points)
polydata.SetVerts(vertsCellArray)
polydata.Modified()

myModel = getNode('deformed_Right_tibia')

distanceFilter = vtk.vtkDistancePolyDataFilter()
distanceFilter.SetInputData(0, myModel.GetPolyData())
distanceFilter.SetInputData(1, polydata)
distanceFilter.SignedDistanceOff()
distanceFilter.Update()


myModel.SetAndObservePolyData(distanceFilter.GetOutput())
"""

thresholdFilter = vtk.vtkThreshold()
thresholdFilter.SetInputData(distanceFilter.GetOutput())
maximumDistanceToPointsmm = 15
thresholdFilter.ThresholdBetween(0, maximumDistanceToPointsmm)
thresholdFilter.SetInputArrayToProcess(0, 0, 0,
    vtk.vtkDataObject.FIELD_ASSOCIATION_CELLS,
    "Distance")
#thresholdFilter.SetInputArrayToProcess(0, 0, 0,
#   vtk.vtkDataObject.
#    FIELD_ASSOCIATION_POINTS, "Distance")
thresholdFilter.Update()
geometryFilter = vtk.vtkGeometryFilter()
geometryFilter.SetInputData(thresholdFilter.GetOutput())
geometryFilter.Update()
myModel2 = getNode('deformed_Right_tibia Copy')
myModel2.SetAndObservePolyData(geometryFilter.GetOutput())

Faces need to be well-shaped for it to work well I think. Here is an example:

lassoan · October 17, 2021, 11:49pm

Thanks for the update, this looks nice. Could you turn this into a Dynamic modeler tool? It should be fairly straightforward to implement and would make the feature conveniently available to end users.

We can already use curves and planes to define regions on surfaces, but sometimes using points is easier and can be more robust (especially if you simply add all points within a certain radius instead of marching on the surface). We could also make the output cropping optional: the end result could be the cropped model; or the original model storing selection in point/cell scalars.

mau_igna_06 · October 18, 2021, 6:39pm

Thanks for the update, this looks nice. Could you turn this into a Dynamic modeler tool? It should be fairly straightforward to implement and would make the feature conveniently available to end users.

I’m doing it. I have done a full Slicer build. I have my DynamicModeler changes in my fork of the SlicerSurfaceToolbox how do I make the built Slicer load my DynamicModeler changes. I guess I have to change a variable of the superbuild to point to my repository. Or how do I do to build my new code of the SlicerSurfaceToolbox?

or the original model storing selection in point/cell scalars.

This is done in a script, I have to implement it yet in the new SelectionTool class of the Dynamic Modeler.

This is the script to have selection scalars:

maximumDistancemm = 30
cellScalars = myModel.GetMesh().GetCellData()
distanceArray = cellScalars.GetArray("Distance")
selectionArray = vtk.vtkIntArray()
selectionArray.SetName("Selection")
selectionArray.SetNumberOfValues(myModel.GetMesh().GetNumberOfCells())
for i in range(myModel.GetMesh().GetNumberOfCells()):
    if distanceArray.GetTuple1(i) < maximumDistancemm:
        selectionArray.SetTuple1(i,1)
    else:
        selectionArray.SetTuple1(i,0)

cellScalars.AddArray(selectionArray)

# Set up coloring by selection array
myModel.GetDisplayNode().SetActiveScalar("Selection", vtk.vtkAssignAttribute.CELL_DATA)
myModel.GetDisplayNode().SetAndObserveColorNodeID("vtkMRMLColorTableNodeWarm1")
myModel.GetDisplayNode().SetScalarVisibility(True)

The script would be faster to execute if the selection-fiducials created just points polydata instead of small spheres. @lassoan you can check the commented part of the first script I posted (more above this comment) to see if you can find out why this faster approach doesn’t work (there are no Distance arrays created for the model).

mau_igna_06 · October 19, 2021, 11:01am

Building Slicer (without changing the source code) with the VisualStudio GUI fails and gives like 500 errors, building it from the command line succeeds. Although I don’t like it building from the command line because it takes very long, not as much as a clean build, but like 70% of the time. Building from the GUI usually takes a lot once, but then it only rebuilds the changed files and it is a lot faster.

Now I changed ‘SuperBuild.cmake’ to point to my forked repository of the surfaceToolbox and changed the commit_tag to the last commit of the selectionToolFeatureBranch. And I’m doing a clean build from the command line. I think it should work.

How can I make Slicer compile faster from the command line? Why there is not match in compiling results using the command line and using the GUI of Visual Studio? (I followed the build instructions for windows)

Here is my branch if you want to check it out: GitHub - mauigna06/SlicerSurfaceToolbox at SelectionToolFeature

lassoan · October 20, 2021, 6:02am

You can, but you don’t have to. You can make changes locally, test everything, and send a pull request to the SurfaceToolbox repository to get your changes integrated into Slicer.

Instead of using vtkDistancePolyDataFilter, you can get the list of points within a certain radius using a fast point locator.

You could also add an option to add points based on geodesic distance (distance along the surface) using these filters.

You must set the build type (debug/release) to the same that you used when you configured your project using the command line.

mau_igna_06 · October 20, 2021, 11:04am

Instead of using vtkDistancePolyDataFilter, you can get the list of points within a certain radius using a fast point locator .

Here is the faster (a lot faster) script for selection using what you suggested:


fiducialList = getNode('MarkupsFiducial_2')
myModel = getNode('normalBone')

modelPolydata = vtk.vtkPolyData()
modelPolydata.ShallowCopy(myModel.GetPolyData());

# Initialize the locator
pointTree = vtk.vtkPointLocator();
pointTree.SetDataSet(modelPolydata);
pointTree.BuildLocator();

selectionArray = vtk.vtkIntArray()
selectionArray.SetName("Selection")
selectionArray.SetNumberOfValues(modelPolydata.GetNumberOfPoints())
selectionArray.Fill(0)

maximumDistanceToPointsmm = 5

for i in range(fiducialList.GetNumberOfControlPoints()):
    position = [0,0,0]
    fiducialList.GetNthControlPointPosition(i,position)
    result = vtk.vtkIdList()
    pointTree.FindPointsWithinRadius(maximumDistanceToPointsmm, position,
                                    result);
    
    for j in range(result.GetNumberOfIds()):
        point_ind = result.GetId(j);
        selectionArray.SetTuple1(point_ind, 1);


pointScalars = modelPolydata.GetPointData()
pointScalars.AddArray(selectionArray)

thresholdFilter = vtk.vtkThreshold()
thresholdFilter.SetInputData(modelPolydata)
thresholdFilter.ThresholdBetween(0.9, 1.1)
thresholdFilter.SetInputArrayToProcess(0, 0, 0,
    vtk.vtkDataObject.FIELD_ASSOCIATION_POINTS,
    "Selection")
thresholdFilter.Update()
geometryFilter = vtk.vtkGeometryFilter()
geometryFilter.SetInputData(thresholdFilter.GetOutput())
geometryFilter.Update()

selectedNormalBoneDistal = slicer.mrmlScene.AddNewNodeByClass('vtkMRMLModelNode','selectedNormalBone')
selectedNormalBoneDistal.CreateDefaultDisplayNodes()
displayNode = selectedNormalBoneDistal.GetDisplayNode()
displayNode.SetColor(0,0,1)
selectedNormalBoneDistal.SetAndObservePolyData(geometryFilter.GetOutput())

If you think the script is okey, I’ll update the selectionToolBranch.

I could filter selected points a little more by checking if the normal of the inside-sphere points (given by the pointLocator) has less than 90 degress with the normal of pointNearestToTheFiducial. That would avoid selection though the surface when radius of the sphere is big I think. What do you think about this?

I could try to make the geodesic distance option. With the change above do you think it will still be needed? Do you know if the geodesic-distace algorithm execution is fast? (if it takes as long as vtkDistancePolyDataFilter, it may not be worth it, or does it?)

lassoan · October 21, 2021, 4:04am

That usually does not work very well (normals can be noisy, so the edge can be very rough and it is also hard to find a good angle threshold), so we can leave it for now.

The geodesic distance filter may be faster. The speed may depend on the maximum distance - for small distances the propagation should stop earlier.

mau_igna_06 · October 22, 2021, 11:14am

Hi @lassoan. I added the geodesicDistance algorithm to the SelectionTool of the dynamicModeler (you can checkout my branch). The code should be 95% okey in theory but I wasn’t able to try it because I’m having trouble compiling Slicer. It fails from Visual Studio GUI with 2 errors. Now I’m trying to build from the command line. I don’t know if it will work out. It takes a lot of time.

Meanwhile I was trying to compile the example of Geodesics and I cannot make it work. These are my commands:

cd GeodesicsOnMeshesRelease
"C:\Program Files\CMake\bin\cmake.exe" -G "Visual Studio 16 2019" -A x64 -DVTK_DIR:PATH=C:\D\S4R\VTK-build ..\GeodesicsOnMeshes
"C:\Program Files\CMake\bin\cmake.exe" --build . --config Release

If fails with this log:

Microsoft (R) Build Engine versión 16.9.0+5e4b48a27 para .NET Framework
Copyright (C) Microsoft Corporation. Todos los derechos reservados.

  Checking Build System
  Building Custom Rule C:/Users/mau_i/Downloads/SourceCode3_GeodesicsOnMeshes/GeodesicsOnMeshes/FastMarching/CMakeLists
  .txt
  GW_Config.cpp
  GW_Face.cpp
  GW_FaceIterator.cpp
  GW_Mesh.cpp
  GW_SmartCounter.cpp
  GW_Vertex.cpp
  GW_VertexIterator.cpp
  GW_GeodesicFace.cpp
  GW_GeodesicMesh.cpp
  GW_GeodesicPath.cpp
  GW_GeodesicPoint.cpp
  GW_GeodesicVertex.cpp
  GW_TriangularInterpolation_Linear.cpp
  GW_TriangularInterpolation_Quadratic.cpp
  GW_TriangularInterpolation_Cubic.cpp
  Generando código...
  MeshGeodesics.vcxproj -> C:\Users\mau_i\Downloads\SourceCode3_GeodesicsOnMeshes\GeodesicsOnMeshesRelease\bin\Release\
  MeshGeodesics.lib
  Building Custom Rule C:/Users/mau_i/Downloads/SourceCode3_GeodesicsOnMeshes/GeodesicsOnMeshes/Source/CMakeLists.txt
  vtkPolyDataGeodesicDistance.cxx
C:\Users\mau_i\Downloads\SourceCode3_GeodesicsOnMeshes\GeodesicsOnMeshes\Source\vtkPolyDataGeodesicDistance.h(16,10): f
atal error C1083: No se puede abrir el archivo incluir: 'vtkPolyDataAlgorithm.h': No such file or directory [C:\Users\m
au_i\Downloads\SourceCode3_GeodesicsOnMeshes\GeodesicsOnMeshesRelease\Source\vtkMeshGeodesics.vcxproj]
  vtkFastMarchingGeodesicDistance.cxx
C:\Users\mau_i\Downloads\SourceCode3_GeodesicsOnMeshes\GeodesicsOnMeshes\Source\vtkPolyDataGeodesicDistance.h(16,10): f
atal error C1083: No se puede abrir el archivo incluir: 'vtkPolyDataAlgorithm.h': No such file or directory [C:\Users\m
au_i\Downloads\SourceCode3_GeodesicsOnMeshes\GeodesicsOnMeshesRelease\Source\vtkMeshGeodesics.vcxproj]
  vtkFastMarchingGeodesicPath.cxx
C:\Users\mau_i\Downloads\SourceCode3_GeodesicsOnMeshes\GeodesicsOnMeshes\Source\vtkFastMarchingGeodesicPath.h(40,10): f
atal error C1083: No se puede abrir el archivo incluir: 'vtkGeodesicPath.h': No such file or directory [C:\Users\mau_i\
Downloads\SourceCode3_GeodesicsOnMeshes\GeodesicsOnMeshesRelease\Source\vtkMeshGeodesics.vcxproj]
  vtkPolygonalSurfaceContourLineInterpolator2.cxx
C:\Users\mau_i\Downloads\SourceCode3_GeodesicsOnMeshes\GeodesicsOnMeshes\Source\vtkPolygonalSurfaceContourLineInterpola
tor2.h(48,10): fatal error C1083: No se puede abrir el archivo incluir: 'vtkPolyDataContourLineInterpolator.h': No such
 file or directory [C:\Users\mau_i\Downloads\SourceCode3_GeodesicsOnMeshes\GeodesicsOnMeshesRelease\Source\vtkMeshGeode
sics.vcxproj]
  Generando código...

Do you have any suggestions of how to make it work? (Maybe the error happens because I’m trying to build using VTK-build and I don’t know if it finished compiling in the Slicer superbuild, commandline is processing)

Also I don’t know what I have to do to make the Geodesics library compile inside the dynamicModeler (or how to make vtkFastMarchingGeodesicDistance findable for my SelectionTool)

mau_igna_06 · October 22, 2021, 7:21pm

I could compile Slicer well from the command line. So now the next thing to do is to compile the new Dynamic Modeler and test it. But I don’t know how to add the vtkFastMarchingGeodesicDistance to the DynamicModeler logic, it has a folder full of dependencies. And I don’t know how to configure the CMake.

I tried to follow how Combine Models does it in its CMakeLists.txt but the Geodesics have a lot more dependencies.

I would like the geodesics filter to be available from the python interpreter also.

mau_igna_06 · October 26, 2021, 6:41pm

Ok. I think I fixed most of the bugs. There is only one left:

Gravedad	Código	Descripción	Proyecto	Archivo	Línea	Estado suprimido
Error	LNK1181	no se puede abrir el archivo de entrada 'FastMarching.lib'	vtkSlicerDynamicModelerModuleLogic	C:\D\S4R\Slicer-build\E\SurfaceToolbox\DynamicModeler\Logic\LINK	1

Link to the code

@lassoan maybe you can help with this

mau_igna_06 · October 27, 2021, 1:24pm

All working. I submitted a pull-request here.

Topic		Replies	Views
STL models - Alignment and Volume measurement Support	16	3306	August 26, 2023
Surface Registration (STL model with cbct data) Support	6	2249	January 8, 2020
Surface Scan Registration Support	4	525	February 20, 2020
Fiducial Registration Wizard Problem Support markups-fiducials	2	828	November 14, 2022
Difference Surface Registration vs Model Registration Support registration , stl	1	683	March 14, 2021

Revisiting Rigid Mesh Registration

Related topics