ShapeRotator: An R tool for standardized rigid rotations of articulated three-dimensional structures with application for geometric morphometrics.

The quantification of complex morphological patterns typically involves comprehensive shape and size analyses, usually obtained by gathering morphological data from all the structures that capture the phenotypic diversity of an organism or object. Articulated structures are a critical component of overall phenotypic diversity, but data gathered from these structures are difficult to incorporate into modern analyses because of the complexities associated with jointly quantifying 3D shape in multiple structures. While there are existing methods for analyzing shape variation in articulated structures in two-dimensional (2D) space, these methods do not work in 3D, a rapidly growing area of capability and research. Here, we describe a simple geometric rigid rotation approach that removes the effect of random translation and rotation, enabling the morphological analysis of 3D articulated structures. Our method is based on Cartesian coordinates in 3D space, so it can be applied to any morphometric problem that also uses 3D coordinates (e.g., spherical harmonics). We demonstrate the method by applying it to a landmark-based dataset for analyzing shape variation using geometric morphometrics. We have developed an R tool (ShapeRotator) so that the method can be easily implemented in the commonly used R package geomorph and MorphoJ software. This method will be a valuable tool for 3D morphological analyses in articulated structures by allowing an exhaustive examination of shape and size diversity.

BACKGROUND

Data on shape and size variation are essential in many fields, including evolutionary biology and ecology, engineering, medical science, and anthropology (Loncaric, 1998; McIntyre & Mossey, 2003; Slice, 2006). For most of these studies, the most widely used tools for analyzing morphological variation within or between a group of organisms or objects are based on Cartesian coordinates of landmarks (Bookstein, 1997).

Of the wide array of methods using Cartesian coordinates, geometric morphometrics (GM) is the most common, especially when analyzing shape and size variation and covariation (Adams, Rohlf, & Slice, 2013; Mitteroecker & Gunz, 2009). The first two steps of this GM procedure consist of a landmark approach that (1) gathers (two‐ or three‐dimensional) coordinates of anatomically defined and homologous loci, followed by (2) a generalized Procrustes analysis (GPA) that superimposes configurations of each set of landmarks in all specimens, by removing all effects of size, translation, and rotation, in order to only obtain shape information (Adams et al., 2013; Klingenberg, 2008). Geometric morphometrics, therefore, allows accurate quantitative analyses of shape and size, in either two‐dimensional (2D) or three‐dimensional (3D) space.

3D morphological analyses are the most accurate, as objects and organisms exist in 3D space. The recent growth in X‐ray micro‐CT scanning and surface scanning has seen a rapid increase in the application of 3D geometric morphometric techniques, but progress has been hampered by the lack of a simple method to incorporate data from complex articulated structures.

In evolutionary biology, identifying morphological differences among different groups or taxa is crucial in order to understand evolutionary processes and their relationship to the environment (Losos, 1990; Pagel, 1999; Ricklefs & Miles, 1994). This can be difficult, especially, if traits have coevolved, or if morphological diversification has been hindered by phylogenetic legacy or trade‐offs imposed by the organism's functional habitat (Ghalambor, McKay, Carroll, & Reznick, 2007). Complex body shape patterns require more detailed analyses of shape, obtained by collecting data from several structures that capture the whole gamut of morphological variation in an organism. One example of this is the extraction and assembly of data from articulated structures, such as skeletons, for 3D analyses with geometric morphometric techniques. This is especially important in functional morphological studies, as they usually involve analyzing more than one structure due to mechanical correlations or morphological integration. For example, jointly analyzing the skull and mandible could be crucial to disentangle the relationship between diet and head shape evolution (Cornette, Baylac, Souter, & Herrel, 2013). Similarly, collectively evaluating different modules in the limbs, especially when correlated to locomotion, or considering several structures across the whole body, could improve our understanding of the effect of environmental conditions on morphological evolution (Vidal‐García & Keogh, 2017).

Unfortunately, nonrigid structures, such as articulated structures, will inevitably suffer the effects of natural or free rotation or translation events and be different in each individual and structure (Adams, 1999). These events could obstruct the correct quantification of shape variation by adding rotation artifacts to GM analyses (Adams, Rohlf, & Slice, 2004). Thus, orientation of these structures needs to be corrected and standardized prior to performing shape analyses. Several approaches for shape analysis of landmark data in articulated structures already have previously been proposed (Adams, 1999; Bookstein, 1991). However, even though these approaches can be used in 2D and also in 3D data with further modifications, currently, it has only been implemented in the two‐dimensional space, such as in geomorph's function fixed.angle()(Adams et al., 2017).

Here, we present the R tool ShapeRotator: a simple geometric rigid rotation approach to study three‐dimensional (3D) shape of articulated structures, or independent structures, within an organism. We describe a method that removes shape variation due to the effect of translation between independent structures and rotation generated by movement in an articulation, among others. Thus, our approach translates and rotates articulated (or even independent) structures in order to obtain a comparable dataset with all effects of random movement and rotations removed.

We apply this method to two landmark‐based datasets for analyzing shape variation using geometric morphometrics: (A) two substructures in a single‐point articulation and (B) two substructures in a double‐point articulation (Figures 1 and 2a,b). We also provide the example dataset used in ShapeRotator (available in GitHub) to execute these two kinds of rigid rotations, which then allows geometric morphometric analyses to be performed in the two most commonly used 3D GM analytical software packages: geomorph (Adams et al., 2017) and MorphoJ (Klingenberg, 2011). This method also will allow exporting of the rotated coordinates for posterior analyses in other software platforms, even outside of the field of geometric morphometrics. As the basis of this method lies upon rigidly spinning any structure defined by 3D coordinates, it could be used in any other shape analyses that use coordinate data, such as continuous surface meshes used in spherical harmonics (Shen, Farid, & McPeek, 2009). Our method is a convenient addition to the rapidly evolving tool kit of geometric morphometrics because it allows a more comprehensive exploration of morphological diversity through the gathering of shape data from complex 3D structures.

Figure 1

Two structures in a simple‐point articulation and a double‐point articulation being rotated using simple.rotation() and double.rotation(), respectively. The structures depicted are a frog humerus and radioulna for the single‐point articulation and a skull and a detached mandible for the double‐point articulation

Figure 2

(a) Application of a translation ‐translate()‐ and the 3D rigid rotation method simple.rotation() in a simple‐point articulation (e.g., humerus and radioulna) for three different “specimens,” by rotating articulated structures to a standardized position relative to each other. (b) Application of the translation ‐translate()‐ and the 3D rigid rotation method double.rotation() in a double‐point articulation for three different “specimens” (e.g., skull and mandibles). (c) Rotation method exemplified for simple.rotation() by depicting the plane spanned by the already translated point p 0 and A. Please note that p 0 depicts the origin point (0, 0, 0). The rotated resulting point p , landmarks B, C, D, E, and F, and angle θT (desired angle between the two structures) are also depicted. (d) Rotation method exemplified for double.rotation(). Please note that as both structures are not attached for the double‐point articulation example, they are both translated to p 0, which depicts the origin point (0, 0, 0). The rotated resulting point p , landmarks A, B, C, D, E, F, G, and H, and angles θ1, θ2, and θT (desired angle between the two structures) are also depicted

METHODOLOGY

We begin with a set of points . which represents a 3D object and are ordered so that represents the base point and represents the end point, by which we mean that this object has an axis starting from and ending at . Our goal was to rotate these points via a rigid motion so that the axis on which these two points sit is either on the X, Y, or Z‐axis in . Rotation of vectors in is a well‐known and easily resolved problem, and various formalisms exist in geometry. Thus, we translate our set of points so that maps to the origin (0, 0, 0). This is a simple transformation T defined by:

Note that the axis X = span {(1, 0, 0)}, Y = span {(0, 1, 0)}, Z = span {(0, 0, 1)}, where each of the generating vectors is unit. Let us fix our desired axis to which we rotate the object to be A = span {a} where a = (1, 0, 0), or a = (0, 1, 0), or a = (0, 0, 1). As we have translated points {p } and vectors correspond to positions, we are simply looking to rotate the vector p to A, and each other vector as a rigid motion with respect to this rotation. There are a number of ways to do this, but the simplest way is to consider the plane spanned by p and A and then to rotate by the angle between p and A within this plane (Figure 2c,d). Such a rotation is performed via rotating on the axis to the plane, which is determined by a normal vector to this plane.

Let us describe this setup slightly more generally. For two vectors , the axis to the plane spanned by these two vectors is determined by a unit normal to the plane (there are two choices due to orientation), which we denote by N (u, v):where × is the cross product. The angle between these vectors is then given by :where · is the dot (scalar) product between vectors. The rotation matrix about an axis , where w = (w 1, w 2, w 3) is a unit vector, of angle θ radians is given by the well‐known matrix:

Thus, to obtain a rotation matrix which is the rigid motion rotating the vector u to v in the plane spanned by u and v, we obtain the expression:

Getting back to our original problem, we set v = p and u = a, and then we have the rotated points:where R (p , a) p is the action of the matrix R (p , a) on the vector p .

It may be necessary to introduce a further constraint in the rotation. For instance, suppose a =  (0, 1, 0) and there is a point p , now rotated to r via the method we describe, which should lie in the Y‐axis. That is, we need to further rotate r to a point . To do this, we simply rotate in the axis a, by an angle θ(r ) = arctan ((r )3/(r )1), where r  = ((r )1, (r )2, (r )3). That is,

IMPLEMENTATION

Overview of the ShapeRotator package

Here, we illustrate the functions available within the ShapeRotator R tool and the basic steps required in order to successfully implement the rotation on a dataset of 3D coordinates. ShapeRotator allows the rigid rotation of sets of both landmarks and semilandmarks used in geometric morphometric analyses, enabling morphometric analyses of complex objects, articulated structures, or multiple parts within an object or specimen. This tutorial uses two example datasets: (A) two neighboring bones of the arm in a frog (humerus and radioulna), representing a single‐point articulation example and (B) a skull and a mandible, representing a two‐point articulation example (Figure 1). The main steps required are as follows: (1) importing the data and fixating the rotation axes, (2) translating the whole dataset of coordinates or points so that the main selected point , and (3) rotating the two structures to the desired angle (as outlined on Figure 2a,b).

Importing a dataset

In the first example dataset (example A), we use two geometric morphometric datasets containing both landmarks and semilandmarks for two neighboring and articulated bones (humerus and radioulna) from a group of several species of frogs (details in Appendix S1), in tps format. We first import the datasets using the R package geomorph (Adams et al., 2017): library(devtools) install_github(“marta‐vidalgarcia/ShapeRotator”) library(geomorph) radioulna <‐ readland.tps(“radioulna.tps”, specID = “ID”, readcurves = F) humerus <‐ readland.tps(“humerus.tps”, specID = “ID”, readcurves = F)

These two GM datasets (radioulna and humerus) will be rotated on different rotation axes in order to conform the aimed angle between them. This process is not exclusive to two neighboring structures, and thus, it could be performed for as many independent subunits as desired by choosing the different angles among different rotation axes and all the translation processes. For more help on importing the GM datasets, please see Adams et al., (2017) and the associated help files. Please note that this method also works for semilandmarks as long as they have been equidistantly positioned prior to the translation and the rigid rotation. The same process is needed for the second example dataset (example B).

Translating

During this step, each structure will be translated to the point of origin so that ; thus, the distance from the coordinates of landmark_a (A , A , A ) is substracted from all the landmarks in all specimens, for example (N  − A , N  – A , N  – A ) for landmark N. This translation is made with the function translate(), as it follows for the one‐point articulation example (example A): data.1_t <‐ translate (T=data_1, landmark=landmarkA) data.2_t <‐ translate (T=data_2, landmark=landmarkD)

And for the double‐point articulation (example B): data.3_t <‐ translate (T=data_3, landmark=landmarkA) data.4_t <‐ translate (T=data_4, landmark=landmarkE)

Please note that the default is to set the origin point to, but this can be changed to another origin point. For example: data.1_t <‐ translate (T=data_1, landmark=landmarkA, origin = c(1,3,5)) data.2_t <‐ translate (T=data_2, landmark=landmarkD, origin = c(1,3,5))

Fixing the rotation axis

In order to fix a rotation axis in a structure, we first need to select in our dataset two suitable landmarks for each structure which the rotation axis will go through landmarks A and B (for the first structure), and landmarks D and E (for the other structure). In the radioulna dataset from example A, landmark A is the landmark in the 1st position and landmark B is in the 10th position. Similarly, for the humerus dataset, landmark D would be the landmark on the 52nd position, and landmark E would be in the 19th position. Please note that in example A, we would perform a simple rotation between two structures, based on the angle formed by three points (one on each structure and the articulation point itself. However, as we aim to remove all articulation‐related variation, including torsional rotations and mirroring artefacts, we need an extra landmark in each structure that shares at least one of their coordinates (but ideally two) with the distal landmark (e.g., y A = y B for the first structure, and z D = z E for the second). This extra landmark is needed for the simple reason that there is not enough information about the orientation of the structure with only two landmarks per structure, so even though the rigid rotation will work properly, it could position this structure in the wrong “mirroring” orientation (Figure 2c). This orientation issue is corrected by adding landmark C. In this example, landmark C is the 17th landmark, while landmark F is the 107th landmark in the humerus dataset. We need to know which landmarks will be selected in both structures prior to the rotation process in order to insure that the rotation process will work properly. Setting up landmarks that will be used in the rotation process of example B is very similar to example A, but it needs four landmarks for each structure instead of three, as it is a more complex articulation that requires two separate rotations per structure (Figure 2d). Even though each of these rotations are calculated internally (only the four landmarks per structure and the desired angle between them need to be provided), it will be beneficial to choose landmarks that are spatially arranged in a way that facilitates the rotation process, and results in the rotating multistructure being placed in a biologically relevant angle between each substructure. ShapeRotator will give a warning message if the landmarks chosen are not optimal.

Rotating

In the rotation step, we will use the function simple.rotation() in order to rigidly rotate the two structures of example A to the desired angle (in degrees), as it follows: joined_dataset <‐ simple.rotation(data.1 = data.1_t, data.2 = data.2_t, land.a = 10, land.b=1, land.c=17, land.d=52, land.e=19, land.f=107, angle = 90)

The input datasets data.1 and data.2 correspond to the two translated datasets (in this case translated_radioulna and translated_humerus. We then use the selected landmarks as explained in the previous section. Finally, we include the angle (in degrees) that we would like to use to position the two structures to one another. One of the options of the function simple.rotation() is to select the desired angle between the two structures so that we can perform the rigid rotation of each structure positioning them in the selected angle in relation to each other. In order to do so, we use the internal function vector.angle(), by providing the desired angle in degrees (from 0° to 360°). The function vector.angle() will return angle.v, a vector that forms that angle with the vector (1, 0, 0). In the example dataset in ShapeRotator, we rotate the coordinates from the two bones so that they form an angle of 90° degrees to each other: New_vector <‐ vector.angle(90)

So that New_vector = c(0, 1, 0). Thus, we could check the vector that the function simple.rotation() will use, based on the input angle. The output from the function simple.rotation() is a 3D array with the two joined datasets (data.1 and data.2). Please note that the two datasets are joined based on their associated specimen names (dimnames). Thus, the order of the specimens in each dataset is not important, as long as all the cases match perfectly between the two datasets in the same specimen. If there are extra specimens for one of the datasets or the names do not match properly, simple.rotation() will not include them in the output rotated joined dataset.

As explained in the previous section, the function used for the rotation step for example B needs four landmarks for each structure instead of three, as it is a more complex articulation that will need to get landmarks from both structures joined (or very close) to form the double‐articulation hinges, as well as internally running separate rotations per structure on top of the final rotation. We will use double.rotation() function. For this function, we will need to match the dimnames: matched = match.datasets(data.3_t, data.4_t)

The datasets that will need to be rotated are the following: data.3_rotate = matched$matched1 data.4_rotate = matched$matched2

Similarly to simple.rotation(), the double.rotation() will be as follows: dr45l = double.rotation(data.1= data.3_rotate, data.2= data.4_rotate, land.a=1, land.b=2, land.c=3, land.d=4, land.e=1, land.f=2, land.g=3, land.h=4, angle=45) joined_dataset <‐ double.rotation(data.1 = translated_skull, data.2 = translated_mandible, land.a = 1, land.b=2, land.c=3, land.d=4, land.e=1, land.f=2, land.g=3, land.h=4, angle = 45)

The function double.rotation() will return a list that will need to be joined as it follows: dr45 = join.arrays (dr45l$rotated1, dr45l$rotated2)

Double.rotation() allows the user to translate one of the objects after the rotation, in the case of not wanting them in contact to one another. For example skull_translate = c(1.7,0.1, 0) dr45_st = join.arrays (dr45l$rotated1, translate(dr45l$rotated2,land.e, skull_translate))

Plotting

Finally, after the rotation process, we can visualize a 3D plot for each specimen with the two rotated structures for both a single‐point and a double‐point articulations using the function plot.rotation.3D(): plot.rotation.3D(joined.data= joined_dataset, data.1=data.1, data.2=data.2, specimen.num =1) plot.rotation.3D(joined.data=dr45_st, data.1=data.3, data.2=data.4, specimen.num =1)

The default plotting colors for the landmarks of the joined rotated structure (joined.data) is black and red for each original structure (data.1 and data.2).

Exporting

After the rotation process, we could either use the joined GM array in further analyses, visualize the resulting joined structure through geomorph (Adams et al., 2017), or we could also export it and save it in order to use it in another software, such as MorphoJ (Klingenberg, 2011). In this step, we will be using the function writeland.tps() in the R package geomorph (Adams et al., 2017) in order to save a tps file from the joined GM array: writeland.tps(A=“joined_arm”, file = “joined_arm.tps“, scale = NULL)

Other applications

Our method is an important addition to the tool kit of the geometric morphometrics field. It will facilitate the analyses of compound 3D morphological datasets in geometric morphometrics analyses but will also be useful outside of this field as it can be applied to any method that uses 3D coordinates. The examples of applications are numerous in different fields of study, such as biology, anthropology, paleontology, medical sciences, archeology, and engineering. For example, in evolutionary biology, ShapeRotator would allow analyses of multiple or articulated hard structures (such as different segments of an exoskeleton, different articulated bones, or neighboring plant structures, among others), different structures from the same object or organism (e.g., different and not adjacent body parts), or pieces from damaged specimens. In medicine and veterinary science, it could be used to examine shape and size variation in different organisms’ growth due to different nutritional treatments or to examine how different structures respond to injuries or surgery. It would be useful in paleontology or archeology when trying to quantify shape of different objects or organisms that might have been preserved in disarticulated pieces. Finally, we would like to include a word of caution on how the angle chosen between different substructures could produce different results. As common sense would suggest, the angle chosen should always be within the range of angles possible for that complex object, but in some cases, it might be necessary to perform the GM analysis on more than one angle to be sure that the results do not differ. In other cases, it might even be better to analyse substructures separately.

CONFLICT OF INTEREST

None declared.

Click here for additional data file.