## X is a matrix of k landmarks by 3 dimensions of one specimen;
# "midline", "right", and "left" are vectors indicating the row indices (or the names if appear 
# as row names of X) of the midline, right and left landmarks. Make sure the left and right landmarks are 
# specified in the same order so that the first element in the "left" vector specifies the left equivalent 
# of the first element in the "right" vector, etc.
# The function aligns the specimen along the midline plane and then fills in missing data by 
# reflecting from the other side ("reflected"), averages the left and right ("averaged"), and 
# provides also the deviation of the midline landmarks from the midline plane ("errors") and the 
# vector length between the right and left after reflection ("asymmetry").
# The resulting aligned specimen will have the midline axis as z, the longest non-midline axis 
# as x, and the second non-midline axis as y. No further rotation and reflection is done, therefore the specimen may be facing any direction along each of the axes.

# NOTE: the resulting configuration will be re-organized to have all the midline landmarks first, then the right lateral ones, and then the left lateral ones.
# an example for an input file and a protocol that implements this function (as well as the unifyVD) can be downloaded from http://home.uchicago.edu/~annat/
# Please email me with comments and questions annat22@gmail.com
# last updated March 25th 2009

AMP <- function(X, midline, right, left) { 
		Xm <- X[midline,]
		Mm <- matrix(apply(Xm, 2, mean, na.rm=TRUE), byrow=TRUE, nr=nrow(X), nc=ncol(X))
		Xc <- X-Mm # translating all the LM's based on the centroid of the midline
		W <- Xc[midline,] 
		i <- which(is.na(W[,1]))
		if (length(i)>0) W<-W[-i,] # discarding missing midline LM's
		V<-svd(W)$v
		Xa <- Xc%*%V # rotating the centered specimen to align along the midline plane
		Mr <- Xa[right,]
		Ml <- Xa[left,]
		i <- which(is.na(Mr),arr.ind=TRUE)
		Mr[i[which(i[,2]==3),1],3] <- -Ml[i[which(i[,2]==3),1],3]
		Mr[which(is.na(Mr))] <- Ml[which(is.na(Mr))] # filling NA's in the right side by reflecting the left
		i <- which(is.na(Ml),arr.ind=TRUE)
		Ml[i[which(i[,2]==3),1],3] <- -Mr[i[which(i[,2]==3),1],3]
		Ml[which(is.na(Ml))] <- Mr[which(is.na(Ml))] # filling NA's in the left side by reflecting the right
		Xaf <- rbind(Xa[midline,], Mr, Ml)
		Mr[,3] <- -(Mr[,3]) # reflecting the right side onto the left for averaging
		Mav <- (Mr+Ml)/2 # averaged
		list(reflected=Xaf, averaged=rbind(Xa[midline,], Mav), errors=Xa[midline,3], asymmetry=sqrt(rowSums(Mr-Ml)^2))}