In cornelmpop/Lithics3D: A toolbox for 3D analysis of archaeological lithics

library(rgl)
library(Lithics3D)
knitr::opts_chunk$set(
  collapse = TRUE,
  rgl.newwindow = TRUE,
  comment = "#>"
)
#knitr::knit_hooks$set(webgl = hook_webgl)

Introduction

The purpose of the Lithics3D package is to provide a flexible R toolbox for working with 3D scans of archaeological lithics. The included functions are meant to be as modular as possible so that they may be combined in different ways to achieve a wide range of goals. This vignette demonstrates some potential uses of the Lithics3D library, which, at the moment, provides the following core components (for details please see the help pages):

Functions to solve basic geometric problems such as line/line intersections or projecting points onto a plane
Functions for orienting 3D scans
Functions for sectioning 3D scans into meaningful components (e.g. ventral surfaces of flakes)
Functions for measuring edge angles and other object attributes (e.g. volume)
Three scans (and landmarks) used for test purposes

About this document

Currently this is but a quick and dirty way of showcasing some of the functionality provided by this early version of the Lithics3D package. Hopefully it will be considerably more polished in the near future.

Feature tours:

1. Measuring edge angles on a flake

Automatic angle measurements can be made along an artifact's edge using the edgeAngles function and some landmarks (can be as few as two!) to define the edge. The following demonstrates how this may be accomplished using one of the 3D scans and landmark sets included in the package.

First, let's load the data. We will use demoFlake1 here and we will use the first three included landmarks to define the edge we're interested in:

library(Lithics3D)
library(rgl)

data(demoFlake1)
lms = demoFlake1$lms[c(1:3),]

shade3d(demoFlake1$mesh, color="green")
points3d(lms, color="red", size=10)

The landmarks we picked to identify an edge don't do a very good job of describing said edge. If we wanted to measure edge angles at 30 locations along this edge we would use the pathResample function to identify 30 equidistant points along the path, as shown below, but with so few landmarks we would be sampling points along two lines, which don't even run along the surface of the flake:

#edge.curve = pathResample(as.matrix(lms), 30, method="npts")
#points3d(edge.curve, color="blue", size=6)
lines3d(lms, color="blue", lwd=3)

This obviously does not work, because these resampled points do not follow the actual edge of the object. What we really want to do is to use the available landmarks to trace the actual edge of the object in 3D space. This can be accomplished using the sPathConnect connect function to trace the shortest path between the given landmarks on the surface of the flake and following ridges along the way. This function returns the IDs of the mesh vertices rather than actual coordinates, so an intermediary step is required. Let's try again:

edge.curve.ids = sPathConnect(lms, demoFlake1$mesh, path.choice="ridges")
edge.curve.path = t(demoFlake1$mesh$vb)[edge.curve.ids,1:3]
lines3d(edge.curve.path, color="blue", size=6, lwd=4) # Let's see the actual path:
edge.curve = pathResample(edge.curve.path, 30, method="npts")
points3d(edge.curve, color="purple", size=10)

That certainly looks better! So, now we have sampling locations for edge angles, and all that's left to do is to actually measure them using the high-level edgeAngles function, at a distance of 3mm:

close3d()
res = edgeAngles(demoFlake1$mesh, edge.curve, 3)

Done! But how? Let's see using the edge_angles_vis3d and edge_angles_vis2d functions. The first of these functions shows where the measurements were taken in 3D, while the second produces 2D plots of the mesh slices (technically, the endpoints of the mesh triangle edges intersected by the plane on which a given angle was measured) and the angle measuring points.

# Visualize the measurements in 3D:
window_id <- edge_angles_vis3d(res, demoFlake1$mesh, edge.curve)

# Generate figures with 2D plots of the angle measurements, with 6 plots per
# figure, arranged in 2 columns and 3 rows:
angle_plots <- edge_angles_vis2d(res, demoFlake1$mesh, ncol = 2, nrow = 3)

# Show the first figure:
angle_plots[[1]]

# Close the 3D window once done visualizing:
close3d(window_id)

As you notice, angles were not measured at the first and last point along the edge curve, because at the ends of the curve it is impossible to determine where the measurements should be performed.

2. Isolating the ventral surface and taking a closer look:

3D scans can be segmented into multiple components with the mesh_segment_by_path function and a set of landmarks defining a closed path to slice the mesh with. Here we will use a set of landmarks set around the perimeter of a test flake scan (demoFlake2) to isolate the ventral surface of the latter, orienting it so that it's aligned with the XY plane, and examining the thickness/size of the bulb of percussion:

First, let's load and see what the data looks like:

data(demoFlake2)
shade3d(demoFlake2$mesh, color="green")
points3d(demoFlake2$lms, color="red", size=10)

Now, let's split the flake connecting the available landmarks by following surface ridges:

library(Morpho)
rgl.close()
m.seg = mesh_segment_by_path(demoFlake2$mesh, demoFlake2$lms)
shade3d(m.seg[[1]], color="green")
shade3d(m.seg[[2]], color="blue")

Done! But how? The mesh_segment_by_path function will first connect the landmarks using the sPathConnect function internally, and then it will use that path to section the mesh. Let's plot that path using the underlying function:

res = sPathConnect(demoFlake2$lms, demoFlake2$mesh, path.choice="ridges")
res = c(res, sPathConnect(demoFlake2$lms[c(nrow(demoFlake2$lms), 1), ],
                          demoFlake2$mesh, path.choice="ridges"))
lines3d(t(demoFlake2$mesh$vb)[res, ], color="yellow", lwd=6)

Now we can take a closer look at the isolated ventral surface, but it makes sense to orient it first so that we can plot it in 2D as well. The surface can be oriented in many different ways, but here we will use the mesh.orient_by_contour_pca function with landmarks along the edges of the flake (i.e. without the landmark located by the point of percussion).

Let's take a look at the data:

rgl.close()
new.path = sPathConnect(demoFlake2$lms[-nrow(demoFlake2$lms),],
                        demoFlake2$mesh, path.choice="ridges")
shade3d(m.seg[[2]], color="blue")
points3d(t(demoFlake2$mesh$vb)[new.path, 1:3], color="red", size=10)

Now we orient the isolated surface by performing a PCA on the path, resampling it at 100 equidistant points:

vent.o = mesh_orient_by_contour_pca(t(demoFlake2$mesh$vb)[new.path, 1:3],
                                    npts=100, m.seg[[2]])

Because this surface is now aligned with the XY plane, we can plot it in 2D! Let's see what it looks like, using Z values to colour the surface by thickness, after ensuring the mean elevation or thickness (Z) value of the contour is zero.

library(ggplot2)
tp = data.frame(t(vent.o$mesh$vb)[,1:3])
tp[,3] = tp[,3] - mean(vent.o$contour.res[,3])
sp3<-ggplot(tp, aes(x=PC1, y=PC2, z=PC3, color=PC3)) + geom_point() + coord_fixed() + labs(col="Thickness (mm)")
sp3+scale_color_gradientn(colours=rainbow(5)) + xlab("X-axis (mm)") + ylab("Y-axis (mm)")

3. Map the thickness of a flake across its suface:

Thickness maps of 3D scans can be easily created and visualized using the mesh_tmap and mesh_tmap_plot functions, as demonstrated here with the included demoFlake2 data.

For thickness maps to be meaningful we must first orient the 3D scan in a meaningful way. What "meaningful" is will depend on the application. Here, we will orient the artifact using the Morpho::pcAlign function, although we could just as easily orient it as we did above with mesh.orient_by_contour_pca, with alignAxis, etc. Morpho::pcAlign provides a good way to automatically orient flat-ish and somewhat elongated objects such as blades, points, etc.

mesh.rot = Morpho::pcAlign(demoFlake2$mesh)

After this step, which requires some thought in real applications, we can produce a thickness map of the oriented object. The mesh_tmap function checks the thickness of the mesh by first overlying a grid, and then measuring the maximum distance between mesh vertices located in individual grid cells. If we request a high resolution grid on a low resolution mesh we may have no valid values, so the function will automatically lower the resolution of the grid until only a small number of cells, defined by the ld.cutoff parameter, have less than 10 vertices mapped to them. The starting resolution is defined by the base.res parameter.

OK, so with that clarification out of the way, let's proceed to create a thickness. We will try to obtain a resolution of 0.2mm, but we will specify that only 0.3% of the underlying grid cells may have less than 5 values. The actual resolution will be in the output plot.

library(ggplot2)

ld.cutoff = 0.003 # 0.3% of grid cells with less than 5 values.
base.res = 0.1 # Baseline resolution for thickness maps (0.2mm)
tmap <- mesh_tmap(mesh.rot, base.res, ld.cutoff)
g <- mesh_tmap_plot(tmap$tmap, 12, y.label="y-axis",
               x.label="x-axis", annot=paste("res:", tmap$res))
g + labs(title="Thickness map!")