Creating an intensity model using tetrahedral meshes

volumemodel.scala

package local

import scalismo.common.DiscreteField
import scalismo.common.interpolation.BSplineImageInterpolator3D
import scalismo.geometry._3D
import scalismo.image.DiscreteImage
import scalismo.io.StatisticalModelIO
import scalismo.mesh.{ScalarVolumeMeshField, TetrahedralMesh}
import scalismo.numerics.PivotedCholesky
import scalismo.statisticalmodel.DiscreteLowRankGaussianProcess
import scalismo.statisticalmodel.dataset.DataCollection

object BuildIntensityModel {

  def main(args: Array[String]): Unit = {
    scalismo.initialize()
    implicit val rng = scalismo.utils.Random(42)

    val referenceMesh: TetrahedralMesh[_3D] = ???

    // the meshes corresponding to the individual cases - in correspondence
    // with the reference
    val tetrahedralMeshes: Seq[TetrahedralMesh[_3D]] = ???

    // sequence of images, volume(i) should be image corresponding to tetrahedralMeshes(i)
    val volumes: Seq[DiscreteImage[_3D, Short]] = ???

    val scalarVolumeMeshFields = for ((tetraMesh, volume) <- tetrahedralMeshes.zip(volumes)) yield {
      val volumeInterpolated = volume.interpolate(BSplineImageInterpolator3D(degree = 3))
      val scalars = tetraMesh.pointSet.points.map(volumeInterpolated).map(_.toFloat)

      // The domain of each mesh is the reference mesh
      ScalarVolumeMeshField(referenceMesh, scalars.toIndexedSeq)
    }
    val dataCollection = DataCollection.fromScalarVolumeMesh3DSequence(scalarVolumeMeshFields)
    val intensityModel =
      DiscreteLowRankGaussianProcess.createUsingPCA(dataCollection, PivotedCholesky.RelativeTolerance(1e-10))

    StatisticalModelIO.writeVolumeMeshIntensityModel3D(intensityModel, new java.io.File("model.h5")).get

    // sampling from the model
    val sample: DiscreteField[_3D, TetrahedralMesh, Float] = intensityModel.sample()
  }
}