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).
## 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"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.
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 2 important arguments for users:
outliersThe 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.
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.
BoundaryDiscardThe second 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 ‘Vertex’ more reliably eliminates problematic faces. ‘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).
## Total Surface DNE = 183.9593DNE3d()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.
Note that the color scale in this plot is set relative to the intrinsic energy values of this surface. When comparing multiple surfaces, setting the scale manually will ensure it is the same for all. This is done with the setRange parameter.
The reported “Total Surface DNE” excludes boundary faces and the faces with the highest 0.1% energy densities (according to the value supplied for outliers). Both sets of faces can be accessed through indexing the object created by the DNE function.
The faces identified as being on the boundary or having the highest localized DNE values are stored in the list headers under ‘Edge_Values’ and ‘Outliers’ respectively. You can view their localized DNE values, areas, and locations.
## List of 5
## $ Surface_DNE : num 184
## $ Face_Values :'data.frame': 9999 obs. of 2 variables:
## ..$ Dirichlet_Energy_Densities: num [1:9999] 9.62 15.03 7.66 19.87 20.11 ...
## ..$ Face_Areas : num [1:9999] 0.00392 0.00372 0.00427 0.00352 0.00344 ...
## $ Boundary_Values:'data.frame': 233 obs. of 2 variables:
## ..$ Dirichlet_Energy_Densities: num [1:233] 1.181 1.395 1.23 0.422 0.466 ...
## ..$ Face_Areas : num [1:233] 0.00417 0.00329 0.0036 0.00492 0.00365 ...
## $ Outliers :'data.frame': 10 obs. of 2 variables:
## ..$ Dirichlet_Energy_Densities: num [1:10] 79.6 60.8 85.9 64.8 84.4 ...
## ..$ Face_Areas : num [1:10] 0.00386 0.00331 0.00383 0.00237 0.00439 ...
## $ plyFile :List of 4
## ..$ vb : num [1:4, 1:5054] 0.168 1.904 3.218 0.69 0.198 ...
## .. ..- attr(*, "dimnames")=List of 2
## .. .. ..$ : chr [1:4] "" "" "" "nor"
## .. .. ..$ : NULL
## ..$ 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"## Dirichlet_Energy_Densities Face_Areas
## 1384 1.1810403 0.004165885
## 1463 1.3953744 0.003292230
## 1465 1.2295707 0.003600112
## 1466 0.4222771 0.004920895
## 1468 0.4663526 0.003650777
## 1469 0.4285510 0.003733139## Dirichlet_Energy_Densities Face_Areas
## 4885 79.59439 0.003856765
## 4982 60.84890 0.003305719
## 4984 85.89515 0.003831460
## 5064 64.79186 0.002366240
## 5098 84.37851 0.004388282
## 5155 68.28464 0.003952371For users trying to better understand DNE, consider exploring its properties on the ‘Hills’ object, which is a simple sine-cosine undulating plane.
The RFI() function will measure the three dimensional surface area of the tooth crown and the two dimensional area of the tooth’s planimetric footprint, then use these values to calculate the relief index. The molaR RFI() function relies on an alpha-shape convex hull algorithm which sometimes requires adjustment to properly measure the 2D footprint area. Thus RFI() has one argument, alpha which by default is set to the value 0.06.
alphaTo calculate the 2D footprint area, all mesh vertices are compressed into a single X-Y plane. The ahull() algorithm from the package alphahull then traces the boundaries of the 2D footprint through successively linking points together in a step-wise manner. The upper left-most vertex is typically the starting point, and a radius equal to a percentage of the distance between the 2 furthest points then pivots around the starting point in a counter-clockwise fashion until it contacts another vertex. The newly contacted vertex is the next point in the convex hull succession and identified as part of the footprint boundary.
The size of the radius used to seek the next point in the succession is determined by the alpha argument. Too small an alpha could result in stranding the succession on a point where the next is too far to reach. Too small an alpha sometimes results in finding a point which is interior to the boundary because the next boundary point is further away than the length of the radius. Alternatively, too large an alpha can result in ‘short-cutting’ an infolding, causing the footprint to be erroneously large.
Experimentation has shown, that for the particular example tooth provided in molaR, the best alpha value is 0.08.
## RFI = 0.5909791
## 3D Area = 69.54166
## 2D Area = 21.32687Check2D()To ensure the proper alpha value has been chosen for your surface, it is often best to check the performance with the Check2D() function.
Check2D() plots the points identified as being on the boundary, and then plots slice-shaped triangles originating from the footprint center to each successive pair of the identified boundary points. The colorized footprint is the one calculated and used by the RFI() function. When using Check2D(), the user should inspect the boundary points to make sure the colorized footprint accurately and closely traces the boundary points. Additionally, users should rotate the surface looking for glinting rays originating from the central point or vertices (other than the center) located inside the apparent footprint boundary. These are indicative of extra slice-shaped triangles layered onto the 2D footprint by errant recycling of boundary points by the ahull() function. When it appears that infoldings are ‘short-cut’ by the ahull() function, users should reduce the alpha-value with the alpha argument. When there are extra pie-slice shaped triangles (see above), or the identified boundary clearly infiltrates the actual boundary of the surface, the alpha-value should be increased.
###RFI3d()
The three dimensional surface and its two dimensional footprint can be plotted adjacently using the RFI3d() function. Users can adjust the opacity and color of the tooth mesh, as well as the color of the footprint with the arguments SurfaceColor, Opacity, and FootColor.
(plot not shown to save on data volume)
The OPC() function bins each triangular face on a mesh surface into one of 8 groups based on the X-Y orientation of the face normal, then determines the number of resulting contiguous “patches” composed of adjacent faces sharing the same orientation. Once all the patches have been identified, they are summed to get the OPC value.
The OPC function has 3 arguments: rotation is how many degrees in the X-Y plane the surface will be rotated prior to assessing face orientations, default is 0. minimum_faces sets the minimum number of faces a patch must possess to be counted for the total OPC value. Default for minimum_faces is 3, following the original 2.5D implementation of OPC. Alternatively, users can disable the minimum_faces argument by supplying a positive value to the minimum_area argument, which stipulates a minimum proportion of the total surface area of the tooth each patch must meet to be counted in the total OPC.
## Total Number of Patches = 68
## Number of Patches per Directional Bin =
## Bin 1: 10
## Bin 2: 9
## Bin 3: 7
## Bin 4: 7
## Bin 5: 9
## Bin 6: 6
## Bin 7: 11
## Bin 8: 9As noted above, the default for the OPC() function counts any patch consisting of three or more faces. This can be changed using the minimum_faces parameter or overridden by the minimum_area parameter, which sets the minimum proportion of total surface area a patch must contain to be counted.
## Total Number of Patches = 42
## Number of Patches per Directional Bin =
## Bin 1: 5
## Bin 2: 6
## Bin 3: 6
## Bin 4: 4
## Bin 5: 4
## Bin 6: 5
## Bin 7: 7
## Bin 8: 5## Total Number of Patches = 19
## Number of Patches per Directional Bin =
## Bin 1: 3
## Bin 2: 2
## Bin 3: 2
## Bin 4: 3
## Bin 5: 1
## Bin 6: 3
## Bin 7: 3
## Bin 8: 2Due to the somewhat arbitrary boundaries of bins, differences in specimen orientation during analysis can result in minor variations of OPC. The OPCr() function attempts to account for this by iteratively rotating the tooth (default is 8 iterative rotations spanning a total of 45°) and calculating the OPC of each iteration. A mean OPC is reported. Users can alter the number of rotations with the Steps argument, and the size of each rotation (in degrees) with the stepSize argument.
OPCr_Example1 <- OPCr(Tooth)
OPCr_Example2 <- OPCr(Tooth, Steps = 5, stepSize = 9, minimum_faces = 2) #minimum_faces & minimum_area are passed to each iteration of OPC## $OPCR
## [1] 68.75
##
## $Each_Run
## Degrees_Rotated Calculated_OPC
## [1,] 0.000 68
## [2,] 5.625 72
## [3,] 11.250 69
## [4,] 16.875 69
## [5,] 22.500 66
## [6,] 28.125 66
## [7,] 33.750 76
## [8,] 39.375 64## $OPCR
## [1] 79.4
##
## $Each_Run
## Degrees_Rotated Calculated_OPC
## [1,] 0 76
## [2,] 9 81
## [3,] 18 78
## [4,] 27 79
## [5,] 36 83The object returned by OPCr() also contains the OPC values and degrees of rotation for each iteration:
## Degrees_Rotated Calculated_OPC
## [1,] 0.000 68
## [2,] 5.625 72
## [3,] 11.250 69
## [4,] 16.875 69
## [5,] 22.500 66
## [6,] 28.125 66
## [7,] 33.750 76
## [8,] 39.375 64## Degrees_Rotated Calculated_OPC
## [1,] 0 76
## [2,] 9 81
## [3,] 18 78
## [4,] 27 79
## [5,] 36 83OPC3d()With an object created by the OPC() function, (but importantly, NOT the OPCr() function), users can plot OPC onto their surfaces with the function OPC3d().
OPC3d() has several arguments that allow the user to customize their plots. binColors is a string of colors users can customize to achieve the desired appearance of their plots; default is set to a rainbow array. patchOutline is a logical argument that when enabled traces patches with an outline (default FALSE), while outlineColor sets the tracing color (default is black).
OPC plots are highly customizable with color palettes and patch outlines:
colors <- c('firebrick', 'whitesmoke', 'deeppink', 'darkorchid', 'cornflowerblue', 'cyan', 'skyblue', 'turquoise')(plot not shown to save data)
As of ver. 4.0, molaR contains the Slope() function, which measures the average slope of the surface. It works by assessing the slope of each individual face, weights each face by its relative area, then averages the weighted slopes. Like RFI() and OPC(), the Slope() function is highly reliant on the proper orientation of the analyzed surface. In the case of a tooth, it is important that the normal to the occlusal plane is parallel to the positive Z-direction. Deviations from this aligment can cause significant errors in calculating slope.
Faces located on the tooth’s undersurface or downward-facing overhangs will have negative slopes and are discarded from the analysis.
## Average Surface Slope= 75.48422GuessGuess is the only argument in the Slope() function: a logical (default FALSE) through which the function will try to guess whether the tooth is oriented as described above. When activated, it will flip the tooth into what it assumes is the proper orientation. This guessing is inconsistent and not guaranteed, and should only be used when making figures-not performing actual analyses.
Slope3d()The slope of each face can be plotted with the Slope3d() function. Color control is achieved with a colorramp assembled from a sequential input string of colors, which can be adjusted with the colors argument. The default is a string of colors as follows: colors=c('blue', 'cornflowerblue', 'green', 'yellowgreen', 'yellow', 'orangered', 'red'). Negative slope faces that do not contribute to the overall surface Slope value are by default masked in black, which users can adjust with the maskNegatives parameter.