SSM Using Scalismo (1): Rigid alignment

Hereby, I am going to use Scalismo to rigidly register/align/fit objects represented by triangular meshes onto a particular object. The alignment is based on partial ordinary Procrustes analysis method (POPA). POPA involves translation and rotation (not scaling) of an object with respect to another in order to align/fit the transforming object to the stationary one as close as possible. Read about POPA here. This initial alignment makes the non-rigid registration easier and faster later.

There are 3 geometries of femoral heads in the form of STL and I select one as the reference object (reference mesh) and I will make the other 2 aligned/registered on the reference one.

Fig.1. 3 femoral heads.

For each shape, I specified 3 anatomically-chosen landmarks (LMs) and had their coordinates (x,y,z) listed in separate CSV files. For example, one file contains:
51.5378,-65.1603,-33.8293
44.1069,-66.5612,-26.1846
56.3146,-67.5174,-30.7776

Remark: Scalismo recognizes/identifies a landmark by its ID (name). This means that landmarks pointing the same anatomical point/feature on different objects should have the same ID (although different coordinates).

I already set up Scala enviroment. A piece of code can be run by lines in Scala console (Ctrl+Shift+X in IntelliJ); in this way, use Ctrl+Enter for I/O in the console environment. The whole object can be rum by extending the App trait.

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.

object InitialRR extends App {
  //To send the code selection to Scala Console, select the code in the Editor, press Ctrl + Shift + X
  //imports
  import scalismo.geometry._
  import scalismo.common._
  import scalismo.ui.api._
  import scalismo.mesh._
  import scalismo.registration.LandmarkRegistration
  import scalismo.io.{MeshIO, LandmarkIO}
  import scalismo.numerics.UniformMeshSampler3D
  import breeze.linalg.{DenseMatrix, DenseVector, csvread}
  import scalismo.registration.{Transformation, RotationTransform, TranslationTransform, RigidTransformation}
  import scala.io.StdIn.{readLine, readInt}

  scalismo.initialize()
  implicit val rng = scalismo.utils.Random(42)
  // ui
  val ui = ScalismoUI()
  //------------------------------------------------------------------------------------------------------------------//
  // This is a function to remove a group view by asking the order from the user
  def removeGroupView(groupView:scalismo.ui.api.Group) : Unit = {
    val s = readLine(s"remove ${groupView.name}  [y/n]? ")
    if (s == "y") groupView.remove()
  }

  // [Loading stl and LM files]
  val stlDir =  "I:/registration/SampleData/Meshes/"
  val LMDir =   "I:/registration/SampleData/LMs/"

  val meshFiles = new java.io.File(stlDir).listFiles()
  meshFiles.foreach(file => println(s"${file}\n"))

  val numOfMeshes = meshFiles.size
  println(s"number of files is ${numOfMeshes}")

  val originalData = ui.createGroup("Original Data")
  var i: Int = 1
  val (meshes, meshViews) = meshFiles.map(meshFile => {
    val mesh = MeshIO.readMesh(meshFile).get
    val meshView = ui.show(originalData, mesh, "mesh" + i)
    i = i + 1
    (mesh, meshView) // return a tuple of the mesh and the associated view
  }).unzip // take the tuples apart, to get a sequence of meshes and meshViews
  //meshViews(0).remove()
  removeGroupView(originalData)
   //------------------------------------------------------------------------------------------------------------------//
  // [removing centroid of each mesh i.e moving the centroid to (0,0,0), and hence the mesh by the same translation]

  val MeshPoints = meshes.map(mesh => (mesh.pointSet.points).toArray) //iterator to array

  def centroid(TriMeshPoints: Array[scalismo.geometry.Point[scalismo.geometry._3D]]): EuclideanVector[_3D] = {
    val vecrdPoints = TriMeshPoints.map { point => point.toVector }
    (vecrdPoints.reduce((v1, v2) => v1 + v2)) / TriMeshPoints.size //centroid
  }

  //finding centroids
  val MeshCentroids = MeshPoints.map(mesh => centroid(mesh))
  // translating meshes to the origin i.e making the centroid equal to the origin
  val translations = MeshCentroids.map(centVec => TranslationTransform[_3D](centVec.*(-1)))
  val translatedMeshes = (0 until numOfMeshes).map(i => meshes(i).transform(translations(i)))
  val zeroCentData = ui.createGroup("ZeroCentData")
  val translatedMeshViews = (0 until numOfMeshes).map(i => ui.show(zeroCentData, translatedMeshes(i), "mesh" + i))

  //selecting a RefMesh
  val refMesh = translatedMeshes(0)

  //------------------------------------------------------------------------------------------------------------------//
  //[generating landmark Objcs of each mesh]
  //reading landmark coordinates from files associated with data set files
  val LmFiles = new java.io.File(LMDir).listFiles
  LmFiles.foreach(file => println(s"${file}\n"))

  val LmCrdnts: Array[DenseMatrix[Double]] = LmFiles.map(file => csvread(file, separator=',', skipLines=0))

  // - creating point objects from the lm crdnts.A vector of vectors; each vector contains translated pointObjs of each set of Lms in a file

  val LmPointObjs = (0 until numOfMeshes).map(i =>{
    (0 until LmCrdnts(i).rows).map(row =>
      translations(i)(Point( LmCrdnts(i)(row,0), LmCrdnts(i)(row,1),LmCrdnts(i)(row,2))))  // translating the pointObjs by the translations that move the meshes to their centroids
  })

  val LMs = LmPointObjs.map(eachSet => {
    var n = 0
    eachSet.map(pointObj => {n = n+1; Landmark(s"lm-${n}", pointObj)}) //Ids' names are important for registeration
  }) // a vector of vectors of landmarks. Each vector contains LmObjs of each set of Lms in a file
  val LMviews = LMs.map (eachSet => {
    eachSet.map(lm => ui.show(zeroCentData,lm,lm.id))
  })
  removeGroupView(zeroCentData)
  // -----------------------------------------------------------------------------------------------------------------//
  // [Rigid registration based on landmarks. Rigid (Procrustes) alignment of each mesh to the RefMesh]
  // RefLMs = LMs(0)
  val bestTransform  = (1 until numOfMeshes).map (i => LandmarkRegistration.rigid3DLandmarkRegistration(LMs(i) ,LMs(0), center = Point(0, 0, 0)))
  //rigid3DLandmarkRegistration pares the corresponding landmarks based on their Ids (names)

  //------------------------------------------------------------------------------------------------------------------//
  // [Transforming the landmarks and the meshes onto the RefMesh based on the bestTransform]
  val alignedData = ui.createGroup("alignedData")
  //transforming LMs
  val rigidTransformedLms = (0 until  numOfMeshes-1).map (i => {
    LMs(i+1).map (lm => lm.transform(bestTransform(i)))  //(i+1) as the RefLms stay where they are
  })

  val rigidTransformedLmsviews = rigidTransformedLms.map (eachSet => {
    eachSet.map(lm => ui.show(alignedData,lm,lm.id))
  })
  val RefLmsView = LMs(0).map(lm => ui.show(alignedData,lm,lm.id))

  //transforming Meshes
  val alignedMeshes = (0 until numOfMeshes-1).map(i => translatedMeshes(i+1).transform(bestTransform(i)))
  // showing the meshes
  val alignedMeshesViews = (0 until numOfMeshes-1).map(i => ui.show(alignedData, alignedMeshes(i), "Amesh" + i))
  val RefMeshesView = ui.show(alignedData, refMesh, "RefMesh")

  //------------------------------------------------------------------------------------------------------------------//
  //Saving aligned mesh and LMs to new files
  val outputMeshDir = drive + "I:/registeration/AlignedData/Meshes/"
  val outputLmDir = drive + "I:/registeration/AlignedData/LMs/"

  MeshIO.writeMesh(refMesh, new java.io.File(outputMeshDir + "mesh0.stl"))
  (0 until numOfMeshes-1).foreach(i => MeshIO.writeMesh(alignedMeshes(i), new java.io.File(outputMeshDir + s"mesh${i+1}.stl")))


  LandmarkIO.writeLandmarksJson[_3D](LMs(0), new java.io.File(outputLmDir + "LM0.json"))
  (0 until numOfMeshes-1).foreach(i => LandmarkIO.writeLandmarksJson[_3D](rigidTransformedLms(i), new java.io.File(outputLmDir + s"LM${i+1}.json")))
  // asking for closeing the ui
  val finito = readLine("finito [y/n]? ")
  if (finito == "y") ui.close()
}

Here is what the code does, POPA has rigidly registered two femoral heads onto the one selected as the reference.