morphomap: An R package for long bone landmarking, cortical thickness, and cross-sectional geometry mapping

Objectives: This study describes and demonstrates the functionalities and application of a new R package, morphomap , designed to extract shape information as semi-landmarks in multiple sections, build cortical thickness maps, and calculate biomechanical parameters on long bones. Methods: morphomap creates, from a single input (an oriented 3D mesh representing the long bone surface), multiple evenly spaced virtual sections. morphomap then directly and rapidly computes morphometric and biomechanical parameters on each of these sections. The R package comprises three modules: (a) to place semilandmarks on the inner and outer outlines of each section, (b) to extract cortical thicknesses for 2D and 3D morphometric mapping, and (c) to compute cross-sectional geometry. Results: In this article, we apply morphomap to femora from Homo sapiens and Pan troglodytes to demonstrate its utility and show its typical outputs. morphomap greatly facilitates rapid analysis and functional interpretation of long bone form and should prove a valuable addition to the osteoarcheological analysis software toolkit. Conclusions: Long bone loading history is commonly retrodicted by calculating bio-mechanical parameters such as area moments of inertia, analyzing external shape and measuring cortical thickness. morphomap is a software written in the open source R environment, it integrates the main methodological approaches (geometric morpho-metrics, cortical morphometric maps, and cross-sectional geometry) used to parame-trize long bones.

The methods of VA are applicable to skeletal studies in all species.
While landmarks and MMs measure the form of the diaphysis, CSG predicts resistance to particular loading scenarios (e.g., compression, bending, etc.) and has proven effective in retrodicting bone loading history (Endo & Takahashi, 1982;Huiskes, 1982;Piziali et al., 1976;Uhthoff & Jaworski, 1978). Analyses of CSG have been widely applied to studies of locomotion, handedness, body mass, and other aspects of habitual loading.
Here, we present morphomap, an R package specifically designed to extract cross-sections from long bone meshes at specified intervals along the diaphysis and to site semilandmarks on the periosteal and endosteal outlines of each cross-section and to calculate 2D and 3D MM and CSG parameters. We demonstrate the validity of this compu-

| MATERIALS AND METHODS
We tested the functionality of the morphomap R package on a human and a chimpanzee femur. The human femur (OdN 244) is from a Northern Italian Bronze Age necropolis of Olmo di Nogara (Canci, Tafuri, Fornaciari, Cupito, & Salzani, 2011). It was imaged at the hospital of the University of Pisa using a GE LightSpeed RT16 medical CT scanner (slice thickness 0.625 mm. slice increment 0.625 mm, voltage 120 kV, X-ray tube current 300 mA, standard reconstruction algorithm, pixel size 479 μm). The chimpanzee femur (Kupri 1022) was obtained from the Digital Morphology Museum, KUPRI (http://dmm. pri.kyoto-u.ac.jp/dmm/WebGallery/index.html). Its acquisition was performed on a Toshiba Asteio medical CT scanner (slice thickness 1.0 mm, slice increment 0.4 mm, voltage 120 kV, X-ray tube current 150 mA, standard reconstruction algorithm, pixel size 399 μm). The external and internal surfaces extracted from these are included in the morphomap R package.

| The morphomap R package
morphomap is designed to produce by default 61 cross-sections along the diaphysis of the long bone defined at increments of 1% between 20 and 80% of its biomechanical length as it is commonly defined in studies of diaphyseal CSG (Ruff, 2002;Marchi, 2005Marchi, , 2007  At each cross-section, morphomap also calculates the periosteal (ExtP, mm) and endosteal perimeters (MedP), the cortical thickness measured at the four anatomical quadrants (anterior, posterior, medial, and lateral, in mm) and the CSG variables: the total area (TotA, mm 2 ), the medullary area (MedA, mm 2 ), the cortical area (CA, mm 2 ), the area moments of inertia around the y and x axis, (I x and I y , mm 4 ), the maximum and minimum area moments of inertia (I max and I min , mm 4 ), the angle between the major axis and the mediolateral axis of the crosssection, theta (θ, radians), the polar moment of inertia (J, mm 4 ), the section moduli (Z x ,Z y , Z max , Z min , and Z pol , mm 4 ) and the maximum chord lengths from y and x axes (d x and d y , mm; see Supporting Information). Table 1 summarizes the main functions of morphomap. The morphomap R package is available for all major platforms (Windows, Linux/Unix and Mac OS X).

| Automatic segmentation of long bone meshes
morphomap requires an oriented long bone mesh as the input ( Figure S1). The orientation of the mesh can be checked using the function morphomapCheck which is described in Supporting Information. It may be necessary to mirror long bones to avoid issues of laterality. This can be achieved using the morphomapMirror function (Table 1). The first step of the measurement procedure is the automatic separation of the periosteal and endosteal surfaces ( Figure 2) by applying the CA-LSE R-tool (Computer Assisted Laser Scanner Emulator; Profico, Schlager, et al., 2018). This is carried out by the function morphomapSegm (Table 1).

| Extraction of cross-sections
Sectioning is carried out using the morphomapCore function. First, the user specifies the biomechanical length of the long mesh.   In functions morphomap2Dmap or morphomap3Dmap (Table 1) the cortical thicknesses are used to draw topographic maps in 2D or 3D. In these, colors represent thickness, standardized by default (but selectable) from 0 (minimum) to 1 (maximum), on a 100-step color (the color palette is fully customizable by the user) scale from dark red

| Exporting images of cross-sections
The morphomapPic function (Table 1)  The morphomapRaster function (Table 1) allows the pixel/mm ratio to be adjusted (0.1 by default) and exports the resulting bitmap (raster image) in TIFF format. This is useful to prepare images without a scale bar for importing into other morphometric or reconstruction software.

| Cross-sectional geometry module and its validation
The morphomapCSG function calculates in each section the CSG parameters most commonly used in biomechanical studies: the area, the area moments of inertia (I x , I y ; Figure S4), the principal area moments of inertia (I max and I min ; Figure S5), theta (θ), the polar moment of inertia (J) (Figure S6), and the section moduli (Z x , Z y , Z max , Z min , and Z pol ). The formulae used to calculate each CSG parameter implemented in morphomap are given in Table S1. For the calculation of these parameters morphomap converts the cortical area defined by periosteal and endosteal contours into picture elements of adjustable side length (delta) ( Figure S2). In this way, all the parameters are calculated without interpolation. By default, delta has length 0.1 mm, this is customizable by the user ( Figure S3).
We tested the calculation of CSG parameters embedded in morphomap by extracting two cross-sections from two surfaces artificially created for this purpose; a hollow circular cross-section and a hollow rectangular cross-section (Figure 6a,b). Table 2 shows the comparison between the values obtained via morphomap and those obtained by applying the formulae reported in Supporting Information.
In addition, we compared the values calculated in morphomap with the corresponding figures obtained when using BoneJ (Doube et al., 2010) for diaphyseal cross-sections at 50% of the total biomechanical length of a human femur (OdN 244) and a Pan troglodytes (Kupri 1022) femur ( Figure 6). Table 3 shows the parameters calculated via morphomap and BoneJ.

| Case study-The human and chimpanzee femora
In this section, we show the application of the main functions embedded in the morphomap R functions to a human and a chimpanzee  Ruff (2002). The biomechanical lengths of the bones are 380.23 mm and 277.13 mm for the human and the chimpanzee, respectively. First, the mesh was uploaded into R by using the file2mesh function belonging to the Morpho R package (Schlager, 2017). Second, the morphomapSegm function was used to separate the periosteal and endosteal contours. The morphomapSegm function was used to produce two suitable surfaces for the extraction of the cross-sections. The morphomapCore function extracted sixtyone cross-sections from 20 to 80% of the biomechanical length of the two femora (see Figure 1a,b and Figures S7-S8).
Third, the morphomapShape function was used to extract from each cross-section 24 equiangular semi-landmarks on the periosteal and endosteal outlines (see Figure 1c). Fourth, the morphomap2Dmap function was used to create 2D MMs ( Figure S9) and fifth, the morphomap3Dmap function was used to create 3D MMs (Figure 7).
In agreement with a previous study (Puymerail, 2013) performed on a sample of 12 modern humans and 10 chimpanzees, the human diaphysis shows greater cortical thickness, resisting posterior bending along the linea aspera, while the chimpanzee femur is characterized by medial and lateral bony reinforcements in the proximal portion of the bone. We have reported in Tables S2 and S3 the CSG parameters computed for the human and chimpanzee femora, respectively. These were calculated for 13 cross-sections taken a at 5% increments (for simplicity) from 20 to 80% of the biomechanical length of the femur.
Finally, morphomap offers the prospect of rapid processing of large samples. Table 4 shows the CPU time required to process the chimpanzee femur.

| DISCUSSION
Cortical bone topographic thickness variation in long bones reflects biomechanical loading history (Kivell, Davenport, Hublin, Thackeray, & T A B L E 2 Cross-sectional geometry of the sections in Figure 6a,b Skinner, 2018; Niinimäki et al., 2017;Ruff, Holt, & Trinkaus, 2006;Shaw & Stock, 2009). The R package, morphomap is intended to facilitate the work of functional morphologists in understanding how cortical bone distribution relates to loading history or other factors.
Morphomap is designed to semi-automatically: (a) place semilandmarks on the inner and outer outlines of each section, (b) extract cortical thicknesses for 2D and 3D morphometric mapping, and (c) to compute cross-sectional geometry.
In this article, we detail the algorithms and tools incorporated in morphomap and show that its calculations are valid. The tests performed concerning the use of morphomap to calculate cross-sectional geometric properties are in agreement with both theoretical expectations and the results obtained using the most popular software for such purpose (i.e., BoneJ) (Doube et al., 2010). Further, we compare 2D and 3D MMs from human and chimpanzee femora and show they are in agreement with those from previous work on larger samples (Puymerail, 2013).
Thus, morphomap is a new tool for the study of cortical variation in long bones. It is written in the open source R environment (R Development Core Team, 2016) and is compatible with other R packages commonly used in mesh manipulation, such as Morpho (Schlager, 2017) and Arothron (Profico, Schlager, et al., 2018), morphological studies (Adams & Otárola-Castillo, 2013) and MMs (Adler, Nenadic, & Zucchini, 2003;Sarkar, 2008). morphomap semiautomatically processes the entire bone, beginning with an input mesh that is separated into outer and inner cortical meshes using the CAL-SE method from the Arothron R package . The R programming language allows for easy editing of the R code for alternative implementations and