SSM Using Scalismo (4): Generalized Procrustes analyses

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)
Fig. 1. Meshes and views of their cross-sectional contours.

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)
Fig. 2. Mean mesh (red) and views of cross-sectional contours of the mean and the original meshes.

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"))
Fig. 3. The original and the partial-GPA fits.

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:

Fig. 4. full-GPA Mean mesh (red) and views of cross-sectional contours of the mean and the original meshes.
Fig. 5. The original meshes and their full-GPA fits.

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 \hat\mu_F. 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.