angularity

Description

Compute the angularity (delta slope).

Usage

angularity(mesh, ratio = FALSE)

Arguments

Value

If ratio = FALSE, a numeric vector of angularity values i.e. delta slope of each polygon compared to their adjacent polygons, for all the polygons of the mesh. If ratio = TRUE, a numeric vector of angularity ratio values. Ratio is computed by dividing polygon slope by 90, replacing vertex elevation by the average ratio of faces adjacent to the vertex, then dividing the slope of polygons from this new mesh by 90. Although it is a non-standard variable, it was implemented because it better discriminates sharp edges than basic angularity. Warning: both options can take up a significant amount of time for large meshes.

References

Ungar and Williamson (2000)

Examples

delta_slope <- angularity(dkmodel$complex, ratio = FALSE)
summary(delta_slope)
#angularity ratio:
delta_slope_ratio <- angularity(dkmodel$complex, ratio = TRUE)
summary(delta_slope_ratio)

#render on a map:
dkmap(dkmodel$complex, delta_slope, col = "slope",
legend.lab = "Angularity (°)")
#angularity ratio:
dkmap(dkmodel$complex, delta_slope_ratio, col = "slope",
legend.lab = "Angularity ratio")

Average-Relative Curvature (ARC)

Description

Compute a scale-free estimate of mean curvature.

Usage

arc(mesh, range = c(0.01, 0.99))

Arguments

Details

Area-relative curvature (ARC) is obtained by dividing the mean curvature of each triangle by the mean curvature of an hemisphere, the surface area of which is the same as the surface area of the total mesh object. Coincidentally, the surface area of a hemisphere is linked to its mean curvature by the following function: 2.4481 * 1 / square root(surface area) As a result, ARC is a scale-free estimate of surface curvature. It can be used to estimate the sharpness of a mesh.

Value

A numeric vector of area-relative curvature values for all the polygons of the mesh.

Examples

curvature <- arc(dkmodel$complex)
summary(curvature)

#There is a default truncature of extreme values below 1% or above 99%.
#Without truncature:
curvature <- arc(dkmodel$complex, range = c(0, 1))
summary(curvature)

#mean positive ARC:
parc <- mean(curvature[curvature >= 0])
#mean negative ARC:
narc <- mean(curvature[curvature < 0])

#render on a map:
dkmap(dkmodel$complex, curvature, col = "arc",
min.range = -20, max.range = 20)
#absolute truncature of the values above 20 or below -20:
dkmap(dkmodel$complex, curvature, col = "arc", min.range = -20, max.range = 20)

2D surface area

Description

Compute the area of a 2d projection on the (xy) plane.

Usage

area2d(mesh, method = "concave")

Arguments

Details

This function can assess 2D surface area of the projection of a mesh on the (xy) plane. The projection is assimilated to a hull, which can be a convex hull (see chull) or a concave hull (see concaveman). Note that if your mesh is a combination of two or more discontinuous surfaces (e.g., isolated cusps), you should not use either approach. As of version 1.42.2, the concave hull fails intermittently on Mac systems, so the function defaults to convex hulls (on other systems, it defaults to concave hulls)

Value

A single 2D surface area value, corresponding to the footprint of the mesh.

See Also

rfi

Examples

#The following objects should have the exact same footprints:
area2d(dkmodel$basin)
area2d(dkmodel$complex)
area2d(dkmodel$cusp)
area2d(dkmodel$flat)

#Graphical rendering of convex hull
x <- dkmodel$cusp
FootprintVerts <- t(x$vb[1:2, ])
Hull <- grDevices::chull(x = FootprintVerts[, 1], y = FootprintVerts[, 2])
plot(x$vb[1, ], x$vb[2, ], xlab = "x", ylab = "y")
points(x$vb[1, Hull], x$vb[2, Hull], col = "orange1")

#Graphical rendering of concave hull
x <- dkmodel$cusp
FootprintVerts <- t(x$vb[1:2, ])
FootprintVerts[, 1] <- FootprintVerts[, 1] - min(FootprintVerts[, 1])
FootprintVerts[, 2] <- FootprintVerts[, 2] - min(FootprintVerts[, 2])
Hull <- concaveman::concaveman(points = FootprintVerts)
plot(x$vb[1, ] - min(x$vb[1, ]), x$vb[2, ] - min(x$vb[2, ]), xlab = "x", ylab = "y")
points(Hull, col = "green1")

dkborder

Description

Selects the border of a 3d triangle mesh.

Usage

dkborder(mesh)

Arguments

Value

A vector of indices corresponding to the triangles with at least one vertex on the border of the mesh.

Examples

border <- dkborder(dkmodel$cusp)

# Map the border in orange:
is_border <- rep(1, Rvcg::nfaces(dkmodel$cusp))
is_border[border] <- 2
dkmap(dkmodel$cusp, is_border, col = c("white", "#E69F00"), col.levels = 2, legend = FALSE,
scalebar = FALSE, smooth = FALSE)

# Compare with vcgBorder from the R package Rvcg, in blue:
vcgborder <- which(Rvcg::vcgBorder(dkmodel$cusp)$borderit == TRUE)
is_border[vcgborder] <- 3
dkmap(dkmodel$cusp, is_border, col = c("white", "#E69F00", "#56B4E9"), col.levels = 3,
legend = FALSE, scalebar = FALSE, smooth = FALSE)
#As you can see, it all depends on what you want to select!

crop a mesh

Description

Crop a 3d triangle mesh.

Usage

dkcrop(mesh, y)

Arguments

Value

A new object of class mesh3d for which all polygons out of y have been removed.

Examples

#Crop above a certain threshold:
mythreshold <- quantile(elev(dkmodel$basin), 0.5)
mypolynetwork <- poly.network(dkmodel$basin, elev(dkmodel$basin),
lwr.limit = mythreshold)
mynewmesh <- dkcrop(dkmodel$basin, mypolynetwork)
dkmap(mynewmesh, elev(mynewmesh))

#Crop the sharpest dental elements:
sharpmesh <- dkcrop(dkmodel$basin, poly.network(dkmodel$basin,
Rvcg::vcgCurve(dkmodel$basin)$meanitmax,
lwr.limit = quantile(Rvcg::vcgCurve(dkmodel$basin)$meanitmax, 0.8),
min.size = 50))
dkmap(sharpmesh, arc(sharpmesh), col = "arc", col.levels = 20,
min.range = -20, max.range = 20)
#Map of the sharpest elements' elevation, slope and orientation;
dkmap(sharpmesh, elev(sharpmesh), col = "elev")
dkmap(sharpmesh, slope(sharpmesh), col = "slope", col.levels = 9,
min.range = 0, max.range = 90)
dkmap(sharpmesh, orient(sharpmesh), col = "orient", col.levels = 8,
min.range = 0, max.range = 360)

3d topographic map

Description

Map topographic variables on a 3d triangle mesh.

Usage

dkmap(
  mesh,
  y,
  alpha = 1,
  alpha.above = TRUE,
  alpha.faces = NULL,
  alpha.thresh = NULL,
  bg = "white",
  col = "slope",
  col.levels = 100,
  col.main = "black",
  col.lab = "black",
  col.sub = "black",
  col.axis = "black",
  max.range = NULL,
  min.range = NULL,
  lit = TRUE,
  cex = 2,
  cex.axis = 2,
  cex.main = 4,
  cex.sub = 3,
  cex.lab = 2,
  family = "sans",
  font.axis = 1,
  font.lab = 2,
  font.main = 3,
  font.sub = 2,
  main = "",
  sub = "",
  legend = TRUE,
  legend.lab = "y",
  legend.type = "stack",
  windowRect = NULL,
  orient = "current",
  bbox = FALSE,
  origin = TRUE,
  scalebar = FALSE,
  smooth = TRUE
)

Arguments

Value

An rgl window displaying the topography of a variable over a 3d mesh.

See Also

rgl

Examples

#Map of orientation:
orient <- orient(dkmodel$complex)
dkmap(dkmodel$complex, orient, col.levels = 8, col = "orient",
legend.lab = "Orientation (degrees)",legend.type = "pie", min.range = 0,
max.range = 360, orient = "occlusal")

#Map of area-relative curvature:
arc <- arc(dkmodel$complex)
dkmap(dkmodel$complex, arc, col = "arc", legend.lab = "ARC",
min.range = -20, max.range = 20, col.levels = 15, orient = "occlusal")

#Map of Dirichlet normal energy:
dne <- dne(dkmodel$complex)
dkmap(dkmodel$complex, dne, col = "dne", legend.lab = "DNE",
legend.type = "log", orient = "occlusal")

#Map of 3d-Area of polygons (for surface checking):
dkmap(dkmodel$complex, Rvcg::vcgArea(dkmodel$complex, perface = TRUE)$pertriangle,
legend.lab = "3d Area (mm\U00B2)", orient = "occlusal")

dkmodel

Description

A list containing theoretical model surfaces corresponding to a flat surface, a single cusp, a shallow basin and a complex surface.

Usage

dkmodel

Format

An object of class list of length 4.

Source

https://github.com/nialsiG/A-comparison-of-relief-estimates-used-in-3d-dental-topography

Examples

Flat_surface <- dkmodel$flat
Single_cusp <- dkmodel$cusp
Shallow_basin <- dkmodel$basin
Complex_surface <- dkmodel$complex

dkorigin

Description

Sets the lowermost point of the mesh to 0 on the Z-axis

Usage

dkorigin(mesh)

Arguments

Value

An object of class mesh3d.

Examples

#Map of elevation before using dkorigin:
dkmap(dkpongo$OES, doolkit::elev(dkpongo$OES), col = "elev", legend.lab = "Elevation (mm)")

#Map of elevation after dkorigin:
leveled <- dkorigin(dkpongo$OES)
dkmap(leveled, doolkit::elev(leveled), col = "elev", legend.lab = "Elevation (mm)")

dkpongo

Description

A dataset containing the OES and the EDJ surfaces of a Pongo pygmaeus upper second molar (SMF-1117)

Usage

dkpongo

Format

An object of class list of length 2.

Source

https://www.morphosource.org/Detail/MediaDetail/Show/media_id/42357

Examples

Pongo_OES <- dkpongo$OES
Pongo_EDJ <- dkpongo$EDJ

cumulative profile, its slope and the area under its curve

Description

A function for drawing the cumulative profile of a variable, computing the area under the curve and the slope of the profile at the arithmetic mean of the variable.

Usage

dkprofile(
  x,
  type = "cartesian",
  xlab = paste("cumulated frequency (%)"),
  ylab = "",
  main = "",
  col = "red",
  alpha = 1,
  size = 1,
  linetype = "solid"
)

Arguments

Value

A list containing (1) the area under the curve of the profile, (2) the profile to be drawn, and (3) the slope of the profile at the mean of the variable.

References

doi:10.3389/fphys.2017.00524Thiery et al. (2017)

Examples

#Elevation (hypsometric) profile (see Thiery et al., 2017):
dkprofile(elev(dkpongo$OES), main = "Elevation profile - Pongo pygmaeus",
ylab = "Elevation (%)", col = "#0072B2", linetype = "solid")

#Enamel-dentine distance (pachymetric) profile:
dkprofile(oedist(dkpongo$OES, dkpongo$EDJ),
main = "Elevation profile - Pongo pygmaeus", ylab = "Distance (%)",
col = "#F0E442", linetype = "dashed")

#Curvature (kurtometric) profile:
dkprofile(Rvcg::vcgCurve(dkpongo$OES)$meanitmax,
main = "Curvature profile - Pongo pygmaeus", ylab = "Curvature (%)",
col = "#D55E00", linetype = "dotted")

preset orientations

Description

A function to orient 3d topographical maps using preset values.

Usage

dksetview(orient = "occlusal")

Arguments

Value

sets the orientation of the 'rgl' window.

See Also

dkmap

Examples

inclinCusp <- inclin(dkmodel$cusp)
dkmap(dkmodel$cusp, inclinCusp, col = "inclin", min.range = 0, max.range = 180)
dksetview()
#possible orientations are "distal", "left", "occlusal", "mesial" and "right"

Dirichlet normal energy

Description

Compute the Dirichlet normal energy.

Usage

dne(mesh, range = 0.999, total = FALSE)

Arguments

Details

The current algorithm gives a different estimate of DNE when compared with the values obtained using the molaR package. Albeit small, this difference likely comes from different methods of border selection. The doolkit package uses the function dkborder, which accurately selects every polygon sharing a vertex (or more) with the border.

Value

If total = FALSE, a numeric vector of dne values for all the polygons of the mesh. If total = TRUE, a single DNE value, calculated as the sum of the products of polygon normal energies * polygon areas.

References

doi:10.1002/ajpa.21489Bunn et al. (2011)

doi:10.1007/s10914-016-9326-0Pampush et al. (2016)

See Also

DNE

dkborder

Examples

dne <- dne(dkmodel$complex)
summary(dkmodel$complex)

#total DNE value corresponds to the sum of products Dne * Area:
round(sum(dne*Rvcg::vcgArea(dkmodel$complex, perface = TRUE)$pertriangle), 3)
#can be directly computed using \code{dne}:
dne(dkmodel$complex, total = TRUE)

#render on a map:
dkmap(dkmodel$complex, dne, legend.type = "log", col = "dne")

elevation

Description

Compute the elevation (z component of triangle barycenter).

Usage

elev(mesh, origin = TRUE)

Arguments

Value

A numeric vector of elevation values for all the polygons of the mesh.

See Also

inclin

rfi

slope

Examples

elev <- elev(dkmodel$cusp)
summary(elev)

#render on a map:
dkmap(dkmodel$cusp, elev)

hypso

Description

Compute the maximum height, length, width and corresponding hypsodonty index (ratio of the maximum height over the maximum length)

Usage

hypso(mesh, origin = TRUE)

Arguments

Value

A list of values for hypsodonty index, height, length and width of the mesh. The hypsodonty index is not expressed relative to 100 but to 1. Note: the tooth surface is expected to be oriented such as the X-axis is the bucco-lingual axis, the Y-axis is the mesio-distal axis, and the occlusal plane is parallel to the (XY) plane and faces upward.

Examples

hypso(dkmodel$cusp)

inclin

Description

Compute inclination i.e. the angle between triangles and the vertical plane in degrees, comprised between 0 and 180.

Usage

inclin(mesh)

Arguments

Value

A numeric vector of inclination values for all the polygons of the mesh.

References

doi:10.1371/journal.pone.0066142Guy et al. (2013)

See Also

slope

Examples

inclin <- inclin(dkmodel$cusp)
summary(inclin)

#render on a map:
dkmap(dkmodel$cusp, inclin, col = "inclin",
min.range = 0, max.range = 180, legend = TRUE)

Distance from outer enamel surface to enamel dentine junction

Description

Compute the distance from enamel vertices to dentine mesh.

Usage

oedist(oes, edj, ray = FALSE)

Arguments

Value

A numeric vector of vertex-to-mesh distance values for all the polygons of the x mesh.

References

doi:10.1371/journal.pone.0066142Guy et al. (2013)

doi:10.1371/journal.pone.0138802Guy et al. (2015)

doi:10.3389/fphys.2017.00524Thiery et al. (2017)

doi:10.1098/rsbl.2019.0671Schwartz et al. (2020)

See Also

meshDist

Examples

edd <- oedist(dkmodel$cusp, dkmodel$flat)
summary(edd)
AETgeom <- mean(edd)
#Geometric relative enamel thickness, obtained by dividing AETgeom by the
#square root of EDJ area
#Note: it is different from classic RET which requires the volume of the
#dentine inside the enamel cap (see Thiery et al., 2017)
AETgeom/sqrt(Rvcg::vcgArea(dkmodel$flat))
#Absolute crown strength:
edj_radius <- max(dist(cbind(dkmodel$flat$vb[1,], dkmodel$flat$vb[2,])))/2
sqrt(mean(edd) * edj_radius)

#render on a map:
oedist <- doolkit::oedist(dkmodel$cusp, dkmodel$flat)
dkmap(dkmodel$cusp, oedist)
#distance map can also be rendered on EDJ surface:
eodist <- oedist(dkmodel$flat, dkmodel$cusp)
dkmap(dkmodel$flat, eodist)

orientation patch count

Description

Count the number of orientation patches using poly.network.

Usage

opc(mesh, bins = 8, min.size = 3, rotation = 0)

Arguments

Value

A data.frame displaying the number of patches and their size (number of triangles) for each orientation bin. Note: if you want the surface area of each patch, seepoly.network

References

doi:10.1038/nature05433Evans et al. (2007)

See Also

orient

opcr

Examples

#8 bins (default):
opc <- opc(dkmodel$complex)
#8 bins starting from mesial, as in Evans et al. 2007:
opc <- opc(dkmodel$complex, rotation = -(360/16))
#4 bins (mesial, buccal, distal and lingual):
opc <- opc(dkmodel$complex, bins = 4, rotation = -(360/8))

orientation patch count rotated

Description

Compute the orientation patch count rotated of a triangle mesh.

Usage

opcr(mesh, bins = 8, min.size = 3)

Arguments

Value

A data.frame displaying the number of patches and their size (number of triangles) for each orientation bin.

References

doi:10.1038/nature10880Wilson et al. (2012)

See Also

opc

orient

OPCr

Examples

#8bins (default):
opcr <- opcr(dkmodel$complex)
#16 bins (computation time increase exponentially):
opcr <- opcr(dkmodel$complex, bins = 16)

orientation of polygons

Description

Returns the occlusal orientation (exposure in GIS)

Usage

orient(mesh)

Arguments

Value

A numeric vector of occlusal orientation values in degrees for all the polygons of the mesh. Let the orientation from above be depicted as a trigonomical circle, then for a tooth positioned as in Guy et al. (2015) an orientation of 0 (mesial) would be located at an angle of pi/2, and an orientation of 90° (buccal) would be located at an angle of 2*pi.

See Also

opc

opcr

Examples

orient <- orient(dkmodel$complex)
summary(orient)

Identify polygon networks

Description

From a selected variable y, identifies patches of adjacent polygons that share a given range of y values. These patches are called ’polygon networks’.

Usage

poly.network(
  mesh,
  y,
  lwr.limit = stats::quantile(y, 0.75),
  upr.limit = stats::quantile(y, 1),
  min.size = 3
)

Arguments

Value

An object of class "polygon.network" composed of the face index and the membership of each triangle answering the set conditions. The function makes patches of contiguous triangles, and each patch is indexed with a unique number corresponding to its membership.

Examples

#Isolate cusps using elevation:
mythreshold <- quantile(elev(dkmodel$cusp), 0.65)
cusps <- poly.network(dkmodel$cusp, elev(dkmodel$cusp), lwr.limit = mythreshold,
min.size = 100)
myvector <- rep(0, Rvcg::nfaces(dkmodel$cusp))
myvector[cusps@faces] <- cusps@membership[]
myvector <- as.factor(myvector)
ncusps <- length(levels(myvector)) - 1
levels(myvector) <- c(0:ncusps + 1)
dkmap(dkmodel$cusp, as.numeric(myvector), col = cbPalette <- c("#000000", "#E69F00",
"#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7"),
col.levels = ncusps + 1, legend.lab = "Elevation (mm)")

#Any other variables could be used to define the clusters
#Mean curvature:
crests <- poly.network(dkmodel$complex, Rvcg::vcgCurve(dkmodel$complex)$meanitmax,
lwr.limit = quantile(Rvcg::vcgCurve(dkmodel$complex)$meanitmax, 0.8), min.size = 10)
doolkit::dkmap(mesh = dkmodel$complex, y = doolkit::arc(dkmodel$complex,
range = c(-20, 20)), col = "arc", col.levels = 256, min.range = -20,
max.range = 20, orient = "occlusal", legend.lab = "ARC",
alpha.thresh = quantile(doolkit::arc(dkmodel$complex), 0.8), alpha = 0.3,
alpha.above = FALSE)

#Orientation and surface of patches:
patch_orient <- data.frame(bin = NULL, patch = NULL, size = NULL, surface = NULL)
for (i in 1:8) {
  Cluster <- poly.network(dkmodel$complex, orient(dkmodel$complex),
  lwr.limit = 45 * (i - 1), upr.limit = 45 * i)
  Patches <- levels(as.factor(Cluster@membership))
  Bins <- rep(paste(45 * (i - 1), "-", 45 * i), length(Patches))
  Areas <- rep(0, length(Patches))
  for (j in 1:length(Patches)) {
    test <- Cluster@faces[Cluster@membership == Patches[j]]
    Areas[j] <- round(sum(Rvcg::vcgArea(dkmodel$complex,
    perface = TRUE)$pertriangle[test]), 3)
  }
  patch_orient <- data.frame(rbind(patch_orient,
  cbind.data.frame(Bins, Patches, Areas)))
}
#Since patches made of 3 or less polygons are discarded,
#sum of patch areas < total surface area:
sum(patch_orient$Areas)
Rvcg::vcgArea(dkmodel$complex)

S4 class for polygon networks

Description

Polygon networks are subgraphs made of polygons (i) sharing topographic features and (ii) in contact with the rest of the subgraph by at least 1 polygon edge. Objects of S4 class polygon.network are typically made using the function poly.network

Slots

membership

a vector of numeric values specifying, for each triangle, the index number of the patch to which the triangle is assigned

faces

a vector of numeric values indicating the mesh triangle indexes

relief index

Description

Compute the relief index of a 3d triangle mesh.

Usage

rfi(mesh, method = "Ungar", hull = "concave")

Arguments

Details

As of version 1.42.2, the concave hull fails intermittently on Mac systems, so the function defaults to convex hulls (on other systems, it defaults to concave hulls).

Value

A single relief index value.

References

doi:10.1016/j.jhevol.2008.08.002Boyer (2008) doi:10.1371/journal.pone.0066142Guy et al. (2013) Ungar and Williamson (2000)

See Also

area2d

RFI

Examples

rfi <- rfi(dkmodel$cusp, method = "Ungar", hull = "convex")
lrfi <- rfi(dkmodel$cusp, method = "Boyer", hull = "convex")
gamma <- rfi(dkmodel$cusp, method = "Guy")

relief rate

Description

Compute the relief rate from a sub-sample of a 3d triangle mesh. For instance, the relief rate could be computed from the portion of a molar above the lowermost point of its central basin, compared to the whole tooth.

Usage

rrate(uncropped, cropped, origin = TRUE)

Arguments

Value

A single relief rate value.

References

doi:10.1002/ajpa.23916Thiery et al. (2019)

Examples

medelev <- median(elev(dkmodel$cusp))
basins <- dkcrop(dkmodel$cusp, which(elev(dkmodel$cusp) < medelev))
cusps <- dkcrop(dkmodel$cusp, which(elev(dkmodel$cusp) > medelev))

rrate(dkmodel$cusp, basins)
rrate(dkmodel$cusp, cusps)

shape.index

Description

Compute various shape indices.

Usage

shape.index(mesh, origin = TRUE)

Arguments

Details

A handful of indices have been developed to characterize the shape of natural landscapes, including drainage basins. While some indices are very scale-sensitive (e.g., Gravelius' compactness coefficient), others are dimensionless. Horton (1932) introduced a form factor computed as the quotient of the basin's surface area over the square of the maximum basin length. Schumm (1956) developed a basin elongation index computed as the quotient of twice the square root of surface area over the product of basin length and the squareroot of pi. Lastly, Chorley et al. (1957) developed a lemniscate ratio which corresponds to the ratio between the surface of a lemniscate of same length over the basin area,and computed as (pi*(Length^2))/(4*Area).

Value

A list of indices:

• Form factor (Horton, 1932)

• Basin elongation (Schum, 1956)

• Lemniscate ratio 'K' (Chorley et al., 1957)

References

doi:10.1029/TR013i001p00350Horton (1932)

doi:10.1130/0016-7606(1956)67[597:EODSAS]2.0.CO;2Schumm (1956)

doi:10.2475/ajs.255.2.138Chorley et al. (1957)

Examples

ShapInd <- shape.index(dkmodel$basin)
ShapInd$FormFactor
ShapInd$Elongation
ShapInd$K

slope

Description

Compute slope i.e. the angle between triangles and the horizontal plane in degrees, comprised between 0 and 90.

Usage

slope(mesh)

Arguments

Value

A numeric vector of slope values for all the polygons of the mesh.

References

Ungar and Williamson (2000)

See Also

inclin

Examples

slope <- slope(dkmodel$cusp)
summary(slope)

#render on a map:
dkmap(dkmodel$cusp, slope, col.levels = 9, col = "slope",
min.range = 0, max.range = 90, legend = TRUE)