Using molaR
James D. Pampush & Paul E. Morse
January 26, 2023
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.
Introduction
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.
library(molaR)
str(Tooth)## 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.
Tooth$vb[,1:10]## [,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.000000The 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:
Tooth$it[,1:10]## [,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 5The 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.
Tooth$normals[,1:10]## [,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.0000000The 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:
outliers
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.
kappa
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.
BoundaryDiscard
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).
oex
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.35434Running 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).
DNE3d()
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')Note that the color scale in the plot above is set relative to the intrinsic energy values of this particular surface. When comparing multiple surfaces, setting the scale manually will ensure it is the same for all, making comparisons much easier. This is done with the setMax parameter. setMax can be set to any positive value. It will set the top of the scale to the supplied value and color all faces exceeding the setMax value the max range color. Alternatively, a setMax value dramatically higher than any surface value will have the effect of muting the color of the surface. This may be necessary to make some surface comparisons. When a user supplies a setMax dramatically smaller than most of the values on the surface, all values exceeding the setMax will be colored with the most extreme reds and pinks. In the example below, the new max range value is set much higher than the original max value of 0.26318. Note how much ‘duller’ and ‘cooler’ the generated surface looks.
DNE3d(DNE1, setMax = 1.3, main='Vervet Tooth')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 ‘Boundary_Values’ and ‘Outliers’ respectively. You can view their localized DNE values, areas, and locations.
str(DNE1)## List of 10
## $ Surface_DNE : num 184
## $ Convex_DNE : num 137
## $ Concave_DNE : num 47.4
## $ Convex_Area : num 52.8
## $ Concave_Area : num 16.8
## $ Kappa : num 0
## $ Face_Values :'data.frame': 9999 obs. of 3 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 ...
## ..$ Kappa_Values : num [1:9999] 2.52 3.64 1.98 4.03 4.39 ...
## $ Boundary_Values:'data.frame': 233 obs. of 3 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 ...
## ..$ kappas : num [1:233] 0.778 0.848 0.973 0.398 0.603 ...
## $ Outliers :'data.frame': 10 obs. of 3 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 ...
## ..$ kappas : num [1:10] 0.678 0.8099 0.2535 -5.8669 0.0125 ...
## $ 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"
## - attr(*, "class")= chr "DNE_Object"head(DNE1$Boundary_Values)## Dirichlet_Energy_Densities Face_Areas kappas
## 1384 1.1810403 0.004165885 0.7775519
## 1463 1.3953744 0.003292230 0.8477520
## 1465 1.2295707 0.003600112 0.9733688
## 1466 0.4222771 0.004920895 0.3978485
## 1468 0.4663526 0.003650777 0.6034634
## 1469 0.4285510 0.003733139 0.6125445head(DNE1$Outliers)## Dirichlet_Energy_Densities Face_Areas kappas
## 4885 79.59439 0.003856765 0.6779536
## 4982 60.84890 0.003305719 0.8098859
## 4984 85.89515 0.003831460 0.2535490
## 5064 64.79186 0.002366240 -5.8668604
## 5098 84.37851 0.004388282 0.0124725
## 5155 68.28464 0.003952371 -4.7117739Users can also index the DED (Dirichlet Energy Density) values, areas, and kappa values on each individual face under the $Face_Values portion of the DNE() output. Users may decide that some faces should be reincorporated into their analyses, and this can be quickly done by adjusting the outliers and BoundaryDiscard arguments. For users trying to better understand DNE, consider exploring its properties on the ‘Hills’ object, which is a simple sine-cosine undulating plane.
DNE is a very complex measurement and researchers should be careful with its application. As of molaR 5.0 there are some additional advanced plotting functions aimed at allowing users to examine their results in finer detail so they can be more confident about their conclusions. Below is a brief summary of some of these new functions. Additionally, remember that the DNE() function produces an object which can be indexed for many more specialized analyses.
DNE3dDiscard()
This function plots the analyzed surface highlighting the faces removed from the analysis either because they were classified as outliers, or boundary faces. Users can customize the color of the surface with the arguments boundCol and outlierCol. Additionally, the user can chose to highlight the concave versus convex aspects of the surface by naming new colors in either the baseCol, concaveCol, or both arguments. Like DNE3d(), users can title their plot, fine adjust the view, and export a 3D ply file. This function can help users quickly identify exactly which faces are being excluded from their DNE summation, but may also offer useful ways of visualizing the concave versus convex areas of the surface.
DNE3dDiscard(DNE1, concaveCol='pink', main='Vervet Tooth')DNEpie()
DNEpie() works on an object made by the DNE() function and plots the proportions of either area or DNE contained within the concave and convex partitions. The partitioning is determined by the user inputted kappa value from the original DNE() analysis. Surface area within each partition is the default plot, but users can plot DNE totals with type=DNE. Users can customize the colors in the plot with the arguments convexCol and concaveCol.
DNEpie(DNE1, convexCol='seashell', concaveCol='purple', main='Vervet Tooth')
DNEpie(DNE1, convexCol='seashell', concaveCol='purple', main='Vervet Tooth', type='DNE')
The domination of the surface by the convex area is pretty typical of Cercopithecoid teeth. Note that the concave areas of the tooth account for slightly more of the total DNE value than is expected based on the percentage of concave surface area. This is a very different arrangement as compared to the Hills object.
Hills1 <- DNE(Hills)## Total Surface DNE = 136.5057
## Convex DNE= 53.50109
## Concave DNE= 83.00464DNEpie(Hills1, convexCol='seashell', concaveCol='purple', main='Vervet Tooth Area')
DNEpie(Hills1, convexCol='seashell', concaveCol='purple', main='Vervet Tooth DNE', type='DNE')Note how much more of the surface is contained within the concave portions of the Hills surface. And as before, the total DNE contained within the concave portion of the surface is slightly higher than would be predicted based purely on the amount of area contained within each partition.
DNEDensities()
DNEDensities() works on an object made by the DNE() function and produces density plots of face DNEs (Dirichlet normal energy per face, i.e. Dirichlet normal energy density) or face areas within the concave and convex partitions. The partitioning is determined by the user supplied kappa value from the original DNE() analysis. DNE densities is the default plot, but users can plot face area densities with type=area. The default location of the legend is topright, but users can adjust this with the legendPos argument, which accepts keywords (i.e., ‘topleft’, ‘bottomright’, etc). Users can customize the colors in the plot with the arguments convexCol and concaveCol. Titles can be added with the main argument.
DNEDensities(DNE1, convexCol='seashell', concaveCol='purple', main='Vervet Tooth DNEs')
DNEDensities(DNE1, convexCol='seashell', concaveCol='purple', main='Vervet Tooth \nFace Areas', type='area')
The density plots can be very useful for examining whether outliers have been appropriately excluded from the surface, and to determine whether or not the faces are of relatively equal size with a normal(ish) distribution. As we can see here, there is only a small leftward skew in face DNE values (very small numbers) and not much of a rightward skew towards very large face values. This indicates that the Total Surface DNE value is well representative of the surface gestalt. Furthermore, the face areas plot shows a normal distribution, suggesting that the surface is well processed with faces of roughly equal size throughout.
Keen users should be able to associate specific features on the mesh surface with particular aspects of the density plots.
DNEbar()
The function DNEbar() is made to quickly compare DNE values across multiple surfaces. DNEbar() plots bar charts of surface DNE. Charts can display any of three options, convex DNE total, concave DNE total, and a stacked barplot of both convex and concave DNE surface totals. Users can toggle between these options using the type argument which accepts ‘Concave’, ‘Convex’, or ‘both’ (default). Users can title their plot with main, relocate the position of the legend with legendPos and legendInset as well as customize the plot colors with convexCol and concaveCol. This function works best when applied to a list of DNE() produced objects, but will also plot a single surface if applied to an output of the DNE() function directly. When using a list, the names of each item in the list will be used for the labels below each bar in the plot. Users can alter the orientation of the names with the las argument, text size with the cex.names argument, and can optionally decide to supply their own names for the bar, which will need to be inserted with the arugment names.arg=c(). All of these customizing arguments are passed to and works as it does in the barplot() base R function.
DNEs <- list()
DNE1 <- DNE(Tooth)
DNE2 <- DNE(Hills)
DNEs$Tooth <- DNE1
DNEs$Hills <- DNE2
DNEbar(DNEs, convexCol='seashell', concaveCol='purple')## Total Surface DNE = 183.9593
## Convex DNE= 136.605
## Concave DNE= 47.35434## Total Surface DNE = 136.5057
## Convex DNE= 53.50109
## Concave DNE= 83.00464
Relief Index (RFI)
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 two arguments both centered around the alpha value. The first is alpha which by default is set to the value 0.06 (more below). The second is the logical findAlpha which uses an iterative search to identify the optimal alpha value. Default is FALSE for findAlpha.
alpha
To 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 (i.e. shorter than the search radius). Too small an alpha may also result 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. For these reasons, surfaces with very short vertical sidewalls tend to have issues calculating RFI because the boundaries of the 2D planimetric footprint are relatively bereft of vertices. Provided there are no significant infoldings, it is often best to increase the size of alpha when working with such surfaces.
Calculating RFI with molaR’s RFI() function is a bit different than the other calculations in that it can be performed on surfaces of any vertex and face count. This is not the case for current iteration of DNE() and OPC() because the calculations are based around individual faces, and very high face counts encounter computational limits.
Experimentation has shown, that for the particular example tooth provided in molaR, the best alpha value is 0.08, but on most surfaces a value of 0.06 seems to work best.
RFI1 <- RFI(Tooth, alpha=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(RFI1)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.
RFI3d(RFI1)Orientation Patch Count (OPC)
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 counted 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.
OPC1 <- OPC(Tooth)## 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: 9OPC patch parameters
As 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.
OPC2 <- OPC(Tooth, minimum_faces = 20)## 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: 5OPC3 <- OPC(Tooth, minimum_area = 0.01)## 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: 2Orientation Patch Count Rotated (OPCr)
Due 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:
OPCr_Example1$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 64OPCr_Example2$Each_Run## 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. If users supply fewer colors than orientation bins, then the additional bins will be colored white. patchOutline is a logical argument that when enabled traces patches with an outline (default is FALSE), while outlineColor sets the tracing color (default is ‘black’).
Other plotting options include the logical maskDiscard which blacks out faces excluded from the analysis for being part of an undersized patch. The logical legend enables users to disable the legend. Alternatively when the logical scaleLegend is engaged, the sector radii of the legend pie plot will be scaled to the relative surface area contained in each orientation bin (colored correspondingly, of course, default is false).
OPC3d(OPC1, scaleLegend=TRUE, main='Vervet Tooth')OPC plots are highly customizable with color palettes and patch outlines:
colors <- c('firebrick', 'whitesmoke', 'deeppink', 'darkorchid', 'cornflowerblue', 'cyan', 'skyblue', 'turquoise')OPC3d(OPC1, binColors=colors, patchOutline=T, outlineColor='yellow')OPCbinareas()
The function OPCbinareas() produces either a barplot or a pie chart of the surface areas for each bin in the OPC analysis. Users can choose a pie chart by including type=pie. Titles can be included using the main argument, and users can modify the colors of the plot by entering a string of color key words in the order of the bins.
OPCbinareas(OPC1, main='Vervet Tooth')
OPCbinareas(OPC1, main='Vervet Tooth', type='pie')
Slope
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.
Slope1 <- Slope(Tooth)## Average Surface Slope= 75.48422Guess
Guess 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.
Slope3d(Slope1)Area Relative Curvature (ARC)
Area Relative Curvature (ARC) can be calculated on a surface using the ARC() function. ARC is a measure of curviness or possibly sharpness in the vein of DNE. Originally developed by Guy et al. and briefly presented in their 2013 PLoS ONE paper “Prospective in (Primate) Dental Analysis through Tooth 3D Topographical Quantification” ARC is a surface average of the mean measures of principal curvature of every surface face. This measure strongly correlates with DNE and Convex DNE, though both measures are integrals of the surface rather than surface-wide averages as ARC is.
arc1 <- ARC(Tooth)## Mean ARC = 1.593289
## Mean Positve ARC= 2.841305
## Mean Negative ARC= -2.699952As of the writing of this vignette, there are and have been multiple ways Area Relative Curvature (ARC) has been reported. The function here automatically returns the surface-wide average, the concave-only average, and the convex only average values of relative curvature are all printed and easily found in the resulting ARC-object.
The results of the ARC analysis can also be mapped onto the PLY-file surface, much like othe other plotting functions in this package. To display the map, you must first create the ARC-object using the ARC() function.
ARC3d(arc1)Note for Mac users
If you use a Mac, you may run into some unique problems when trying to install and use molaR. Two common problems are discussed here.
In order to produce the three-D visual maps of topographic surfaces in a Mac OS 10.14 or later, users will need to download and keep updated XQuartz.
Some users may run into problems trying to install the molaR dependency alphahull or its dependencies. If the resulting error includes the term xcrun this indicates that the Mac OS the user is running doesn’t include Xcode, which offers many features but importantly for alphahull Xcode translates commands between code languages. Something necessary for some of the alphahull calculations. To install Xcode on a Mac, open Terminal and run: xcode-select --install. This will download and install Xcode. If that does not work you may need to reset Xcode. Do that by opening Terminal and running: sudo xcode-select --reset. The easiest way to open Terminal on Mac is to hold the command key and hit the space bar.
Tips and Concluding Thoughts
Do you need to process a lot of surfaces all at once? Do you want a simple way to output them all together? Then you want to explore using molaR_Batch().
Basically, molaR_Batch() will apply your desired dental topography analyses to a folder full of .ply files and output a .csv file of the numbers. Its quick and dirty and really meant as a kind of exploratory tool for a quick numbers. He’s an example of how you could code it:
molaR_Batch(pathName=~PathToYourFolder, filename=DNEBatch.csv, DNE=TRUE, RFI=FALSE, OPCr=FALSE, Slope=FALSE)This chunk of code would produce a .csv file of DNE values from all the *.ply files in the folder at the end of ~PathToYourFolder. I have a much more useful bit of code further down in this section.
Are you running into a surface which is giving you some problems? After importing it as an object, try overwriting the original file after running it through the molaR_Clean() function. For example:
Tooth.cleaned <- molaR_Clean(Tooth)## Removed 0 faces with area = 0
## Removed 0 unreferenced vertices from meshThis is a poor example because the Vervet Tooth file we provide is a high quality mesh and the molaR_Clean() function cannot identify any problematic faces or vertecies. A face which has no area (i.e. two legs of the face are on top of one another), or vertices which are not associated with any faces (‘floating’ vertices) need to be removed or they can foul up many of molaR‘s functions. These types of defects can be removed with molaR_Clean(), which has two arguments (cleanType, and verbose). cleanType is a string with three arguments defining what to clean: ’Vertices’, ‘Faces’, or ‘Both’. verbose is a logical which toggles on and off the printing to the console of the alterations made to the ply file. Often this is a good check to make sure that was the problem with your analysis.
The quality of the surface is everything. The tools packaged here in molaR cannot work around problems related to the surfaces themselves.
If you’re like me, and you want to really explore your data when you analyze it so that you are certain you understand the story it is telling, then consider using some list() objects in combination with the tools in molaR.
For example. Here is a little piece of code you can use if you have a folder full of *.ply files and you want to play around with them in R with molaR.
Lets assume you have a folder on your desktop called ‘teeth’ and its full of *.ply files. You can read them into R with the following bit of code.
files.teeth <- list.files(path='/Desktop/teeth', pattern='.ply')
files <- list()
l.files <- length(files.teeth)
for(i in 1:l.files) {
file <- paste('~/Desktop/teeth/', files.teeth[i], sep='')
temp <- vcgPlyRead(file)
files[[i]] <- molaR_Clean(temp)
names(files)[i] <- files.teeth[i]
}That bit of code just created a list object files which is a set of *.ply files all read in together into a single object. This can be real handy now. For example, you can now make similar objects of this one, but with the completed dental topography analyses from molaR.
Here is a useful set of examples using DNE and its related functions.
DNEs <- list()
for(i in 1:length(names(files))) {
DNEs[[i]] <- DNE(files[[i]])
names(DNEs)[i] <- names(files)[i]
}DNEs is now an object similar to files in which all the surfaces have been analyzed for DNE and is as list of DNE_Objects.
Try running:
DNEbar(DNEs)
DNEpie(DNEs[[1]], main=as.character(names(DNEs)[1]))
DNE3d(DNEs[[1]], main=as.character(names(DNEs)[1]))
DNE3dDiscard(DNEs[[1]], main=as.character(names(DNEs)[1]))
DNEDensities(DNEs[[1]], main=as.character(names(DNEs)[1]))If you ran these commands and properly produced the DNEs list object, then you should see how easy it is to index the DNE_Objects out of that list. This should make keeping track of your data, anlyzinging, plotting it, and understanding it, much easier (well, maybe not that last part…).
There are problems, for sure. We are amateur coders; anatomists and anthropologists by original training and self taught developers, sorry. But still, we believe this is an excellent set of dental topography tools, which may as yet find many new applications. We have seen so many creative applications since molaR was first published. Additionally, please, feel free to find creative new extensions of this work, so that it too may evolve as a tool.