Automatic generation of objective footprint outlines

Model import, orientation, and contour line generation

Import of 3D-models is achieved using the vcgImport function of the Rvcg package (Schlager, 2017), which supports commonly used formats including the widely used PLY. The supplied 3D-model should contain only a single complete footprint as well as a margin of surrounding surface. The exact extent of the surrounding margin will not alter results except when an increased margin size includes additional large-scale continuous contours that can be mistakenly recognized as forming part of the footprint margin. After import, the script rotates the point cloud of the extracted vertex point xyz-coordinates to fit the horizontal plane (i.e., the tracking surface on which the animal walked) by employing principal component analysis (PCA) on the three variables (x,y,z). The PCA fits three orthogonal axes (PC1–3) to the point cloud. PC1 is defined as the axis of greatest variation, with PC2 and PC3 capturing successively less variation. In most situations, PC1 and PC2 will represent the horizontal plane (the plane of greatest variation), and PC3 the relief (i.e., the deviation from the horizontal plane). Problematically, the point cloud can get mirror-inverted during PCA fitting as the signs of the columns of the rotation matrix are arbitrary, a problem also occurring in respective implementations in 3D-mesh software like Meshlab (tested with version v2016.12) or CloudCompare (tested with version 2.9). The present script calculates the Procrustes distance (a measure of shape difference) of a subset of points of the model before and after the PCA fit, and will mirror back when detecting a significant difference. If the footprint is a cast (convex hyporelief) rather than a mold (concave epirelief), it will be automatically mirrored into a mold. The scale of the input 3D model is preserved throughout the process; measurements of the generated outline can be taken from both the plots and coordinates returned by the script.

All subsequent computations are based on contour lines of 30 equally spaced height levels extracted from the point cloud using the “getContourLines” function of the contoureR R package (Hamilton, 2015, Fig. 1B). Contour lines reduce the complex three-dimensional problem to a simpler and easier-to-handle two-dimensional one, and form the natural basis for 2D footprint outlines. Before the objective outline can be extracted, a number of additional steps are required, including (1) the removal of contours not representative for the footprint wall and (2) the establishment of homology between the points of separate contours.

Contours not representative for the footprint wall are excluded based on simple criteria. First, all open contours are removed, eliminating structures that continue beyond the margins of the model. Second, only the longest contour of each height level is selected and kept, respectively, with all others removed. This results in a stack of continuous contours, with one contour per height level. Third, all contours less than 50% of the length of the longest contour are removed, while assuring that no gaps within the stack are being created. This approach eliminates smaller structures within the footprint that are unlikely to contain relevant information on the footprint wall (Fig. 1C). An option is implemented that allows for processing multiple impression per model, which is useful in cases where the footprint is not defined by a single outline. For each additional impression, the contour selecting procedure is repeated with the contours selected for the previous stacks excluded.

The resulting stack of contours may still include a number of contours that convey little or no information on the footprint wall, including roundish contours around the actual footprint. To eliminate these contours as well, and to limit the height range under consideration, areas of all contours are computed as a measure of form difference. Starting from the middle contour of the stack, the differences in area of each contour with its next lowest (or highest) neighbor are compared; if the difference in area between two contours exceeds a pre-defined threshold-value, the upper (or lower) of this contour and all following contours are removed (Fig. 1D). Different threshold-values are defined for the lower and the upper half of the stack. An option for adjusting these values is available, allowing to influence how many lower or upper contours are to be included, possibly changing the height level of the resulting calculated outline. All footprints presented in this work were calculated using the default parameters.

Homologization of contours

Even if the starting point would correspond between all contours and if each contour would contain an equal number of equidistant points (requirements not fulfilled a priori), the individual points of the separate contours would tend to deviate from each other when far from the starting point, as the shapes of the contours are not identical. For this reason, when producing a simple mean shape, points would be averaged obliquely rather than perpendicularly to the footprint wall, leading to erroneous results. The implemented solution detects a number of “bottlenecks”—pairs of points with minimum distance between the innermost and outermost contour. Points forming the bottleneck will be considered homologous (define a line that is assumed to be perpendicular to the footprint wall), and the points in-between the bottlenecks will be interpolated by resampling.

First, all contours are resampled to the same number of equidistant points, using 500 points per default. The resulting contours can be variously oriented clockwise or counter-clockwise; contours are reversed accordingly to achieve uniform orientations. Second, Euclidean distances between all possible pairs of the inner and outer contour of the stack are calculated and stored in a matrix with the dimensions n x n. The pair with the minimum distance, the first bottleneck, is then extracted, and those points of the intermediate contours are detected that are closest to a line defined by the bottleneck points. The resulting set of homologous points is then defined as the starting point of the contours. Third, additional bottlenecks are detected to establish homology. The implemented algorithm first detects a second bottleneck on the side of the footprint opposite to the first bottleneck; two additional bottlenecks are then found on each side between the first and second bottleneck. More bottlenecks are detected within the intermediate sections if the latter are long enough. In all cases, bottlenecks in proximity to existing bottlenecks are prevented, assuring a roughly equal distribution of bottlenecks along the outline (Fig. 1D). Finally, the individual sections between the bottlenecks are resampled to equal numbers of equidistant points, which can now be considered homologous.

Tracing of the steepest slope

The objective outline will be traced along the steepest slope of the track wall. For each point within each set of homologous points, the minimum distance between the neighboring contours is measured. A set of homologous points is not always fully perpendicular to the footprint wall, especially when the section between the bottleneck points is long and contours differ much in orientation. For this reason, the algorithm does not simply compute the distances within the set of homologous points, but the distances between each of the homologous points and all nearby points within and outside of the set. The steepest slope computed this way (Fig. 1E) is seldom continuous along the whole outline, but rather tends to fade out and continue on a different height level, frequently leading to abrupt steps in the outline that are obviously incompatible with human interpretations. For this reason, the algorithm does not only detect the steepest point, but takes into account the steepness at all other points. Then, the final coordinate of the steepest slope is computed as the weighted arithmetic mean of all these points. Points will receive equal weight only when the steepness is equal; the lower the steepness compared to the steepest slope, the less weight is given.

The continuous and smooth final outline is produced by applying an approximating function to the resulting set of points. Of various tested options, including B-splines and Bézier curves, elliptic Fourier transforms were found to produce results most similar to those expected from a human interpreter (Figs. 1F, 2). Elliptic Fourier transforms are performed using the “efourier” function of the R package Momocs (Bonhomme et al., 2014) using 25 harmonics and 10 smoothing iterations.

As the resulting objective outline and contour stack will be rotated arbitrarily, an algorithm attempts to rotate both outline and stack into an upright orientation. This algorithm is based on the assumption that digit impressions are facing upwards and require a longer contour segment (i.e., more equidistant points) to be described. In a first step, the centroid, or center of mass, of the outline is computed. Then, the length of the part of the outline above the centroid is calculated for all rotation angles (1–360°), and the outline rotated according to the angle that maximizes this length.

Automatic landmark placement

The script, as described above, can be applied to any kind of depression, as no a priori assumptions on trackmaker anatomy are introduced. Such assumptions are required when linear and angular measurements are to be extracted or when different outlines are to be aligned for quantitative shape analysis. Therefore, a function is included that returns landmark coordinates for tridactyl footprints, which are the most common dinosaur footprints and one of the most common footprint types in general. The algorithm first detects and separates the three digit impressions, assuming that the outline was correctly rotated in a more-or-less upright orientation in the previous step. Then, the rotation angle of the outline is further refined by considering the central digit impression only, which is commonly considered to approximate the mid-axis of the footprint (Leonardi et al., 1987). A total of six landmarks are defined as in Lallensack, Van Heteren & Wings (2016), on the tips of the three digit impressions, on the two hypex points, and on the heel. The landmarks on the tips of the digit impressions are defined as the distal ends of the digital axes; this definition reduces the influence of claw impressions, which may vary strongly depending on preservation and behavior (Lallensack, Van Heteren & Wings, 2016). The digital axes were computed by rotating the respective digit into the upright position and taking the mean of the x-coordinates; the intersection with the outline was found using the rgeos R package (Bivand & Rundel, 2018). The two hypex points were computed by finding the lowest point between the enclosing digit impressions relative to a line connecting the landmarks at the tips of these digit impressions. The landmark on the heel region is the intersection of the axis of the central digit impression and the proximal margin of the outline.

Case study and sensitivity analysis

The applicability of the presented approach for the quantitative analysis of footprints was tested on 13 footprints pertaining of a single theropod trackway (T3) from the Lower Cretaceous of Münchehagen, Germany (3D-models are provided as data S2). Geometric morphometric analysis of the same footprints was conducted in earlier studies based on interpretive drawings (Lallensack, Van Heteren & Wings, 2016; Wings, Lallensack & Mallison, 2016), suggesting that (1) the landmarks on the hypex positions and on the heel are more variable than the landmarks on the digital tips, and that (2) the lateral hypex is more variable than the medial hypex. In the present case study, the calculated objective outlines were aligned using Generalized Procrustes Analysis (GPA) with the geomorph R package (Adams, Collyer & Kaliontzopoulou, 2018). The mean shape of the aligned shapes was then compared with the meanshape published by Lallensack, Van Heteren & Wings (2016). This is expected to reveal potential systematic differences between the traditional interpretive approach and the automated approach presented herein. Furthermore, landmarks are automatically extracted from all 13 generated outlines as described above, and their variability in y-direction (parallel to the footprint mid-axis) computed, in the same way as has been done in Lallensack, Van Heteren & Wings (2016) based on interpretive outlines.

Limitations

Human-made interpretational footprint drawings aim to capture as much information about the trackmaker anatomy as preserved. Although slope steepness is the most important criterion, the steepest slope will frequently fade out along the footprint wall to continue at a different height level, making interpolation unavoidable. Furthermore, humans tend to make a priori assumptions about trackmaker anatomy when producing the outlines, which allows them to take into account extramorphological (unrelated to the foot anatomy) features and include anatomical features of interest such as digital pad impressions and claw marks.

The present algorithm is successful in detecting and interpolating outlines even when the steepest parts of the slope are indistinct (Fig. 2). It does currently not include a priori assumptions that would emphasize anatomically important details and account for extramorphological features, which keeps the algorithm simple and predictable, and applicable to a wide range of footprint types. The broad applicability is demonstrated in Fig. 2C with a hominin footprint from the famous Laetoli tracksite of Tanzania (Leakey & Hay, 1979). In the absence of a priori assumptions, however, the outlines expectedly tend to provide less information on the presumed foot anatomy than interpretative drawings, and may be unusable in cases where track morphology is obscured by extramorphological features. Artifacts caused by a crack running transversally through the central digit impression are shown in Fig. 2A for a medium-sized theropod footprint from the Late Jurassic Courtedoux–Bois de Sylleux tracksite of Switzerland (Castanera et al., 2018).

Additional limitations currently arise from the necessity to limit the vertical extent of the contour stack, which may exclude relevant anatomical features not captured by the selected contours. Furthermore, the described approach requires that impressions can be described by single contours, which is not the case in all cases, especially when the footprint is very shallow and indistinct. However, a preliminary option is implemented to process footprints that are composed of more than one impression, as demonstrated with an ornithopod footprint attributed to Caririchnium kyoungsookimi from the Jindong Formation in Goseong County, South Korea (Figs. 2B, 2D data provided by Anthony Romilio). This footprint is relatively shallow (maximum depth is 3.8% of maximum length of model) and comprises three separate impressions.

Qualitative comparisons with human-made outlines

The similarity of the generated outlines with human-made interpretations is demonstrated on five fossil footprints that have been previously published in the literature (Figs. 2D–2I). All five examples are compared with published interpretative drawings that had been produced using the same 3D-models; all examples, as is the case for examples used elsewhere in this work, were processed using the same script version and parameters. Footprints T3/47, T3/37, and I1-35 (Figs. 2D–2F, respectively) come from the Lower Cretaceous of Münchehagen, Germany (Lallensack, Van Heteren & Wings, 2016; Wings, Lallensack & Mallison, 2016). All three footprints were left by the right foot. T3/47 and T3/37 were part of a larger theropod and I1/35 (Fig. 2F) of an ornithopod trackway. The objective outlines (red continuous lines) are generally in accordance with the interpretational drawings published in Lallensack, Van Heteren & Wings (2016) (dotted blue lines). However, the sediment bars between the digital impressions tend to be less extensive (e.g., Figs. 2E–2F), and digital impression IV in T3/47 is somewhat abbreviated in the objective outline due to sediment infilling in the distal tip of the impression. The generally good match between the generated outlines presented herein and the interpretive ones of Lallensack, Van Heteren & Wings (2016) may be partly due to the use of similar criteria for the definition of the track margin, which are discussed in detail in the latter study, most importantly the criterion of the steepest slope.

Tracks 5 and 6 of QM F10322 (Figs. 2G–2H) are the left and right pedal impressions of a large tridactyl trackway from the Upper Cretaceous of Lark Quarry in Queensland, Australia (Thulborn & Wade, 1984; Romilio & Salisbury, 2011; Romilio & Salisbury, 2014). The 3D-models, based on photographs of the in situ specimen taken in 2013, were provided by Anthony Romilio. The trackway became famous after it was suggested to have been left by a large theropod causing a dinosaur stampede (Thulborn & Wade, 1984), a hypothesis that has been discussed controversially in recent years (Romilio & Salisbury, 2011; Romilio & Salisbury, 2014; Thulborn, 2013; Thulborn, 2017; Romilio, Tucker & Salisbury, 2013), with the identification of the large tridactyl trackmaker as either a theropod or an ornithopod constituting a major point of disagreement. This discussion is instructive in showing how much interpretative outlines can differ when produced by separate researchers with different preconceptions about the responsible trackmaker taxon. Well aware of the subjectivity problem, Romilio & Salisbury (2014, fig. 7F, 8F) did not produce traditional outline drawings but selected a single contour line they considered representative. Still, this approach is not completely objective, as separate contours can differ greatly in shape and dimensions (Falkingham, 2016). The calculated outlines presented herein correspond well with those of Romilio & Salisbury (2014) (Figs. 2G–2H). However, they differ considerably from outlines of the same imprints presented by Thulborn & Wade (1984, plate 17), and especially from a more recent interpretation of track 5 by Thulborn (2017, fig. 5 (3)), which is redrawn here for comparison (Fig. 2I). Track 5 (Figs. 2H–2I) is a prime example of how different interpretational drawings can be when based on fundamentally different assumptions. Given their interpretational nature, it is, on principle, not possible to discard one of these disparate interpretations as incorrect (Falkingham, 2016). Furthermore, the good correspondence of the generated outline with the outline presented by Romilio & Salisbury (2014) could be partly due to the fact that both outlines were produced based on contour lines, while the differing outline of Thulborn (2017) was traced on a photograph. Nevertheless, the presented approach may offer a standard to which interpretive outlines can be compared. It furthermore allows for an objective qualitative comparison of different footprints, as generated outlines would be free of preconceived assumptions on the trackmaker taxon and computed using the same parameters. Although not a replacement for interpretive outlines given the limitations outlined above, generated outlines may be used in combination with interpretive ones, and may form the objective basis for the production of the latter.

Uses in quantitative shape analysis

All 13 footprints of the T3 trackway included in the original geometric morphometric analysis of Lallensack, Van Heteren & Wings (2016), Wings, Lallensack & Mallison (2016) were processed using the same program version and parameters (Fig. S3), and analyzed following the protocol provided by the mentioned studies. The mean shape derived from the Procrustes-aligned generated outlines closely matches the previously published mean shape that was based on interpretive outlines of the same tracks (Fig. 3C). This indicates that, at least in this case, systematic differences in generated and interpretive outlines are minor. Analysis of the variability of the six landmarks (as indicated in Fig. 3A), however, reveals more substantial differences (Fig. 3D). The objective approach confirms that the lateral and medial hypex (landmarks 2 and 4) are more variable than the remaining landmarks. However, the previously published observation that the heel is highly variable is not supported. Furthermore, the central digit (landmark 3) is markedly more variable according to the objective approach. This is partly due to sediment infilling of the distal end of this digit impression (an extramorphological feature) in one of the footprints, which is not taken into account by the objective approach (Fig. 3B). This case study demonstrates that generated outlines can be used to objectively reproduce results derived from interpretive outlines, although the possible influence of extramorphological features and other potential errors in the objective outlines need to be taken into account.

To test for repeatability of results generated by the present approach, one model (track 6 of QM F10322) was saved at successively lower resolutions using MeshLab, and the resulting outlines were compared (Figs. 3E–3F). Observed differences resulting from the separate model resolutions are found to be negligible, even when the highest resolution (Fig. 3E; 196,236 faces, 9.4 MiB) is compared with the lowest tested resolution (Fig. 3F; 5,000 faces, 246.2 KiB).