library(gdi)This vignette explores some other applications for the functions provided in the gdi package.
Estimating airspace proportion (or bone compactness indices) The image in question is a greyscale png, showing a black bone with a white medullary cavity on a grey background. We hence analyze the color values in channel 1. The imghist() function will provide a histogram plot and a table with the number of pixels for different color values
fdir <- system.file(package="gdi")
imghist(file.path(fdir,"exdata","femur.png"), channel=1, unique=TRUE)
#> count proportion
#> 0 86236 0.3923385
#> 0.501960784313725 73463 0.3342266
#> 1 60101 0.2734349The setting unique=TRUE will cause unique color values to be listed in the output. If unique==FALSE, the output will instead be listed for the bins, as defined by the parameter “breaks”.
Alternatively, cscorr() can be used for measuring areas of differently colored pixels:
cscorr(file.path(fdir,"exdata","femur.png"), channel=1, method="less", threshold=0.2, return="area")->compacta #measure compact bone (black pixels)
cscorr(file.path(fdir,"exdata","femur.png"), channel=1, method="greater", threshold=0.8, return="area")->medulla #measure medulla (white pixels)
medulla/(compacta+medulla) #calculate the airspace proportion
#> [1] 0.4107027Estimating second moment of area for cross-sections For a cross-section provided as an image file (or matrix), the second moment of area and polar moment of inertia can be calculated using the function csI().
Consider an example of a square cross-section, here given as a 10×10 matrix representing the pixels of an image:
matrix(rep(1,100),nrow=10,ncol=10)->m #generate sample matrix representing the pixel color values of an image
csI(m) #return the second moments of area (Ix and Iy) and the polar moment of inertia (Iz) of the cross-section
#> I_x I_y I_z
#> 833.3333 833.3333 1666.6667
#> attr(,"centroid")
#> x y
#> 5.5 5.5
#> attr(,"bounding_box")
#> min max
#> x 1 10
#> y 1 10The result here equal the moments of inertia of a square with a side length of 10 unity: 833 = a^4/12 = 10^4/12 = 833
The csI() function is applicable to cross-sections of any geometry, which are approximated by separately calculating moments of inertia contributions for each pixel of the image and summing up the results. For example, the moments for the bone cross-section analyzed above can be calculated as:
csI(file.path(fdir,"exdata","femur.png"), channel=1, method="less", threshold=0.2)
#> I_x I_y I_z
#> 613791415 3393241735 4007033150
#> attr(,"centroid")
#> x y
#> 344.3405 148.8596
#> attr(,"bounding_box")
#> min max
#> x 19 694
#> y 13 303To get results in real-world dimensions, a scale can be provided in the function call. In this case, the bone measures 11.6 cm wide, equating to 676 px in the image, so the appropriate scale to get results in cm⁴ is 676/11.6:
csI(file.path(fdir,"exdata","femur.png"), channel=1, method="less", threshold=0.2, scale=676/11.6)
#> I_x I_y I_z
#> 53.21891 294.21172 347.43063
#> attr(,"centroid")
#> x y
#> 344.3405 148.8596
#> attr(,"bounding_box")
#> min max
#> x 19 694
#> y 13 303Precision of the results depends on the resolution of the images in question. In all cases the orientation of I_x and I_y is based on the orientation of the image given as input.
Estimating moment of (rotational) inertia for volumetric models For this example, we will again use the sample images provided with the package to perform a volumetric estimation (see other vignettes):
fdir <- system.file(package="gdi")
measurements_lateral <- measuresil(file.path(fdir,"exdata","lat.png"), return="all") # using the "return" parameter, we can make measuresil output a data.frame containing the centers on the y axis
measurements_dorsal <- measuresil(file.path(fdir,"exdata","dors.png")) #because our organism is bilaterally symmetric, we are not interested in the location of the COM on the transverse axis, and can therefore disregard recording centers for the dorsal view
gdi(measurements_lateral, measurements_dorsal, scale=100, return="all")->gdiresults # the parameter setting for "return" makes gdi() return a data.frame with all individual segment volumes, rather than only a single volume for the entire shapeIn order to estimate rotational inertia for this volumetric model, we can now use the rotI() function:
rotI(gdiresults, scale=1000, axis="yaw", densities=1)
#> total_mass I_point_masses I_elliptical_cs I_rectangular_cs
#> 40.256774 4.748666 4.851946 4.886372
#> I_corrected_cs
#> 4.851946Note that since desired results will most likely be in kg*m², the scale for gdi() and rotI() will have to be different (in this case, 1000 px = 1 m, but 100 px = 1 dm). The setting for parameter “axis” defaults to “yaw” (or “y”), but “pitch” (or “z”) and “roll” (or “x”) are also possible settings. The parameter “densities” can be used to supply a unit density for the entire model (defaults to 1), or alternatively a vector of segment densities, if these vary throughout the model (e.g. due to the presence of pneumatic structures in some segments).
In addition, the parameter “axis_coord” allows the specification of the coordinate of the axis around which rotation is desired. If not specified, the default is for rotation around the center of mass of the shape.
The output provides three different approximations of the rotational inertia; treating the segments as point masses (definite underestimate), treating the segments as elliptical cylinders, and treating the segments as cuboidal (likely overestimate).
For cross-sections deviating from these shapes, a correction factor to take into account different mass distributions in a cross-section can be applied. These correction factors can be estimated using cscorr(), if a representation of the cross-section as an image is available:
cscorr(file.path(fdir,"exdata","cross_section.png"),return="rotI")->c
c
#> I_corr_x_pitch I_corr_y_yaw I_corr_z_roll
#> 1.114281 1.091465 1.106729The values represent the ratios between the moments of inertia of a shape with the given cross-section, and a hypothetical ellipse with the same horizontal and vertical diameters. Note that it is assumed that mass is constant for this comparison, so accurate segment volumes and densities will still need to be supplied to rotI().
rotI(gdiresults, scale=1000, corr=c[1])
#> total_mass I_point_masses I_elliptical_cs I_rectangular_cs
#> 40.256774 4.748666 4.851946 4.886372
#> I_corrected_cs
#> 4.863749Using rotI() with the corr-argument will correct for the difference in cross-sectional geometry by multiplying the segment moments of inertia calculated using the cylindrical approximation by the correction factor supplied. As we can see in the above example, the output under “I_corrected_cs” is now slightly different from the elliptical approximation. Plausible real-world moments of inertia should be expected to always fall somewhere in between the elliptical and the cuboidal approximations.
Note that for most animal body shapes, the effect of cross-sectional geometry on rotational inertia is likely to be minor, as the anteroposterior distance from the center of rotation is the dominant component.