Gaussian process sampling and marginalization

The goal in this tutorial is to experiment with the notions of Gaussian process sampling and marginalization, and also learn how to evaluate the probability of shapes according to a given model.

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.

Sampling from discrete and continuous Gaussian processes:

As we have seen previously, Gaussian processes allow us to model normal distributions over both discrete and continuous functions.

Let us load a GP from a statistical face model:

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

and sample from it:

val f : DiscreteVectorField[_3D,_3D] = model.gp.sample 
show(f, "discreteSample")

As you can see here the sampled function is a discrete function, defined over a discrete set of points. This is due to the fact that our Gaussian Process is stored in a file and was therefore discretized over the points of the reference face.

Getting deformations over continuous domains

As we saw in a previous tutorial, it is possible to interpolate a Discrete Vector Field and turn it into a function defined over a continuous domain:

val contF : VectorField[_3D,_3D] = f.interpolateNearestNeighbor

This allows us then to evaluate the deformation field on any point we choose in the 3-dimensional Real space. In this particular type of interpolation, when evaluating the function at a point p, the nearest point to p from the discrete domain of f is identified and its deformation value returned as the value of contF at p.

What if we wanted to sample such deformations defined over continuous domains from our GP?

Of course, one could sample a discrete deformation and interpolate it as we just did above. The same functionality is accessible a much handier way in Scalismo by interpolating the GP itself:

val contGP = model.gp.interpolateNearestNeighbor

Now, when sampling from the continuous GP, we obtain a vector-valued function that we can evaluate anywhere:

val contSample: VectorField[_3D,_3D] = contGP.sample

From continuous to discrete: marginalization

A very interesting property of continuous GPs, differentiating them from Multivariate Normal Distributions (MVN) is that they allow us to evaluate the sampled function values at any set of chosen points.

Now that we have a continuous GP, we can choose to sample a deformation function either on the entire face:

val fullSample = contGP.sampleAtPoints(model.mean)
show(fullSample, "fullSample")

or on a single point:

remove("fullSample")
val singlePointDomain = UnstructuredPointsDomain(IndexedSeq(model.mean.point(PointId(8156))))
val singlePointSample = contGP.sampleAtPoints(singlePointDomain) 
show(singlePointSample, "singlePointSample")

(should be a vector at the tip of the nose, can be behind the face)

Here we started by creating a discrete domain by using the UnstructuredPointsDomain constructor, that when given a finite list of 3D points, defines a discrete domain in Scalismo consisting only of the indicated points.

We then restricted our attention on this discrete domain and sampled a discrete deformation field defined only on this domain, by using the sampleAtPoints method.

This notion of restricting the points of interest to a particular set of points (i.e. a discrete domain), is not only possible when sampling, but also when modeling the variations of the functions over the points of interest. Hence, it is not only possible to get instance functions, that are samples from the distributions, but the actual distribution over functions evaluated at the specified points.

This process is the previously seen marginalization of Gaussian processes:

val rightEyePt: Point[_3D] = model.mean.point(PointId(4281))
val leftEyePt: Point[_3D] = model.mean.point(PointId(11937))
val dom = UnstructuredPointsDomain(IndexedSeq(rightEyePt,leftEyePt))
val marginal : DiscreteGaussianProcess[_3D, _3D] = contGP.marginal(dom)

In the code above, we marginalized all points of the Real Space, except for two eye points on the mean mesh, and obtained a marginal Gaussian process that is defined only on those two points.

As Gaussian Processes are a generalization of Multivariate Normal distributions (MVN), they inherit an interesting property of those, when it comes to marginalization.

A marginal of an MVN yields another MVN.

Similarly, in the case of Gaussian processes, marginalization over a set of variables, yields a Gaussian Process. This time however, as we are interested in modelling variations over a discrete set of points, the Gaussian process is itself discrete as indicated by the returned type above.

Notice that a DiscreteGaussianProcess is in our case nothing more than a multivariate normal distribution, allowing to model deformations defined on top of the two indicated points.

Let's now sample from this marginal GP:

val sample : DiscreteVectorField[_3D, _3D] = marginal.sample 
show(sample, "marginal_sample") 

(deformations might appear behind the eyes)

As you can see, the sample is now a discrete deformation field containing only 2 deformations over the desired points.

Exercise: build a marginal Gaussian process defined over three points, the two eye points above and the tip of the nose and display a few samples from it.

val dom3 = UnstructuredPointsDomain(IndexedSeq(rightEyePt,leftEyePt, model.mean.point(PointId(8156))))
// dom3: scalismo.common.UnstructuredPointsDomain[scalismo.geometry._3D] = scalismo.common.UnstructuredPointsDomain3D@3a139012

val marginal3 : DiscreteGaussianProcess[_3D, _3D] = contGP.marginal(dom3)
// marginal3: scalismo.statisticalmodel.DiscreteGaussianProcess[scalismo.geometry._3D,scalismo.geometry._3D] = scalismo.statisticalmodel.DiscreteGaussianProcess@593f78e2

show(marginal3.sample, "sample3") 

Marginal of a statistical mesh model

Given that a StatisticalMeshModel is in reality just a wrapper around a GP, it naturally allows for marginalization as well:

val noseTipModel : StatisticalMeshModel = model.marginal(IndexedSeq(PointId(8156)))

Notice in this case, how the passed argument to the marginal function is an indexed sequence of point identifiers instead of point coordinates. This is due to the fact that we are marginalizing a discrete Gaussian process. Since the domain of the GP is already discrete, marginalization in this case is done by selecting a subset of the discrete domain. Hence the use of point identifiers instead of 3D coordinates.

Let us now sample from this nose tip model:

val tipSample : TriangleMesh = noseTipModel.sample
println("nb mesh points " + tipSample.numberOfPoints)

Given that the marginal model is a StatisticalMeshModel, sampling from it returns a TriangleMesh. When inspecting the points of the returned sample, we see that it contains only one point, the nose tip.

Exercise: Complete the code below to sample and show 100 nose tips. What do you notice about the dimensionality of the dataset? Why do you think this is the case?

val noseTips : IndexedSeq[Point[_3D]]= ???

val noseTips : IndexedSeq[Point[_3D]]= (0 until 100).map{ i => noseTipModel.sample.point(PointId(0))}
// noseTips: IndexedSeq[scalismo.geometry.Point[scalismo.geometry._3D]] = Vector(Point3D(113.188965,36.937553,247.20032), Point3D(118.66289,34.828724,248.26424), Point3D(99.580864,33.443222,248.4863), Point3D(121.92694,29.234987,248.24617), Point3D(110.34369,32.049225,248.94214), Point3D(110.2979,32.256817,249.3414), Point3D(112.78375,29.229954,248.66107), Point3D(111.52891,33.37366,247.87367), Point3D(112.450356,32.733086,248.46556), Point3D(107.151436,35.95451,248.69998), Point3D(112.7719,35.470657,249.14215), Point3D(121.6468,31.257717,247.71683), Point3D(114.66894,29.63125,249.17279), Point3D(97.787056,36.619854,249.98549), Point3D(117.02041,35.743923,247.68755), Point3D(115.09602,33.09323,247.97017), Point3D(116.22282,32.364777,248.05087), Point3D(111.22469,33.47232,248.38594), Point3...

// Even though the point cloud of nose tips is in 3D, it has a rather planar distribution. This is due to the fact that the tip of the nose was used in the rigid alignment of the face meshes before computing the statistical shape model. As a result, the shape variations at that point were strongly reduced, resulting in this particular case in 2-dimensional variations only.

show(noseTips, "noseTipSamples")

Nose marginal

Let us now generate a marginal StatisticalMeshModel that is defined on more points. In this case, we are interested in modelling shape variations for points of the nose only.

To do so, we first need to extract the list of point identifiers of the nose points:

def isCloseToNoseMiddle(pt: Point[_3D]) : Boolean =  {
  val middleNose = model.referenceMesh.point(PointId(8152))
  (pt - middleNose).norm > 40
}

val refNose : TriangleMesh = Mesh.clipMesh(model.referenceMesh, isCloseToNoseMiddle)
show(refNose,"refNose")

Here we simply clipped away all the face points that are more than 40 mm away from a point on the middle of the nose. As a result, we obtain a mesh of the nose of the reference face.

To retrieve the identifier of each of the nose points on the reference mesh, we can now use the TriangleMesh's findClosestPoint method:

val nosePtIDs = refNose.points.map(p => model.referenceMesh.findClosestPoint(p).id)

Finally, we use the retrieved point identifiers to marginalize our shape model:

val noseModel = model.marginal(nosePtIDs.toIndexedSeq)
show(noseModel, "noseModel")

Exercise: Use the graphical user interface to sample a few random noses.

Probability of shapes and deformations:

Given that GPs and shape models are probability distributions, a typical functionality one expects from such entities is:

Given a shape or a deformation, can I know its probability?

This is implemented in Scalismo via the pdf (Probability Density Function) method:

val defSample : DiscreteVectorField[_3D, _3D] = noseModel.gp.sample
val probDefSample : Double = noseModel.gp.pdf(defSample) 

Here, after sampling a nose deformation field, we use the pdf method to retrieve the probability density value for this particular instance.

Due to the fact that the returned value is a density, one is generally more interested in the ordinal relation between the returned PDFs rather than their actual value. For this reason, it is more common to use the logpdf method, that simply returns the natural logarithm of the density value (pdf) which gives a numerically more stable value than the pdf. The natural logarithm being a monotonous function, it preserves the ordinal relation of the two values, thus allowing to determine which is more probable according to the model :

val defSample1 = noseModel.gp.sample
val defSample2 = noseModel.gp.sample

val logPDF1 = noseModel.gp.logpdf(defSample1)
val logPDF2 = noseModel.gp.logpdf(defSample2)

println("Sampled nose deformation 1 is more probable than sample 2 is : " + (logPDF1 > logPDF2) )

Exercise: Sample a set of 20 nose deformations and compare their logpdf to that of the marginal mean. How predictable was the result ? :)

val meanLogPDF = noseModel.gp.logpdf(noseModel.gp.mean)
// meanLogPDF: Double = -9.189385332096727

(0 until 20).foreach { i => 
  println(noseModel.gp.logpdf(noseModel.gp.sample) <= meanLogPDF)
}
// true
// true
// true
// true
// true
// true
// true
// true
// true
// true
// true
// true
// true
// true
// true
// true
// true
// true
// true
// true

// The mean of a normal distribution being the sample with the highest probability, we could have priorly guessed that the logpdf of any sample will be lower or equal to of the mean.