Version: Next

Shape modelling with Gaussian processes and kernels

In this tutorial we learn how to define our own Gaussian processes using analytically defined kernels. Further, we experiment with different kernels that are useful in shape modelling.

To enhance your understanding of this tutorial, we recommend the following resources from our online course:

Covariance functions (Video)
Constructing kernels for shape modelling (Article)
Enlarging the flexibility of statistical shape models (Article)

To run the code from this tutorial, download the following Scala file:

Tutorial07.scala

Preparation

As in the previous tutorials, we start by importing some commonly used objects and initializing the system.

import scalismo.geometry.*
import scalismo.common.*
import scalismo.common.interpolation.TriangleMeshInterpolator3D
import scalismo.mesh.*
import scalismo.io.{StatisticalModelIO, MeshIO}
import scalismo.statisticalmodel.*
import scalismo.numerics.UniformMeshSampler3D
import scalismo.kernels.*

import scalismo.ui.api.*

import breeze.linalg.{DenseMatrix, DenseVector}

import java.io.File

import scalismo.utils.Random.FixedSeed.randBasis

  val ui = ScalismoUI()

In the following we will visualize the effect different Gaussian process models have, on a given reference mesh. We start by loading this reference mesh and visualizing it in a separate group.

  val referenceMesh = MeshIO.readMesh(File("datasets/lowResPaola.ply")).get

  val modelGroup = ui.createGroup("gp-model")
  val referenceView = ui.show(modelGroup, referenceMesh, "reference")

Modelling deformations using Gaussian processes:

A Gaussian Process is defined by two components: the mean function and the covariance function.

The mean:

As our goal is to model deformation fields, the mean of the Gaussian process will naturally be a deformation field. It is the deformation field that deforms our reference mesh into the mean shape.

In the case where the reference shape we select roughly corresponds to an average shape, and without any additional knowledge about our shape space, a zero mean is a reasonable choice.

  val zeroMean = Field(EuclideanSpace3D, (pt:Point[_3D]) => EuclideanVector3D(0,0,0))

The covariance function:

The covariance function, or kernel function, defines the properties that characterize likely deformations in our model. Recall, that a covariance function is a symmetric, positive semi-definite function, $k: X \times X \to R^{ d \times d}$ , which defines the covariance between the values at any pair of points $x, x'$ of the domain. Since we are modelling deformation fields (I.e. vector-valued functions), the covariance function is matrix-valued.

To define a kernel in Scalismo, we need to implement the following methods of the abstract class MatrixValuedPDKernel, which is defined in Scalismo:

abstract class MatrixValuedPDKernel[D]() {

    def outputDim: Int;
    def domain: Domain[D];
    def k(x: Point[D], y: Point[D]): DenseMatrix[Double];
  }

Note: This class is already defined as part of Scalismo. Don't paste it into your code.

The outputDim field sets the dimensionality of the values we model. As we model 3D vectors, outputDim should be 3. The domain indicates the points where our kernel is defined, often set to the entire Euclidean space RealSpace3D. Lastly, k denotes the covariance function.

A frequently used kernel is the Gaussian kernel. Here's how it's defined in Scalismo:

  case class GaussianKernel[D](sigma: Double, scaleFactor: Double = 1) extends PDKernel[D]:
    val sigma2 = sigma * sigma
    val s2 = scaleFactor * scaleFactor

    override def domain = EuclideanSpace[D]

    override def k(x: Point[D], y: Point[D]): Double = 
      val r = x - y
      scala.math.exp(-r.norm2 / sigma2) * s2    

A Gaussian kernel is defined by two parameters, named sigma and scaleFactor. The parameter sigma determines how strongly neighboring points are correlated (i.e. how smooth likely functions are) and scaleFactor determines the amplitude of the function values.

Note that GaussianKernel is not a MatrixValued kernel, but a PDKernel. This means, it outputs only a scalar value. It can be turned into Matrix-valued kernels, using the class DiagonalKernel:

  val scalarValuedGaussianKernel : PDKernel[_3D]= GaussianKernel3D(sigma = 100.0, scaleFactor = 10)
  val matrixValuedGaussianKernel = DiagonalKernel3D(scalarValuedGaussianKernel, scalarValuedGaussianKernel, scalarValuedGaussianKernel)

Note that the DiagonalKernel assumes, that every space dimension is independent. In above example, where we used the same kernel in every space dimension, there is an even shorter way of wrting the same:

  DiagonalKernel3D(scalarValuedGaussianKernel, 3)

Building the GP :

Now that we defined the mean and covariance functions, we can build a Gaussian process as follows:

  val gp = GaussianProcess3D[EuclideanVector[_3D]](zeroMean, matrixValuedGaussianKernel)

We now have the capability to generate deformations from our Gaussian process at any set of points we choose. In the following instance, we select the points of the reference mesh:

  val sampleGroup = ui.createGroup("samples")
  val sample = gp.sampleAtPoints(referenceMesh)
  ui.show(sampleGroup, sample, "gaussianKernelGP_sample")

The result is an instance from the Gaussian Process evaluated at the points we indicated; in this case on the points of the reference mesh.

We can visualize the result of this process by interpolating the deformation field and then applying it to deform the reference mesh:

  val interpolatedSample = sample.interpolate(TriangleMeshInterpolator3D())
  val deformedMesh = referenceMesh.transform((p : Point[_3D]) => p + interpolatedSample(p))
  ui.show(sampleGroup, deformedMesh, "deformed mesh")

Low-rank approximation

Whenever we create a sample using the sampleAtPoints method of the Gaussian process, internally a matrix of dimensionality $nd \times nd$ ,where $n$ denotes the number of points and $d$ the dimensionality of the output space, is created. For a large number of points, we will quickly run out of memory.

We can get around this problem by computing a low-rank approximation of the Gaussian process. To obtain such a representation we use the method approximateGPCholesky of the LowRankGaussianProcess object:

  val lowRankGP = LowRankGaussianProcess.approximateGPCholesky(
    referenceMesh,
    gp,
    relativeTolerance = 0.01,
    interpolator = TriangleMeshInterpolator3D[EuclideanVector[_3D]]()
  )

This call computes a finite-rank approximation of the Gaussian Process using a Pivoted Cholesky approximation. The procedure automatically chooses the rank (I.e. the number of basis functions of the Gaussian process), such that the given relative error is achieved. (The error is measures in terms of the variance of the Gaussian process, approximated on the points of the reference Mesh).

Once we've established this low-rank Gaussian process, we are now equipped to directly sample continuous deformation fields. We can use this deformation field to warp the reference mesh:

  val  defField : Field[_3D, EuclideanVector[_3D]]= lowRankGP.sample()
  referenceMesh.transform((p : Point[_3D]) => p + defField(p))

A more convenient way to sample meshes is, to build a Point Distribution Model from the lowrank Gaussian process. This allows us to directly sample meshes that following the distribution represented by the Gaussian process.

  val pdm = PointDistributionModel3D(referenceMesh, lowRankGP)
  val sampleFromPdm : TriangleMesh[_3D] = pdm.sample()

This model can also be visualized directly in ScalismoUI.

  val pdmView = ui.show(modelGroup, pdm, "group")

Building more interesting kernels

In the following we show a few more examples of kernels that are of interest in the shape modelling context.

Kernels from Statistical Shape Models

As discussed previously, a Point Distribution Model in Scalismo is a discrete Gaussian process. We have seen how to interpolate it to obtain a continuously defined Gaussian Process. As any Gaussian process is completely defined by its mean and covariance function, it follows that this is also true for the GP that representes the PDM.

This allows us to use this sample covariance kernel in combination with other kernels. For example, we often want to slightly enlarge the flexibility of our statistical models. In the following, we show how this can be achieved.

In a first step, we get the Gaussian process from the model an interpolate it.

  val pcaModel = StatisticalModelIO.readStatisticalTriangleMeshModel3D(File("datasets/lowresModel.h5")).get
  val gpSSM = pcaModel.gp.interpolate(TriangleMeshInterpolator3D())

We can then access its covariance function and add to it additional variance using a Gaussian kernel:

  val covSSM : MatrixValuedPDKernel[_3D] = gpSSM.cov
  val augmentedCov = covSSM + DiagonalKernel(GaussianKernel[_3D](100.0, 1.0), 3)

Finally, we build the Gaussian process with the new kernel.

  val augmentedGP = GaussianProcess(gpSSM.mean, augmentedCov)

From here on, we follow the steps outlined above to obtain the augmented SSM.

  val lowRankAugmentedGP = LowRankGaussianProcess.approximateGPCholesky(
    referenceMesh,
    augmentedGP,
    relativeTolerance = 0.01,
    interpolator = TriangleMeshInterpolator3D[EuclideanVector[_3D]]()
  )
  val augmentedSSM = PointDistributionModel3D(pcaModel.reference, lowRankAugmentedGP)

Changepoint kernel:

Another very useful kernel is the changepoint kernel. A changepoint kernel is a combination of different kernels, where each kernel is active only in a certain region of the space.

Here we show how we can define a kernel, which has different behavior in two different regions. The function $s$ acts as a (soft) switch, which determines which kernel is active.

  case class ChangePointKernel(kernel1 : MatrixValuedPDKernel[_3D], kernel2 : MatrixValuedPDKernel[_3D])
    extends MatrixValuedPDKernel[_3D]():

    override def domain = EuclideanSpace[_3D]
    val outputDim = 3

    def s(p: Point[_3D]) =  1.0 / (1.0 + math.exp(-p(0)))
    def k(x: Point[_3D], y: Point[_3D]) = 
      val sx = s(x)
      val sy = s(y)
      kernel1(x,y) * sx * sy + kernel2(x,y) * (1-sx) * (1-sy)
    
  

Let's visualize its effect with two different Gaussian Kernels

  val gk1 = DiagonalKernel3D(GaussianKernel3D(100.0, 1.0), 3)
  val gk2 = DiagonalKernel3D(GaussianKernel3D(10.0, 1.0), 3)
  val changePointKernel = ChangePointKernel(gk1, gk2)
  val gpCP = GaussianProcess3D(zeroMean, changePointKernel)
  val sampleCP =  gpCP.sampleAtPoints(referenceMesh)
  ui.show(sampleGroup, sampleCP, "ChangePointKernelGP_sample")

As you can see each kernel is now active only on one half of the face.

Symmetrizing a kernel

Quite often, the shapes that we aim to model exhibit a symmetry. This is particularly valid in the case of faces. It turns out that we can symmetrize deformations for any given covariance function using the following formula below:

where x_m is the symmetric point to x around the YZ plane.

The resulting kernel will preserve the same smoothness properties of the deformation fields while adding the symmetry around the YZ plane.

Let's turn it into code:

  case class xMirroredKernel(kernel : PDKernel[_3D]) extends PDKernel[_3D]:
    override def domain = kernel.domain
    override def k(x: Point[_3D], y: Point[_3D]) = kernel(Point(x(0) * -1.0 ,x(1), x(2)), y)


  def symmetrizeKernel(kernel : PDKernel[_3D]) : MatrixValuedPDKernel[_3D] = 
    val xmirrored = xMirroredKernel(kernel)
    val k1 = DiagonalKernel(kernel, 3)
    val k2 = DiagonalKernel(xmirrored * -1f, xmirrored, xmirrored)
    k1 + k2
  

  val symmetrizedGaussian = symmetrizeKernel(GaussianKernel[_3D](100, 1.0))

  val gpSym = GaussianProcess3D(zeroMean, symmetrizedGaussian)
  val sampleGpSym =  gpSym.sampleAtPoints(referenceMesh)

  ui.show(sampleGroup, sampleGpSym, "ChangePointKernelGP_sample")

Related resources​

Preparation​

Modelling deformations using Gaussian processes:​

The mean:​

The covariance function:​

Building the GP :​

Low-rank approximation​

Building more interesting kernels​

Kernels from Statistical Shape Models​

Changepoint kernel:​

Symmetrizing a kernel​

Related resources