gdi

library(gdi)

Basic Workflow for Volume Estimation

This vignette demonstrates the basic workflow of the gdi package for estimating the volume of an (extinct) animal based on multiview photos or reconstructions.

We start with a set of image files showing orthogonal views of an animal for which a volumetric estimate is desired. These image files should be aligned and with the long axis horizontal, and complex features that result in a major violation of the assumption of elliptical or superelliptical cross-sections (e.g. limbs or fins that are protruding from the silhouette) should be removed and estimated separately.

Example images are provided with the package in the folder exdata. To import them, execute:

fdir <- system.file(package="gdi")

measurements_lateral <- measuresil(file.path(fdir,"exdata","lat.png"))
measurements_dorsal <- measuresil(file.path(fdir,"exdata","dors.png"))

In this case the images are saved in rgba format and with transparent background, which is recommended. We therefore do not need to change the default options for threshold (0.5), method (“greater”) and channel (4) provided with the function. The default selection will automatically count all pixels with an opacity (alpha channel value) of more than 0.5 as part of the silhouette and treat everything else as background.

Verify that our measurement worked and resulted in two vectors of equal length:

length(measurements_lateral)
#> [1] 2341
length(measurements_dorsal)
#> [1] 2341

We can now perform our graphic double integration to estimate the volume:

gdi(measurements_lateral, measurements_dorsal, scale=100, method="raw")
#> x_dim_23.35_units 
#>          40.25677
gdi(measurements_lateral, measurements_dorsal, scale=100, method="smooth")
#> x_dim_23.35_units 
#>          40.25668

The result of this function is the estimated volume in the chosen unit of measurements. The silhouette is at a scale of 100px:10cm, so a scale value of 100 will report the volume in liters (cubic decimeters), namely 40.2566754 for this example. The principal difference between “raw” and “smooth” for the parameter “method” is that “raw” approximates the body as a stack of elliptical cylinders (“classic” gdi), whereas “smooth” approximates it as a sequence of elliptical frusta (with base areas based on the diameters of segments i and i+1).

To illustrate the difference, consider the following example:

sil <- c(0,1)
gdi(sil, sil, method="raw", scale=1)
#> x_dim_1_units 
#>     0.7853982
gdi(sil, sil, method="smooth", scale=1)
#> x_dim_1_units 
#>      1.047198

Details of the gdi() function As we can see, method “raw” gives the volume of a cylinder with a diameter of 1 and a length (because scale = 1) of 1. On the other hand, method “smooth” (any other string will also work), estimates the volume as a series of frusta, with base areas based on each segment and the one succeeding it. In this case, because segment 1 has a diameter of 0, this results in the estimated volume being that of a cone with length 1 and base diameter 1 and a cylinder for segment two.

This difference in methods has very minor effects in analyses using pixel-precise measurements of digital images with reasonable levels of resolutions (in the case of the first example, the difference is only ca. 9.9^{-5} l), but can have a notable effect at lower resolutions (i.e. with fewer, longer segments, such as occur when performing gdi on manual measurements).

Cross-sectional corrections In the event that cross-sectional geometry of our shape deviates markedly from an ellipse, the package provides two approaches to account for this. The first is the specification of a superellipse exponent as part of the function call to gdi():

gdi(measurements_lateral,measurements_dorsal,k=2.3, scale=100)
#> x_dim_23.35_units 
#>          42.26999

In this case we specified an exponent of k=2.3, which results in a cross-sectional area, and thus a volume, approximately 5% greater than a standard ellipse. For a review of cross-sectional geometry in vertebrates, see Motani (2001).

The function sellipse.coo() used in conjunction with the standard plotting functions can be used to visualize possible cross-sectional shapes based on their superellipse exponents:

sellipse.coo(2.0)->ellipse
plot(ellipse$x,ellipse$y,col="grey", type="l")
polygon(ellipse$x,ellipse$y,col="grey", border="grey")

sellipse.coo(2.3)->se2.3
lines(se2.3$x,se2.3$y, col="blue") #plot a superellipse with exponent 2.3

sellipse.coo(3)->se3
lines(se3$x,se3$y, col="red") #plot a superellipse with exponent 3

The second option is to use a graphical reconstruction of the cross-sectional shape to calculate a correction factor determined from the ratio in areas between the actual cross-section, and an ellipse matching its vertical and horizontal maximum diameters (which correspond to the diameters that would be measured from silhouettes for the GDI). Functionality for determining this correction factor from a silhouette image of the cross-section is implemented in the cscorr()-function:

fdir <- system.file(package="gdi")
correction_factor <- cscorr(file.path(fdir,"exdata","cross_section.png"))
print(correction_factor)
#> [1] 1.092215

The resulting correction factor can then be multiplied by the volume estimated by the uncorrected gdi()-results, or directly supplied to the gdi()-function using the ‘corr’ parameter.

Both the k and corr parameters can either be a single value applicable to the entire volume, or a numeric vector containing a correction factor for every segment (this vector needs to be of the same length as the diameter measurements used with the gdi()-function, otherwise an error will be returned). If desired, interpolation between different cross-sectional geometries over the length of the model (e.g. using approx()) might be useful in order to further increase accuracy of the estimation.

The volume of more complex shapes can be estimated using this technique, by manually splitting the silhouettes into several pairs of image files, each containing two orthogonal views of each body part. Total volume can then be estimated as the sum of all body parts. Alternatively, if orthogonal views are not available for every body part, simplified assumptions can be made, e.g. based on known height-width ratios. For example:

gdi(measurements_lateral, measurements_lateral*0.9, scale=100)
#> x_dim_23.35_units 
#>          43.85454

Here it is assumed that the body width is, on average, 0.9 times the body depth. Such an assumption can be error-prone, but, if used appropriately, and grounded in a solid understanding of the morphology of the taxon under study, might provide reasonable approximations of body volumes even where no two orthogonal views are available.

The cscorr()-function provides a way of estimating an appropriate factor to use here when specifying the “aspect_ratio” setting for parameter ‘return’:

aspect_ratio <- cscorr(file.path(fdir,"exdata","cross_section.png"))
print(aspect_ratio)
#> [1] 1.092215

Here the estimated depth/width-ratio is 1.09, so an appropriate width depth ratio would be 1/1.09=0.916.

Complex Shapes

Complex protruding structures that violate the assumption of roughly elliptical cross-sections are best estimated separately. Notably this includes limbs.

hindlimb_lateral <- measuresil(file.path(fdir,"exdata","hl.png"),align="v")
forelimb_lateral <- measuresil(file.path(fdir,"exdata","fl.png"), align="v")

By setting the parameter “align” to “v”, we can specify that the silhouettes in question are aligned vertically, i.e. that the elliptical slices should be taken horizontally. In this case, we only have lateral views of the limbs, so we will make a simplified assumption of the transverse diameters being approximately 0.7 times the anteroposterior diameters.

gdi(hindlimb_lateral,0.7*hindlimb_lateral,scale=100)->hindlimb
gdi(forelimb_lateral, 0.7*forelimb_lateral, scale=100)->forelimb
gdi(measurements_lateral, measurements_dorsal, scale=100, method="raw")->axial_total
knitr::kable(data.frame(axial_total, forelimb, hindlimb))

To get the total body volume, we can now simply sum up our volume (taking care to multiply paired structures by two):

axial_total+2*forelimb+2*hindlimb
#> x_dim_23.35_units 
#>          42.42419