, Silvester Bartsch
[aut], Nicole Grunstra [aut], Philipp Mitteroecker [aut]| Title: | Warping Landmark Configurations |
|---|---|
| Description: | Compute bending energies, principal warps, partial warp scores, and the non-affine component of shape variation for 2D landmark configurations, as well as Mardia-Dryden distributions and self-similar distributions of landmarks, as described in Mitteroecker et al. (2020) <doi:10.1093/sysbio/syaa007>. Working examples to decompose shape variation into small-scale and large-scale components, and to decompose the total shape variation into outline and residual shape components are provided. Two landmark datasets are provided, that quantify skull morphology in humans and papionin primates, respectively from Mitteroecker et al. (2020) <doi:10.5061/dryad.j6q573n8s> and Grunstra et al. (2020) <doi:10.5061/dryad.zkh189373>. |
| Authors: | Anne Le Maitre [aut, cre]
, Silvester Bartsch
[aut], Nicole Grunstra [aut], Philipp Mitteroecker [aut] |
| Maintainer: | Anne Le Maitre <[email protected]> |
| License: | GPL-3 |
| Version: | 1.0.1 |
| Built: | 2024-12-16 07:02:33 UTC |
| Source: | CRAN |
Convert a three-dimensional array of landmark coordinates into a two-dimensional matrix
array.to.xxyy(A)array.to.xxyy(A)
A |
A 3D array (p x k x n) containing landmark coordinates for a set of specimens |
Function returns a two-dimensional matrix of dimension (n x [p x k]), where rows represent specimens and columns represent variables. The p first columns correspond to X coordinates, etc.
A <- array(rnorm(40), c(5, 2, 4)) # random 2D coordinates of 5 landmarks for 4 specimens array.to.xxyy(A)A <- array(rnorm(40), c(5, 2, 4)) # random 2D coordinates of 5 landmarks for 4 specimens array.to.xxyy(A)
Computes the principal warps and the bending energy of a reference shape configuration, as well as the variance of the partial waprs, the partial warp scores and the non-affine component of shape variation for 2D landmark coordinates (3D not implemented). Small-scale and large-scale components of shape variation can also be computed.
create.pw.be(A, M_ref, d = NULL)create.pw.be(A, M_ref, d = NULL)
A |
a k x 2 x n array, where k is the number of 2D landmarks, and n is the sample size. |
M_ref |
a k x 2 refence matrix (usually the sample mean shape), where k is the number of 2D landmarks |
d |
(optional) an integer value comprised between 1 and (k-3) to compute small-scale shape components (between 1 and d) and large-scale shape components (between d+1 and k-3) |
A list containing the following named components:
bendingEnergy |
the bending energy (the (k-3) eigenvalues of the bending energy matrix) |
principalWarps |
the k x (k-3) matrix of principal warps (the k eigenvectors of the bending energy matrix) |
partialWarpScores |
the n x (2k-6) matrix of partial warp (the projection of the vectors of shape coordinates, expressed as deviations from the reference shape, onto the principal warps) |
variancePW |
the variance of the (k-3) partial warps |
Xnonaf |
the n x 2k matrix of the non-affine component of shape variation |
Xsmall |
the n x 2k matrix of the small-scale shape variation (if d is provided) |
Xlarge |
the n x 2k matrix of the large-scale shape variation (if d is provided) |
Bookstein FL. (1989). Principal Warps: Thin-plate splines and the decomposition of deformations. IEEE Transactions on pattern analysis and machine intelligence 11(6): 567–585. https://ieeexplore.ieee.org/abstract/document/24792
See CreateL for the creation of the bending energy matrix
# 2D landmark coordinates library("geomorph") data("HomoMidSag") # dataset n_spec <- dim(HomoMidSag)[1] # number of specimens k <- dim(HomoMidSag)[2] / 2 # number of landmarks homo_ar <- arrayspecs(HomoMidSag, k, 2) # create an array # Procrustes registration homo_gpa <- Morpho::procSym(homo_ar) m_overall <- homo_gpa$rotated # Procrustes coordinates m_mshape <- homo_gpa$mshape # average shape # Computation of bending energy, partial warp scores etc. homo_be_pw <- create.pw.be(m_overall, m_mshape) # Partial warp variance as a function of bending energy logInvBE <- log((homo_be_pw$bendingEnergy)^(-1)) # inverse log bending energy logPWvar <- log(homo_be_pw$variancePW) # log variance of partial warps mod <- lm(logPWvar ~ logInvBE) # linear regression # Plot log PW variance on log BE^-1 with regression line plot(logInvBE, logPWvar, col = "white", asp = 1, main = "PW variance against inverse BE", xlab = "log 1/BE", ylab = "log PW variance") text(logInvBE, logPWvar, labels = names(logPWvar), cex = 0.5) abline(mod, col = "blue")# 2D landmark coordinates library("geomorph") data("HomoMidSag") # dataset n_spec <- dim(HomoMidSag)[1] # number of specimens k <- dim(HomoMidSag)[2] / 2 # number of landmarks homo_ar <- arrayspecs(HomoMidSag, k, 2) # create an array # Procrustes registration homo_gpa <- Morpho::procSym(homo_ar) m_overall <- homo_gpa$rotated # Procrustes coordinates m_mshape <- homo_gpa$mshape # average shape # Computation of bending energy, partial warp scores etc. homo_be_pw <- create.pw.be(m_overall, m_mshape) # Partial warp variance as a function of bending energy logInvBE <- log((homo_be_pw$bendingEnergy)^(-1)) # inverse log bending energy logPWvar <- log(homo_be_pw$variancePW) # log variance of partial warps mod <- lm(logPWvar ~ logInvBE) # linear regression # Plot log PW variance on log BE^-1 with regression line plot(logInvBE, logPWvar, col = "white", asp = 1, main = "PW variance against inverse BE", xlab = "log 1/BE", ylab = "log PW variance") text(logInvBE, logPWvar, labels = names(logPWvar), cex = 0.5) abline(mod, col = "blue")
2D Cartesian coordinates of 87 landmarks quantifying the skull morphology along the midsagittal plane for 24 adult modern humans.
data(HomoMidSag)data(HomoMidSag)
A data frame with 24 rows and 184 variables
Bartsch, Silvester (2019) The ontogeny of hominid cranial form: A geometric morphometric analysis of coordinated and compensatory processes. Master's thesis, University of Vienna.
Mitteroecker, Philipp et al. (2020) Morphometric variation at different spatial scales: coordination and compensation in the emergence of organismal form. Systematic Biology, 69(5): 913–926. doi: 10.1093/sysbio/syaa007
Mitteroecker, Philipp et al. (2020) Data form: Morphometric variation at different spatial scales: coordination and compensation in the emergence of organismal form. Dryad Digital Repository. doi: 10.5061/dryad.j6q573n8s
Create a matrix of 2D shape coordinates drawn from a Mardia-Dryden distribution (3D not implemented)
md.distri(M_ref, n, sd = 0.05)md.distri(M_ref, n, sd = 0.05)
M_ref |
a k x 2 refence matrix (usually the sample mean shape), where k is the number of 2D landmarks |
n |
the number of observations |
sd |
the standard deviation of the distribution (default = 0.02) |
the n x 2k matrix of shape coordinates drawn from a Mardia-Dryden distribution
# 2D landmark coordinates library("geomorph") data("HomoMidSag") # dataset n_spec <- dim(HomoMidSag)[1] # number of specimens k <- dim(HomoMidSag)[2] / 2 # number of landmarks homo_ar <- arrayspecs(HomoMidSag, k, 2) # create an array # Procrustes registration homo_gpa <- Morpho::procSym(homo_ar) m_mshape <- homo_gpa$mshape # average shape # Mardia-Dryden distribution Xmd <- md.distri(m_mshape, n = n_spec, sd = 0.005) # Visualization plot(Xmd[, 1:k], Xmd[, (k+1):(2*k)], asp = 1, las = 1, cex = 0.5, main = "Mardia-Dryden distribution", xlab = "X", ylab = "Y")# 2D landmark coordinates library("geomorph") data("HomoMidSag") # dataset n_spec <- dim(HomoMidSag)[1] # number of specimens k <- dim(HomoMidSag)[2] / 2 # number of landmarks homo_ar <- arrayspecs(HomoMidSag, k, 2) # create an array # Procrustes registration homo_gpa <- Morpho::procSym(homo_ar) m_mshape <- homo_gpa$mshape # average shape # Mardia-Dryden distribution Xmd <- md.distri(m_mshape, n = n_spec, sd = 0.005) # Visualization plot(Xmd[, 1:k], Xmd[, (k+1):(2*k)], asp = 1, las = 1, cex = 0.5, main = "Mardia-Dryden distribution", xlab = "X", ylab = "Y")
2D Cartesian coordinates of 70 landmarks quantifying the skull morphology along the midsagittal plane for 67 adult modern primates (mostly papionins). The data correspond to a list with the 6 following elements:
data(papionin)data(papionin)
A list of 6 elements.
coords The 3D array of landmark coordinates
species The vector of species names
semi_lm The vector of semilandmark numbers of the full dataset
curves The list of curves for sliding semilandmarks of the full dataset
links The matrix of links between landmarks for the full dataset
outline A list of 4 elements for the analysis of the outline shape: subset, the landmark numbers for the susbet; semi_lm, the vector of semilandmark numbers; curves, the list of curves for sliding semilandmarks; links, the matrix of links between landmarks for the full dataset.
Grunstra, Nicole D. S. et al. (2021) Detecting phylogenetic signal and adaptation in papionin cranial shape by decomposing variation at different spatial scales. Systematic Biology, 70(4): 694–706. doi: 10.1093/sysbio/syaa093
Grunstra, Nicole D. S. et al. (2020) Data form: Detecting phylogenetic signal and adaptation in papionin cranial shape by decomposing variation at different spatial scales. Dryad Digital Repository. doi: 10.5061/dryad.zkh189373
Create a matrix of 2D shape coordinates drawn from a self-similar distribution (3D not implemented)
ssim.distri(M_ref, n, sd = 0.02, f = 1)ssim.distri(M_ref, n, sd = 0.02, f = 1)
M_ref |
a k x 2 refence matrix (usually the sample mean shape), where k is the number of 2D landmarks |
n |
the number of observations |
sd |
the standard deviation of the distribution (default = 0.02) |
f |
a scaling factor |
the n x 2k matrix of shape coordinates drawn from a self-similar distribution
# 2D landmark coordinates library("geomorph") data("HomoMidSag") # dataset n_spec <- dim(HomoMidSag)[1] # number of specimens k <- dim(HomoMidSag)[2] / 2 # number of landmarks homo_ar <- arrayspecs(HomoMidSag, k, 2) # create an array # Procrustes registration homo_gpa <- Morpho::procSym(homo_ar) m_mshape <- homo_gpa$mshape # average shape # Self-similar distribution Xdefl <- ssim.distri(m_mshape, n = n_spec, sd = 0.05, f = 1) # Visualization plot(Xdefl[, 1:k], Xdefl[, (k+1):(2*k)], asp = 1, las = 1, cex = 0.5, main = "Self-similar distribution", xlab = "X", ylab = "Y")# 2D landmark coordinates library("geomorph") data("HomoMidSag") # dataset n_spec <- dim(HomoMidSag)[1] # number of specimens k <- dim(HomoMidSag)[2] / 2 # number of landmarks homo_ar <- arrayspecs(HomoMidSag, k, 2) # create an array # Procrustes registration homo_gpa <- Morpho::procSym(homo_ar) m_mshape <- homo_gpa$mshape # average shape # Self-similar distribution Xdefl <- ssim.distri(m_mshape, n = n_spec, sd = 0.05, f = 1) # Visualization plot(Xdefl[, 1:k], Xdefl[, (k+1):(2*k)], asp = 1, las = 1, cex = 0.5, main = "Self-similar distribution", xlab = "X", ylab = "Y")
Maps landmarks via thin plate spline based on a reference and a target configuration in 2D or in 3D. This function is an extension of the tps3d function for a set of specimens.
tps.all(X_array, REF_array, TAR_matrix)tps.all(X_array, REF_array, TAR_matrix)
X_array |
original coordinates - a 3D array (p x k x n) containing original landmark coordinates for a set of specimens |
REF_array |
reference coordinates (e.g., outline landmarks for all specimens) - a 3D array (p x k x n) containing reference landmark coordinates for a set of specimens |
TAR_matrix |
target coordinates (e.g., average outline landmarks) - a matrix (n x k) containing target landmark coordinates |
p is the number of landmark points, k is the number of landmark dimensions (2 or 3), and n is the number of specimens.
Function returns a 3D array (p x k x n) containing the deformed input (original landmark set warped onto the target matrix).
Bookstein FL. (1989). Principal Warps: Thin-plate splines and the decomposition of deformations. IEEE Transactions on pattern analysis and machine intelligence, 11(6): 567–585. https://ieeexplore.ieee.org/abstract/document/24792
See tps3d
data("papionin") # load dataset # Full dataset: 70 landmarks papionin_ar <- papionin$coords # Outline dataset: subset of 54 landmarks outline_ar <- papionin_ar[papionin$outline$subset, , ] # Subset: Macaca only mac <- grep("Macaca", papionin$species) # genus Macaca papionin_macaca <- papionin_ar[, , mac] outline_macaca <- outline_ar[, , mac] # Landmark sliding by minimizing bending energy + superimposition (GPA) library("Morpho") papionin_gpa <- procSym(papionin_macaca, SMvector = papionin$semi_lm, outlines = papionin$curves) outline_gpa <- procSym(outline_macaca, SMvector = papionin$outline$semi_lm, outlines = papionin$outline$curves) # Warping the slid landmarks of the full landmark dataset to the average outline shape residual_shape <- tps.all(X_array = papionin_gpa$dataslide, REF_array = outline_gpa$dataslide, TAR_matrix = outline_gpa$mshape)data("papionin") # load dataset # Full dataset: 70 landmarks papionin_ar <- papionin$coords # Outline dataset: subset of 54 landmarks outline_ar <- papionin_ar[papionin$outline$subset, , ] # Subset: Macaca only mac <- grep("Macaca", papionin$species) # genus Macaca papionin_macaca <- papionin_ar[, , mac] outline_macaca <- outline_ar[, , mac] # Landmark sliding by minimizing bending energy + superimposition (GPA) library("Morpho") papionin_gpa <- procSym(papionin_macaca, SMvector = papionin$semi_lm, outlines = papionin$curves) outline_gpa <- procSym(outline_macaca, SMvector = papionin$outline$semi_lm, outlines = papionin$outline$curves) # Warping the slid landmarks of the full landmark dataset to the average outline shape residual_shape <- tps.all(X_array = papionin_gpa$dataslide, REF_array = outline_gpa$dataslide, TAR_matrix = outline_gpa$mshape)
Convert a matrix of landmark coordinates into a three-dimensional array
xxyy.to.array(M, p, k = 2)xxyy.to.array(M, p, k = 2)
M |
A matrix of dimension n x [p x k] containing landmark coordinates for a set of specimens. Each row contains all landmark coordinates for a single specimen. The first columns correspond to the X coordinates for all landmarks, etc. |
p |
Number of landmarks |
k |
Number of dimensions (2 or 3) |
Function returns a 3D array (p x k x n), where p is the number of landmark points, k is the number of landmark dimensions (2 or 3), and n is the number of specimens. The third dimension of this array contains names for each specimen if specified in the original input matrix.
X <- matrix(rnorm(40), nrow = 4) # Random 2D coordinates of 5 landmarks for 4 specimens xxyy.to.array(X, 5, 2)X <- matrix(rnorm(40), nrow = 4) # Random 2D coordinates of 5 landmarks for 4 specimens xxyy.to.array(X, 5, 2)