# 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")
Run the code above in your browser using DataLab