Building Blocks for Sampling

Welcome to the second day of the coding afternoons. Today it is all about sampling. We start with a short, general introduction of the software used for sampling. Then we address the problem of drawing samples of our generative face model given observed landmarks.

Basic sampling

First we setup the basic stuff. Create a new scala application file (e.g. Day02.scala). The basic skeleton for the application is the same as provided yesterday:

object Day02 extends App {}

We first initialize scalismo and seed a random number generator. The random number generator is implicitly passed to several parts of the framework where randomness is needed. Hence, setting the seed in the beginning makes the results reproducible.

import scalismo.utils.Random

scalismo.initialize()
implicit val rng = Random(1024L)

Sampling Building Blocks

From the lecture you know what sampling is about. In the framework this maps to pieces of code. In the following A is a type parameter to introduce the software building blocks. You do not need to execute or copy the following code as we will revisit them step-by-step later on:

// An initial state, i.e. a given set of parameters to start with
val intialState: A = ???

// A way to change a set of parameters to propose a new set of parameters
val proposal: ProposalGenerator[A] = ???

// A point wise evaluation of the posterior given a set of parameters
val evaluator: DistributionEvaluator[A] = ???

// A sampling algorithm
val algorithm: Metropolis[A] = Metropolis[A](proposal,evaluator)

// A sampling chain which is lazy in the software
val chainIterator: Iterator[A] = algorithm.iterator(intialState)

// A sequence of realized samples
val states: IndexedSeq[A] = chainIterator.take(5000).toIndexedSeq

Next, we present the example you should know from the lecture.

Toy example

We reuse the example of a 2d Gaussian distribution with strong correlation among the dimensions as the target distribution. We will generate samples from an isotropic Gaussian distribution used as proposal distribution. Then we transform these proposed samples into samples from the target distribution by only point wise evaluations of the target distribution.

The parameters for this example are represented as a 2d vector:

import breeze.linalg.DenseVector

val initialState: DenseVector[Double] = DenseVector(0.0,0.0)

Next we create the proposal distribution. We use the generic trait ProposalGenerator[A] and override the function propose which proposes a new state given the current state. We have chosen to update either of the components with a Gaussian perturbation update. The Gaussian perturbation has a symmetric transition ratio. This means the transition from one state to another is the same as the other way around. This is a nescessary condition for the Metropolis algorithm.

import scalismo.sampling.ProposalGenerator
import scalismo.sampling.SymmetricTransitionRatio
val proposal = new ProposalGenerator[DenseVector[Double]] with SymmetricTransitionRatio[DenseVector[Double]] {

  override def propose(current: DenseVector[Double]): DenseVector[Double] = {
    // selection of the component
    val componentId = rng.scalaRandom.nextInt(2)
    
    // random update
    val delta = rng.scalaRandom.nextGaussian()
    
    val next = current.copy
    next(componentId) += delta
    next
  }
  
}

Then we define the evaluator. We use directly the target density we try to infer. The basic class for the evaluator is the DistributionEvaluator[A] where we override the function logValue(sample: A).

import breeze.linalg.{DenseMatrix, DenseVector}
import scalismo.statisticalmodel.MultivariateNormalDistribution
import scalismo.sampling.DistributionEvaluator

val evaluator = new DistributionEvaluator[DenseVector[Double]] {

  val posterior = MultivariateNormalDistribution.apply(
    mean = DenseVector(1.5,2.0),
    cov = DenseMatrix(
      (2.0,1.9),
      (1.9,2.0)
    )
  )

  override def logValue(sample: DenseVector[Double]): Double = {
    posterior.logpdf(sample)
  }
}

The proposal and the evaluator are then combined in the Metropolis algorithm. We get a lazy iterator over the future samples of the markov chain by providing an initial state:

import scalismo.sampling.algorithms.Metropolis

val algorithm = Metropolis[DenseVector[Double]](proposal,evaluator)
val chainIterator = algorithm.iterator(initialState)

We can draw the number of samples we want from the iterator. Transforming the samples to an IndexedSeq forces the computer to calculate them.

val samples: IndexedSeq[DenseVector[Double]] = chainIterator.take(5000).toIndexedSeq

To compare our samples we draw the same number of points directly from the posterior distribution.

val posteriorSamples = for (i <- (0 until 5000)) yield evaluator.posterior.sample()

Then we visualize the two sets of samples:

import breeze.plot._
val f2d = Figure("Scatter plot - distribution")
val p2d = f2d.subplot(0)
p2d.legend = true

// data to plot
p2d += plot(samples.map(_(0)),samples.map(_(1)), style = '.',name = "samples")
p2d += plot(posteriorSamples.map(_(0)),posteriorSamples.map(_(1)), style = '.',name = "posterior")
p2d.xlim = (-3,7)
p2d.ylim = (-3,7)

Do the two sets of sample follow the same distribution? Change now the scale of the update (delta) in the proposal by an order of magnitude, e.g. multiply the update once with 0.1 and once with 10. Can you observe a change in the result?

While the initial configuration leads to a more or less good result, increasing the updates by a factor of ten will result in much less points. Many samples will be rejected. Nevertheless the distribution still seems to be approximated nicely. Decreasing the update step size by a factor of ten leads to high correlation between the samples. Nearly all samples are accepted, but the 5000 samples we draw are not enough to cover the full range of the distribution. We have to increase the number of samples to better exploit the lengthy structure of the distribution.

Sampling using Landmarks

We will now develop a sampling algorithm to infere the posterior distribution of render parameters given observed landmarks. For this we load the model and some parameters provided for today:

import java.io.File
import scalismo.faces.io._

val modelFile = new File("data/model2009-bfm.h5")
val model = MoMoIO.read(modelFile).get

val targetParameters = RenderParameterIO.read(new File("data/day2/exampleParameters.rps")).get

The above parameters are the underlying ground-truth parameters. We will infer a distribution of parameters from a set of landmarks we observe. The observed landmarks are assumed to be noisy observations. So we will not aim to get a single best estimate.

But first let us create a set of initial parameters for our algorithm to start with:

import scalismo.faces.parameters.{MoMoInstance, RenderParameter}

val modelInstance = MoMoInstance.fromCoefficients(model.zeroCoefficients, modelFile.toURI)
val initialParameters = RenderParameter.defaultSquare.
  withMoMo(modelInstance).
  forImageSize(512, 512)

Next, we visualize the images for the loaded and the initial parameters. For this we create a renderer, render an image for each set of parameters and then display the generated images in a GUI:

import scalismo.faces.sampling.face.MoMoRenderer

val renderer = MoMoRenderer(model)
val target = renderer.renderImage(targetParameters)
val image = renderer.renderImage(initialParameters)


import scalismo.faces.gui.ImagePanel
import scalismo.faces.gui.GUIBlock._
val imagePanel = ImagePanel(image)
val targetPanel = ImagePanel(target)
shelf(targetPanel,imagePanel).displayIn("Sampling")

Next we have a look at landmarks.

Generating the landmarks

The model you have loaded has some predefined landmarks.

val allLandmarkNames = model.landmarks.keys

We can print all names of the landmarks by:

allLandmarkNames.foreach(lmName => println(lmName))

There are mainly two sets of landmarks. The first set of landmarks is known under the name Farkas. We will use our internal set, the second set provided with the model. To select this subset of landmarks we can filter the names of the landmarks and keep those with at least one '.' in the name:

val landmarkNames = allLandmarkNames.filter(lmName => lmName.contains(".")).toSeq

Let us generate the location of the landmarks in the image. We use the renderer to get the 2d position in the image. Using the LandmarksDrawer.drawLandmarks function we can visualize the positions of the landmarks on top of an image. Finally we display the so obtained image in the UI:

import scalismo.faces.landmarks.LandmarksDrawer
import scalismo.faces.color.RGBA

val landmarks = landmarkNames.map(tag => renderer.renderLandmark(tag,targetParameters).get)
val targetLM = LandmarksDrawer.drawLandmarks(target,landmarks,RGBA(1.0,0,0,1),5)
targetPanel.updateImage(targetLM)

Proposals Distribution

From the two images we can see that we need to rotate and translate the model in order to match the target face. Hence we need to sample more than a single parameter. We will adapt the three rotation angles and a 2d translation. We use the GaussianRotationProposal class for generating rotations. The generated rotations have a specified rotation axis and an angle sampled from a zero-mean Gaussian distribution. Further we use a GaussianTranslationProposal which change the translation parallel to the camera, along the x-, and y-axis.

import scalismo.faces.sampling.face.proposals._
import scalismo.geometry.{Vector, Vector3D}

val rollProposal = GaussianRotationProposal(axis = Vector3D.unitZ, sdev = 0.1)
val yawProposal = GaussianRotationProposal(axis = Vector3D.unitY, sdev = 0.1)
val pitchProposal = GaussianRotationProposal(axis = Vector3D.unitX, sdev = 0.1)

val translationProposal = GaussianTranslationProposal(sdev = Vector(10.0,10.0))

To combine the proposals we can use the MixtureProposal. Note that we can specify the weights for the mixture explicitly, otherwise the weights are equally set to 1.0 and then normalized over the full mixture. Important is that we can build a mixture proposal only from the same type of proposals, here all proposals are of type SymmetricProposalGenerator[Pose].

import scalismo.sampling.proposals.MixtureProposal
import MixtureProposal.implicits._

val rotationProposal = MixtureProposal(rollProposal+yawProposal+pitchProposal)
val poseProposal = MixtureProposal(0.4*:rotationProposal + 0.6*:translationProposal)

When at the end of this tutorial block we also adapt the shape we need to generate samples of type RenderParameter. Also it is easier to develop code when as many blocks of our sampling application use RenderParameter as there type parameter.

We can uplift a SymmetricProposalGenerator[Pose] to be SymmetricProposalGenerator[RenderParameter] using the function toParameterProposal. The function applies an implicit type conversion contained in the package ParameterProposal.implicits.

import ParameterProposals.implicits._
val poseParameterProposal = poseProposal.toParameterProposal

Landmark Likelihood

We define now the landmark likelihood. We assume that the individual observed landmarks positions are independent given the model. Additionally we assume that the observations are affected by additive isotropic Gaussian noise. We can create an evaluator for this likelihood using the LandmarkPointEvaluator.isotropicGaussian function:

import scalismo.faces.sampling.face.evaluators._

val landmarksEvaluator = LandmarkPointEvaluator.isotropicGaussian(
  landmarks,    // observed target landmarks
  4.0,          // isotropic uncertainty of each landmark
  renderer      // renderer used to generate the landmark positions
)

Sampling Chain

Now we plug the parts together. First we create the chain using the new proposal and evaluator:

val poseAlgorithm = Metropolis[RenderParameter](poseParameterProposal,landmarksEvaluator)
val chainIterator = poseAlgorithm.iterator(initialParameters)

Then we draw the samples from the chain by:

val landmarkBasedSamples = chainIterator.take(5000).toIndexedSeq

We extract from the full set of parameters the roll angle only:

val rollAngles = landmarkBasedSamples.map(_.pose.roll)

We can visualize the histogram of the sampled roll angles together with the underlying ground truth value from the targetParameters:

val f = Figure("Roll angle - Histogram")
val p = f.subplot(0)
p.legend = true

// data to plot
p += hist(rollAngles,100,name = "roll")

// ugly groundtruth hack
p += plot(
  IndexedSeq.fill(2)(targetParameters.pose.roll),
  IndexedSeq(0.0,-p.ylim._2*0.02),
  name = "ground truth"
)

The logger for Visualization during sampling

We can visualize the samples while they are produced. Like this we can get a feeling for the produced samples. For this we use the ChainStateLogger[A] class. As the name suggest, the current state of a chain is processed further. We can overwrite the function logState(sample: A).

We will write a logger that visualizes the face and the location of the landmarks. For this we first render the image with the current sample. Then we draw the target and the current landmark positions. At the end we update one of the images in the UI.

import scalismo.sampling.loggers.ChainStateLogger

val lmLogger = new ChainStateLogger[RenderParameter] {

  override def logState(sample: RenderParameter): Unit = {
    val currentImage = renderer.renderImage(sample)
    val imageWithTargetLMs = LandmarksDrawer.drawLandmarks(currentImage,landmarks,RGBA(0,0,1,0.33),8)
    val generatedLandmarks = landmarkNames.map(tag => renderer.renderLandmark(tag, sample).get)
    val imageWithBothLMs = LandmarksDrawer.drawLandmarks(
      imageWithTargetLMs,
      generatedLandmarks,
      RGBA(1.0, 0.0, 0.0),
      3
    )
    imagePanel.updateImage(imageWithBothLMs)
  }
  
}

To use the logger in a sampling run we attach it to the iterator using the function loggedWith. This function uses again an implicit conversion as for the function toParameterProposal introduced earlier. We need to import from ChainStateLogger.implicits to make the implicit conversion available. Additionally we used the method subSampled(n: Int). This reduces the computational overhead as only every n^th sample triggers the logger.

import scalismo.sampling.loggers.ChainStateLogger.implicits._

val poseSamples = chainIterator.take(5000).
  loggedWith( // attaches a logger
    lmLogger.subSampled(100) // the logger logs no longer all samples
  ).toIndexedSeq
println("Done!")

When you rerun the code, you should see that until the end the pose wobbles around. Additionally not all landmarks match closely the target landmarks. To correct for that we will next adapt the shape too.

Sampling the shape too

For adapting the shape it is enough to change the shape parameters while sampling. For this we can add another term to our mixture distribution we use as proposal. The framework provides the proposal GaussianMoMoShapeProposal:

val shapeProposal = GaussianMoMoShapeProposal(0.1)

You can combine this proposals with the former pose proposals in a further mixture. Make sure that you uplift both proposals to RenderParameter proposals using the function toParameterProposal. Then run the sampling again. What can you see?

While the landmarks match now better, the overall impression of the faces will be more unnatural, as we evaluate only the landmark likelihood. We do not incorporate yet our prior, the distribution of the parameters of the shape model.

Adding the prior

Following bayes rule we will evaluate the posterior up to the normalization by multiplying the likelihood with the prior. We can do this by first using the class GaussianShapePrior to define an evaluator for the prior. Then you can use the ProductEvaluator to combine the prior with the landmark evaluator which we already defined. As guidance you can use the following snippet.

import scalismo.faces.sampling.face.evaluators.PriorEvaluators.GaussianShapePrior

val shapePrior = GaussianShapePrior(0.0,1.0)

import scalismo.sampling.evaluators.ProductEvaluator
import ProductEvaluator.implicits._
val productEvaluator = ProductEvaluator(landmarksEvaluator * shapePrior)

Change the sampling now to use the product as evaluator. The faces should now look more natural while the landmarks are well matched. However the shape is still not constrained very much. Can you confirm this result?

This concludes the basic building blocks. Next we will introduce the challenge you work on the rest of the week. For the challenge you will start to incorporate also the image color. Further we provide you with information of more advanced topics. Using all the tools you should then be able to come-up with an algorithm for 3d reconstruction of a face from a single image.