-------- Original Message --------
Subject: Re: GPA possible based on subset of landmarks?
Date: Mon, 27 Feb 2012 13:58:17 -0500
From: annat haber <[email protected]>
To: [email protected]
Hi Dominique,
If you have only two configurations, you can use my R function
unifyVD.R, pasted in the end and also downloadable from Rohlf's website
http://life.bio.sunysb.edu/morph/
It's designed to unify the ventral and dorsal aspects of the skull based
on a subset of landmarks that are taken from both aspects. Your problem
sounds exactly the same, only your "ventral" and "dorsal" configurations
are two different specimens, if I understand correctly.
But, since you probably want to simultaneously superimpose more than two
configurations you need to modify it and combine it with the gpa
function - a version of which is also pasted in the end of this email
("PROC.R"). If you're familiar with R and can do this yourself, then the
two functions should give you a good starting point. Otherwise, I'll try
to do it when I get a chance.
Basically, the easiest procedure, I think, would be to translate every
complete configuration based on the centroid of the subset, then use the
PROC function (or the procGPA function in package shapes) on all the
subsets. Then find the rotation matrix between each original subset and
its superimposed coordinates and use this rotation matrix to rotate the
complete configuration(after it has been already translated onto the
subset centroid) using this bit of code (repeat this part separately for
each subset):
M<- t(na.omit(S)) %*% na.omit(O) # S is the superimposed coordinates and
O is the original translated coordinates of a specific subset.
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
Xsi<- Xc %*% RM# rotate the translated complete configuration Xc
You should of course also scale everything to centroid size before doing
the gpa, but I'm actually not sure how exactly (input anyone?). I
suppose you should scale every configuration by the centroid of its subset.
I hope that helps (and correct...). Please feel free to contact me
directly with any questions or comments.
Cheers,
Annat
On Mon, Feb 27, 2012 at 7:21 PM, morphmet
<[email protected]
<mailto:[email protected]>> wrote:
-------- Original Message --------
Subject: GPA possible based on subset of landmarks?
Date: Thu, 16 Feb 2012 18:36:59 -0500
From: Dominique Adriaens <[email protected]
<mailto:[email protected]>>
To: <[email protected] <mailto:[email protected]>>
Dear all,
Does anyone know how it would be possible to perform GPA on a set of
landmark configurations, where the actual GPA is performed based on a
subset of landmarks in the configurations but where the other landmarks
are modified accordingly (but are not involved in the GPA itself). We
want to combine different sets of digitised landmark configurations,
being a set of landmarks on bony elements (which comprise a subset of
landmarks shared in all configurations) and different sets of landmarks
on muscles (a sequence of images are used, where muscles are
subsequently peeled off). The final goal is then to pool all the muscle
landmarks with that shared set of bony landmarks into one dataset,
after
they have been superimposed based on the coordinates of the shared
‘bony
landmarks’.
Any suggestion would be welcome!
Thanks and best regards
Dominique
*Prof. Dr. Dominique Adriaens*
Ghent University
Evolutionary Morphology of Vertebrates & Zoology Museum
K.L. Ledeganckstraat 35, B-9000 Gent
BELGIUM
tel: +32 9 264.52.19 <tel:%2B32%209%20264.52.19>, fax: +32 9
264.53.44 <tel:%2B32%209%20264.53.44>
E-mail: [email protected]
<outbind://5-__000000002B873A8E03CAD411AE9E00__5004995AACC4103900/[email protected]
<mailto:[email protected]>>
URL: http://www.fun-morph.ugent.be/
http://www.zoologymuseum.__ugent.be/
<http://www.zoologymuseum.ugent.be/>
http://www.evolutietheorie.be/
Logo labo klein
### unifyVD
## This function works with both 2D and 3D; no need to specify which it is
# X is a configuration of one specimen in the format of k landmarks by m
dimensions
# "ventral" and "dorsal" are vectors indicating the names (NOT indices)
# of the ventral and dorsal landmarks, respectively.
# comV and comD are vectors indicating the common landmarks' names (NOT
indices) for the ventral and dorsal sides respectively (in the same order);
# the function attaches the ventral side to the dorsal side
# (or any other two aspects of the same structure) based on the common
landmarks.
# 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 landmarks 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
landmarks are averaged
# and listed only for the dorsal, so the output will have fewer rows
than the input.
# NOTE: the resulting configuration will be re-organized to have all the
dorsal
# landmarks first, then all the ventral ones.
# It is therefore HIGHLY RECOMMENDED that landmarks will be designated
by names in X
# (i.e., rownames(X) is a charcter vector rather than NULL)
# The output includes the unified data ($unified) and the unification
errors ($errors)
# an example for an input file and a protocol that implements this
function can be downloaded from
# http://home.uchicago.edu/~annat/
# Please email me with comments and questions [email protected]
<mailto:[email protected]>
# last updated Feb 28 2011
unifyVD<- function(X, dorsal, ventral, comD, comV, average=FALSE) {
V <- X[comV,]
D <- X[comD,]
Xv<- X[ventral,]
Xd <- X[dorsal,]
V[which(is.na <http://is.na>(D))] <- NA
D[which(is.na <http://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)
dv[comV[which(is.na <http://is.na>(dv[comV,1]))],] <-
dv[comD[which(is.na <http://is.na>(dv[comV,1]))],]
dv[comD[which(is.na <http://is.na>(dv[comD,1]))],] <-
dv[comV[which(is.na <http://is.na>(dv[comD,1]))],]
if (average==TRUE) {
dv[comD,] <- (dv[comD,]+dv[comV,])/2
dv<- dv[-match(comV, rownames(dv)),]}
list(unified=dv, errors=sqrt(rowSums((Xvr[comV,]-Dc)^2)))}
#### PROC
# AA is an array of the k x m x N format; k landmakrs, m dimensions, N
specimens; partial procrustes superimposition;
scale.centroid<- function(X) {X/sqrt(sum(X^2, na.rm=TRUE))} # X is
centered already
fitX <- function(X, Ref) {
X[which(is.na <http://is.na>(Ref))] <- NA
Ref[which(is.na <http://is.na>(X))] <- NA
M<- t(na.omit(Ref)) %*% na.omit(X)
SVD <- svd(M)
V<- SVD$v# eigenvectors scaled to 1 unit length
U <- SVD$u
L <- diag(SVD$d)
S <- ifelse(L<0, -1, L)
S <- ifelse(L>0, 1, L)
RM<- V%*% S%*% t(U) # the rotation matrix
Xr<- X%*% RM} # rotate X onto the reference
PROC<- function(AA, sizescaled=TRUE) {
Xc<- array(apply(AA, 3, scale, scale=FALSE), dim=dim(AA)) # centered
if(sizescaled==TRUE) {Xs<- array(apply(Xc, 3, scale.centroid),
dim=dim(AA))} else{Xs<- Xc}
R <- Xs[,,1]
d <- 1
while (d>0.001) {
Xf<- array(apply(Xs, 3, fitX, Ref=R), dim=dim(AA), dimnames=dimnames(AA))
MC<- matrix(nr=nrow(Xf), nc=ncol(Xf)) # mean configuration
for (i in 1:nrow(Xf)) {for (j in 1:ncol(Xf)) {MC[i,j] <- mean(Xf[i,j,],
na.rm=T)}}
for (i in 1:nrow(R)) { if (is.na <http://is.na>(R[i,1])) MC[i,]<-NA}
d <- sqrt(sum((MC-R)^2, na.rm=T))
R <- MC
Xs <- Xf}
list(rotated=Xf, opt.stat=d)}
--
Annat Haber, PhD
Postdoc fellow
Tel-Aviv University
Department of Zoology
Tel Aviv 69978, Israel
c: +972 50 2517 901
http://home.uchicago.edu/~annat/ <http://home.uchicago.edu/%7Eannat/>
