In mannfred/curvr: Calculate Total Curvature From Landmarked Specimens

knitr::opts_chunk$set(echo=FALSE, warning=FALSE, message=FALSE)

library(curvr)
library(dplyr)

Abstract

Flower-pollinator curvature as viewed from the side (dorsiventral plane) has been a trait of special interest through the post-Darwinian era of pollination ecology. However, within pollination ecology there are at least four metrics in use with few references to their origins or the meaning of the associated units. [@boehm_2020] recently suggested clarifying the definition of curvature by using concepts from differential geometry.
Here, we introduce an R package curvr for computing curvature from lines fitted to landmarked specimens. curvr is designed to take inputs from Momocs, an open-source toolkit for 2D morphometrics in R. We also outline possible advances to the curve-fitting protocol that may subsequently improve estimations of total curvature.
Example dataset.

\newpage

1. Introduction

The curvature of flowers and pollinator mouthparts (e.g. bird bills) is a widespread trait contributing to morphological diversity in the angiosperms and their pollinators. The evolutionary ecology and history of measuring curvature in pollination ecology was recently reviewed in [@boehm_2021]. To briefly summarize; the curvature of floral tubes as viewed from the dorsiventral plane (hereafter 'curvature') has long been recognized as morphological feature that contributes to floral isolation [sensu @grant_1949]. As curvature represents one of three spatial dimensions through which floral morphology can evolve (in addition to width and length), it is likely to play a role in floral diversification [@kay_2009; @armbruster_2009; @vamosi_2018]. Similarly, pollinators that specialize on curved flowers evolve curved mouthparts, contributing to morphological diversification of these clades [@weinstein_2017; @maruyama_2018]. Thus, flower-pollinator curvature has been of special interest to pollination ecologists, resulting in over 40 studies since the late 1970's [@boehm_2021]. Among these studies, Boehm et al. [-@boehm_2021] identified four metrics for defining and quantifying curvature, in addition to the more general approach of quantifying shape via geometric morphometrics. Despite the many approaches to measuring curvature there is a lack of consensus as to what curvature means, how to interpret their units, and how varying methods could or should be compared. To this end, [@boehm_2021] proposed a revised definition of curvature that directly borrowed from the mathematical literature, including an interpretation of accompanying units.

Curvature, as defined in differential geometry, is the rate of change of the tangent, $\mathbf{T}$, with respect to arc length, $s$, along the curve:

$$\kappa = \frac{d \mathbf{T}}{ds}$$

When the angle, $\phi$ of $\mathbf{T}$ is related to an axis of reference (e.g. the $x$ axis), curvature can be expressed as:

$$\kappa = \frac{d\phi}{ds}$$

For a straight line, the direction of the tangent is constant and therefore, curvature is zero. For a curved line, the tangent vector $\mathbf{T}$ changes direction from point to point, resulting in some non-zero value of curvature. Integrating point-wise curvature gives the total curvature [@milnor_1954] of a curve:

$$\kappa_{total} = \int_{0}^{s_{max}} \kappa \; ds$$

At some arbitrary point $i$, point-curvature ($\kappa_i$) is reported as $radians$ per unit length. For example, if measurements have been made at every millimeter along the curve, point-wise curvature is reported as $radians \cdot mm^{-1}$. Total curvature, the sum of total tangent rotations across the length of the curve, is expressed as $radians$. Converting curvature measurements to $degrees$ can be done by multiplying the radian measure by $180\ degrees/\pi$.

The R package curvr was created to allow this definition of curvature to be applied in pollination ecology, and any other system that requires the measurement of 2-D curvature. This package is designed to work alongside common protocols in geometric morphometrics (GM), the study of shape and its covariance with associated or causal variables [@bookstein_1991; @macleod_2002]. While relevant concepts from GM are outlined here, we recommend the concise introduction to the topic by Webster and Sheets [-@webster_2010], and for detailed reading, Zelditch et al [-@zelditch_2012].

Geometric morphometrics aims at quantifying shape, i.e. the geometric information remaining after after specimen size, location, and orientation have been accounted for [@kendall_1977; @bookstein_1991]. A typical GM protocol begins with assigning landmarks to the specimen. Landmarks are points on a specimen that are topologically or biologically homologous [@bookstein_1991]. Generally, the more landmarks assigned to a specimen, the more information about its shape that has been captured. In a comparative study, landmarked specimens are scaled, translated, and rotated to remove features not contributing to shape [@rohlf_1990_b]. In outline-based GM, shape functions are fit to the landmarks, for example using the tangent angle function [@zahn_1972]. In landmark-based GM, the set of landmarks are taken together as a 'landmark configuration' and exist in a multidimensional shape space, defined by the number of landmarks used, and the spatial dimensions in which they were assigned [2D or 3D, @kendall_1984]. In both approaches, shape has been summarised as a single variable so that useful statistical procedures can be applied; most GM protocols involve principal components analysis (PCA) of shape variables, enabling the major axes of shape variation to be calculated and plotted in PCA 'shape space' [@macleod_2002].

Inspection of deformation grids [@bookstein_1991] or end-member specimens [@macleod_2002] plotted onto these axes are suggestive of the geometric features contributing to shape variation (e.g. curvature). The role of curvr is to re-purpose the landmark assignments used in the GM analysis to calculate total curvature from each specimen. By explicitly quantifying curvature, the investigator can explicitly test the contribution of curvature to shape variation along some principal axis of interest. In the following section the algorithms for computing total curvature are described, followed by a demonstration of its utility alongside a typical GM analysis.

\newpage

2. The `curvr` package

Computing total curvature necessitates fitting a curve to the specimen of interest. Fortunately, the landmarks used in a GM protocol can be re-purposed for regressing curves to the specimen. However, in GM, landmarks are placed throughout the entire specimen to quantify shape. Quantifying curvature typically implies a segment of the shape that is of particular interest [but see: @silk_1978], therefore the landmarks used for curvature analysis will be a subset of those used in the GM protocol.

The curvr package can be used alongside a typical GM protocol where specimens have been digitally landmarked, e.g. using tpsDig [@rohlf_2015], Momocs [@bonhomme_2014], or geomorph [@adams_2013]. Once a set of specimens have been landmarked, the corresponding files can be imported into R using e.g. geomorph::readland.tps(). For open curves, we recommend rotating so that the longest axis of the shape is parallel to the x-axis, so that single-valued functions (SVF) can be fit &emdash; this can be done with Momocs::coo_alignxax(). Total curvature can then be computed using curvr::curvature_spline(). This function takes three arguments: the first is the matrix object containing the landmark data, the second is the upper and lower bounds of the curve from which to estimate curvature, and the third is the type of spline to fit to the landmarks.

Since GM practitioners will be well-versed in handling landmark matricies, we will focus on the latter two arguments. The x_range argument will, in most cases, be a numeric vector of the first and last x-axis landmark of the specimen. If handling multiple specimens, we recommend extracting these boundaries and storing in a list. Deciding on the type argument will vary between systems. Two options are available: smooth (smoothing) and ip (interpolating). Perperoglou et al. [-@perperoglou_2019] recently reviewed spline fitting procedures in R, and we recommend this primer to assess which is appropriate for the study system at hand.

The curvature_spline() function first takes the landmark matrix and fits a spline via stats::smooth.spline() or stats::spline(). From these splines, the first and second derivatives are computed and stored as functions using stats::splinefun() with the appropriate deriv argument.

A function for point-wise curvature $\kappa_i$ is then constructed as:

$$\kappa_i = \frac{\lvert{y_i''}\rvert} {(1 + y_i'^2)^{3/2}}$$ Which is curvature at $x=i$ for a single-valued function $y=f(x)$ [@huber_2008; @rutter_2000, pp. 133-139]. The arc-length segment $ds_i$ is simply:

$$ds_i = \sqrt{1^2 + {y'_i}^2}$$

These two equations are assigned to functions using the Hessian and gradient functions constructed above, and passed to stats::integrate() to calculate total curvature:

$$\kappa_{total} = \int_{i=0}^{i=s_{max}} \frac{\lvert{y_i''}\rvert} {(1 + y_i'^2)^{3/2}} \ \ \sqrt{1^2 + {y'i}^2}= \int{0}^{s_{max}} \kappa \; ds$$

3. Applications

In this section, we illustrate the utility of quantifying curvature along a typical analysis of morphometric data. Here, we test whether bill curvature is correlated with variation along a principal component of shape space. For this demonstration, landmark was kindly provided by Navalon et al [-@navalon_2020] — their study tested whether skull and bill morphology are integrated or decoupled phenotypes. One of the results from Navalon et al (2020, p.271) states:

"We found that each of the major clades of landbirds diverged to unique cranial morphologies. Parrots (Psittaciformes) are characterized by a single ancestral shift towards very high rates of skull shape evolution, resulting in a characteristic cranial anatomy with short, curved beaks and expanded braincases."

short: https://besjournals.onlinelibrary.wiley.com/doi/epdf/10.1111/2041-210X.13385 https://besjournals.onlinelibrary.wiley.com/doi/epdf/10.1111/2041-210X.13293

pgs 184, 186o https://ebookcentral.proquest.com/lib/ubc/reader.action?docID=422384&ppg=18

Before point-wise curvature is estimated from the derivatives, the polynomial function is parameterized by arc length. The advantage of this step is that curvature can be calculated at even increments of arc length, $s$. The polynomial is parameterized by arc length by first splitting the curve into sub-curves of incrementally increasing proportions of total length, from 0 to 1. For example, a curve split into 10 parts would have 10 sub-curves with arc length calculated as $0.1(s)$, $0.2(s)$, ..., $1.0(s)$. Each sub-curve has a lower $x$ bound of known value (the lower bound of the original polynomial). The upper $x$ bound (of each sub-surve) is unknown. The root-finding function stats::uniroot solves for the unknown upper $x_i$ value, $u_i$, for each $i$ sub-curve:

f(u) = pracma::arclength(polynomial, lower x boundary, u)$length

The result is a vector of $x$ values that correspond to even increments of arc length along the polynomial. We then use the gradient (first derivatives, $y'$) and Hessian matricies (second derivatives, $y''$) to calculate point-wise $\kappa_i$ at every position $x_i$:

$$\kappa = \frac{|y''|}{(1+y'^2)^{\frac{3}{2}}}$$