Hello Scalismo!

The goal in this tutorial is to familiarize with a few important data structures in our library and the visualization capabilities of Scalismo Lab.

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.

Meshes (surface data)

Let's start by first reading a triangle mesh from a file:

val mesh : TriangleMesh = MeshIO.readMesh(new File("datasets/Paola.stl")).get

Now, we add the mesh to the 3D scene and assign it a name there. In this case we chose the name Paola to identify the mesh in the 3D scene:

show(mesh, "Paola")

Now that the mesh is displayed in the "Scalismo Viewer's 3D view", you can interact with it as follows:

Note: if you are a Mac user, please find out how to emulate these events using your mouse or trackpad

Note also that you can use the RC, X, Y and Z buttons in the 3D view to recenter the camera on the displayed object.

A 3D triangle mesh in our library is simply a list of 3D points and a list of triangle cells. Here we just print out the first point (the printed String should appear on the result pane in the Scalismo viewer window)

println("first point " + mesh.point(PointId(0)))

and the first cell. 

println("first cell " + mesh.cells(0))

The first cell is a triangle between the first, second and third points of the mesh. Notice here that the cell indicates the identifiers of the points (their index in the point sequence) instead of the geometric position of the points.

Let us now display the points forming the mesh. 

show(mesh.points.toSeq, "pointCloud")

This should add a new point cloud element to the scene with the name "pointCloud".

Note: depending on your computer, visualizing the full point cloud may slow down the visualization performance.

Note that to clean up the 3D scene, you can delete the objects either from the user interface (by right-clicking on the object's name), or programmatically as such :

remove("pointCloud")

Points and Vectors

Different operations exist on 3D points in Scalismo:

The difference between two points yields a vector:

val v : Vector[_3D] = Point(4f,5f,6f) - Point(1f,2f,3f)

The sum of a point with a vector yields a new point:

val p : Point[_3D] = mesh.point(PointId(0)) + v

Points can be converted to vectors:

val v2 : Vector[_3D]= p.toVector

and vice versa:

val p2 : Point[_3D] = v.toPoint
Exercise: Given a list of points, return their center (solution does not have to be a one-liner). Hint: this could be a good occasion to try out the map and reduce functions presented in the cheat sheet.
val pointList : List[Point[_3D]] = List(Point(4f,5f,6f), Point(1f,2f,3f), Point(14f,15f,16f), Point(7f,8f,9f), Point(
10f,11f,12f))
val center : Point[_3D] = ???
val vectors = pointList.map{p : Point[_3D] => p.toVector}  // use map to turn points into vectors
// vectors: List[scalismo.geometry.Vector[scalismo.geometry._3D]] = List(Vector3D(4.0,5.0,6.0), Vector3D(1.0,2.0,3.0), Vector3D(14.0,15.0,16.0), Vector3D(7.0,8.0,9.0), Vector3D(10.0,11.0,12.0))

val vectorSum = vectors.reduce{ (v1,v2) => v1+v2} // reduce the collection of vectors by summing it
// vectorSum: scalismo.geometry.Vector[scalismo.geometry._3D] = Vector3D(36.0,41.0,46.0)

val centerV: Vector[_3D] = vectorSum * (1f / pointList.length ) // divide the sum by the number of points  
// centerV: scalismo.geometry.Vector[scalismo.geometry._3D] = Vector3D(7.2000003,8.2,9.2)

val center = centerV.toPoint
// center: scalismo.geometry.Point[scalismo.geometry._3D] = Point3D(7.2000003,8.2,9.2)

Scalar Images

A discrete scalar image (e.g. gray level image) in Scalismo is simply a function from a discrete domain of points to a scalar value.

Let's read and display a 3D image (MRI of a human):

val image = ImageIO.read3DScalarImage[Short](new File("datasets/PaolaMRI.vtk")).get
show(image, "mri")
Note: depending on your view on the scene, it could appear as if the image is not displayed. In this case, make sure to rotate the scene and change the position of the slices as indicated below.

To visualize the different image slices in the viewer, select "Scene" (the upper node in the scene tree graph) and use the X,Y,Z sliders.

You can also change the way of visualizing the 3D scene under the

View -> Perspective menu.

Scalar Image domain

Let's inspect the domain of the image :

val origin : Point[_3D] = image.domain.origin
val spacing : Vector[_3D] = image.domain.spacing
val size : IntVector[_3D] = image.domain.size

The discrete image domain is a 3-dimensional regular grid of points originating at point (92.5485, -121.926, 135.267), with regular spacing of 1.5 mm in each dimension and containing 171, 171, 139 grid slots in the x, y and z directions respectively.

To better see this, let's display the first 172 points of the image domain 

val imagePoints : Iterator[Point[_3D]] = image.domain.points.take(172)
show(imagePoints.toSeq, "imagePoints")

Scalar image values

The other important part of a discrete image are the values associated with the domain points

val values : Iterator[Short] = image.values

This is an iterator of scalar values of type Short as encoded in the read image.

Let's check the first value, that is the value associated with the origin : 

image.values.next

Given that the origin point has index (0,0,0), this same value can be obtained by accessing the image at this index :

image(IntVector(0,0,0))

Naturally, the number of scalar values should be equal to the number of point domains 

image.values.size == image.domain.points.size

Notice that you can check the intensity value at a particular point position in the image, by maintaining the Ctrl key pressed and hovering over the image. The intensity value will then be displayed in the lower left corner of the Scalismo viewer window.

Creating scalar images

Given that discrete scalar images are simply a mapping between points and values, we can very easily create such images programmatically.

Here we create a new image defined on the same domain of points with artificially created values: we have a value of 1000 at every 50 th point, otherwise 0 : 

val values = (0 until image.domain.numberOfPoints).map{ i:Int => if (i % 50 == 0) 1000 else 0}
val image2 = DiscreteScalarImage(image.domain, ScalarArray(values.toArray))
show(image2,"pattern")
Exercise: create and display a thresholded version of the MRI image, where all intensity values above 300 are replaced with 0 (make sure to use 0.toShort to have the correct type)
val threshValues : Array[Short] = image.values.map{v : Short => ???}.toArray

val thresholded : DiscreteScalarImage[_3D,Short] = ???  
val threshValues : Array[Short] = image.values.map{v :Short => if(v<= 300) v else 0.toShort}.toArray
// threshValues: Array[Short] = Array(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 7, 7, 5, 5, 6, 6, 7, 7, 8, 9, 8, 8, 8, 9, 8, 8, 10, 9, 9, 7, 8, 12, 10, 7, 8, 9, 9, 8, 8, 9, 12, 9, 9, 8, 9, 11, 8, 9, 8, 10, 10, 10, 10, 8, 13, 10, 10, 10, 10, 11, 10, 13, 12, 11, 12, 7, 10, 10, 11, 10, 9, 8, 11, 9, 12, 11, 13, 12, 11, 11, 13, 1...

val thresholded : DiscreteScalarImage[_3D,Short] = DiscreteScalarImage[_3D, Short](image.domain, ScalarArray(threshValues))
// thresholded: scalismo.image.DiscreteScalarImage[scalismo.geometry._3D,Short] = <function1>

show(thresholded, "thresh")

Note that the same result can be achieved much easier, by directly mapping the thresholding function to the image as such :

val thresholded2 : DiscreteScalarImage[_3D,Short] = image.map{v :Short => if(v<=300) v else 0.toShort}

Statistical Mesh Models

Now let's load and show a statistical shape model of faces.

val faceModel = StatismoIO.readStatismoMeshModel(new File("datasets/bfm.h5")).get
show(faceModel, "faceModel")

(suggestion: previously created objects are not needed in the following and can be removed from the scene. You can do this by right-clicking on the Static Objects element and selecting "Remove all" )

Sampling in the UI

Exercise: Sample random instances of faces by using the graphical tools in the scene pane : click on the "Instance 1" tree node and then the "Random" button
Exercise: click a landmark on a position of the face model, e.g. chin or eye corner.. (use the toggle button "LM" in the toolbar to activate landmark clicking). Now continue sampling from the model. What happens to the selected point?
Note: in case landmark clicking did not work for you, you can also add landmarks by hovering over the position and pressing the X key, once the LM toggle button activated. Otherwise, you can add a landmark at the tip of the nose programmatically by executing:
addLandmarksTo(Seq(Landmark("A", faceModel.mean.point(PointId(8156)))), "faceModel")

As you can see, a new instance of the face model is displayed each time along with the corresponding landmark point. Notice how the position of the landmark point changes in space while it keeps the same "meaning" on the face (eye corner, tip of nose ..)

Retrieving Landmarks from the scene

Once you clicked a Landmark point on an object in the scene, the landmark information can be retrieved using the getLandmarksOf method : 

addLandmarksTo(IndexedSeq(Landmark("B", faceModel.mean.point(PointId(86)))),"faceModel")
val landmarks : Seq[Landmark[_3D]] = getLandmarksOf("faceModel").get
val landmarkId : String = landmarks(0).id
val landmarkPosition : Point[_3D] = landmarks(0).point

Sampling programmatically

Let's now sample from the model programmatically. Execute the code below several times (remember that you can select it and press Shift+Enter in the code pane)

remove("sampledFace")
val sampledFace : TriangleMesh = faceModel.sample
show(sampledFace, "sampledFace")
addLandmarksTo(IndexedSeq(Landmark("B", landmarkPosition + Vector(2f,2f,2f))),"faceModel")
Exercise: sample a 100 faces from the model and trace the position of point with id 610 in each sample. Show all the positions taken by this point as a point cloud. What does the point distribution look like? How could one possibly represent it?
val pc610 :IndexedSeq[Point[_3D]] = ???
val pc610 :IndexedSeq[Point[_3D]] = (0 until 100).map(i => faceModel.sample.point(PointId(610)))
// pc610: IndexedSeq[scalismo.geometry.Point[scalismo.geometry._3D]] = Vector(Point3D(176.07132,51.22436,306.5037), Point3D(175.84587,46.11668,299.1928), Point3D(176.52913,47.31555,296.9788), Point3D(178.95473,47.983837,301.95932), Point3D(174.49808,47.562614,298.26947), Point3D(175.89851,52.443333,310.99637), Point3D(174.27841,49.18123,304.5494), Point3D(177.27287,45.63107,296.2007), Point3D(173.76054,49.382957,301.4943), Point3D(173.803,53.141193,306.78156), Point3D(175.2414,48.565052,303.42282), Point3D(174.92194,46.97611,298.46335), Point3D(175.87683,48.28586,299.4449), Point3D(178.40048,46.032158,297.63226), Point3D(175.25412,49.745464,302.80408), Point3D(178.01096,47.822166,304.83115), Point3D(175.4864,49.557175,305.47015), Point3D(179.29689,48.250854,304.80005), Point3D(172.19485,46...
show(pc610, "pc610")
Exercise: compute the mean of the point positions pc610 above and compare it to the position of the same point id on the mean mesh of the face model (either by displaying both, or computing the distance between the two). What do you notice?

Finally, let's sample several faces from the model and compare the number of points in each sample

(0 until 10).foreach{ i :Int => println(faceModel.sample.numberOfPoints) }
Why is the number of points always the same?

We will find out the answer in the coming week.