About this document

This vignette was composed using rmarkdown within RStudio ver. 1.2.5033. It contains WebGL figures that were produced with rglwidget() from the rgl package. To properly view WebGL figures, your browser must be running Javascript and WebGL. Visit https://get.webgl.org for further information. 3D plots in this vignette can be rotated using the left mouse button and zoomed with the mouse wheel.


The R package molaR provides functions that allow the user to quantitatively measure and graphically represent dental surface topography. The following is a demonstration of the primary functions in molaR, as well as some recommended best practices.

molaR analyzes three-dimensional embedded triangular mesh files (*.ply files). These files can be imported into R with the function vcgPlyRead() from the package Rvcg, which can also clean meshes for users. Two sample mesh data files are provided with the molaR package for function demonstration and for users to experiment with: ‘Tooth’ is a scanned M1 of Chlorocebus sabaeus (a vervet monkey: USNM 112176). ‘Hills’ is an undulating plane produced with the formula: z = 3cos(x/2) + 3sin(y/2). Please note that many of embedded graphics have been left out of this document to save on data volume. CRAN only allows packages to be 5mb in size and the inclusion of these large complex rgl surfaces required too much data, and so many have been removed.

## List of 4
##  $ vb      : num [1:4, 1:5054] 0.168 1.904 3.218 1 0.198 ...
##  $ it      : num [1:3, 1:9999] 1 6 5 9 8 1 8 6 1 4 ...
##  $ normals : num [1:4, 1:5054] -0.946 -0.269 0.181 1 -0.843 ...
##  $ material: list()
##  - attr(*, "class")= chr "mesh3d"

Above demonstrates the data of a typical *.ply file. The $vb component is a data frame containing the locations in x-, y-, z- space (in that order) for each of the vertices of the surface. The rows represent from top to bottom the x-, y-, and z- cooridnates for each vertex. The first set of cooridnates corresponds to the first vertex. The second to the second and so on.

##          [,1]     [,2]     [,3]    [,4]     [,5]     [,6]     [,7]     [,8]
## [1,] 0.168007 0.198457 0.187041 0.21913 0.168292 0.145556 0.137716 0.157722
## [2,] 1.904370 1.816000 2.022110 1.83946 1.876530 1.986120 2.068540 2.011160
## [3,] 3.217730 3.182750 3.349860 3.29128 3.111460 3.131930 3.057850 3.254230
## [4,] 1.000000 1.000000 1.000000 1.00000 1.000000 1.000000 1.000000 1.000000
##          [,9]    [,10]
## [1,] 0.191781 0.173828
## [2,] 1.933730 2.113470
## [3,] 3.323510 3.312430
## [4,] 1.000000 1.000000

The command above calls the first 10 vertices. The first vertex is located at: x=0.168007, y=1.90430, and z=3.217730. The last row doesn’t have any meaning in this particular *.ply file but is used as a category for expressing additional information about each vertex.

The $it component stores the face information. Each face is defined through a list of three vertices defining the three legs of the triangle. The convention is to use right-hand rule to indicate which side of the face is the outer vs inner side (counter-clockwise rotation around the vertices). View the first ten faces with:

##      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
## [1,]    1    9    8    4    4    3    9    6   28     6
## [2,]    6    8    6    9    1    8   11    7    5    15
## [3,]    5    1    1    1    2    9    3   15   15     5

The last important part of the *.ply file is the data frame containing the vertex normals ($normals). Each column entry contains information about the x-, y-, z- orientation of the vertex normal and the legnth of the normal. The top 3 rows are ordered x-, y-, z- like the vertex information, and the bottom row contains the length of the normal. Typically the normals are all of unit length.

##            [,1]       [,2]       [,3]       [,4]        [,5]        [,6]
## [1,] -0.9459373 -0.8428705 -0.8828807 -0.7634751 -0.93698084 -0.99140865
## [2,] -0.2692401 -0.5223141 -0.0696271 -0.5095993 -0.34885481 -0.11502183
## [3,]  0.1808654  0.1294501  0.4644069  0.3967549 -0.01916402  0.06227988
## [4,]  1.0000000  1.0000000  1.0000000  1.0000000  1.00000000  1.00000000
##             [,7]        [,8]       [,9]      [,10]
## [1,] -0.99833763 -0.97690153 -0.8743064 -0.9358898
## [2,] -0.04310859 -0.06837288 -0.2620185  0.0668773
## [3,]  0.03825735  0.20245624  0.4085763  0.3458870
## [4,]  1.00000000  1.00000000  1.0000000  1.0000000

The other two components describe any extra information about the surface and the class of the object (a 3d-mesh, i.e., mesh3d).

Dirichlet normal energy (DNE)

Dirichlet normal energy can be calculated on a surface with the DNE() function. Face energy density values (i.e. the measures of localized curvature) can then be rendered onto a three dimensional surface plot using DNE3d(). DNE() has 4 important arguments for users:


The first argument, outliers, sets the percentage range of outlier faces to be excluded from the DNE summation. Default for outlier exclusion is the top one tenth of one percent (0.001). Some surfaces will require a larger outlier exclusion value to account for irregularities on the surface.

Typically, outlier faces are associated with dimples, cracks, spikes, or other imperfections on the mesh which are not representative of the overall curvature of the surface. These imperfections can arise due to the molding, casting, scanning, or downstream digital processing of teeth, but may also be ‘real’ surface features. In the case of these imperfections arising from the original specimen, users should exercise their best judgment when incorporating these features into their analyses. Artifacts arising from the production of the surfaces should ideally be eliminated prior to importation of the surface into R. The DNE measurement is extremely sensitive to these imperfections and uncritical application to surfaces can lead to considerable mismeasurement.

Outlier exclusion is very important for DNE calculation because DNE is a geometrically-based summation. As a curve tightens on a surface, the localized calculation of DNE increases exponentially—not linearly. In some cases, outlier faces will have localized DNE values larger than the rest of the surface combined. Including outliers in the DNE summation is therefore often unrepresentative of the gestalt surface curvature, however they may represent real features or even the main focus of the investigation. Users will have to make informed decisions while considering the goals and hypotheses of each project.


The kappa value defines the partitioning of the surface into concave vs convex portions. In addition to the intensity of the curvature, the DNE function also measures surface curvature orientation. The kappa value dictates the division between putative convex and concave faces. A flat face will have a kappa value of 0, negative values are concave. Positive values are convex. The default for kappa is 0 but users can adjust it to ensure a large portion of the surface is considered either concave or convex, depending on their analysis goals. This value DOES NOT affect the overall calculation of DNE, it simply just determines the partitioning of the values into the concave and convex portions. Users can of course, ignore this subsetting of the DNE output and continue to use the conventional summed measurement which is returned in a form unchanged from prior version of molaR.


The third argument, BoundaryDiscard, sets the criteria for excluding surface faces on the boundary of the mesh. Excluding faces on the boundary of the mesh is important because often these faces have highly irregular and inconsistent vertex normals—due to the lack of an adjacent face with which to calibrate the orientation of the vertex normal. Because the orientations of the vertex normals are included in the DNE calculation, it is important that they be accurate, however this is often not the case with vertices on the mesh boundary.

There are three options for BoundaryDiscard: ‘Vertex’ (default) will exclude faces with at least 1 vertex on the boundary from the DNE summation. ‘Leg’ will exclude faces which have a leg (i.e., 2 vertices) on the boundary; this setting was formerly the default and many previous publications have used this boundary exclusion criterion. However, recent studies have shown that ‘Leg’ often fails to eliminate problematic faces, so as of molaR 4.5 the default has been to exclude faces with any vertex on the boundary. ‘None’ will not discard any boundary faces, this option should be used when working on a closed surface (i.e., a mesh with no boundary).


Outlier exclusion partitioning, oex is a string with two options. c (default) excludes outliers off the surface by considering the entire surface as a whole and then discarding the user-defined outlier percentage. s splits the surface into the concave and convex portions and then discards the user-defined percentage from each of the two portions. The two options eliminate the same number of faces but the percentage of each orientation removed will be defined by this parameter. In the case of most dental surfaces, the majority of faces identified as outliers tend to be included in the concave portions of the surface. Yet in many taxa the convex portions of the surface account for the majority of the total surface area. Because of these factors, users should carefully evaluate how they want their outliers excluded. It is the recommendation of the authors not to adjust the exclusion criteria and leave it as removing the top 0.001 of faces, regardless of concavity and convexity.

DNE1 <- DNE(Tooth)
## Total Surface DNE = 183.9593 
## Convex DNE= 136.605 
## Concave DNE= 47.35434

Running the DNE() function returns the Convex, Concave, and Total Surface DNE. Users uninteresed in the paritioning of the surface can ignore the Convex and Concave reportings and focus on the Total Surface DNE, which in this tooth mesh is measured at 183.9593. For more details users can explore the components of the DNE1 object (see below).


The energy densities calculated across the surface can be plotted using the DNE3d() function. Due to the skewed distribution of exponentially-increasing energy densities, with relatively few mesh faces typically contributing large values to the total surface DNE, DNE plots by default display log-transformed surface energy, which users can disable with the logColors parameter. Additionally, as of molaR 5.0 DNE3d() also applies two different color patterns—inspired from precipitation maps showing combinations of rain and snow—for the concave and convex portions. Users can disable this feature with the argument signColor=FALSE, restoring the conventional plotting color scheme.

Users can title the plot with the main argument. Font sizes are adjustable with the cex and cex.main arguments. Users can make additional tweaks to their plots with the other arguments leftOffset and legend. If the user wishes to export their plot (without the scale, however) into a surface mesh PLY file, they can by specifying a name after the fileName argument. This will print a DNE-colored *.ply file of the supplied name to the working directory. These files can then be imported into other visualizing software packages with the DNE coloration intact.

DNE3d(DNE1, main='Vervet Tooth')