Gaussian processes and Point Distribution Models

The goal in this tutorial is to understand the relation between Point Distribution Models (PDM) or Statistical Shape Models, as they are known in Scalismo, and Gaussian Processes.

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.

GPs and PDMs:

Let's start by loading a shape model (or PDM) of faces : 

val faceModel = StatismoIO.readStatismoMeshModel(new File("datasets/bfm.h5")).get
show(faceModel, "faceModel")

As this models a normal distribution of face meshes, we can for example randomly sample from it:

val sampledFace : TriangleMesh = faceModel.sample
show(sampledFace, "sampledFace")

or inspect its mean:

val meanFace : TriangleMesh = faceModel.mean
show(meanFace, "meanFace")

But what happens when we do that?

Well, as we saw previously, when given a set of meshes in correspondence, such as the face meshes used to build this face model, we can see this dataset as a set of deformation fields. This is done by selecting a reference face, and computing the deformation field from this reference to every other face in the dataset.

We also saw that Gaussian Processes allow us to model normal distributions over functions.

Deformation fields, being functions mapping a given D-dimensional point (3D in our face example) into a D-dimensional deformation vector (also 3D here), they are therefore perfect candidates to be modelled by a GP.

In fact, the face model we just loaded is in reality a Gaussian process model of deformation fields.

In the rest of the document, we will see how, when generating a mesh using our PDM, we are in fact generating a deformation field using the GP associated with the PDM and using it to warp the reference mesh of the PDM.

The GP behind the PDM:

The face model we just loaded contains a Gaussian Process : 

val faceGP : DiscreteLowRankGaussianProcess[_3D,_3D] = faceModel.gp

(we will see the meaning of Discrete and Low Rank a bit further in the course. For now, all that matters is that faceGP is a Gaussian Process).

A GP being a normal distribution, it has a mean : 

val meanField : DiscreteVectorField[_3D,_3D] = faceGP.mean

As you see here, the mean of this GP is a deformation field. Let's visualize it :

show(meanField, "meanField")
Exercise : make everything invisible in the 3D scene, except for "meanField" and "meanFace". Now zoom in (right click and drag) on the vector field. Where are the tips of the vectors ending?

As you hopefully saw for yourself, all the tips of the mean deformation vectors end on points of the mean face.

To find out where they start from, let's display the face model's reference mesh : 

val referenceFace : TriangleMesh = faceModel.referenceMesh
show(referenceFace, "referenceFace")
Exercise : Zoom in on the scene and observe the deformation field. Where are the vectors starting this time?

As you can see, the mean deformation field of the Gaussian Process contained in our face model is a deformation from the reference mesh of the model into the mean face mesh.

Hence when calling faceModel.mean, what is really happening is

  1. the mean deformation field is obtained (by calling faceModel.gp.mean)
  2. the mean deformation field is then used to deform the reference mesh (faceModel.referenceMesh) into the triangle Mesh representing the mean face

The same is happening when randomly sampling from the face model :

  1. a random deformation field is sampled (faceModel.gp.sample)
  2. the deformation field is applied to the reference mesh to obtain a random face mesh
Exercise : Perform the 2 steps above in order to sample a random face (that is sample a random deformation first, then use it to warp the reference mesh).
val sampleDef :DiscreteVectorField[_3D, _3D] = faceModel.gp.sample
// sampleDef: scalismo.common.DiscreteVectorField[scalismo.geometry._3D,scalismo.geometry._3D] = <function1>

val interp = sampleDef.interpolateNearestNeighbor 
// interp: scalismo.common.VectorField[scalismo.geometry._3D,scalismo.geometry._3D] = <function1>

def trans(pt: Point[_3D]) = pt + interp(pt) 
// trans: (pt: scalismo.geometry.Point[scalismo.geometry._3D])scalismo.geometry.Point[scalismo.geometry._3D]

val sampleFace = referenceFace.transform(trans)
// sampleFace: scalismo.mesh.TriangleMesh = TriangleMesh(Vector(Point3D(159.80699,-16.824137,302.87454), Point3D(159.79189,-16.54446,302.8275), Point3D(159.50099,-16.553707,302.73486), Point3D(159.5667,-16.900185,302.80557), Point3D(159.7949,-16.265621,302.7584), Point3D(159.45605,-16.208763,302.61874), Point3D(159.86273,-15.928817,302.71658), Point3D(159.48698,-15.848027,302.5246), Point3D(159.92737,-15.604397,302.66724), Point3D(159.5283,-15.483684,302.4357), Point3D(160.01692,-15.207039,302.58197), Point3D(159.59595,-15.065459,302.31705), Point3D(160.11647,-14.81326,302.49014), Point3D(159.66661,-14.649347,302.2188), Point3D(160.25197,-14.413159,302.40152), Point3D(159.76935,-14.237692,302.1105), Point3D(160.40154,-14.015871,302.32104), Point3D(159.8743,-13.827415,301.99094), Point3D(16...

show(sampleFace, "sampleFace")