Finally, I arrive at statistical shape models. I am going to create a point distribution model (PDM) of shapes. Reading some theory on the point distribution models, I need **superimposed shapes** and their **mean shape**. Superimposed shapes and their mean shapes are available from my previous tutorial on GPA with Scalismo. The outcomes of the full GPA is going to be utilized for creating a PDM. The PDM in this tutorial uses shapes (not forms), and it is based on the multivariate normal distribution of the points’ coordinates collected in a shape vector (Eq. 16 in this post).

This post contains code blocks; if your browser does not automatically display the code blocks, please click on the displayed dark narrow rectangular areas to display the containing code.

I start with initializing Scalismo and loading the superimposed configurations, i.e. shapes. Following the previous tutorial, I use the femoral head shapes.

object PDM extends App { import scalismo.geometry._ import scalismo.common._ import scalismo.ui.api._ import scalismo.mesh._ import scalismo.io.{StatisticalModelIO, MeshIO} import scalismo.statisticalmodel._ import scalismo.registration._ import scalismo.statisticalmodel.dataset._ import scala.io.StdIn.{readLine, readInt} scalismo.initialize() implicit val rng = scalismo.utils.Random(42) // ui val ui = ScalismoUI() // Loading data val drive = "I:/" val meshDir = drive + "registration/Full GPA/" val unsortedMeshFiles = new java.io.File(meshDir).listFiles() unsortedMeshFiles.foreach(file => println(s"\n$file")) println(unsortedMeshFiles) // specifying the file name of the mean mesh generated through the full GPA val meanMeshName = "meanMesh.stl" val meshFileNames = unsortedMeshFiles.map(file => file.getName()) val indexOfMeanMesh = meshFileNames.indexOf(meanMeshName) // bringing the meanMesh to the head (first entry) of the meshFiles sequence val meshFiles = unsortedMeshFiles.updated(0, unsortedMeshFiles(indexOfMeanMesh)).updated(0, unsortedMeshFiles(indexOfMeanMesh))

I have to choose a reference mesh and some training meshes. The full GPA mean mesh is chosen as the reference mes and the three (registered and superimposed) meshes as the training meshes.

// specifying the reference mesh and the training meshes; the full GPA mean mesh is set to be the reference mesh val refMesh = MeshIO.readMesh(meshFiles.head).get val trainingMeshes = meshFiles.tail.map(file => MeshIO.readMesh(file).get)//the other 3 meshes

Once a reference mesh and training meshes are chosen, the deformation vector from each point of the reference mesh to its corresponding point of each training mesh should be calculated. Thereby, the deformation fields morphing the reference mesh to the training meshes are obtained.

// [The meshes are already aligned using the full GPA, and in correspondence, so we can create a discrete GP] // creating the deformation vector field val deformationFields = trainingMeshes.map(mesh => { val deformationVectors = refMesh.pointSet.pointIds.map { id => mesh.pointSet.point(id) - refMesh.pointSet.point(id) }.toIndexedSeq // a discrete deformation field needs a domain and the deformation vectors at each point of the domain DiscreteField[_3D, UnstructuredPointsDomain[_3D], EuclideanVector[_3D]](refMesh.pointSet, deformationVectors) })

The calculated deformation fields are samples (realizations) of a discrete *random field*. The domain of the random field is the reference mesh points. The random field of deformations is assumed to have a multivariate normal distribution. Here, the discrete random field is expressed by a random vector. Thereby, a shape, described by a shape vector, can be produced by adding the shape vector of the reference shape to the normally distributed random vector of deformations expanded by the deformation modes (this is explained by Eq. 16 in this post).

Technically, Scalismo works with the concept of Gaussian process (GP) which is a continuous random field whose any finite sub-collection of random variables (here deformations) has a multivariate normal (Gaussian) distribution. Therefore, continuous sample deformation fields should be approximated based on the known discrete realizations/samples of the random field. By interpolating a discrete field, Scalismo creates a continuous field. Different interpolation approaches are possible; Scalismo uses the Nearest Neighbor Interpolation method. The following code block defines a discrete GP (of the deformations) which is in fact a normally distributed random vector with elements indexed by the reference mesh points (a discrete domain). The discrete GP (random vector of deformations) is expressed in terms of its principal axes (modes). This is done using the PCA as described in this post (see Eq. 16).

val interpolator = NearestNeighborInterpolator[_3D, EuclideanVector[_3D]]() val continuousDefFields = deformationFields.map(field => field.interpolate(interpolator)) val gp = DiscreteLowRankGaussianProcess.createUsingPCA(refMesh.pointSet, continuousDefFields)

As I have 3 training meshes and therefore 3 training deformation fields, the rank of the sample covariance matrix of deformations is at most 3. This means that at most 3 principal axes are used to express the random vector of deformations.

I can sample the GP (random field) of deformations. By adding a sample deformation field to the reference mesh, I can generate a shape vector. This shape vector is a sample mesh from meshes whose points’ coordinates are normally distributed. Sampling some meshes from the *mesh distribution*, I’ll have the following realization:

val PDMGroup = ui.createGroup("refMeshGroup") val refMeshView = ui.show(PDMGroup, refMesh, "refMesh") refMeshView.color = java.awt.Color.GREEN val gpDefView = ui.addTransformation(PDMGroup, gp, "RefMeshDefrmdByGp")

The horizontal bar-like controls in the figures are proportional to the coefficients of the deformation modes (). Each coefficient varies within three standard deviations (Eq. 4 in this post).

Scalismo assembles the reference mesh and the discrete GP as a statistical mesh model:

// [creating a statistical shape model (PDM) val model = StatisticalMeshModel(refMesh, gp.interpolate(interpolator))

This is the statistical PDM.