Shape completion using Gaussian process regression
The goal in this tutorial is to learn how to use GP regression to predict missing parts of a shape.
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.
Enlarging the flexibility of a shape model
Let's start by loading an incomplete face that we need to reconstruct:
val noseless = MeshIO.readMesh(new File("datasets/noseless.stl")).get
show(noseless,"noseless")As you can see, the nose is missing. In the remainder of the tutorial, we will use a simple face model built from 10 face scans to reconstruct the missing nose:
val littleModel = StatismoIO.readStatismoMeshModel(new File("datasets/model.h5")).get
show(littleModel, "littleModel")As this model was built from a very little dataset, chances that it manages to reconstruct the missing nose properly are rather slim.
To increase the shape variability of the model, we can combine it with a symmetric Gaussian kernel:
val zeroMean = VectorField(RealSpace[_3D], (pt:Point[_3D]) => Vector(0,0,0))
val scalarValuedKernel = GaussianKernel[_3D](30) * 10
case class XmirroredKernel(ker : PDKernel[_3D]) extends PDKernel[_3D] {
override def domain = RealSpace[_3D]
override def k(x: Point[_3D], y: Point[_3D]) = ker(Point(x(0) * -1f ,x(1), x(2)), y)
}
def SymmetrizeKernel(ker : PDKernel[_3D]) : MatrixValuedPDKernel[_3D,_3D] = {
val xmirrored = XmirroredKernel(ker)
val k1 = DiagonalKernel(ker)
val k2 = DiagonalKernel(xmirrored * -1f, xmirrored, xmirrored)
k1 + k2
}
val sim = SymmetrizeKernel(scalarValuedKernel)
val gp = GaussianProcess(zeroMean, sim)
val model = StatisticalMeshModel.augmentModel(littleModel, gp, 50)
show(model, "model")
Here, we started by creating a Gaussian process yielding symmetric smooth deformations as we did in a previous tutorial. We then used this GP to enlarge the flexibility of our little face model by calling the augmentModel method of the StatisticalMeshModel object where we indicated the shape model to enlarge, the Gaussian process we just built, and the desired number of eigenfunctions to be kept.
The augment method used above is a utility function that combines two steps we have seen in previous tutorials:
It starts by building a new continuous Gaussian process where the covariance function is now the sum of two kernels. In our case, these are the covariances learned from data and the kernel function of the indicated GP (the symmetric Gaussian GP in our case).
It then performs a low-rank decomposition of the combined GP (by means of a KL-expansion) and retains only the indicated number of eigenfunctions (50 in our case).
As a result, we obtain a more flexible shape model containing variations from both sample data and a generic smooth symmetric shape model.
Exercise: compare the old and new model by sampling random faces from both.
Model fitting when given correspondence
We will now use our model to perform the reconstruction as follows:
We will try to fit the face model to the given partial face using Gaussian process regression. This means that we seek to find a model instance that resembles well enough the given partial shape.
We will then take the nose part of the fit as a reconstruction proposal for the missing nose.
As we saw previously, to perform GP regression we need observations of the deformation vectors at some points. We can obtain such observations by indicating some correspondences manually:
remove("littleModel")
val lms = LandmarkIO.readLandmarksJson[_3D](new File("datasets/noseFittingLandmarks.json")).get
addLandmarksTo(lms, "noseless")Exercise: click landmarks on the model mean that correspond to the landmarks we just added to our noseless mesh. Attention: it is *very* important that the order of the landmarks in the two lists is the same. Therefore, make sure to click the model landmarks according to the alphabetical order in the figure below:
Note: in case you were unable to click the landmarks, you can add them programmatically by executing the code below
// execute this only if you were unable to click the landmarks
val lm = LandmarkIO.readLandmarksJson[_3D](new File("datasets/noseFittingModelLms.json")).get
addLandmarksTo(lm,"model")val modelPts : Seq[Point[_3D]] = getLandmarksOf("model").get.map{lm => lm.point}
val noselessPts = getLandmarksOf("noseless").get.map{lm => lm.point}By indicating these correspondences above, we just indicated how each selected point of the model should be deformed to its corresponding point on the target mesh. In other words, we observed a few deformation vectors at the selected model points.
Exercise: visualize the partial deformation field that we obtain by clicking these correspondences.
val modelPtIds : Seq[PointId] = modelPts.map(p => model.mean.findClosestPoint(p).id )
// modelPtIds: Seq[scalismo.common.PointId] = Vector(PointId(2215), PointId(6216), PointId(10129), PointId(14000), PointId(8206), PointId(8119), PointId(8173))
val referencePoints = modelPtIds.map(id => model.referenceMesh.point(id)) // deformations are defined on the reference points
// referencePoints: Seq[scalismo.geometry.Point[scalismo.geometry._3D]] = Vector(Point3D(-44.588898,33.49212,91.77969), Point3D(-22.792572,34.274506,94.471886), Point3D(20.83371,37.082195,96.12834), Point3D(43.95854,35.464653,92.0858), Point3D(1.4750519,-54.254044,112.0056), Point3D(-1.4502869,42.292107,114.24812), Point3D(2.4739685,-18.920044,116.90587))
val domain = UnstructuredPointsDomain(referencePoints.toIndexedSeq)
// domain: scalismo.common.UnstructuredPointsDomain[scalismo.geometry._3D] = scalismo.common.UnstructuredPointsDomain3D@7bac9856
val deformations = (0 until referencePoints.size).map(i => noselessPts(i) - referencePoints(i) )
// deformations: scala.collection.immutable.IndexedSeq[scalismo.geometry.Vector[scalismo.geometry._3D]] = Vector(Vector3D(0.36021423,2.733448,-3.9664154), Vector3D(0.64364624,2.6426086,-3.0349808), Vector3D(1.0522346,1.9592819,-3.2956848), Vector3D(0.6336441,2.504036,-2.7499847), Vector3D(-1.2590135,2.7793884,-1.8868179), Vector3D(1.1124456,0.9776001,-1.7943115), Vector3D(-2.7564373,3.9540195,-1.9024124))
val defField = DiscreteVectorField(domain, deformations)
// defField: scalismo.common.DiscreteVectorField[scalismo.geometry._3D,scalismo.geometry._3D] = <function1>
show(defField, "partial_Field")We can now perform GP regression and retrieve the rest of the deformations fitting our observations:
val littleNoise = NDimensionalNormalDistribution(Vector(0,0,0), SquareMatrix((0.5f,0,0), (0,0.5f,0), (0,0,0.5f)))
val trainingData = (modelPts zip noselessPts).map{ case (mPt, nPt) =>
(model.mean.findClosestPoint(mPt).id, nPt, littleNoise)
}
val posterior = model.posterior(trainingData.toIndexedSeq)
show(posterior, "posterior")In this case, we performed the regression directly at the StatisticalMeshModel level and retrieved the posterior mesh model fitting our observations.
With this posterior model, we get a normal distribution of faces satisfying our observations by having the selected characteristic points at the indicated positions.
Exercise: sample a few random faces from the posterior model and compare the position of the selected landmarks points to the corresponding points on the target mesh. How well do the sampled faces resemble our target?
What if we had more correspondences?
As you could hopefully see for yourself, even with the little amount of indicated correspondences, we start to obtain face model instances that resemble our target mesh. To obtain more similar face instances, we need more correspondences:
val modelLMs = LandmarkIO.readLandmarksJson[_3D](new File("datasets/modelLandmarks.json")).get
addLandmarksTo( modelLMs, "model")
val noselessLMs = LandmarkIO.readLandmarksJson[_3D](new File("datasets/noselessLandmarks.json")).get
addLandmarksTo( noselessLMs, "noseless")Here, we just loaded 200 pre-indicated corresponding landmarks from file and added them to both the model and the target mesh.
We can now reiterate our posterior model computation, this time however, with more landmarks:
val modelLandmarks = getLandmarksOf("model").get
val noselessLandmarks = getLandmarksOf("noseless").get
val trainingData = (modelLandmarks zip noselessLandmarks).map{ case (mLm, nLm) =>
(model.mean.findClosestPoint(mLm.point).id, nLm.point, littleNoise)
}
val betterPosterior = model.posterior(trainingData.toIndexedSeq)
show(betterPosterior, "betterPosterior")Exercise: sample a few random faces from the new posterior model. How much variability is left in the model? In which region are the sampled shapes varying the most? (Hint: you can check the first principal components of the shape model to visualize the biggest variations)
Answer: hopefully you managed to see that most of the remaining variability is now at the nose region. This is very noticeable when varying the first principal component of the posterior model.
Finally, as we are interested in the nose region only, let us marginalize our posterior to obtain a posterior nose model as we did in a previous tutorial:
val nosePtIDs = model.referenceMesh.pointIds.filter { id =>
(model.referenceMesh.point(id) - model.referenceMesh.point(PointId(8152))).norm <= 42
}
val posteriorNoseModel = betterPosterior.marginal(nosePtIDs.toIndexedSeq)
show(posteriorNoseModel, "posteriorNoseModel")Exercise: hide all the objects in the 3D scene, except for the noseless face and the nose model. Sample a few instances from the model by changing the coefficients of the principal components and also by drawing random samples. How do the reconstructions look?
Exercise: the borders of the computed nose model instances do not always perfectly match the incomplete mesh. How would you amend this problem?
Answer: one possible solution could be to click more landmarks on the edges that need to be stitched together.
As you could see, using more corresponding points leads to a better nose reconstruction. The question however remains: How did we obtain the 200 correspondences we used for the second posterior model?