Rigid alignment
The goal in this tutorial is to learn how to perform rigid alignment of shapes in Scalismo.
In case you would like to learn more about the syntax of the Scala programming language, have a look at the Scala Cheat Sheet. This document is accessible under:
Documents -> Scala Cheat Sheet
Note: for optimal performance and to keep the memory footprint low, we recommend to restart Scalismo Lab at the beginning of every tutorial.
Quick view on Transformations
Let's start by loading and showing Paola's mesh again:
val mesh : TriangleMesh = MeshIO.readMesh(new File("datasets/Paola.stl")).get
show(mesh, "Paola")Scalismo allows us to perform geometric transformations on meshes.
Transformations are functions that map a given 2D or 3D point, into a new transformed point.
For example, such a transformation could be a Translation:
val translation = TranslationTransform[_3D](Vector(100,0,0))When given a point, the defined transform will then simply return a new point translated by a vector (100,0,0)
val pt : Point[_3D] = translation(Point(1,1,1))Another example of classical transformations are Rotations:
val rotation : RotationTransform[_3D] = RotationTransform(0f,3.14f,0f)Here we created a rotation transform around the space origin, by specifying the 3 Euler angles defining the rotation. In this particular case, we rotate every point with approximately 180 degrees around the Y axis (centered at the origin of the space).
val pt2 : Point[_3D] = rotation(Point(1,1,1))In Scalismo, such transformations can be applied either to single points individually or to collections of points such as triangle meshes:
val translatedPaola : TriangleMesh = mesh.transform(translation)
show(translatedPaola, "translatedPaola")Here, we used the transform method of the TriangleMesh class. This method takes as input a transformation function that maps a 3D point into another 3D point and applies it to all the points of the mesh while maintaining the same triangulation. As a result, we obtain a new translated version of Paola.
Exercise: Apply the rotation transform to the original mesh of Paola and show the result
val rotated : TriangleMesh = ???val rotated = mesh.transform(rotation)
// rotated: scalismo.mesh.TriangleMesh = TriangleMesh(Vector(Point3D(-161.78987,-11.115057,-301.44525), Point3D(-161.78427,-10.857458,-301.40735), Point3D(-161.5174,-10.88125,-301.31372), Point3D(-161.57983,-11.205217,-301.39404), Point3D(-161.78093,-10.600695,-301.3632), Point3D(-161.46243,-10.561865,-301.21378), Point3D(-161.80467,-10.281149,-301.32828), Point3D(-161.46033,-10.202814,-301.13956), Point3D(-161.84917,-9.97047,-301.24756), Point3D(-161.48087,-9.865271,-301.00626), Point3D(-161.90215,-9.595915,-301.15878), Point3D(-161.51779,-9.481844,-300.90262), Point3D(-161.9766,-9.239666,-301.0281), Point3D(-161.56274,-9.10297,-300.78204), Point3D(-162.05594,-8.832981,-300.93127), Point3D(-161.62268,-8.678461,-300.66003), Point3D(-162.13518,-8.426069,-300.83377), Point3D(-161.6905,-8.257...show(rotated, "rotatedPaola") (In case you do not see the rotated mesh, make sure to right click and drag down in the 3D scene to zoom out. You can also click the "RC" button on top of the 3D view to Reset the Camera to view all objects in the scene)
Rigid transformations
A rigid transformation is simply a combination of a rotation and a translation transform.
Let's now compute a rigid transformation of Paola, by rotating the translated mesh:
val rigidPaola = translatedPaola.transform(rotation)
show(rigidPaola, "rigidPaola")Note: since the rotation is around the origin, you might have to zoom out (hold right click and drag down) to see the result.
rigidPaola is now the result of a rigid transformation combining both a translation and a rotation on Paola's mesh.
Exercise: Build a composed transformation from the individual translation and rotation transforms (translation and rotation) that, when applied to Paola's original mesh, would give the same result as rigidPaola
def composedTransformation(p: Point[_3D]) : Point[_3D] = ??? def composedTransformation(p: Point[_3D]) : Point[_3D] = rotation(translation(p))
// composedTransformation: (p: scalismo.geometry.Point[scalismo.geometry._3D])scalismo.geometry.Point[scalismo.geometry._3D]val rigid2 = mesh.transform(composedTransformation)
show(rigid2, "composedRigid")Rigid alignment
Let's clean up the 3D scene a bit:
remove("translatedPaola")
remove("rotatedPaola")
remove("composedRigid")Suppose now that we are left with the original mesh and the rigidly transformed one (rigidPaola in the scene).
Suppose also that we do not know the parameters of the translation and rotation that led to rigidPaola.
How can we retrieve those parameters and obtain a transformation from the original mesh to rigidPaola ?
Well, as you saw previously, this can be solved by indicating a few corresponding points, usually referred to as landmarks, and performing a Procrustes anaylsis to retrieve the best transform mapping a set of the landmarks on the first shape into the set of landmarks on the second shape.
Exercise: Click the following landmarks, in the indicated order on both meshes : 1) right corner of the right eye, 2) left corner of the left eye, 3) tip of the nose, 4) middle of chin
Note: In case you have trouble clicking the landmarks, you can read them from file and add them to the scene as such:
// Execute this only if you had trouble clikcing landmarks
val paolaFS = LandmarkIO.readLandmarksJson[_3D](new File("datasets/paolaFailsafe.json")).get
addLandmarksTo(paolaFS, "Paola")
val rigidPaolaFS = LandmarkIO.readLandmarksJson[_3D](new File("datasets/rigidPaolaFailsafe.json")).get
addLandmarksTo(rigidPaolaFS, "rigidPaola")Let's now programmatically retrieve the landmarks:
val originalLms = getLandmarksOf("Paola").get
val rigidLms = getLandmarksOf("rigidPaola").getNow that we have the lists of landmarks, we can apply Procrustes analysis to retrieve the best rigid transformation from the original set of landmarks to the new one as follows:
val bestTransform = LandmarkRegistration.rigid3DLandmarkRegistration(originalLms, rigidLms)Note: in case this resulted in an exception(Empty set of landmarks), please make sure that you clicked the indicated set of landmarks in the exercise above
Here we used the rigid3DLandmarkRegistration method that takes as input the list of original landmarks to be transformed as well as the destination landmarks that should be matched at best.
As you can see in the result pane, the returned type of the retrieved transform is RigidTransformation[_3D]
Let's now apply it to the original set of landmarks, to see how well they are transformed :
val transLms = originalLms.map{l => bestTransform(l.point)}
show(transLms, "transLMs")Question: Are the transformed landmarks matching perfectly the set of rigidLms? If not, why not?
Answer: Most probably not, as it is practically impossible to manually click truly corresponding points. Therefore, when aligning the two sets of clicked points, chances are high that they don't overlap exactly.
Let's now also apply the transform to the entire mesh :
val alignedPaola = mesh.transform(bestTransform)
show(alignedPaola, "alignedPaola") (colour the mesh to visualize it better)
Question: Is the transformed mesh matching perfectly rigidPaola? If not, why not?
Answer: For the same reason as the answer above, the rigid transform computed based on non-truly corresponding landmarks will contain slight translation and rotation effects therefore resulting in a slight mismatch between the meshes.
Exercise: Load the mesh file named "datasets/323.stl" and perform a similar alignment of Paola to this mesh. Feel free to pick any set of landmarks you wish.
Solution: no new code needed here. Simply read the new mesh and reiterate the steps above.
Exercise: Did we really have to click landmarks? Both the original mesh of Paola and rigidPaola are in correspondence. Using this knowledge, can you come up with 2 sets of corresponding landmarks, without clicking them? Use this set to perform rigid alignment. How does the result compare to the case where we clicked landmarks?
val ptIds = Seq(PointId(100), PointId(400), PointId(500), PointId(600)) // any 4 point ids would do
// ptIds: Seq[scalismo.common.PointId] = List(PointId(100), PointId(400), PointId(500), PointId(600))
val originalLms = ptIds.map(id => Landmark("L"+ id.id, mesh.point(id)))
// originalLms: Seq[scalismo.geometry.Landmark[scalismo.geometry._3D]] = List(Landmark(L100,Point3D(167.59204,14.908397,302.99237),None,None), Landmark(L400,Point3D(161.62909,-5.67869,298.76025),None,None), Landmark(L500,Point3D(177.61209,58.655556,296.2613),None,None), Landmark(L600,Point3D(170.1445,42.28475,300.32037),None,None))
val rigidLms = ptIds.map(id => Landmark("L"+ id.id, rigidPaola.point(id)))
// rigidLms: Seq[scalismo.geometry.Landmark[scalismo.geometry._3D]] = List(Landmark(L100,Point3D(-267.10916,14.908397,-303.41815),None,None), Landmark(L400,Point3D(-261.15295,-5.67869,-299.17654),None,None), Landmark(L500,Point3D(-277.13995,58.655556,-296.70303),None,None), Landmark(L600,Point3D(-269.6659,42.28475,-300.7502),None,None))
val bestTransform = LandmarkRegistration.rigid3DLandmarkRegistration(originalLms, rigidLms)
// bestTransform: scalismo.registration.RigidTransformation[scalismo.geometry._3D] = <function1>