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