In this post, I am going to implement the partial and the full generalized Procrustes analyses (partial GPA and full GPA) with Scalismo. Before getting into the implementation, I may want read this post explaining the theory of GPA’s.

## Partial GPA

Through partial GPA, more than 2 configurations/meshes get fitted/superimposed onto a common (and initially unknown) mean configuration. The procedure is iterative and involves translation and rotation of the meshes.

GPA’s need meshes in correspondence; here, I load the femoral head meshes that I previously had their correspondence established through the non-rigid registration (this post). Initializing and loading the meshes:

This post contains code blocks; if your browser does not automatically display the code blocks, please click on the displayed dark narrow rectangular areas to display the containing code.

object GPA extends App { import scalismo.geometry._ import scalismo.common._ import scalismo.ui.api._ import scalismo.mesh._ import scalismo.io.{StatisticalModelIO, MeshIO} import scalismo.statisticalmodel._ import scalismo.registration._ import scalismo.statisticalmodel.dataset._ import java.awt.Color scalismo.initialize() implicit val rng = scalismo.utils.Random(42) val ui = ScalismoUI() // Loading data val drive = "I:/" val meshDir = drive + "registration/Registered/" val meshFiles = new java.io.File(meshDir).listFiles() meshFiles.foreach(file => println(s"\n$file")) println(meshFiles)

Visualizing the meshes:

val (meshes, meshViews) = meshFiles.map(meshFile => { val mesh = MeshIO.readMesh(meshFile).get val meshView = ui.show(dsGroup, mesh, "mesh") (mesh, meshView) // returns a tuple of the mesh and the associated view }) .unzip //unzip the tuple into two separate lists // setting some colors meshViews(0).color = new Color(85,235,40) meshViews(1).color = new Color(20,40,250) meshViews(2).color = new Color(250,150,20)

Here, methods are defined to calculate the centroid (center of mass) of a mesh and centering the mesh. This is not necessary for GPA, as it involves translation. Also, GPA method in Scalismo, calculates the centroid of a given reference mesh, and uses it as the center of rotation of the meshes through the process. It does not harm if all the meshes are initially centered. Also, it is consistent with the GPA algorithm explained in the related post.

def CM (mesh:TriangleMesh[_3D]): EuclideanVector3D = { val numOfPoints = mesh.pointSet.numberOfPoints var total = EuclideanVector3D(0,0,0) val meshPoints = mesh.pointSet.points //iterator while (meshPoints.hasNext){ //iterator total = total + (meshPoints.next()).toVector } total*(1d/numOfPoints) } val centeredMeshes = meshes.map(mesh => { val translationVec = CM(mesh)*(-1) val translation = TranslationTransform[_3D](translationVec) mesh.transform(translation) })

To do GPA’s with Scalismo, a DataCollection should be defined. A DataCollection needs a reference mesh and a sequence of meshes; all in correspondence. A GPA will be performed on the sequence of the meshes. The reference mesh is used as the domain of transformations. A mean mesh, in Scalismo, is defined as the mean deformations (vectors) of the meshes from the reference mesh, added to the reference mesh; it should be remembered that all the meshes are in correspondence. It does not mater which, but, one of the meshes should be selected as the reference mesh.

// General Procrustes Analysis: Aligns a set of meshes with respect to some unknown “mean” shape (independent of ordering of shapes) // setting a reference mesh. it is needed for data collection but no effect on the GPA val refMesh = centeredMeshes(2) val allMeshes = centeredMeshes val dataClctn = ((DataCollection.fromMeshSequence(refMesh,allMeshes))._1).get // to check how many meshes are in the data collection println(dataClctn.size)

Performing the partial GPA (number of iterations = 10, precision of *average point-wise distance* between meshes (halt distance) = 1e-5):

// doing the partial GPA val gpaDataClctn = DataCollection.gpa(dataClctn ,10,1e-5)

Retrieving the mean mesh and visualization:

// retrieving the mean mesh val meanMesh = gpaDataClctn.meanSurface //CM(meanMesh) // to check whether the mean mesh is centered or not val meanGroup = ui.createGroup("meanMesh") val meanMeshViewg = ui.show( meanGroup, meanMesh,"meanMesh") meanMeshViewg.color = new Color(245,20,50)

The transformed versions of the original meshes onto the mean mesh are called the **partial Procrustes fits**. They can be retrieved and viewed:

//aligning by translation and rotation all the meshes onto the mean mesh val meanMeshPoints = meanMesh.pointSet.points.toIndexedSeq val gpaTransforms = for (mesh <- allMeshes) yield { val meshPoints = mesh.pointSet.points.toIndexedSeq LandmarkRegistration.rigid3DLandmarkRegistration(meshPoints.zip(meanMeshPoints), center = Point3D(0, 0, 0)) } val gpaGroup = ui.createGroup("partial Gpa meshes") //partial GPA fits val transformedMeshes = (allMeshes.zip(gpaTransforms)).map(z => { val mesh=z._1 val transToMean = z._2 mesh.transform(transToMean) }) val transfrmdMeshViews = transformedMeshes.map(mesh => ui.show(gpaGroup, mesh, "trMesh"))

The centroid size (see this post) of each mesh in this example can be calculated by this method:

def centroidSize (mesh:TriangleMesh3D): Double ={ var total = 0d val meshPoints = mesh.pointSet.points //iterator val centroid = CM(mesh) // is a vector here while (meshPoints.hasNext){ //iterator total = total + scala.math.pow(((meshPoints.next().toVector - centroid).norm),2) } scala.math.sqrt(total) }

For example:

centroidSize(TransformedMeshes(0)) centroidSize(allMeshes(0)) centroidSize(MeanMesh)

The centroid sizes are (unit of length):

original mesh 1: 487.222, mesh 1 partial-GPA fit: 487.222

original mesh 2: 692.725, mesh 2 partial- GPA fit : 692.725

original mesh 3: 566.668, mesh 3 partial- GPA fit : 566.668

mean mesh: 577.732

As expected, the original and their partial-GPA fits have similar sizes due to no-scaling in the partial GPA. However, the mean mesh must have a different size; as it has.

## Full GPA

The full GPA entails scaling, translation and rotation of the meshes in order to superimpose each mesh onto the (unknown) mean mesh. Scalismo does not currently have a method performing the full GPA, but its partial *gpa* can be manipulated to do the full GPA. To this end, we should access the package *scalismo.statisticalmodel.dataset* in the file *DataCollection.scala* to copy the needed pieces of codes regarding the partial GPA method. This is already done as below. The copied lines of code are used to define new functions methods. The only difference between the method doing the full GPA and the partial one, is that *similarity3DLandmarkRegistration* method in the former one instead of *rigid3DLandmarkRegistration* in the latter one.

import scalismo.utils.Random def fGpa(dc: DataCollection, maxIteration: Int = 3, haltDistance: Double = 1e-5)(implicit rng: Random): DataCollection = { fGpaComputation(dc, dc.meanSurface, maxIteration, haltDistance) } def fGpaComputation(dc: DataCollection, meanShape: TriangleMesh[_3D], maxIteration: Int, haltDistance: Double)(implicit rng: Random): DataCollection = { if (maxIteration == 0) return dc val referencePoints = dc.reference.pointSet.points.toIndexedSeq val numberOfPoints = referencePoints.size val referenceCenterOfMass = referencePoints.foldLeft(Point3D(0, 0, 0))((acc, pt) => acc + (pt.toVector / numberOfPoints)) val meanShapePoints = meanShape.pointSet.points.toIndexedSeq // align all shape to it and create a transformation from the mean to the aligned shape val dataItemsWithAlignedTransform = dc.dataItems.par.map { dataItem => val surface = dc.reference.transform(dataItem.transformation) val transform = LandmarkRegistration.similarity3DLandmarkRegistration(surface.pointSet.points.toIndexedSeq.zip(meanShapePoints), referenceCenterOfMass) DataItem("gpa -> " + dataItem.info, Transformation(transform.compose(dataItem.transformation))) } val newdc = DataCollection(dc.reference, dataItemsWithAlignedTransform.toIndexedSeq) val newMean = newdc.meanSurface if (MeshMetrics.procrustesDistance(meanShape, newMean) < haltDistance) { newdc } else { fGpaComputation(newdc, newMean, maxIteration - 1, haltDistance) } }

Let’s compute the mean mesh and then superimpose (by the similarity transformations) the meshes onto the mean mesh; the code is:

val fGpaDataClctn = fGpa(dataClctn ,20,1e-5) val fGpaMeanMesh = fGpaDataClctn.meanSurface //CM(meanMesh) // to check whether the mean mesh is centered or not val fGpaMeanGroup = ui.createGroup("fullGpaMeanMesh") val fGpaMeanMeshView = ui.show( fGpaMeanGroup,fGpaMeanMesh ,"fullGpaMeanMesh") fGpaMeanMeshView.color = new Color(245,20,50) val fGpaGroup = ui.createGroup("Full Gpa Meshes") val fGpaMeanMeshPoint = fGpaMeanMesh.pointSet.points.toIndexedSeq val fullGpaTransforms = for (mesh <-allMeshes) yield { val meshPoints = mesh.pointSet.points.toIndexedSeq // transformations include scaling LandmarkRegistration.similarity3DLandmarkRegistration(meshPoints.zip(fGpaMeanMeshPoint), center = Point3D(0, 0, 0)) } //full GPA fits val fGpaTransformedMeshes = allMeshes.zip(fullGpaTransforms).map (z=> z._1.transform(z._2)) val fGpaTransformedMesheViews = fGpaTransformedMeshes.map(mesh => ui.show(fGpaGroup, mesh, "trMesh")) fGpaTransformedMesheViews(0).color = new Color(85,235,40) fGpaTransformedMesheViews(1).color = new Color(20,40,250) fGpaTransformedMesheViews(2).color = new Color(250,150,20)

Here are the views:

Again, the centroid sizes can be calculated as:

original mesh 1: 487.222, mesh 1 full-GPA-transformed: 418.329

original mesh 2: 692.725, mesh 2 full-GPA-transformed: 417.958

original mesh 3: 566.668, mesh 3 full-GPA-transformed: 419.082

mean mesh: 421.707

As it can be observed, the original meshes are scaled in the case of the full GPA analysis to fit the mean mesh. The sizes of all the transformed and the mean mesh are close to each other, because, scaling is part of the transformations fitting the meshes onto a common mean. Also it should be noted that the (relative) scales of the meshes are not preserved.

Saving the mean meshes and the fits into files and closing the object:

MeshIO.writeMesh(meanMesh, new java.io.File(drive + "registration/Partial GPA/" + "meanMesh.stl")) (0 until 3).foreach(i => MeshIO.writeMesh(transformedMeshes(i), new java.io.File(drive + "registration/Partial GPA/" + s"mesh${i}.stl"))) MeshIO.writeMesh(fGpaMeanMesh, new java.io.File(drive + "registration/Full GPA/" + "meanMesh.stl")) (0 until 3).foreach(i => MeshIO.writeMesh(fGpaTransformedMeshes(i), new java.io.File(drive + "registration/Full GPA/" + s"mesh${i}.stl"))) }

## Full GPA and Shapes

Until this moment, I did not refer to a mesh/surface as a *shape*, because a shape has its own technical definition. The shape of an object/configuration/mesh is all its geometrical information modulo its scale, location, and rotational effects. In other words, **shape** is an **equivalence class of geometrical objects** modulo the similarity transformation. The full GPA removes the scale, location, and rotational effects of objects/meshes (relative to a mean object/mesh) in order to fit them all onto the mean mesh, which is the **sample full Procrustes mean shape** . Therefore, whatever geometrical variations that remain between the mean shape and each of the full Procrustes fits (full-GPA transformed meshes) reflect the variations due to the differences/dissimilarities in their **shapes**.