# X is a configuration of one specimen in the format of k landmarks by 3 dimensions
# "ventral" and "dorsal" are vectors indicating the indices or names of the ventral and dorsal landmarks, respectively. 
# comV and comD are vectors indicating the common LM's names (NOT indices) for the ventral and dorsal sides respectively (in the same order); 
# the function attaches the ventral side to the dorsal side based on the common LM's. 
# It finds the rotation matrix that minimizes the sum of squared deviations between the ventral and the dorsal (Rohlf 1990. Proceedings of the Michigan Morphometrics Workshop), 
# and calculates the errors as the vector length between the ventral and the dorsal equivalent LM's after fitting.
# If either the dorsal or the ventral copy of a common landmark is missing, that landmark is ignored for the unification.
# Any other missing data does not affect the procedure at all and can be safely included.
# If average=TRUE then the dorsal and ventral copies of the common LM's are averaged and listed only for the dorsal, so the output will have fewr rows than the input (length(comV) fewer).

# NOTE: the resulting configuration will be re-organized to have all the dorsal landmarks first, then the ventral ones.
# an example for an input file and a protocol that implements this function (as well as the AMP) can be downloaded from http://home.uchicago.edu/~annat/
# Please email me with comments and questions annat22@gmail.com
# last updated March 25th 2009
	
unifyVD <- function(X, ventral, dorsal, comV, comD, average=FALSE) {
		V <- X[comV,]
		D <- X[comD,]
		Xv <- X[ventral,]
		Xd <- X[dorsal,]
		V[which(is.na(D))] <- NA
		D[which(is.na(V))] <- NA # making sure both ventral and dorsal of the same LM are NA's whenever one of them is
		mV <- matrix(apply(V, 2, mean, na.rm=TRUE), byrow=TRUE, nr=nrow(Xv), nc=ncol(Xv))
		Xvc <- Xv-mV # translating all the ventral LM's based on the centroid of the common ones
		Vc <- scale(V, scale=F)
		mD <- matrix(apply(D, 2, mean, na.rm=TRUE), byrow=TRUE, nr=nrow(Xd), nc=ncol(Xd))
		Xdc <- Xd-mD # translating all the dorsal LM's based on the centroid of the common ones
		Dc <- scale(D, scale=F)
		M <- t(na.omit(Dc)) %*% na.omit(Vc) # arbitrarily choosing the dorsal as the reference
		SVD <- svd(M)
		L <- diag(SVD$d) 
		S <- ifelse(L<0, -1, L)
		S <- ifelse(L>0, 1, L)
		RM <- SVD$v %*% S %*% t(SVD$u) # the rotation matrix
		Xvr <- Xvc %*% RM # rotate all the translated ventral LM's
		dv <- rbind(Xdc, Xvr)
		if (average==TRUE) {
			dv[comV[which(is.na(dv[comV,1]))],] <- dv[comD[which(is.na(dv[comV,1]))],]
			dv[comD[which(is.na(dv[comD,1]))],] <- dv[comV[which(is.na(dv[comD,1]))],]
			dv[comD,] <- (dv[comD,]+dv[comV,])/2
			dv <- dv[-which(!is.na(match(rownames(dv),comV))),]}
		list(unified=dv, errors=sqrt(rowSums((Xvr[comV,]-Dc)^2)))}