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.

Related resources

The following resources from our online course may provide some helpful context for this tutorial:

• 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}

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

val ui = ScalismoUI()

In the following we will always visualize the effect of different Gaussian process models, by applying the deformations to a reference mesh. We therefore start by loading the mesh and visualizing it in a separate group.

val referenceMesh = MeshIO.readMesh(new java.io.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 we are modelling deformation fields, the mean of the Gaussian process will, of course, itself be a deformation field. In terms of shape models, we can think of the mean function as the deformation field that deforms our reference mesh into the mean shape.

If the reference shape that we choose corresponds approximately to an average shape, and we do not have any further knowledge about our shape space, it is entirely reasonable to use a zero mean; I.e. a deformation field which applies to every point a zero deformation.

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

The covariance function:

The covariance function, which is also referred to as a kernel, defines the properties that characterize likely deformations in our model. Formally, it 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 field outputDim determines the dimensionality of the values we model. In our case, we model 3D vectors, and hence outputDimshould be 3. The domain indicates the set of points on which our kernel is defined. Most often, we set this to the entire Euclidean space RealSpace3D. Finally k denotes the covariance function.

The most often used kernel is the Gaussian kernel. Recall that the scalar-valued Gaussian kernel, is defined by the following formula:

$$k_g(x,x') = \exp^{-\frac{\left\lVert x-x'\right\rVert^2}{\sigma^2} }.$$

A corresponding matrix-valued kernel can be obtained by multiplying the value with an identity matrix (which implies, that we treat each space dimension as independent). In Scalismo, this is defined as follows:

case class MatrixValuedGaussianKernel3D(sigma2 : Double) extends MatrixValuedPDKernel[_3D]() {

  override def outputDim: Int = 3
  override def domain: Domain[_3D] = EuclideanSpace3D;

  override def k(x: Point[_3D], y: Point[_3D]): DenseMatrix[Double] = {
    DenseMatrix.eye[Double](outputDim) * Math.exp(- (x - y).norm2 / sigma2)
  }
}

This constructions allows us to define any kernel. For the most commonly used ones, such as the Gaussian kernel, there is, however, an easier way in Scalismo. First, the scalar-valued Gaussian kernel is already implemented in Scalismo:

val scalarValuedGaussianKernel : PDKernel[_3D]= GaussianKernel3D(sigma = 100.0)

Further, the class DiagonalKernelallows us to turn any scalar-valued kernel into a matrix-valued kernel, by specifying for each dimension of the output-space a kernel and assuming them to be independent. To obtain the same kernel as defined above, we can write:

val matrixValuedGaussianKernel = DiagonalKernel3D(scalarValuedGaussianKernel, scalarValuedGaussianKernel, scalarValuedGaussianKernel)

In this case, since we are using the same kernel in every space dimension, we can write this even more succinct as:

DiagonalKernel3D(scalarValuedGaussianKernel, 3)

Building the GP :

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

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

We can now sample deformations from our Gaussian process at any desired set of points. Below we choose the points to be those 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 its effect by interpolating the deformation field, which we then use 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. Hence if we want to sample from many points we 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 in Scalismo, we can 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). Using this low rank Gaussian process, we can now directly sample continuous deformation fields:

val  defField : Field[_3D, EuclideanVector[_3D]]= lowRankGP.sample()

These in turn, can be used to warp a reference mesh, as discussed above:

referenceMesh.transform((p : Point[_3D]) => p + defField(p))

More conveniently, we can visualize the sampled meshes by building again a Point Distribution Model:

val pdm = PointDistributionModel3D(referenceMesh, lowRankGP)

This model can 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, which are interesting for shape modelling.

Kernels from Statistical Shape Models

As discussed previously, a Statistical Shape Model (SSM) 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 in our statistical shape model.

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(new java.io.File("datasets/lowresModel.h5")).get
val gpSSM = pcaModel.gp.interpolate(TriangleMeshInterpolator3D())

We can then access its covariance function, which is a kernel:

val covSSM : MatrixValuedPDKernel[_3D] = gpSSM.cov

In the next step, we model the additional variance using a Gaussian kernel and add it to the sample covariance kernel.

val augmentedCov = covSSM + DiagonalKernel(GaussianKernel[_3D](100.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.

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), 3)
val gk2 = DiagonalKernel3D(GaussianKernel3D(10.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. Therefore when modelling over such shapes, one would want deformation fields that yield symmetric shapes.

Once we obtained a kernel yielding the type of deformations we desire, it is possible to symmetrize the resulting deformation fields by applying the 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))

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

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