SSM Using Scalismo (5): Creating a point distribution model

Finally, I arrive at statistical shape models. I am going to create a point distribution model (PDM) of shapes. Reading some theory on the point distribution models, I need superimposed shapes and their mean shape. Superimposed shapes and their mean shapes are available from my previous tutorial on GPA with Scalismo. The outcomes of the full GPA is going to be utilized for creating a PDM. The PDM in this tutorial uses shapes (not forms), and it is based on the multivariate normal distribution of the points’ coordinates collected in a shape vector (Eq. 16 in this post).

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.

I start with initializing Scalismo and loading the superimposed configurations, i.e. shapes. Following the previous tutorial, I use the femoral head shapes.

object PDM 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 scala.io.StdIn.{readLine, readInt}

  scalismo.initialize()
  implicit val rng = scalismo.utils.Random(42)
  // ui
  val ui = ScalismoUI()

  // Loading data
  val drive = "I:/"
  val meshDir = drive + "registration/Full GPA/"
  val unsortedMeshFiles = new java.io.File(meshDir).listFiles()
  unsortedMeshFiles.foreach(file => println(s"\n$file"))
  println(unsortedMeshFiles)
  // specifying the file name of the mean mesh generated through the full GPA
  val meanMeshName = "meanMesh.stl"
  val meshFileNames = unsortedMeshFiles.map(file => file.getName())
  val indexOfMeanMesh = meshFileNames.indexOf(meanMeshName)
  // bringing the meanMesh to the head (first entry) of the meshFiles sequence
  val meshFiles = unsortedMeshFiles.updated(0, unsortedMeshFiles(indexOfMeanMesh)).updated(0, unsortedMeshFiles(indexOfMeanMesh))

I have to choose a reference mesh and some training meshes. The full GPA mean mesh is chosen as the reference mes and the three (registered and superimposed) meshes as the training meshes.

// specifying the reference mesh and the training meshes; the full GPA mean mesh is set to be the reference mesh
  val refMesh = MeshIO.readMesh(meshFiles.head).get
  val trainingMeshes = meshFiles.tail.map(file => MeshIO.readMesh(file).get)//the other 3 meshes

Once a reference mesh and training meshes are chosen, the deformation vector from each point of the reference mesh to its corresponding point of each training mesh should be calculated. Thereby, the deformation fields morphing the reference mesh to the training meshes are obtained.

 // [The meshes are already aligned using the full GPA, and in correspondence, so we can create a discrete GP]
  // creating the deformation vector field
  val deformationFields = trainingMeshes.map(mesh => {
    val deformationVectors = refMesh.pointSet.pointIds.map { id =>
      mesh.pointSet.point(id) - refMesh.pointSet.point(id)
    }.toIndexedSeq
    // a discrete deformation field needs a domain and the deformation vectors at each point of the domain
    DiscreteField[_3D, UnstructuredPointsDomain[_3D], EuclideanVector[_3D]](refMesh.pointSet, deformationVectors)
  })

The calculated deformation fields are samples (realizations) of a discrete random field. The domain of the random field is the reference mesh points. The random field of deformations is assumed to have a multivariate normal distribution. Here, the discrete random field is expressed by a random vector. Thereby, a shape, described by a shape vector, can be produced by adding the shape vector of the reference shape to the normally distributed random vector of deformations expanded by the deformation modes (this is explained by Eq. 16 in this post).

Technically, Scalismo works with the concept of Gaussian process (GP) which is a continuous random field whose any finite sub-collection of random variables (here deformations) has a multivariate normal (Gaussian) distribution. Therefore, continuous sample deformation fields should be approximated based on the known discrete realizations/samples of the random field. By interpolating a discrete field, Scalismo creates a continuous field. Different interpolation approaches are possible; Scalismo uses the Nearest Neighbor Interpolation method. The following code block defines a discrete GP (of the deformations) which is in fact a normally distributed random vector with elements indexed by the reference mesh points (a discrete domain). The discrete GP (random vector of deformations) is expressed in terms of its principal axes (modes). This is done using the PCA as described in this post (see Eq. 16).

  val interpolator = NearestNeighborInterpolator[_3D, EuclideanVector[_3D]]()
  val continuousDefFields = deformationFields.map(field => field.interpolate(interpolator))
  val gp = DiscreteLowRankGaussianProcess.createUsingPCA(refMesh.pointSet, continuousDefFields)

As I have 3 training meshes and therefore 3 training deformation fields, the rank of the sample covariance matrix of deformations is at most 3. This means that at most 3 principal axes are used to express the random vector of deformations.

I can sample the GP (random field) of deformations. By adding a sample deformation field to the reference mesh, I can generate a shape vector. This shape vector is a sample mesh from meshes whose points’ coordinates are normally distributed. Sampling some meshes from the mesh distribution, I’ll have the following realization:

  val PDMGroup = ui.createGroup("refMeshGroup")
  val refMeshView = ui.show(PDMGroup, refMesh, "refMesh")
  refMeshView.color = java.awt.Color.GREEN
  val gpDefView = ui.addTransformation(PDMGroup, gp, "RefMeshDefrmdByGp")
Fig. 1. some samples from the statistical point distribution model.
Fig. 1. some samples from the statistical point distribution model. One of the coefficients of the shape modes attains a relatively large value (close to its three standard deviation).

The horizontal bar-like controls in the figures are proportional to the coefficients of the deformation modes (b_i=\alpha_i\sqrt{\lambda_i}). Each coefficient varies within three standard deviations (Eq. 4 in this post).

Scalismo assembles the reference mesh and the discrete GP as a statistical mesh model:

// [creating a statistical shape model (PDM)
  val model = StatisticalMeshModel(refMesh, gp.interpolate(interpolator))

This is the statistical PDM.

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.

SSM Using Scalismo (3): Non-rigid registration

Non-rigid registration

The rigid registration preserves the shape and size (form) of objects by merely transforming them through rotation and translation. Unlike the rigid registration, non-rigid registration (NRR) deforms (morphs) an object/configuration to have it fitted onto another object. Here, an object is a valid surface that its geometry can be discretely represented by a triangular mesh. A configuration is the set (or subset) of points/nodes of the mesh of the object. In a word, the geometry of an object is closely approximated with a finite set of points and triangles connecting the points/nodes; all together a triangular mesh. The triangles linearly interpolate the actual geometry between the points/nodes. I am going to refer to a points of a triangular mesh as a node.

But why we want to morph an object onto another object. One application is to establish one-to-one correspondence between the configurations’ points. For statistical shape modelling/analysis correspondence is essential. Another difference between the rigid and NRR is that the rigid registration entails known landmark correspondence; however, NRR seeks correspondence. Briefly, the registration problem (NRR) is the problem of establishing correspondence between two objects (their points).

Imagine two objects having configurations/meshes X and Y with different number of points; one with m and the other n points. A 2D example is depicted below. Definitely, a bijective correspondence cannot be established between every pair of points of X and Y. So, we cannot think about some point-to-point solution. The correct way of establishing the correspondence between the geometry of two objects is to morph the mesh of one of them onto the other one. In the figure below, the configuration/mesh X is morphed onto the configuration/mesh Y. Thereby, X’ closely matches Y and also its points are in correspondence with X. We have to pay attention to the difference between each two of these three terms, and their individual role in the problem: reference mesh, target mesh, fitted mesh; here, X, Y, and X’ respectively.

Terms: By deformation of an object we mean the deformation vectors or field. By transformation of an object we mean the transformed geometry.

Here, two questions bug my brain (there must be some brain there..!): 1- how to morph?, 2- how much morph? Devising and utilizing a model simulating plausible deformations is the answer to the first question, and defining a distance function to measure the distance between two meshes answers the second one.

Fig. 1. Two configurations with different number of points; one morphs to the other one.

A reference mesh \Gamma_R also acting as the domain of the functions, a target mesh \Gamma_T, a model of deformation fields u, a distance function D, and a weighted regularization term on the deformations formulate the registration problem as the following optimization problem [1]:

\hat u=\argmin_u (D[\Gamma_R,\Gamma_T,u]+\eta R[u])\tag{1}

The deformation field is a function u:\mathbb{R}^3\to\mathbb{R}^3. The distance function eats two (constant) configurations, the reference and the target ones, and also the deformation field, then it computes the distance between the target mesh and the deformed reference mesh, i.e. \Gamma = \Gamma_R+u(\Gamma_R). The regularization term R involves the prior assumption. It penalises deformations that do not fit the the prior assumption [1]. This method is called parametric NRR.

The theory of parametric NRR is the content of another post. Here, I try to explain the practical steps toward putting the theory into practice with Scalismo.

1- Initialising Scalismo

In addition to the main code, I defined some functions and organised them in a package object. The package code will be at the end of the post. How to create a package, click here.

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 NRR extends App {
  //To send the code selection to Scala Console, select the code in the Editor, press Ctrl + Shift + X
  //imports
  import scalismo.geometry._
  import scalismo.common._
  import scalismo.ui.api._
  import scalismo.mesh._
  import scalismo.registration.LandmarkRegistration
  import scalismo.io.{MeshIO, LandmarkIO}
  import scalismo.numerics.UniformMeshSampler3D
  import breeze.linalg.{DenseMatrix, DenseVector, csvread}
  import scalismo.registration.{Transformation, RotationTransform, TranslationTransform, RigidTransformation}
  import scala.io.StdIn.{readLine, readInt}
  import scalismo.numerics._
  import scalismo.kernels._
  import scalismo.statisticalmodel._
  import scalismo.registration._
  scalismo.initialize()
  implicit val rng = scalismo.utils.Random(42)
  // ui
  val ui = ScalismoUI()

This is to include the package:

import custom.methods._

2- Loading meshes and landmarks

We use our rigidly aligned meshes of the femoral heads and their anatomical landmarks; available in the post SSM Using Scalismo (1): Rigid alignment. The initial rigid alignments of the meshes and some landmarks in correspondence are going to mitigate the registration problem. This time, we intend to fit/morph (NRR) the reference mesh to the target meshes.

 // [Loading aligned meshes and LM files]
  val meshDir = "I:/registeration/AlignedData/Meshes/"
  val LMDir = "I:/registeration/AlignedData/LMs/"
  val registeredMeshDir = "I:/registeration/Registered/"
  // meshes
  val meshFiles = new java.io.File(meshDir).listFiles()
  meshFiles.foreach(file => println(s"${file}\n"))

  val numOfMeshes = meshFiles.size
  println(s"number of files is ${numOfMeshes}")

  // locating the RefMesh file i.e mesh0.stl. RefMesh has to be named mesh0.stl
  val refMeshIndex : Int = (0 until numOfMeshes).map(i => meshFiles(i).getName).indexOf("mesh0.stl")
  if (refMeshIndex != 0) {println("Rename refMesh File"); sys.exit() }
  val refMesh = MeshIO.readMesh(meshFiles(refMeshIndex)).get

  val targetMeshes: Array[TriangleMesh[_3D]] = meshFiles.tail.map(meshFile => MeshIO.readMesh(meshFile).get) //meshes to get fitted to refMesh
  // Loading Landarks. 
  val LMFiles = new java.io.File(LMDir).listFiles()
  // Landmarks of the reference mesh should be in LMs(0)
  val LMs = LMFiles.map(f => LandmarkIO.readLandmarksJson[_3D](f).get)

A reference mesh view if you want:

 val refMeshGroup = ui.createGroup("refMeshGroup")
 val refMeshView = ui.show(refMeshGroup, refMesh, "refMesh")
 refMeshView.color = java.awt.Color.GREEN
Fig. 1. The reference mesh (dimensions are added manually)

3- Defining a Gaussian process

The deformation field is going to be modelled by a Gaussian process (GP). GP outputs plausible/permissible deformation components at each point of the reference mesh. To define a GP we need a mean deformation field and a kernel (this post). The GP should produce deformation fields capturing the deformations needed to morph the reference mesh onto each of the target meshes. Looking at the reference mesh, we see a combination of global and local geometrical features; a dome-like shape with an indentation. So, maybe it is better to define a kernel taking different scale of deformations into account. The mean field is set to be a zero field.

val mean = VectorField(RealSpace[_3D], (_ : Point[_3D]) => EuclideanVector.zeros[_3D]) //
  val kernel1 = DiagonalKernel[_3D](GaussianKernel(sigma = 3) * 0.1, outputDim = 3)
  val kernel2 = DiagonalKernel[_3D](GaussianKernel(sigma = 7) * 2.0, outputDim = 3)
  val kernel = kernel1 + kernel2
  val GP = GaussianProcess(mean, kernel)
  val lowRankGP = LowRankGaussianProcess.approximateGPCholesky(refMesh.pointSet,GP, 0.05, NearestNeighborInterpolator())

We can view samples of deformations of the reference mesh by the GP to check whether the defined GP kernel is suitable of capturing the needed deformations or not.

 val gpSampleGroup =  ui.createGroup("gpSampleGroup")
 val refMeshView2 = ui.show(gpSampleGroup, refMesh, "refMesh")
 // Transforming whatever mesh in the gpSampleGroup by the GP-sampled deformation fields
  val gpDefView = ui.addTransformation(gpSampleGroup, lowRankGP, "RefMeshDefrmdByGp")
  removeGroupView(gpSampleGroup) // see the definition in the package
Fig.2. GP-Sampled deformations on the reference mesh.

The deformations seems ok; nothing weird like intersecting triangles occurred. Let’s hope that the kernel is going to work and continue.

4- Defining regression data for posterior GPs

Now, it is the moment that we use the known landmark correspondence information. On each mesh there are 3 (inter-mesh) corresponding landmarks. Therefore, there are known deformation vectors from the landmarks on the reference mesh to their corresponding landmarks on each of the target mesh. We want to incorporate these known deformations into our GP deformation model. This can be done through GP regression, and the resultant model, with the incorporation of the known deformations, is called the posterior GP model. The theory is touched in this post.

The landmarks have IDs and coordinates. The coordinates are not associated with the nodes. Therefore, we first have to associate/find a node closest to each landmark of the reference mesh. Here is how:

//[Defining posterior lowrankGp given the landmark transformations from reference to a target mesh]
  // finding points (nodes) corresponding/associated to/with landmarks of the refMesh
  // The nodes on the mesh closest to the landmarks are the nodes associated with the landmarks
  val refLMs = LMs(0)
  val refLMNodeIds = refLMs.map(eachLM => refMesh.pointSet.findClosestPoint(eachLM.point).id)

Generating the regression data, i.e. known deformation vectors (from the reference mesh to each target mesh).

 // generating regression data i.e known deformations from the LMs' corresponding points on the RefMesh to corresponding landmark points of the target mesh
  //landmarks are already sorted in the order of correspondence
  val knownDeformations =  (0 until numOfMeshes-1).map(i => {
    (0 until LMs(0).size).map(j => (LMs(i+1)(j).point - refMesh.pointSet.point(refLMNodeIds(j))))
  })
  // obtaining the coordinates of the refMesh points associated with the refMesh landmarks
  val refLMpoints = refLMNodeIds.toIndexedSeq.map(eachId =>refMesh.pointSet.point(eachId))

The known deformation fields (regarding each target mesh) at the reference nodes (associated with the reference landmarks) can be viewed:

// visualizing known deformation fields
  // generating a domain of points needed for a vector field definition/construction
  val refLMdomain =  UnstructuredPointsDomain(refLMpoints)
  val knownDeformationFields = knownDeformations.map(eachDefSet => {
    DiscreteField[_3D, UnstructuredPointsDomain[_3D], EuclideanVector[_3D]](refLMdomain, eachDefSet)
  })
  val knownDeformationGroup = ui.createGroup("knownDeformations")
  knownDeformationFields.map(eachfield => ui.show(knownDeformationGroup, eachfield, "knownDefs"))
  removeGroupView(knownDeformationGroup)
Fig. 3. Deformation vectors (fields) on the nodes (the landmarks) of the reference mesh. There are two known deformation vectors at each node (associated with a landmark) because there are two target meshes.

Because the known/observed deformations is usually polluted with sort of uncertainty, a noise term should be added onto the known deformations. We add a little bit of noise as:

val noise = MultivariateNormalDistribution(DenseVector.zeros[Double](3),0.001* DenseMatrix.eye[Double](3))

Each set of regression data (known/observed deformation fields plus noise) regarding each target mesh contains sequences of data in the tuple form of (data point, deformation, noise). Data point is the a point of the reference mesh at which the known noisy deformations are observed. We define the following function to generate sets of regression data regarding each target mesh. Then, we use a function to generate the sets of regression data regarding the target meshes.

def RegressionDataGenerator(dataPoints:IndexedSeq[Point[_3D]], deformations:IndexedSeq[EuclideanVector[_3D]],noise:MultivariateNormalDistribution)={
    val regressionDataSet = (0 until dataPoints.size).map(i => (dataPoints(i),deformations(i),noise))
    regressionDataSet //.toIndexedSeq
  }
// regression data sets regarding target meshes
  val regressionDataSets = (0 until numOfMeshes - 1).map(i =>  RegressionDataGenerator(refLMpoints,knownDeformations(i),noise))

5- Posterior GP and registration

The following steps are going to be taken for each target mesh, so I use a for loop:

for( n <- (0 until numOfMeshes - 1)) {

5-1- Defining the Posterior GP based on each regression data

 val posteriorGP = lowRankGP.posterior(regressionDataSets(n))

Let’s view some samples of transformation of the reference mesh deformed by the posterior GP based on one of the target meshes regression data:

 val postGpGroup = ui.createGroup("postGpGroup")
 val refMeshView3 = ui.show(postGpGroup, refMesh, "refMesh")
 val PGpView = ui.addTransformation(postGpGroup, posteriorGP, s"TtransformedByPGp $n")
 removeGroupView(postGpGroup)
Fig. 4. Sampled transformation of the reference mesh deformed by the posterior GP.

5-2- Solving the NRR problem

To solve an NRR problem based on Eq. 1, an agent modelling the deformation fields, a distance function, a regularizer, and the weight of the regularizer are needed. Scalismo reduces the NRR problem into the following minimization problem [1]:

\argmin_{\alpha_1,…,\alpha_r}D\Big[\Gamma_R,\Gamma_T,\mu+\sum_{i=1}^{r}\alpha_i\sqrt{\lambda_i}\phi_i\Big]+\eta \sum_{i=1}^{r}\alpha_i^2\tag{2}

where, \alpha_i‘s are registration coefficients minimizing the above expression, \mu is the average of the GP (here zero every where), \lambda_i‘s and \phi_i‘s are Eigen values and Eigen functions of the low-rank approximation of the GP (see this post), r is the rank of the low-rank approximation of the GP .These Eignen values/functions are known as the GP or the posterior GP is known. \eta is the weight of the regularization term \sum_{i=1}^{r}\alpha_i^2. The regularization term is an L2 regularizer.

The distance function measuring the distance between the transformed reference mesh \Gamma and the target mesh can be defined as the mean squared distance. For a moment, let’s assume that the surfaces are continuous. A point-to-point distance function is defined as:

\text{CP}_{\Gamma_T}(x):=\|x-\argmin_{x’ \in \Gamma_T} \|x-x’\|\|

where x\in\Gamma a surface, and \|.\| is the Euclidean vector norm. This function calculates the distance between each point on \Gamma and its closest point on \Gamma_T. The arg min operator finds the closest point on \Gamma_T.

The average of the distance function over the continuous domain \Gamma_R represents the average distance between \Gamma and \Gamma_T:

D[\Gamma_R, \Gamma_T,u]:=\frac{1}{S(\Gamma_R)}\int_{\Gamma_R}\big (\text{CP}_{\Gamma_T}(x+u(x))\big )^2 ds

in which S(\Gamma_R) is the area of \Gamma_R. This distance measure goes in Eq. 2 and is evaluated numerically. It is good to know that numerical evaluation of the above integral considering discrete surfaces, i.e. meshes, gets the following form:

D[\Gamma_R, \Gamma_T,u]\approx \frac{1}{n} \sum_{i=1}^{n}\big (\text{CP}_{\Gamma_T}(x_i+u(x_i))\big )^2

such that x_i\in\Gamma, n is the number of nodes/points of \Gamma_R; here the surfaces are represented by meshes. The integral can be evaluated at  m\lt n number of sampled points/nodes of \Gamma_R; in this case we replaces n with m.

Minimizing Eq. 2 is through iteration. To start, we need initial values of the coefficients \alpha_i, number of sampled points/nodes for numerical evaluation of the distance measure, initial value for the regularization weight, and number of iterations.

The initial values of \alpha_i are set to be zeros by a \mathbb{R}^{r\times 1} vector:

val initialCoefficients = DenseVector.zeros[Double](posteriorGP.rank)

Now, instead of working with fixed values for number of sampled points, regularization weight, and number of iterations, the iteration procedure is divided into several stages in order to fine-tune the minimization problem. Physically, finer geometrical features (up to the capability of the GP) are considered at at each stage. To this end, we define a sequence of registration parameters, which are, the reqularization weight, number of iterations, and number of sampled points/nodes. This is how it looks like:

 val registrationParameters = Seq(
      RegistrationParameters(regularizationWeight = 1e-1, numberOfIterations = 20, numberOfSampledPoints = 1000),
      RegistrationParameters(regularizationWeight = 1e-2, numberOfIterations = 30, numberOfSampledPoints = 1000),
      RegistrationParameters(regularizationWeight = 1e-4, numberOfIterations = 40, numberOfSampledPoints = 2000),
      RegistrationParameters(regularizationWeight = 1e-6, numberOfIterations = 50, numberOfSampledPoints = 4000)
    )

where, “RegistrationParameters” is a case class; see the defining code in the package file at the end of the post.

We consider Eq. 2 to clarify the above approach. When the regularization weight is relatively large, the minimization problem leads to some relatively small coefficients, therefore modes of deformations associated with relatively small Eigen values in \sum_{i=1}^{r}\alpha_i\sqrt{\lambda_i}\phi_i, are ruled out. This means only modes having large effects in the deformations remain. These modes of deformations are more toward global deformations rather than local deformations leading to fine geometrical features. Having said that, it is reasonable to start the minimization procedure with less number of points, uniformly sampled from our discrete domain \Gamma_R, and fewer number of iterations. Thereby, a set of coefficients appropriate for an approximate global fit (of the reference mesh to the target mesh) is calculated. Then, the regularization weight is decreased so that other deformation modes with relatively smaller Eigen values can effectively appear in the minimization. This means taking care of more local deformations hence finer geometrical features. As we want to consider smaller deformations, we have to fine-sample the points, i.e. more number of points. At this level, number of iterations should also increase.

The output of the minimization is the alpha-coefficients. The procedure starts with an initial value and after each (sub) iteration the input coefficients are updated. Also, once each stage of iteration as specified above is done, its output coefficients should travel to the start of the next stage. This can be done by a single line in Scala:

 val finalCoefficients = registrationParameters.foldLeft(initialCoefficients)((modelCoefficients, regParameters) =>
      doRegistration(posteriorGP, refMesh, targetMeshes(n), regParameters, modelCoefficients))

The function doRegistration is elucidated in the package file.

Finally, the minimizing coefficients are calculated and stored in finalCoefficients.

6- The fitted mesh

The coefficients are now used to deform the reference mesh (all points) onto a fitted mesh approximating the target mesh. The following two lines calculates the deformation field based on \mu+\sum_{i=1}^{r}\alpha_i\sqrt{\lambda_i}\phi_i with the calculated \alpha_i‘s.

 val transformationSpace = GaussianProcessTransformationSpace(posteriorGP)
 val registrationTransformation = transformationSpace.transformForParameters(finalCoefficients)

Applying the transformation on the reference mesh to get the fitted mesh:

 val fittedMesh = refMesh.transform(registrationTransformation)

7- Visualising the fitted mesh

Here, the fitted mesh on to each target mesh is visualised.

 val fittedMeshGroup =  ui.createGroup("fittedMeshGroup")
 val fittedMeshView = ui.show(fittedMeshGroup, fittedMesh, "fittedMesh")
 val targetMeshView = ui.show(fittedMeshGroup, targetMeshes(n), s"TargetMesh $n")

The fitted mesh on to the first target mesh:

Fig.5. Fitted mesh tot he first target mesh.

To visually assess the accuracy of the fit, we can scroll through the cross sectional contours (normal to the coordinate system’s axes) of the fitted and the target meshes simultaneously. For example:

Fig. 6. Some cross sectional contours.

Based on the visualization of the cross sectional contours, the fitted mesh closely approximates the target mesh except at some end cross sections where local geometrical variations are relatively steep.

The other target mesh is approximated/fitted as below:

Fig.7. Fitted mesh tot he second target mesh.

The worst mismatches of cross sectional contours for this case.

Fig. 8. Some cross sectional contour mismatches.

To accept or not to accept a fitted mesh; this is the question.

It depends on how much discrepancy is allowed in a particular application. Nonetheless, life goes on regardless of our choice. If we are picky, we can play with the kernel (and the registration parameters) to get a better fit.

Average distance and the Hausdorff distance are available as quantitative measures of fit, in Scalismo:

 val AveDist  = scalismo.mesh.MeshMetrics.avgDistance(refMesh, fittedMesh)
 val HsdrfDist  = scalismo.mesh.MeshMetrics.hausdorffDistance(refMesh, fittedMesh)
 println(s"\nAverage distance is $AveDist")
 println(s"\nHsdrf distance is $HsdrfDist \n")

Distance of first target mesh to its fitted mesh:
Average: 1.2 mm
Hausdorff: 3.9 mm
Distance of second target mesh to its fitted mesh:
Average: 1.4 mm
Hausdorff : 4.9 mm

8- Writing the fitted mesh into a file

Before writing the fitted meshes to files, we can observe that the reference mesh and each fitted mesh have the same number of nodes (here, 4470); indeed all in (one-to-one) correspondence. Use the following to access the attribute containing the number of points of the reference mesh (or any other ones).

 refMesh.pointSet.numberOfPoints

Writing to STL files:

 val fittedMeshFile =  new java.io.File(registeredMeshDir + s"ftdMesh$n.stl")
 MeshIO.writeSTL(fittedMesh, fittedMeshFile)

9- Finishing the code

We close the for-loop and use some lines to exit the run:

// to remove the groupp view after each fitting procedure is done 
 removeGroupView(fittedMeshGroup) 

 } //closes the for-loop

 val finito = readLine("finito [y/n]? ")
 if (finito == "y") ui.close()

} // closes the object

[1] Gaussian Process Morphable Models, Marcel Luthi , Thomas Gerig, Christoph Jud, and Thomas Vetter. IEEE Transactions on Pattern Analysis and Machine Intelligence (Volume: 40 , Issue: 8 , Aug. 1 2018).

Appendix: package

package custom

import scalismo.geometry._
import scalismo.common._
import scalismo.ui.api._
import scalismo.mesh._
import scalismo.registration.LandmarkRegistration
import scalismo.io.{MeshIO, LandmarkIO}
import scalismo.numerics.UniformMeshSampler3D
import breeze.linalg.{DenseMatrix, DenseVector, csvread}
import scalismo.registration.{Transformation, RotationTransform, TranslationTransform, RigidTransformation}
import scala.io.StdIn.{readLine, readInt}
import scalismo.numerics._
import scalismo.kernels._
import scalismo.statisticalmodel._
import scalismo.registration._

package object methods {
  implicit val rng = scalismo.utils.Random(42)
  // removing group view
  def removeGroupView(groupView:scalismo.ui.api.Group) : Unit = {
    val s = readLine(s"remove ${groupView.name} and/orJust continue [y/n]? ")
    if (s == "y") groupView.remove()
  }

  // registration

  case class RegistrationParameters(regularizationWeight : Double, numberOfIterations : Int, numberOfSampledPoints : Int)

  def doRegistration(
                      lowRankGP : LowRankGaussianProcess[_3D, EuclideanVector[_3D]],
                      referenceMesh : TriangleMesh[_3D],
                      targetMesh : TriangleMesh[_3D],
                      registrationParameters : RegistrationParameters,
                      initialCoefficients : DenseVector[Double]
                    ) : DenseVector[Double] =
  {
    val transformationSpace = GaussianProcessTransformationSpace(lowRankGP)
    val fixedImage = referenceMesh.operations.toDistanceImage
    val movingImage = targetMesh.operations.toDistanceImage
    val sampler = FixedPointsUniformMeshSampler3D(referenceMesh, registrationParameters.numberOfSampledPoints)
    val metric = MeanSquaresMetric(fixedImage,movingImage,transformationSpace,sampler)
    val optimizer = LBFGSOptimizer(registrationParameters.numberOfIterations)
    val regularizer = L2Regularizer(transformationSpace)
    val registration = Registration(
      metric,
      regularizer,
      registrationParameters.regularizationWeight,
      optimizer
    )
    val registrationIterator = registration.iterator(initialCoefficients)
    val visualizingRegistrationIterator = for ((it, itnum) <- registrationIterator.zipWithIndex) yield {
      println(s"object value in iteration $itnum is ${it.value}")
      //gpView.coefficients = it.parameters
      it
    }
    val registrationResult = visualizingRegistrationIterator.toSeq.last
    registrationResult.parameters
  }

}

SSM Using Scalismo (2): GP mean and kernel

Statistical shape modeling (SSM) using Scalismo is based on Gaussian process (a random field. A Gaussian process, like normal probability distribution, is completely defined by its mean ( µ ) and covariance function ( k ). The covariance function is called kernel. We write GP(µ , k). The mean and kernel are both functions. In SSM, GP models shape deformations (a vector field) in a 3D (or 2D) space. It means that the GP tells us the probable or plausible deformation vector at each point in the space. The mean function eats a point in space, and output a real vector of deformations; kernel eats two points in the space and output a matrix. If the points in our 3D geometric space are represented by their coordinates, the aforementioned functions are:

\mu:\mathbb{R^3}\to\mathbb{R^3} \newline k:\mathbb{R^3}\times\mathbb{R^3}\to\mathbb{R^{3\times3}}\tag{1}

Once the mean and the kernel functions are determined, we can construct our SSM based on the GP. Statistical shape models based on GPs are called Gaussian Process Morphable Models (GPMM) [1]. A GPMM is composed of a GP describing the deformation field, and a reference shape. A reference shape is a subset of our geometric space. Here, the geometric space is a set of point. Describing the deformations through a GP, we can have the GP produce plausible deformations of the reference shape. The plausibility of the deformations is determined by the kernel. Hence, a shape produced/described by a GPMM is mathematically defined as [1]:

\Gamma_S=\{x+u(x)|x\in\Gamma_R\subset\Omega\text{ , }u:\Omega\to\mathbb{R^3}\text{ and }\Omega\supset\Gamma_R\}\tag{2} \newline \Omega \text{ open in }\mathbb{R^3}

The set of points of any plausible shape \Gamma_S is simply formed through giving each point of the reference shape a new location in the space. The new location of each point is obtained by the map x+u(x). The map eats a point, translate it by (along) the vector u(x), and output a point. Although both points and vectors are analytically recognized by tuples in \mathbb{R^3}, I need to pay attention to the conceptual difference between a point and a vector in the space.

Where is the GP? Here, where the GP comes into the scene. All plausible functions (maps), u(x), are distributed according to a GP distribution:

u(x)\sim\mathnormal{GP}(\mu(x),k(x,x’))\text{ s.t }x\in\Omega\tag{3}

The mean (average) shape of the sets of all plausible shapes is naturally defined as:

\overline{\Gamma}=\{x+\mu(x)|x\in\Gamma_R\}\tag{4}

If, for simplicity, we let \mu(x)=0, we’ll get \overline{\Gamma}=\Gamma_R.

Now, the kernel. The kernel is actually the covariance function measuring the variability or correlation of the components of two deformation vectors at two points. As each deformation vector has 3 components, there will be 9 covariance values organized in a 3 by 3 matrix:

k(x,x’):=cov(u(x),u(x’))\tag{5}

where u(x) is a vector in  \mathbb{R^3} . If we consider u(x) (and \mu(x)) as a column vector in the space of 3 by 1 matrices (  \mathbb{R^{3\times1}} ), we have:

cov(u(x),u(x’))=E[(u(x)-\mu(x))(u(x’)-\mu(x’))^\text{T}]\tag{6}

With (.)^T being the transpose operator. Expanding the above, we’ll get the matrix of k(x,x') with the components:

k_{i,j}(x,x’)=cov(u_i(x),u_j(x’))=E[(u_i(x)-\mu_i(x))(u_j(x’)-\mu_j(x’))] \newline\tag{7}

Or in an expanded form with some (easy to capture) change of notations:

k_{i,j}(x,x’)=\begin{bmatrix} cov(u_x^1,u_{x’}^1)&cov(u_x^1,u_{x’}^2)&cov(u_x^1,u_{x’}^3)\\[0.3em] cov(u_x^2,u_{x’}^1)&cov(u_x^2,u_{x’}^2)&cov(u_x^2,u_{x’}^3)\\[0.3em] cov(u_x^3,u_{x’}^1)&cov(u_x^3,u_{x’}^2)&cov(u_x^3,u_{x’}^3)\end{bmatrix}\tag{8}

We can observe that the function k(x,x') is symmetric in the sense k(x,x') = k(x',x) . Why? Because covariance is a symmetric function, i.e. cov(a,b) = cov(b,a) for two random variables (NOT random vectors), a and b.

Remark: symmetry of k(x,x') does not imply that the matrix [k_{i,j}(x,x')] is symmetric as well. If x \neq x', then it is not necessarily true that cov(u_x^1,u_{x'}^2) \neq  cov(u_x^2,u_{x'}^1) , and so on. Note that u_x^1 , u_{x'}^2 , u_x^2 and u_{x'}^1 are four different random variables.

However, the kernel leads to a symmetric matrix when evaluated at the same point, i.e. k(x,x), which is equal to cov(u_i(x),u_j(x))

What sort of kernel should we use?

Let’s define a scalar Gaussian kernel  k_s:\mathbb{R^3}\times\mathbb{R^3}\to\mathbb{R} as [1]:

k_s(x,x’):=s\ \exp({-\frac{\|x-x’\|^2}{\sigma^2}})\tag{9}

where  \|.\| is the Euclidean norm expressing the distance between two (geometric) points in the space, s is a scale, and \sigma^2 is a parameter defining the distance range over which the random variables are correlated. This kernel only depends on the distance of two points and not their locations. This means the correlation between the deformations at two points decreases as their distance increases.

One choice of generating the kernel matrix regarding the deformations is to assume that the deformation components are not correlated. Since we are working with a GPs, being uncorrelated implies independence. In mathematical notations:

cov(u_x^i,u_{x’}^j)=0\text{ for }i\neq j\text{ s.t }i,j\in\{1,2,3\}\tag{10}

which implies:  \rho(u_x^i|u_{x'}^j)=\rho(u_x^i) for  i\neq j , where \rho(.) denotes the related probability distribution. As a side note  u_x^i  and  u_{x'}^j are random variables, therefore above follows this: \rho_{X|Y}(x|y)=\rho_{X}(x) for any random variables x and y.

Now we use the scalar kernel as the building component of our matrix kernel and define a diagonal kernel with the same parameters in each direction (may be I can call it a diagonal kernel with homogeneous parameters):

k_{i,j}(x,x’)=\begin{bmatrix} k_s(x,x’)&0&0\\[0.3em] 0&k_s(x,x’)&0\\[0.3em] 0&0&k_s(x,x’)\end{bmatrix}\tag{11}

Let’s puts all these into practice. We consider a 1 by 1 (unit of distance) plane/plate meshed with triangular elements (STL). A zero-mean GP on the deformations describes all plausible/admissible deformations of the plane. The diagonal kernel with homogeneous parameters is considered. The following figures demonstrate samples of the deformations (sample path or trajectory of the GP) with different kernel’s parameters.

Fig.1. Sampled deformations of a plane from a GP with a diagonal kernel having different homogeneous parameters.

The kernel determines how much and how far (the radius of effect) the deformations at a point in the space should affect the deformations at its neighboring points. Consider the plots of k_s for different parameters (Fig. 2). Remember that k_s provides the covariance values. The higher the covariance, the higher the dependency of deformations. Keeping \sigma constant and increasing s, increases the dependency while keeping the radius of effect constant (Fig. 2 left image). At a point itself, increasing s, increases the variance of deformations (around its mean) at that point; hence, larger deformations are expected. On the other hand, keeping the scale, s, constant and increasing \sigma leads to a larger radius of effect while a constant variance. This means the possible range of deformations (around the mean) does not change at a point, but its deformations (the deformations at that point), now, affect the deformations at the surrounding points to a larger distance, i.e. larger radius of effect. To summarize it, s determines the variance (possible ranges) of deformation magnitudes at a point, and  \sigma controls how locally/globally the deformations at a point should affect its surrounding points’ deformations. Final word: here, working with GPs, zero kernel/zero covariance implies independency of deformations at two points; I deform this way, you deform the way you like. Non-zero kernel means I deform this way and you follow me as much as I say while deforming.

Fig. 2. Effects of the scalar kernel parameters on the kernel values. Horizontal axis: distance between two points, vertical axis: scalar kernel values.

Another interesting observation that we can made is to use the kernel to eliminate the deformations in specific directions. It can be done by setting the GP mean function to zero, as already done, and also set the variance of the components of the deformations to zero. A random variable vanishes when its mean and its variance are both zero. Let’s eliminate the deformations in x and y directions being the in-plain directions and let the plane deform only out-of-plane, i.e. orthogonal to itself (the z direction). This enforcement of deformations resembles plate bending. Here, is the kernel:

k_{i,j}(x,x’)=\begin{bmatrix} 0&0&0\\[0.3em] 0&0&0\\[0.3em] 0&0&k_s(x,x’)\end{bmatrix}\tag{12}
Fig.3. Sampled deformations of a plane from a GP with a kernel that eliminates in-plane deformations.

Some power of kernel combination:

Sometimes it is needed to model different scale of deformations, i.e. modeling global and local deformations simultaneously. For example, consider a “global” bending of a piece of wiggly wire. The deformations are shown in Fig. 4.

Fig. 4. Bending of a piece of wiggly wire.

To consider different scale of deformation we should let the deformations at a point, say point P, strongly affect the deformations of the nearby neighboring points in order to model the local deformations around that point. Also, the deformations at P, should affect the deformations of farther points to some extent. This is achieved by adding kernels (element-wise in the case of a matrix). For example, the following figure shows how a GP can model a global bending of a plate on which there are localized spiky round bumps.

Fig. 5. Different scales of deformations: global bending and localized round bumps.

The kernel used to model the deformations (Fig. 5) follows Eq. 12. The graphs of the scalar functions and their sum are as follows:

Fig. 6. Scalar functions and their sum.

It may be interesting to observe different deformation modes, .i.e eigen functions. The following figure shows some of the modes:

Fig. 7. Some global and local deformation modes.

Here is the code:

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.

import scalismo.geometry._
import scalismo.common._
import scalismo.ui.api._
import scalismo.mesh._
import scalismo.registration._
import scalismo.io.{MeshIO}
import scalismo.numerics._
import scalismo.kernels._
import scalismo.statisticalmodel._
import breeze.linalg.{DenseVector}
import java.awt.Color

scalismo.initialize()
implicit val rng = scalismo.utils.Random(42)
val ui = ScalismoUI()

// loading our mesh file and let it be named refMesh
val refMesh = MeshIO.readMesh(new java.io.File("I:/Scalismo/plane.stl")).get

//viewing the refMesh
val modelGroup = ui.createGroup("model")
val refMeshView = ui.show(modelGroup, refMesh, "referenceMesh")
val refMeshColor = new Color(115,238,70)
refMeshView.color = refMeshColor
// creating the mean vector fields
val mean = VectorField(RealSpace[_3D], (_ : Point[_3D]) => EuclideanVector.zeros[_3D])

Kernel form regarding Fig. 1:

val kernel = DiagonalKernel[_3D](GaussianKernel(sigma = 0.8) * 0.1,3)

Kernel form regarding Fig. 3:

val scalarValuedGaussianKernel3 : PDKernel[_3D]= GaussianKernel(sigma = 0.3)* 0.1
val scalarValuedGaussianKernel2 : PDKernel[_3D]= GaussianKernel(sigma = 0.1)* 0
val scalarValuedGaussianKernel1 : PDKernel[_3D]= GaussianKernel(sigma = 0.1)* 0
val kernel = DiagonalKernel( scalarValuedGaussianKernel1, scalarValuedGaussianKernel2, scalarValuedGaussianKernel3)

Kernel regarding Fig. 5:

val s1 = 0.009
val s2 = 0.15
val diagDernel1 = DiagonalKernel(GaussianKernel(sigma = 1)*0,GaussianKernel(sigma = 1)* 0,GaussianKernel(sigma = 0.032)* s1)
val diagDernel2 = DiagonalKernel(GaussianKernel(sigma = 1)*0,GaussianKernel(sigma = 1)* 0,GaussianKernel(sigma = 0.8)* s2)
val kernel = diagDernel1 + diagDernel2

The rest of the code:

// defining the GP
val gp = GaussianProcess(mean,  kernel)
val lowRankGP = LowRankGaussianProcess.approximateGPCholesky(refMesh.pointSet,  gp,  0.07,  NearestNeighborInterpolator())

// visualizing and sampling
val deformedGroup = ui.createGroup("deformed")
val refMeshView2 = ui.show(deformedGroup, refMesh, "referenceMesh")
refMeshView2.color = refMeshColor
val gpDefView = ui.addTransformation(deformedGroup, lowRankGP, "RefMeshDefrmdByGp")

References:

[1] Gaussian Process Morphable Models, Marcel Luthi , Thomas Gerig, Christoph Jud, and Thomas Vetter. IEEE Transactions on Pattern Analysis and Machine Intelligence (Volume: 40 , Issue: 8 , Aug. 1 2018).

SSM Using Scalismo (1): Rigid alignment

Hereby, I am going to use Scalismo to rigidly register/align/fit objects represented by triangular meshes onto a particular object. The alignment is based on partial ordinary Procrustes analysis method (POPA). POPA involves translation and rotation (not scaling) of an object with respect to another in order to align/fit the transforming object to the stationary one as close as possible. Read about POPA here. This initial alignment makes the non-rigid registration easier and faster later.

There are 3 geometries of femoral heads in the form of STL and I select one as the reference object (reference mesh) and I will make the other 2 aligned/registered on the reference one.

Fig.1. 3 femoral heads.

For each shape, I specified 3 anatomically-chosen landmarks (LMs) and had their coordinates (x,y,z) listed in separate CSV files. For example, one file contains:
51.5378,-65.1603,-33.8293
44.1069,-66.5612,-26.1846
56.3146,-67.5174,-30.7776

Remark: Scalismo recognizes/identifies a landmark by its ID (name). This means that landmarks pointing the same anatomical point/feature on different objects should have the same ID (although different coordinates).

I already set up Scala enviroment. A piece of code can be run by lines in Scala console (Ctrl+Shift+X in IntelliJ); in this way, use Ctrl+Enter for I/O in the console environment. The whole object can be rum by extending the App trait.

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 InitialRR extends App {
  //To send the code selection to Scala Console, select the code in the Editor, press Ctrl + Shift + X
  //imports
  import scalismo.geometry._
  import scalismo.common._
  import scalismo.ui.api._
  import scalismo.mesh._
  import scalismo.registration.LandmarkRegistration
  import scalismo.io.{MeshIO, LandmarkIO}
  import scalismo.numerics.UniformMeshSampler3D
  import breeze.linalg.{DenseMatrix, DenseVector, csvread}
  import scalismo.registration.{Transformation, RotationTransform, TranslationTransform, RigidTransformation}
  import scala.io.StdIn.{readLine, readInt}

  scalismo.initialize()
  implicit val rng = scalismo.utils.Random(42)
  // ui
  val ui = ScalismoUI()
  //------------------------------------------------------------------------------------------------------------------//
  // This is a function to remove a group view by asking the order from the user
  def removeGroupView(groupView:scalismo.ui.api.Group) : Unit = {
    val s = readLine(s"remove ${groupView.name}  [y/n]? ")
    if (s == "y") groupView.remove()
  }

  // [Loading stl and LM files]
  val stlDir =  "I:/registration/SampleData/Meshes/"
  val LMDir =   "I:/registration/SampleData/LMs/"

  val meshFiles = new java.io.File(stlDir).listFiles()
  meshFiles.foreach(file => println(s"${file}\n"))

  val numOfMeshes = meshFiles.size
  println(s"number of files is ${numOfMeshes}")

  val originalData = ui.createGroup("Original Data")
  var i: Int = 1
  val (meshes, meshViews) = meshFiles.map(meshFile => {
    val mesh = MeshIO.readMesh(meshFile).get
    val meshView = ui.show(originalData, mesh, "mesh" + i)
    i = i + 1
    (mesh, meshView) // return a tuple of the mesh and the associated view
  }).unzip // take the tuples apart, to get a sequence of meshes and meshViews
  //meshViews(0).remove()
  removeGroupView(originalData)
   //------------------------------------------------------------------------------------------------------------------//
  // [removing centroid of each mesh i.e moving the centroid to (0,0,0), and hence the mesh by the same translation]

  val MeshPoints = meshes.map(mesh => (mesh.pointSet.points).toArray) //iterator to array

  def centroid(TriMeshPoints: Array[scalismo.geometry.Point[scalismo.geometry._3D]]): EuclideanVector[_3D] = {
    val vecrdPoints = TriMeshPoints.map { point => point.toVector }
    (vecrdPoints.reduce((v1, v2) => v1 + v2)) / TriMeshPoints.size //centroid
  }

  //finding centroids
  val MeshCentroids = MeshPoints.map(mesh => centroid(mesh))
  // translating meshes to the origin i.e making the centroid equal to the origin
  val translations = MeshCentroids.map(centVec => TranslationTransform[_3D](centVec.*(-1)))
  val translatedMeshes = (0 until numOfMeshes).map(i => meshes(i).transform(translations(i)))
  val zeroCentData = ui.createGroup("ZeroCentData")
  val translatedMeshViews = (0 until numOfMeshes).map(i => ui.show(zeroCentData, translatedMeshes(i), "mesh" + i))

  //selecting a RefMesh
  val refMesh = translatedMeshes(0)

  //------------------------------------------------------------------------------------------------------------------//
  //[generating landmark Objcs of each mesh]
  //reading landmark coordinates from files associated with data set files
  val LmFiles = new java.io.File(LMDir).listFiles
  LmFiles.foreach(file => println(s"${file}\n"))

  val LmCrdnts: Array[DenseMatrix[Double]] = LmFiles.map(file => csvread(file, separator=',', skipLines=0))

  // - creating point objects from the lm crdnts.A vector of vectors; each vector contains translated pointObjs of each set of Lms in a file

  val LmPointObjs = (0 until numOfMeshes).map(i =>{
    (0 until LmCrdnts(i).rows).map(row =>
      translations(i)(Point( LmCrdnts(i)(row,0), LmCrdnts(i)(row,1),LmCrdnts(i)(row,2))))  // translating the pointObjs by the translations that move the meshes to their centroids
  })

  val LMs = LmPointObjs.map(eachSet => {
    var n = 0
    eachSet.map(pointObj => {n = n+1; Landmark(s"lm-${n}", pointObj)}) //Ids' names are important for registeration
  }) // a vector of vectors of landmarks. Each vector contains LmObjs of each set of Lms in a file
  val LMviews = LMs.map (eachSet => {
    eachSet.map(lm => ui.show(zeroCentData,lm,lm.id))
  })
  removeGroupView(zeroCentData)
  // -----------------------------------------------------------------------------------------------------------------//
  // [Rigid registration based on landmarks. Rigid (Procrustes) alignment of each mesh to the RefMesh]
  // RefLMs = LMs(0)
  val bestTransform  = (1 until numOfMeshes).map (i => LandmarkRegistration.rigid3DLandmarkRegistration(LMs(i) ,LMs(0), center = Point(0, 0, 0)))
  //rigid3DLandmarkRegistration pares the corresponding landmarks based on their Ids (names)

  //------------------------------------------------------------------------------------------------------------------//
  // [Transforming the landmarks and the meshes onto the RefMesh based on the bestTransform]
  val alignedData = ui.createGroup("alignedData")
  //transforming LMs
  val rigidTransformedLms = (0 until  numOfMeshes-1).map (i => {
    LMs(i+1).map (lm => lm.transform(bestTransform(i)))  //(i+1) as the RefLms stay where they are
  })

  val rigidTransformedLmsviews = rigidTransformedLms.map (eachSet => {
    eachSet.map(lm => ui.show(alignedData,lm,lm.id))
  })
  val RefLmsView = LMs(0).map(lm => ui.show(alignedData,lm,lm.id))

  //transforming Meshes
  val alignedMeshes = (0 until numOfMeshes-1).map(i => translatedMeshes(i+1).transform(bestTransform(i)))
  // showing the meshes
  val alignedMeshesViews = (0 until numOfMeshes-1).map(i => ui.show(alignedData, alignedMeshes(i), "Amesh" + i))
  val RefMeshesView = ui.show(alignedData, refMesh, "RefMesh")

  //------------------------------------------------------------------------------------------------------------------//
  //Saving aligned mesh and LMs to new files
  val outputMeshDir = drive + "I:/registeration/AlignedData/Meshes/"
  val outputLmDir = drive + "I:/registeration/AlignedData/LMs/"

  MeshIO.writeMesh(refMesh, new java.io.File(outputMeshDir + "mesh0.stl"))
  (0 until numOfMeshes-1).foreach(i => MeshIO.writeMesh(alignedMeshes(i), new java.io.File(outputMeshDir + s"mesh${i+1}.stl")))


  LandmarkIO.writeLandmarksJson[_3D](LMs(0), new java.io.File(outputLmDir + "LM0.json"))
  (0 until numOfMeshes-1).foreach(i => LandmarkIO.writeLandmarksJson[_3D](rigidTransformedLms(i), new java.io.File(outputLmDir + s"LM${i+1}.json")))
  // asking for closeing the ui
  val finito = readLine("finito [y/n]? ")
  if (finito == "y") ui.close()
}

Here is what the code does, POPA has rigidly registered two femoral heads onto the one selected as the reference.