 Methodology article
 Open Access
 Published:
Alpha shapes: determining 3D shape complexity across morphologically diverse structures
BMC Evolutionary Biologyvolume 18, Article number: 184 (2018)
Abstract
Background
Following recent advances in bioimaging, highresolution 3D models of biological structures are now generated rapidly and at lowcost. To use this data to address evolutionary and ecological questions, an array of tools has been developed to conduct shape analysis and quantify topographic complexity. Here we focus particularly on shape techniques applied to irregularshaped objects lacking clear homologous landmarks, and propose a new ‘alphashapes’ method for quantifying 3D shape complexity.
Methods
We apply alphashapes to quantify shape complexity in the mammalian baculum as an example of a morphologically disparate structure. Micro computedtomography (μCT) scans of bacula were conducted. Bacula were binarised and converted into point clouds. Following application of a scaling factor to account for absolute size differences, a suite of alphashapes was fitted per specimen. An alpha shape is formed from a subcomplex of the Delaunay triangulation of a given set of points, and ranges in refinement from a very coarse mesh (approximating convex hulls) to a very fine fit. ‘Optimal’ alpha was defined as the refinement necessary in order for alphashape volume to equal CT voxel volume, and was taken as a metric of overall ‘complexity’.
Results
Our results show that alphashapes can be used to quantify interspecific variation in shape ‘complexity’ within biological structures of disparate geometry. The ‘stepped’ nature of alpha curves is informative with regards to the contribution of specific morphological features to overall ‘complexity’. Alphashapes agrees with other measures of complexity (dissection index, Dirichlet normal energy) in identifying ursid bacula as having low shape complexity. However, alphashapes estimates mustelid bacula as being most complex, contrasting with other shape metrics. 3D fractal dimension is identified as an inappropriate metric of complexity when applied to bacula.
Conclusions
Alphashapes is used to calculate ‘optimal’ alpha refinement as a proxy for shape ‘complexity’ without identifying landmarks. The implementation of alphashapes is straightforward, and is automated to process large datasets quickly. We interpret alphashapes as being particularly sensitive to concavities in surface topology, potentially distinguishing it from other shape complexity metrics. Beyond genital shape, the alphashapes technique holds considerable promise for new applications across evolutionary, ecological and palaeoecological disciplines.
Background
The morphology of an organism is both a function of its evolutionary past and its adaptation to present surroundings. Quantifying morphology is fundamental to the study of ecology and evolution. Organisms are, quite literally, shaped by their evolutionary history. And morphology is often the only source of evidence upon which phylogenetic relationships may be reconstructed through deep time. Morphology also plays important role in linking the phenotype to ecology, by establishing causal relationships between anatomy and performance [1]. Flexible tools for quantifying organismal morphology are therefore highly desirable amongst users spanning the disciplines of ecology and evolutionary biology [2]. More broadly, the comparison of morphological features is of interest to applied scientists from a diverse array of background, including archaeology, chemistry, computer science and medicine.
The morphology of an organism and its component parts can be described in terms of size, shape, structure, colour and patterning. Of these, shape has historically been difficult to consistently and objectively quantify, and this challenge forms the basis of the field of morphometrics. ‘Shape’ can be defined as all the geometric information contained within an object, once the effects of rotation, translation and scale have been removed [3]. Traditionally shape has been quantified as a series of single measurements, including ratios [4] and angles [5]. Such measures clearly ignore a wealth of potential shape data however.
Shape variation vs. shape complexity
Associated with recent advances in specimen digitization, a suite of new techniques has been developed to analyse 3D biological shape data. Chiefly, these methods facilitate the quantification of variation in form (size and shape) among specimens using multivariate methods, allowing for either the study of covariation between shapes or between shapes and extrinsic factors [6]. Analyses of biological shape variation most often fall within the paradigm of geometric morphometrics (GMM).
GMM studies typically proceed via the identification of homologous morphological landmarks across a range of specimens, and subsequent Procrustes superimposition to remove the effects of translation, rotation and scale. The placement of landmarks on sutures, muscle attachment scars and tuberosities is therefore common. This approach has proved effective in ecological and evolutionary studies across a range of biological structures, including vertebrate skulls [7], insect wings [8] and tree leaves [9]. More problematic is the landmarking of less featured objects, such as the diaphyses of long bones [10], the body of ribs [11], otoliths [12], seeds [13] and anthropological artefacts [14].
A class of related outline or surfacebased shape analysis tools exist however, that do not necessarily require homologous landmarks to be defined a priori. ‘Eigenshapes’ [15], ‘Eigensurfaces’ [16, 17], ‘Canonical Sampling’ [18, 19], fully automated landmarking (‘auto3Dgm’) [20, 21], ‘Elliptical Fourier Analysis’ [22,23,24,25] and ‘Spherical Harmonics’ [26,27,28,29] are all important contributions to the morphometricians toolbox of shape analysis techniques. In all such cases, the principle goal of the analysis remains the same: to describe the shape of objects, and the specific ways in which objects differ in shape between themselves and as a function of external factors.
Yet a second suite of morphometric techniques seeks to quantify shape complexity. Within the field of biology, complexity may be broadly defined as the number of ‘parts’ comprising an organism or a landscape (be that genes, cell types, organ systems or habitat patches). In the context of shape analysis, here we focus on topographic complexity. Whilst topographic ‘complexity’ has numerous definitions across the literature (see below), complex shapes can intuitively be thought of as those formed by combining parts, or the entirety, of several simple ‘primitive’ shapes. Complexity indices may differ in the specific ‘aspect’ of shape complexity captured, ranging from the degree of selfsimilarity (fractal architecture) displayed, to simpler metrics of surface rugosity. Shape complexity has found numerous important applications within the disciplines of ecology and evolutionary biology. Root complexity may be indicative of the health of a plant [30], whilst tooth complexity has been used to predict the palaeodiet of fossil vertebrates [31]. The complexity of landscape patches has been linked to habitat quality [32], and invertebrate genital complexity has been interpreted in the context of sexual selection mechanisms [22].
It is important to reiterate that the two suites of morphometric techniques highlighted above measure two very different aspects of form (namely shape variation vs. complexity), such that two objects may occupy very similar GMM morphospace whilst being characterised by different values of shape complexity. Two outwardly similar surface meshes with similar landmark configurations may differ in shape complexity if the surfaces deviate in terms of surface rugosity, for example.
Methods for quantifying shape complexity
Several metrics have been advanced for the quantification of spatial or topographic complexity within biological systems:
Dissection index and relief index
In two dimensions, dissection index (DI) is the ratio of an object’s perimeter to the square root of its area. Dissection index is therefore a dimensionless number, providing an indication of the extent to which a shape is more complex than a circle [33]. In three dimensions, the related Relief Index (RFI) is calculated as the ratio of an object’s surface area to its planar area [34], and thus provides an index of rugosity. Both metrics are simple to calculate and intuitive to understand, and represent singleparameter shape descriptors of complexity. A corollary of this however, is neither metric provides an indication of the distribution of complexity across an object. Furthermore, the value of planar area incorporated into RFI is necessarily orientationdependent (total planar area is dependent upon on the orientation of the object relative to the observer when the plan view is taken). When applied to tooth complexity, the preferred orientation is obvious; planar area is calculated in the occlusal plane when RFI is taken as a proxy for hysodonty [34]. Should RFI be extended beyond tooth crown complexity to other biological structures however, the orientation of planar area will need further consideration. Additionally, the calculation of RFI requires a mesh from which to derive surface area, involving an intermediate processing stage for point clouds or voxel based data.
Fractal dimension
The fractal dimension (FD) is a measure of complexity applicable to objects that are selfsimilar (exhibiting repetitive patterns across scales) [35]. FD metrics have commonly been applied to physical landscapes [36] in addition to biological organisms perceived to display selfsimilarity, including plant roots [37], plant leaves [38], stony corals [39], and brain structures [40]. Simply speaking, the fractal dimension captures the ability of an object to fill the Euclidean space within which it is located. The most common implementation of FD applies a ‘boxcounting’ approach, in which a regular grid of boxes of side length s is overlain across the 2D data and the number of occupied boxes counted as N(s). This process is repeated whilst varying the size of s. log N(s) is subsequently plotted as a function of log (1/s), and the slope of this graph is taken as an estimate of FD [41]. Whilst originally implemented on 2D data, fractal analysis has since been extended to operate in 3D [39, 40].
Fractal analysis is undoubtedly a powerful tool that provides an objective and scaleindependent single metric of shape complexity. However, numerous caveats have been expressed when applying FD to biological datasets (see [41] for a review). Most notably, when an object or pattern is not obviously selfsimilar, the application of fractal dimensions can be problematic [42]. Indeed, rather than being truly selfsimilarity, some authors have gone so far as to suggest that most ‘complex’ structures differ in their extent of selfsimilarity across spatiotemporal scales, and are actually best described as selfdissimilar [43]. Furthermore, the value of the fractal dimension for a given outline is a function of several ‘somewhat arbitrary’ decisions, including the location of the grid starting point and the selected values of minimum and maximum s [44]. Within the ecological literature, occupancy is typically calculated across s values spanning ~two orders of magnitude, yet such a limited scaling relationship cannot be taken as strong evidence of genuine fractality [41].
Dirichlet Normal energy
Dirichlet normal energy (DNE) effectively quantifies the ‘curviness’ of a mesh. Most simply “DNE measures the deviation of a surface from being planar” ([45]; p249), and ranges from zero in the case of a flat plane, to higher values associated with steep crests and troughs. It is calculated as the sum of energy values across all faces of a mesh surface, where the energy value at each face is quantified as changes in the normal map. The process does not require the assignment of landmarks, and is unaffected by scale or orientation. Additionally, energy is calculated for every face of the mesh, facilitating energy variation to be visualized across the surface of the object. In this way, DNE allows for specific regions of ‘high’ and ‘low’ complexity to be identified across a specimen. Thus far, DNE has found extensive use in the mammal tooth literature [45,46,47], but has also been applied to quantify the shape complexity of developing embryos [48]. DNE does however, require a mesh a priori, and has been shown to be sensitive to commonlyused mesh preprocessing operations such as smoothing and decimating [49].
Alphashapes
In this study, our objective is to develop a straightforward method for quantifying threedimensional shape complexity that is orientationindependent, does not require assumptions of selfsimilarity or an intermediate meshing stage, and is capable of quantify topographic complexity across multiple scales. Our approach is based on the concept of ‘alphashapes’.
An alphashape is formed from the boundary of an alphacomplex, which is itself a subcomplex of the Delaunay triangulation for a given set of points [50]. For a given set of points in space, a family of alphashapes may be defined, ranging from a very coarse (a convex hull) to very fine fit around said points (Fig. 1). The parameter ‘alpha’ dictates the level of refinement, with a larger alpha resulting in coarser fits and a smaller alpha in finer fits (see [50] for a comprehensive description of 3dimensional alphashapes). The level of refinement necessary in order for an alphashape’s volume to match that of the original dataset to which it is fitted may be taken as a measure of shape complexity: more complex objects will require a more refined alphashape fit in order for volumes to converge.
The resulting alphashape fit may comprise one volume (larger alphas; Fig. 1ce) or multiple smaller volumes (smaller alphas; Fig. 1f). Hence, as alpha decreases, the refinement of the fit changes from a convex hull (the special case when sphere radius is infinite, Fig. 1b) to finer fits as more regions are removed by smaller spheres (Fig. 1ce). Eventually the radius of the sphere decreases to such an extent that no points are intersected and no alphashape is created.
The convex hull (Fig. 1b) and ‘coarser’ alphashapes (Fig. 1cd) occupy a volume equal to or larger than that of the underlying object (Fig. 1a). In contrast, very fine alphashapes (Fig. 1f) will have a volume smaller than the original structure. At some ‘optimal’ level of refinement, alphashape volume and specimen volume will be equal (Fig. 1e), and it is this ‘optimal’ alpha upon which we base our metric of 3D shape complexity.
Within the biological sciences, alphashapes have previously been used to describe characteristics of protein surface shape [51], to segment forested areas from aerial LiDAR data [52] and to describe the spatial distribution of fish within schools [53]. In a practical sense, alphashapes has been implemented in the freeware ‘Meshlab’ [54] as a means of generating surface meshes from point cloud data. The authors have previously applied an alphashapes approach to the problem of body mass estimation in fossil species [55]. In this implementation, a predictive relationship between alpha shape volume and body mass was derived from a suite of articulated modern mammal skeletons digitised using LiDAR. The predictive model was subsequently applied to extinct mammal taxa and their fossil body mass estimated. To the authors knowledge, alphashapes has not previously been applied to explicitly quantify shape complexity however.
Genital shape complexity
Within the field of evolutionary biology, genital form and function has received considerable attention, albeit with a heavy bias towards invertebrates [56]. Genitals are amongst the most diverse, complex and rapidly evolving structures observed in living organisms [57]. Genital shape, rather than size, is often used by taxonomists as a means of distinguishing between closely related species [58, 59], implying greater divergence in genitalic shape than size [22]. Indeed, numerous experimental evolution studies have found direct evidence for sexual selection acting on genital shape across a range of taxa [23, 60, 61].
There is therefore considerable interest in developing automated methods capable of quantifying shape across such complex and diverse structures as animal genitalia. In some instances, traditional landmarkbased GMM techniques have been applied [60, 62,63,64,65]. Such studies frequently consider genital shape variation intraspecifically, or between morphologically similar sister taxa [66, 67]. Yet elsewhere, GMM methods have been applied to broader interspecific samples of genitalia [65, 68, 69], highlighting the applicability of these techniques to quantify shape change in rapidly evolving structures, or those comprised entirely/predominantly of soft tissue [70].
Here we use the mammalian baculum as a test case for the application of alphashapes to qualifying morphological complexity. In the past, bacula have been used as a taxonomic character to differentiate between otherwise indistinguishable sister taxa, such is their morphological disparity between closely related species. Whilst this is predominantly true for rodents and bats [71], baculum morphometrics have also been developed as a diagnostic tool for differentiating between species of carnivore [72]. As far as the authors are aware, a traditional geometric morphometric analysis of baculum shape has not been attempted between species however, potentially due to difficulties associated with the identification of discrete homologous landmarks. The development of a simple and intuitive method for quantifying ‘complexity’ in the mammal baculum in the absence of homologous landmarks therefore has the potential to reinvigorate the study of mammal genital evolution. In addition, the present study is of significance for both ecologists and evolutionary biologists (and those working more broadly in the fields of archaeology and computer science) who will benefit from a new tool for comparing 3D shape complexity across samples of extreme shape diversity.
Methods
Raw data
Twelve mammalian bacula were scanned as an example dataset using micro computed tomography (μCT). The taxa include three families of modern Carnivora (Mustelidae, Canidae and Ursidae; Table 1) and span a range of shapes (Fig. 2) from simple rodlike bones (Ursidae) to complex curved, grooved and notched structures (Mustelidae).
CT scans were conducted at Manchester Xray Imaging Facility using a Nikon 320/225 kV Custom Bay microCT instrument, and the Natural History Museum London using a Nikon 225 kV microCT instrument. Raw CT scans were converted to binary data in ImageJ by automated thresholding according the histogram of raw CT grayscale values. Binarised CT scans were read into MATLAB R2017a (The MathWorks Inc., Natick, MA, USA) slice by slice, and any internal cavities present were filled using two separate automatic gap and hole filling algorithms, (imclose.m and imfill.m) from MATLAB’s Image Processing toolbox. Imclose.m performs a morphological closing on each binary image slice, using a 2D disc of a given radius. In this instance, 6 pixels was found to be the minimum radius that consistently closed the periosteal contour across the sample. Imfill.m identifies holes as being background pixels that cannot be reached from the edge, and subsequently flood fills them with foreground pixels. The relative ‘hollowness’ of bacula has not previously been described, yet all specimens included in the present study did possess internal void spaces. Here we chose to focus on the shape complexity of the external morphology, and hence filled any internal cavities. Nevertheless, the alphashapes technique will function equally well for instances when the internal geometry is pertinent to the research question.
Having filled internal void spaces, CT data were converted directly to point clouds. This highlights an important advantage of the alphashapes approach as a means of directly calculating shape ‘complexity’ from a CT dataset via the process of surface meshing, rather than requiring a surface mesh beforehand. The process of converting CT volumes to surface meshes necessarily involves some degree of smoothing to avoid faceting and topological artefacts resulting from image artefacts and noise. This process ought to be, but is rarely, documented in the metadata [73] and the effect of smoothing on subsequent data analysis is seldom explored. Raw point clouds were generated by designating the xyz coordinates of every voxel in the CT segmentation associated with the baculum as being a single point in space. That is, unlike surfacebased point clouds generated by other popular digitisation techniques such as LiDAR (light detection and ranging) or photogrammetry, here point clouds also comprise ‘internal’ points representing the solid infilled bone. Raw point clouds were randomly downsampled to 100,000 points each, ensuring all specimens were represented by equally sized datasets (but see ‘Sensitivity Analysis’ below).
Alphashapes
Alpha and reference length
Prior to fitting alphashapes, the issue of scale must be dealt with. Here we are interested in quantifying shape ‘complexity’ in the absence of potential size signals. Alpha radii are calculated in the same units as the underlying point clouds, therefore an alpha radius of 100 mm may entirely enclose one smaller specimen yet only half of another larger specimen, for example. Size normalisation may be achieved in one of two ways: either by scaling all point clouds to the same size, or by scaling alpha radii to the overall size of each specimen.
In this implementation of alphashapes, we choose the later. In doing so, the underlying point cloud remains the same scale throughout. A specimen with a maximum length of 100 mm will remain represented by a ~ 100 mmlong point cloud. The resulting alpha shapes mesh, comprising all the triangles formed when the points contributing to the alpha shape are connected, will likewise remain at this original scale and may then be used in downstream functional analyses (see Discussion). Therefore, if two objects are identical in shape and both comprise an equal number of points, yet one is twice the size, the larger specimen requires an alpha radius (α) twice as large to ensure an equal refinement of fit. To calculate the alpha radius for each specimen we use the following equation:
where α is the alpha radii, k is the refinement coefficient and l_{ref} is the point cloud reference length (as described in the following section). Here we are interested in identifying an ‘optimal’ level of alphashape refinement (see below) and therefore chose 200 values of refinement coefficient k that result in alphashapes ranging from coarse fits (convex hulls) to highly refined shapes. Refinement coefficients ranged from 0.1 to 10,000 and were evenly distributed on a logarithmic scale. At the smallest values of refinement, the alphashape ceases to be one continuous volume and the sphere passes inside the point cloud to create multiple small volumes, hence no longer representing the overall shape of the object.
Scaling reference length
The point cloud reference length l_{ref} is a scaling factor allowing equivalencies to be drawn between alphashapes fitted to specimens of absolute different size, as discussed above. Yet arriving at a single ‘reference’ length that adequately describes the overall size of a point cloud is nontrivial. A simple approach is to use the maximum diagonal of the bounding box as a reference length. Alternatively, an earlier implementation of alphashapes as a mass estimation tool [55] utilised the average distance of all points from the centroid of the point cloud. Here, we investigate a third technique, in which the average distance of each point to its nearest 100 neighbours in the downsampled point cloud is used as a descriptor of overall point cloud size. Ultimately, the nearest neighbour technique was preferred as this resulted in alphashapes ‘breaking down’ (i.e. becoming multiple small volumes) at the same refinement coefficient (Fig. 3  see asterisk *), implying alpha radii is well scaled to the relative distance between points in the point cloud, and therefore the overall size of the specimen.
Optimal refinement coefficient
Having calculated alpha radii using the above equation, alphashapes were fitted to point clouds using the MATLAB ‘alphavol’ function written by Jonas Lundgren (http://www.mathworks.co.uk/matlabcentral/fileexchange/28851alphashapes), which both calculates the fit of the alphashape and its associated volume. Alphashapes were fitted for a range of refinement coefficient across all specimens, and volumes extracted. All analyses were run on a laptop computer with 8GB 1600MHz DDR3 RAM and a 1.1GHz Intel Core M processor.
Each specimen is described by a representative curve of alphashape volume against refinement coefficient. As alphashapes become more refined (smaller refinement coefficients), their associated volumes decrease. However, the profile of this alpha curve is a function of the shape complexity of the bone, from gross overall shape (straight vs curved, for example), to specific morphological features (such as grooves or forked tips) and ultimately surface texture/roughness (i.e. pitted or smooth). The alpha curve often has a stepped appearance, with steep regions corresponding to a sudden reduction in alpha volume when alpha radius becomes sufficiently small so as to represent a particular feature or surface texture.
For each specimen, we therefore identify the ‘optimal’ refinement coefficient best reflecting overall shape complexity by comparing alpha volume against ‘raw’ volume. ‘Raw’ volume is an estimate of the biological volume of the specimen, as calculated from the hole filled CT data by multiplying the number of threshold voxels by scan resolution cubed, prior to point cloud downsampling. The refinement coefficient producing an alpha volume that is closest to raw volume will be taken as the ‘optimal’ refinement. To identify the ‘optimal’ refinement coefficient an optimisation approach was undertaken using the ‘fminsearch’ function of MATLAB’s optimisation toolbox, which applies a ‘NelderMead’ search method. The optimisation routine searches for the refinement coefficient that produces the smallest difference between alpha volume and raw volume. This process continues until two conditions have been satisfied: the difference between volumes (alpha volume vs raw volume) is less than 1e4, and the difference between subsequent values of refinement coefficient is less the 1e4. The final refinement coefficient after both conditions have been satisfied is taken as ‘optimal’.
Using our approach, 3D shape complexity is reduced first into one curve per specimen and ultimately into one refinement value per specimen (see supplementary material for MATLAB code). We predict that simple rodlike structures will require a relatively coarser refinement coefficient (relatively larger alpha radii) to match ‘raw’ volumes compared to complex, curved or grooved specimens that will require a more refined alpha shapes (relatively smaller alpha radii) to accurately represent total volume.
Comparison to other shape complexity measures
Here we compare alphashapes to three additional metrics of topographic complexity commonly applied in the field of evolutionary biology. Firstly, we calculate the orientationindependent 3D ‘dissection index’ (DI) which represents the ratio of the squared root of surface area to the cubed root of volume. 2D dissection indices have previously been applied to quantify shape complexity in invertebrate genitalia [22], and here we modify this technique to work on 3dimensional data. Isosurface meshes (comprising 10,000 faces, see DNE section below) were generated from the binarised .raw CT stack in Horos [74], decimated in Geomagic (3D Systems, North Carolina, USA) and surface areas and volumes calculated using the ‘compute geometric measures’ function in Meshlab [54].
Threedimensional fractal dimension (FD) was estimated using a boxcounting algorithm written in MATLAB by Frederic Moisy (https://uk.mathworks.com/matlabcentral/fileexchange/13063boxcount). The function was applied to the binarized CT data after the internal cavities present were filled, but prior to conversion to a point cloud and downsampling (see ‘Raw data’ section above). The function calculates the number of cubes required to cover the baculum N(s) at sequential sizes of box, where the size of the cube s is the length of one side. The slope of the relationship of log(1/s) to log N(s) is taken as an estimate of the FD of the object, where objects with higher topographic complexity have a higher values FD.
Dirichlet Normal Energy (DNE) was calculated in the R package ‘MolaR’ following the methodology of Pampush et al. [47], using the same surface meshes as produced for the 3D DI calculation above. As per previous applications, meshes comprised 10,000 faces [47]. Higher values of DNE are indicative of higher topographic complexity.
Sensitivity analysis
In theory, alpha shapes can cope with infinitely detailed point clouds, yet practically the number of points comprising an object will be dictated by several factors. The scanning technique used can impact the density of the point cloud, with μCT scans often producing very dense point clouds compared to LiDAR or photogrammetry (although this does strongly depend upon the specifics of a given imaging setup). Larger point clouds necessitate longer computational times, which may problematic for large comparative studies. More importantly, the particular research question ought to have a large bearing on the density of the point cloud. If the question under investigation pertains to ‘gross’ morphology, a less dense point cloud may be justifiable, whereas those focusing on features of surface texture may require more detail. Whilst final point cloud size is ultimately determined by the users’ needs, ensuring that all specimens within a comparative study comprise an equal number of points is necessary in order that one sample is not represented in significantly more detail than another, which may potentially skew the results.
We therefore conducted a sensitivity analysis to examine the effect of point cloud density on calculated values of ‘optimal’ refinement coefficients (and associated computing time). Optimal refinement coefficients were calculated for each specimen comprising points cloud sizes ranging from 10^{4} to 10^{6} points, typical for datasets derived from LiDAR, photogrammetry or CT. Reference lengths (see Eq. 1) of the 10^{4} point clouds were used to scale alpha radii for all point cloud sizes, ensuring consistent alpha radii (at each refinement) between point cloud sizes and that results are equivalent.
Results
The alphashapes methodology described here distils the complexity of three dimensional baculum shape, firstly into a single representative curve and ultimately into a single parameter to facilitate further comparative analysis. We consider the shapefitting protocol to be straightforward and relatively computationally inexpensive when operating on point clouds of ~ 100,000 data points. For a typical specimen (Mustela itatsi,14 MB 8 bit raw file), data import and hole filling took 14 s, the calculation of reference length on the basis of 100 nearest neighbours took 14 min, and the calculation of the optimal refinement coefficient took 2 min.
In specimens appearing outwardly similar, the relationship between alphashape volume and refinement coefficient is characterised by similar profiles. Ursid bacula, for example, share a simple rodlike appearance which is smooth and lacking in features such as grooves, curvature or complex apices (Fig. 2 and Fig. 4), and likewise the four bear bacula share similar alpha curves (Fig. 3). At the highest value refinement coefficients (approaching a convex hull), ursid alphashapes overestimate raw volume by ~ 25–50%, and only the outermost points of the point cloud contribute to shape fitting. As refinement coefficients decrease, divergence between alphashape volume and raw volume is quickly reduced, and ‘optimal’ alpha is reached (occurring at refinement coefficients between 11 and 36; Table 2). Beyond which, alphashape volume decreases with refinement coefficient at a slower rate, until the alphashape fit breaks down to form several disconnected volumes (refinement coefficients below 0.6).
Whilst also lacking distinct curvature or a complex tip, the canid baculum does possess a welldeveloped broad urethral groove on the ventral surface (Fig. 2 and Fig. 4). Optimal refinement coefficients of canid bacula are therefore intermediate between those of ursids and mustelids (Table 2). It follows that these specimens have a more complex relationship between alphashape volume and refinement coefficient, with curves taking on a multistepped appearance (Fig. 3). Steps coincide with refinement coefficient values becoming small enough to allow specific morphological features to be detailed. At high values of refinement, alphashapes overestimate canid baculum raw volume by ~ 200%. As the refinement coefficient is reduced, canids display a very pronounced ‘step’ (Fig. 3 iii to ii) at a refinement coefficient of ~ 5. This coincides with alpha radii falling below ~half urethral groove width (Fig. 5) and the distinctive feature suddenly being resolved. Optimal alpha occurs soon after at refinement coefficients of 2.7–3.6.
Finally, mustelids require low values of refinement coefficient to accurately represent raw volume, as expected due to their complex geometry. In Mustela itatsi (Fig. 4a) for example, alphashape volumes generated by high values of refinement coefficient vastly exceed raw volume (by a factor of ~ 3), due to the highlycurved nature of the bone. As refinement coefficient is reduced, previously unseen morphological features become apparent. Overall dorsoventral curvature is defined at a refinement of ~ 10 (Fig. 4a iii) whereas more detailed morphological features such as the urethral groove and the rugose proximal portion associated with attachment to the corpora cavernosa become apparent at a refinement of 1.7 (Fig. 4a ii). Due to the overall complex shape of this structure, alphashape volume converges upon raw volume to produce an ‘optimal’ refinement fit at low values of refinement.
Both DI and DNE complexity metrics agree with the alpha shapes methodology presented here in finding ursids to possess low complexity bacula (Table 2, Fig. 6). Alphashapes is the only metric tested here in which all taxonomic groups are entirely differentiated from each other on the basis of surface complexity. In contrast, there is very little differentiation between family groupings when baculum complexity is quantified by fractal dimension (Table 2, Fig. 6). Average DI and DNE values of canid bacula exceed those of mustelids, reversing the trend present in optimal alpha.
Sensitivity analysis
The results of the sensitivity analysis indicate that optimal refinement coefficients decrease with increasing point cloud size (Fig. 7a). Less dense point clouds require relatively coarser refinement coefficients in order to produce alpha shapes of equal volume to the original dataset. This phenomenon has previously been documented elsewhere, and has been referred to as the ‘coastline paradox’ ([75], see Discussion). Between 10^{5} and 10^{6} points, the rank order of optimal refinement coefficients remains relatively consistent across taxa (Fig. 7a). At the lowest point cloud densities, canids are considered the most ‘complex’, whilst mustelid bacula would appear most complex at point cloud densities of ~ 10^{5} points (Fig. 7a). This simply reflects the scale at which shape complexity is present. Canid bacula possess ‘gross’ complexity (e.g. presence of a deep urethral groove), whilst mustelid bacula are characterised by a more refined level of complexity (e.g. a shallow urethral groove, curved tip and complex apices) which may only be recovered at higher point cloud densities. The time taken to compute optimal refinement coefficients increases dramatically between 10^{5} points (1–2 min per specimen) and 10^{6} points (15–35 min per specimen) (Fig. 7b).
Discussion
The alphashapes methodology presented here represents an additional tool for quantifying 3D shape complexity across biological samples characterised by high morphological disparity. Alphashapes operates by converting thresholded CT data directly to point clouds, thereby removing the requirement to surface mesh structures beforehand. The alphashapes algorithm does produce a suite of surface meshes as an output however, which may be incorporated into subsequent functional analyses. For example, the impact of the canid urethral groove on the biomechanical performance of the baculum may be quantified by constructing a suite of finite element models, based on coarser (groove absent) alphashapes and finer (groove present) alphashapes. The alphashapes algorithm is implemented in programming languages including MATLAB (‘alphaShape’) and R via the ‘alphahull’ package [68], thereby facilitating greater automatization in the future. Furthermore, alphashapes functionality is also present in the freeware software ‘Meshlab’ [54] for those preferring a graphical user interface.
A recent phylogenetic reconstruction of mammalian baculum presence/absence found support for the independent evolution of the structure on 8–9 occasions, with at least two independent gains of baculum within primates [76]. As alphashapes does not require the placement of homologous landmarks, it may therefore be extended to the analysis of potentially analogous structures or used to quantify shape complexity through ontogenetic sequences. We do also urge caution against the a priori assumption of analogous baculum function for mammals however, as no consistent relationship has yet been identified linking features of the baculum to underlying organismal biology across the whole group.
Here we find agreement between alphashapes, DI and DNE techniques in identifying ursid bacula as possessing low topographic complexity (Table 2, Fig. 6). This is perhaps unsurprising, as bear bacula lack both grooves/ridges/curvature at a macro scale and possess a smooth surface texture on a finer scale. In contrast, alpha shapes departs from the other complexity metrics in classifying mustelid bacula as more complex than canids, a pattern that is reversed in DI and DNE values (Table 2, Fig. 6). Disagreement between metrics of shape complexity is not unprecedented [22], and suggests the methods are simply capturing different aspects of complexity.
In this instance, we interpret these differences as being due to the relative sensitivity of each metric to concave versus convex topology. In DI and DNE, any change in topology (concave or convex) will contribute approximately equally to the complexity metric. In contrast, the calculated optimal alpha appears to be more influenced by the presence of concave sections. In Fig. 5, for example, alpha shapes fitted to the convex dorsal surface of the canid baculum change very little across two orders of magnitude in refinement coefficient. In contrast, the form of the alpha shape fitted to the highly concave ventral margin varies substantially alongside refinement coefficient. Thus, for specimens possessing large concave surfaces such as the urethral groove or distal tip curvature, small values of refinement coefficient are necessary for said features to be resolved. We consider alphashapes complexity to therefore be weighted more towards gross concave features than corrugatedlike surface rugosity, in which convex and concave sections occur with approximately equal frequency and magnitude.
Relative to other metrics of topographic complexity considered here, fractal dimension does not distinguish between taxonomic groupings (Table 2, Fig. 6). Indeed, Fig. 8b would suggest carnivore bacula do not exhibit selfsimilarity, and the application of FD to this structure is not justified. In the boxcounting technique applied here, true fractal behaviour would be identified by a ‘plateauing’ in local slope values across several scales of boxsize [77]. As can be seen in Fig. 8b, no such plateaus exist, and bacula cannot be considered to behave in a fractal manner across several orders of magnitude scale.
Whilst the alphashapes method is not heavily user intensive, the process of shapefitting can be computationally costly. To expedite the process, point clouds are downsampled. However, our sensitivity analysis does indicate that optimal refinement coefficients are a function of point cloud density (Fig. 7a). In denser point clouds, surface textural information (such as attachment scars, and small fossae) are preserved and a finer ‘fit’ around such features in necessary in order to recreate the original volume. At lower point cloud densities, only gross morphology is preserved and a coarser ‘optimal’ refinement coefficient is sufficient.
This effect is related to a well know phenomena known as the ‘coastline paradox’ [65], in which the length of a country’s coastline increases as the scale of the measuring unit is decreased. Intuitively, more features of a coastline can be resolved and incorporated into a metric of length when using a shorter ‘measuring stick’. In the case of alphashapes, as point cloud density is downsampled, the likelihood of removing points lying on the outer contour is increased. As the outermost points define the margins of the specimen, downsampling results in an apparent ‘smoothing’ of the object and hence a coarser optimal refinement coefficient. To illustrate this effect, the alphashapes methodology was applied to a 2D point cloud map of Great Britain (Fig. 9). The results mirror our baculum dataset, with denser point cloud maps requiring more refined alphashapes in order to match the original area (Fig. 9). Low density point clouds loose many of the finer features of the coastline (thin peninsulas, bays etc.) and only gross shape is preserved.
Furthermore, the sensitivity analysis highlighted a change in the rank order of species’ optimal refinement coefficients associated with downsampling between 10^{4} and 10^{5} points (Fig. 7a). As discussed above, this pertains to the ‘hierarchy’ of complexity which may be revealed at a given point cloud size. Canids possess ‘gross’ complexity which may be resolved in low resolution point clouds, whilst mustelids are characterised by concave features of microcomplexity which require higher density point clouds to be detected. Beyond 10^{5} points, rank orders are relatively stable yet computational time increases dramatically (Fig. 6b).
Ultimately, point cloud density will be at the discretion of the user. This is not unusual, and similar decisions are made (implicitly or explicitly) whenever selecting the required resolution of a digital photograph or μCT scan. As a rule of thumb in μCT scanning, voxel size must be at most onequarter to onethird of the size of the feature of interest in order to resolve said feature and avoid partial volume effects. Similarly, to guarantee their inclusion in an alphashapes analysis, we recommend the minimum dimension of a given feature (for example, the width of a groove or diameter of a fossae) comprise at least 3–4 data points within the point cloud. Beyond this, the final point cloud density will reflect a compromise between the level of detail required by the user and computer processing time.
That ‘optimal’ refinement coefficients are a function of point cloud density is not problematic for the application of alphashapes within a comparative analysis framework. Minimum point cloud density should be dictated by the smallest feature of interest across the whole sample, and all specimens downsampled to this same degree. Thus, ‘optimal’ refinement coefficients are equivalent across a given dataset. The absolute values of refinement coefficients will be specific to that given dataset however.
In addition, the current implementation of alphashapes is limited in the sense that betweensubject variation in alpha volume can be difficult to ascribe any one particular geometric feature. Figure 10 represents an initial attempt to address this shortcoming, in which data points of the point cloud are coloured according to the coarsest alphashape to which they contribute. The urethral groove of the canid requires a similar level of refinement in order to be resolved as the curved tip of the mustelid (Fig. 10, green), and would therefore be consider equally ‘complex’ in the current implementation of alpha—shapes. In contrast, the urethral groove of the mustelid required a more refined ‘fit’ of alpha in order to be distinguished (Fig. 10, red), contributing to the low values of optimal alpha calculated for all mustelids here. Future implementations of alphashapes will seek to further quantify regional variation in shape complexity within specimens, and will explore means of extracting additional information from alphacurves.
Conclusions
The alphashapes methodology presented here is an important addition to the biologist’s tool kit, providing a metric of topographical complexity that complements and extends preexisting techniques such as Dissection Index, Fractal Index and Dirichlet Normal Energy. In particular, alphashapes appear sensitive to surface concavities, with features such as grooves and pits considered particularly ‘complex’. Alphashapes differs from methods that have previously been applied to genital shape, such as GMM and spherical harmonics, in that it describes the extent to which an object is structural complex, as opposed to how objects differ in the positioning of particular features. We therefore consider alphashapes to be especially useful for measuring the functional properties of shapes, be those animal genitals, corals, or the occlusal surfaces of teeth. Because optimal alpha values reflect the topographical complexity of a surface, rather than the specifics of how that complexity is achieved, it does not require the placement of homologous landmarks and may therefore be used to compare shape complexity across unrelated structures.
Abbreviations
 CT:

Computed tomography
 DI:

Dissection index
 DNE:

Dirichlet normal energy
 FD:

Fractal dimension
 GMM:

Geometric morphometrics
 LiDAR:

Light detection and range
 SPHARM:

Spherical harmonics
References
 1.
Lauder GV. Functional morphology and systematics: studying functional patterns in an historical context. Annu Rev Ecol Syst. 1990;21:317–40.
 2.
Klingenberg CP. MorphoJ: an integrated software package for geometric morphometrics. Mol Ecol Resour. 2011;11.
 3.
Kendall DG. Shape manifolds, procrustean metrics, and complex projective spaces. Bull London Math Soc Oxford University Press. 1984;16:81–121.
 4.
Andermann S. The cnemic index: a critique. Am J Phys Anthropol. 1976;44:369–70.
 5.
Brauer G. Osteometrie. In: Krussmann R, editor. Anthropol I. Stuttgart: Fischer Verlag; 1988. p. 160–232.
 6.
O’Higgins P, Cobb SN, Fitton LC, Gröning F, Phillips R, Liu J, Fagan MJ. Combining geometric morphometrics and functional simulation: an emerging toolkit for virtual functional analyses. J Anat. 2011;218:3–15.
 7.
Dumont M, Wall CE, BottonDivet L, Goswami A, Peigné S, Fabre AC. Do functional demands associated with locomotor habitat, diet, and activity pattern drive skull shape evolution in musteloid carnivorans? Biol J Linn Soc. 2016;117:858–78.
 8.
Ray RP, Nakata T, Henningsson P, Bomphrey RJ. Enhanced flight performance by genetic manipulation of wing shape in drosophila. Nat Commun. 2016;7:10851.
 9.
Klein LL, Caito M, Chapnick C, Kitchen C, O ‘Hanlon R, Chitwood DH, et al. Digital Morphometrics of two north American grapevines (Vitis: Vitaceae) quantifies leaf variation between species, within species, and among Individuals Front Plant Sci 2017;8:373.
 10.
Frelat MA, Katina S, Weber GW, Bookstein FL. Technical note: a novel geometric morphometric approach to the study of long bone shape variation. Am J Phys Anthropol. 2012;149:628–38.
 11.
Weaver AA, Schoell SL, Stitzel JD. Morphometric analysis of variation in the ribs with age and sex. J Anat. 2014;225:246–61.
 12.
Ponton D. Is geometric morphometrics efficient for comparing otolith shape of different fish species? J Morphol. 2006;267:750–7.
 13.
Ros J, Evin A, Bouby L, Ruas MP. Geometric morphometric analysis of grain shape and the identification of tworowed barley (Hordeum vulgare subsp. distichum L.) in southern France. J Archaeol Sci. 2014;41:568–75.
 14.
Buchanan B, O’Brien MJ, Collard M. Continentwide or regionspecific? A geometric morphometricsbased assessment of variation in Clovis point shape. Archaeol Anthropol Sci. 2014;6:145–62.
 15.
Generalizing MN. Extending the Eigenshape method of shape space visualization and analysis. Paleobiology. 1999;25:107–38.
 16.
Polly PD, Macleod N. Characterization and comparison of 3D shapes using eigensurface analysis: locomotion in tertiary carnivora. Palaeontol Electron. 2008;11:1–13.
 17.
Polly PD. Adaptive zones and the pinniped ankle: a threedimensional quantitative analysis of carnivoran tarsal evolution. In: Mammalian evolutionary morphology. Dordrecht: Springer; 2008. p. 167–96.
 18.
Parr WCH, Ruto A, Soligo C, Chatterjee HJ. Allometric shape vector projection: a new method for the identification of allometric shape characters and trajectories applied to the human astragalus (talus). J Theor Biol Elsevier. 2011;272:64–71.
 19.
Parr WCH, Soligo C, Smaers J, Chatterjee HJ, Ruto A, Cornish L, et al. Threedimensional shape variation of talar surface morphology in hominoid primates. J Anat. 2014;225:42–59.
 20.
Boyer D, Lipman Y, Clair ES, Puente J, Funkhouser T, Patel B, et al. Algorithms to automatically quantify the geometric similarity of anatomical surfaces. Proc Natl Acad Sci U S A. 2011;108:18221–6.
 21.
Boyer DM, Puente J, Gladman JT, Glynn C, Mukherjee S, Yapuncich GS, et al. A new fully automated approach for aligning and comparing shapes. Anat Rec. 2015;298:249–76.
 22.
Rowe L, Arnqvist O. Sexual selection and the evolution of genital shape and complexity in water striders. Evolution. 2012;66:40–54.
 23.
Arnqvist G, Danielsson I. Copulatory behavior, genital morphology, and male fertilization success in water striders. Evolution. 1999;53:147–56.
 24.
Arnqvist G. Comparative evidence for the evolution of genitalia by sexual selection. Nature. 1998;393:784–6.
 25.
Holwell GI. Geographic variation in genital morphology of Ciulfina praying mantids. J Zool. 2008;276:108–14.
 26.
Shen L, Farid H, McPeek MA. Modeling threedimensional morphological structures using spherical harmonics. Evolution. 2009;63:1003–16.
 27.
McPeek MA, Shen L, Farid H. The correlated evolution of threedimensional reproductive structures between male and female damselflies. Evolution. 2009;63:73–83.
 28.
McPeek MA, Symes LB, Zong DM, McPeek CL. Species recognition and patterns of population variation in the reproductive. Evolution. 2010;65:419–28.
 29.
McPeek MA, Shen L, Torrey JZ, Farid H. The tempo and mode of threedimensional morphological evolution in male reproductive structures. Am Nat The University of Chicago Press. 2008;171:158–78.
 30.
Wang H, Siopongco J, Wade L, Yamauchi A. Fractal analysis on root Systems of Rice Plants in response to drought stress. Environ Exp Bot. 2009;65:338–44.
 31.
Prufrock KA, Boyer DM, Silcox MT. The first major primate extinction: an evaluation of paleoecological dynamics of north American stem primates using a homology free measure of tooth shape. Am J Phys Anthropol. 2016;159:683–97.
 32.
Imre AR, Bogaert J. The fractal dimension as a measure of the quality of habitats. Acta Biotheor. 2004;52:41–56.
 33.
McLellan T, Endler J. The relative success of some methods for measuring and describing the shape of complex objects. Syst Biol. 2008;47:264–81.
 34.
Boyer DM. Relief index of second mandibular molars is a correlate of diet among prosimian primates and other euarchontan mammals. J Hum Evol. 2008;55:1118–37.
 35.
Mandelbrot BB. The fractal geometry of nature, vol. 460p. New York, NY: W.H. Freeman and Company; 1982.
 36.
AlHamdan M, Cruise J, Rickman D, Quattrochi D. Effects of spatial and spectral resolutions on fractal dimensions in forested landscapes. Remote Sens. 2010;2:611–40.
 37.
Nielsen KL, Lynch JP, Weiss HN. Fractal geometry of bean root systems: correlations between spatial and fractal dimension. Am J Bot. 1997;84:26–33.
 38.
Plotze R, de O, Falvo M, Pádua JG, Bernacci LC, Vieira MLC, Oliveira GCX, Bruno OM. Leaf shape analysis using the multiscale Minkowski fractal dimension, a new morphometric method: a study with Passiflora (Passifloraceae). Can. J Bot. 2005;83:287–301.
 39.
Reichert J, Backes AR, Schubert P, Wilke T. The power of 3D fractal dimensions for comparative shape and structural complexity analyses of irregularly shaped organisms. Methods Ecol Evol. 2017;8:1650–8.
 40.
Liu JZ, Zhang LD, Yue GH. Fractal dimension in human cerebellum measured by magnetic resonance imaging. Biophys J. 2013;85:4041–6.
 41.
Halley JM, Hartley S, Kallimanis AS, Kunin WE, Lennon JJ, Sgardelis SP. Uses and abuses of fractal methodology in ecology. Ecol Lett. 2004;7:254–71.
 42.
Parrott L. Measuring ecological complexity. Ecol Indic. 2010;10:1069–76.
 43.
Wolpert DH, Macready W. Using selfdissimilarity to quantify complexity. Complexity. 2007;12:77–85.
 44.
Slice D. The fractal analysis of shape. In: Contributions to Morphometrics. Madrid: Museo national de Ciencias Naturales; 1993. p. 164–90.
 45.
Bunn JM, Boyer DM, Lipman Y, St. Clair EM, Jernvall J, Daubechies I. Comparing Dirichlet normal surface energy of tooth crowns, a new technique of molar shape quantification for dietary inference, with previous methods in isolation and in combination. Am J Phys Anthropol. 2001;145:247–61.
 46.
Winchester JM, Boyer DM, St. Clair EM, GosselinIldari AD, Cooke SB, Ledogar JA. Dental topography of platyrrhines and prosimians: convergence and contrasts. Am J Phys Anthropol. 2014;153:29–44.
 47.
Pampush JD, Winchester JM, Morse PE, Vining AQ, Boyer DM, Kay RF. Introducing molaR: a new R package for quantitative topographic analysis of teeth (and other topographic surfaces). J Mamm Evol. 2016;23:397–412.
 48.
SalvadorMartínez I, SalazarCiudad I. How complexity increases in development: an analysis of the spatialtemporal dynamics of gene expression in Ciona intestinalis. Mech Dev. 2017;144:113–24.
 49.
Spradley JP, Pampush JD, Morse PE, Kay RF. Smooth operator: the effects of different 3D mesh retriangulation protocols on the computation of Dirichlet normal energy. Am J Phys Anthropol. 2017;163:94–109.
 50.
Edelsbrunner H, Mücke EP. Threedimensional alpha shapes. ACM Trans Graph. 1994;13:43–72.
 51.
Albou LP, Schwarz B, Poch O, Wurtz JM, Moras D. Defining and characterizing protein surface using alpha shapes. Proteins Struct Funct Bioinforma. 2009;76:1–12.
 52.
Eysn L, Hollaus M, Vetter M, Mucke W, Pfeifer N, Regner B. Adapting αshapes for forest delineation using ALS data. 10th Int Conf LiDAR Appl assess for. Ecosystems. 2010:14–7.
 53.
Carette V, Mostafavi M, Devilliers R, Rose G, Hashemi Beni L. Extending marine GIS Capabilities : 3D dynamic and interactive representation of fish aggregations using Delauney Tetrahedralisation and alpha shapes. Geomatica. 2008;62:247–56.
 54.
Cignoni P, Corsini M, Ranzuglia G. Meshlab: an opensource 3d mesh processing system. ERCIM News. 2008;73:45–6.
 55.
Brassey C, Gardiner J. An advanced shapefitting algorithm applied to the quadrupedal mammals: improving volumetric mass estimates. R Soc Open Sci. 2015;2:150302.
 56.
AhKing M, Barron AB, Herberstein ME. Genital evolution: why are females still understudied? PLoS Biol. 2014;12:1–7.
 57.
Eberhard WG. Sexual selection and animal genitalia: Harvard University Press; 1985.
 58.
Song H, Wenzel JW. Mosaic pattern of genital divergence in three populations of Schistocerca lineata scudder, 1899 (Orthoptera: Acrididae: Cyrtacanthacridinae). Biol J Linn Soc. 2008;94:289–391.
 59.
Eberhard WG. Static allometry and animal genitalia. Evolution. 2009;63:48–66.
 60.
Hopwood PE, Head ML, Jordan EJ, Carter MJ, Davey E, Moore AJ, et al. Selection on an antagonistic behavioral trait can drive rapid genital coevolution in the burying beetle, Nicrophorus vespilloides. Evolution. 2016;70:1180–8.
 61.
HeinenKay JL, Langerhans RB. Predationassociated divergence of male genital morphology in a livebearing fish. J Evol Biol. 2013;26:2135–46.
 62.
Arnqvist G, Thornhill R, Rowe L. Evolution of animal genitalia: morphological correlates of fitness components in a water strider. J Evol Biol. 1997;10:613–40.
 63.
Holwell GI, Winnick C, Tregenza T, Herberstein ME. Genital shape correlates with sperm transfer success in the praying mantis Ciulfina klassi (Insecta: Mantodea). Behav Ecol Sociobiol. 2010;64:617–25.
 64.
Simmons LW, GarciaGonzalez F. Experimental coevolution of male and female genital morphology. Nat Commun. 2011;2:374.
 65.
Macagno ALM, Pizzo A, Parzer HF, Palestrini C, Rolando A, Moczek AP. Shape  but not size  Codivergence between male and female copulatory structures in Onthophagus beetles. PLoS One. 2011;6:e28893.
 66.
Pizzo A, Mercurio D, Palestrini C, Roggero A, Rolando A. Male differentiation patterns in two polyphenic sister species of the genus Onthophagus Latreille, 1802 (Coleoptera: Scarabaeidae): a geometric morphometric approach. J Zool Syst Evol Res. 2006;44:54–62.
 67.
Dinca V, Dapporto L, Vila R. A combined geneticmorphometric analysis unravels the complex biogeographical history of Polyommatus icarus and Polyommatus celina common blue butterflies. Mol Ecol. 2011;20:3921–35.
 68.
Mutanen M, Pretorius E. Subjective visual evaluation vs. traditional and geometric morphometrics in species delimitation: a comparison of moth genitalia. Syst Entomol. 2007;32:371–86.
 69.
Parzer HF, Polly PD, Moczek AP. The evolution of relative trait size and shape: insights from the genitalia of dung beetles. Develop Gen Evol. 2018;228:83–93.
 70.
Orbach DN, Hedrick B, Würsig B, Mesnick SL, Brennan PL. The evolution of genital shape variation in female cetaceans. Evolution. 2017:261–73.
 71.
Clark WK. The baculum in the taxonomy of Peromyscus boylei and P. pectoralis. J Mammal. 1953;34:189–92.
 72.
Vercillo F, Ragni B. Morphometric discrimination between Martes martes and Martes foina in Italy: the use of the baculum. Hystrix. 2011;22:325–31.
 73.
Davies TG, Rahman IA, Lautenschlager S, Cunningham JA, Asher RJ, Barrett PM, et al. Open data and digital morphology. Proc R Soc B Biol Sci. 2017;284:20170194.
 74.
Project Horos. DICOM image viewing and measuring. 2015. Available from: http://www.horosproject.org/. Accessed 1 Dec 2018.
 75.
Mandelbrot B. How long is the coast of Britain? Statistical selfsimilarity and fractional dimension. Science. 1967;156:636–8.
 76.
Schultz NG, LoughStevens M, Abreu E, Orr T, Dean MD. The Baculum was gained and lost multiple times during mammalian evolution. Integr Comp Biol. 2016;56:644–56.
 77.
Bouda M, Caplan JS, Saiers JE. Boxcounting dimension revisited: presenting an efficient method of minimizing quantization error and an assessment of the selfsimilarity of structural root systems. Front Plant Sci. 2016;7:1–15.
Acknowledgements
The authors would like to thank Andrew Kitchener and Georg Hantke (National Museum of Scotland, Edinburgh) and Louise Tomsett and Amin Garbout (Natural History Museum, London).
Funding
CB is funded by a BBSRC Future Leader Fellowship (BB/N010957/1). The Manchester Xray Imaging Facility is funded in part by the EPSRC (grants EP/F007906/1, EP/F001452/1 and EP/102249X/1). The funders had no role in the study design, data collection and analysis, decision to publish, or preparation of the manuscript.
Availability of data and materials
All MATLAB code and raw CT data are available on Figshare at https://doi.org/10.6084/m9.figshare.5558557
Author information
Affiliations
Contributions
JG, JB and CB conceived the ideas and designed the methodology; JB and CB collected the data; JG developed the code; CB led the writing of the manuscript. All authors contributed critically to the drafts and gave final approval for publication.
Corresponding author
Correspondence to Charlotte A. Brassey.
Ethics declarations
Ethics approval and consent to participate
Not applicable.
Consent for publication
Not applicable.
Competing interests
The authors declare that they have no competing interests.
Publisher’s Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Rights and permissions
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.
About this article
Received
Accepted
Published
DOI