# Inquiring into a triangular mesh

All codes are written in Python and based on 3D plotting and mesh analysis library PyVista.

## Cells connected to a node

To speed up finding the cells connected to a node, I initiate a dictionary with items {pointId : {cellIds}}. Using a set for cellIds leads to quick update and allows union properties of sets.

numOfPoints = mesh.n_points
meshPointIds = range(numOfPoints)
meshCells = mesh.faces.reshape(-1, 4)[:, 1:4]
pntIdCellIdDict = {pntId: set() for pntId in meshPointIds}
for cellIndx, cell in enumerate(meshCells):
for pntId in cell:
pntIdCellIdDict[pntId].update({cellIndx})

# 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.

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 scalismo.registration.{Transformation, RotationTransform, TranslationTransform, RigidTransformation}

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.