create.pw.be {prWarp}

Principal warps

Description

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.

Usage

create.pw.be(A, M_ref, d = NULL)

Arguments

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)

Value

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)

References

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 Also

See CreateL for the creation of the bending energy matrix

Examples

# 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")