Using molaR
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.
## 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.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:
## [,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.
## [,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.
## 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.
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.
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.
## 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"## 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.6125445## 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.
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.
## 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.
## 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.
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.
## 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.
## 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: 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:
## 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. 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).
OPC plots are highly customizable with color palettes and patch outlines:
colors <- c('firebrick', 'whitesmoke', 'deeppink', 'darkorchid', 'cornflowerblue', 'cyan', 'skyblue', 'turquoise')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.
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.
## 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.
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.
## 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.
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:
## 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.