morphometric functions for R users

[marta rufino]

Dear list members,

I have been following and enjoying Claude's book outline's part
(morphometrics with R), which is great (at least for me... I am an R
fan).
Thus, I have been creating some functions to help me in some tasks. I
though that maybe there are more pleople, working on this subject,
that might be interested. Thus, I am sharing them here. The code is
not the best one (I am not a programer - for example, the 'for' should
be 'apply', etc.), but it is working :) when I have time, I can try to
improve it.
Feel free to use it and improve it :)

:)
Cheers,
M

##############################################################3
##############################################################3
# functions to do morphometric analysis
##############################################################3
# functions done after Claude's book. 'Morphometrics with R'
# Marta M. Rufino
# 15-7-2010
##############################################################3
##############################################################3
##############################################################3
##############################################################3

open.tps.outline=function(fil){
# to test: fil="D:\\MARTA\\Pos Doc\\navalha\\navalha\
\nav_outline2_no path.tps"
# fil is the file name, with respective extension and path
# note the file is better to not include the path for each image.
# X and Y are outline coordinates, ind is the individual number and
facto the fcator (i.e. image name)
# note that facto should then be arranged accordingly

# open the file
kk=scan(fil, what="char", quote="",sep="\n", strip.white=T,
quiet=T)
kk=casefold(kk, upper=F)

#detect outlines, number of points, etc.
sp=grep("points=",kk) # outline points
n=length(sp) # number of configurations
p=as.numeric(sub("points=","",kk[sp])) #number of outline points
in each animal
sp.start=sp+1 #starting location of the
sp.end=sp+p #ending position of each curve
im=kk[grep("image=", kk)]#list of images
im=sub("image=","",im)
im=sub(".tif","",im); im=sub(".jpg","",im)

# extract outline coordinates
# this would be much better, if done with a apply function... :(
newd=array(NA, c(1,4))
colnames(newd)= c("X","Y","ind","facto")
for(ii in 1:n){

kk1=data.frame(X=as.numeric(unlist(strsplit(kk[sp.start[ii]:sp.end[ii]],
" ")))[seq(1,p[ii]*2,2)],

Y=as.numeric(unlist(strsplit(kk[sp.start[ii]:sp.end[ii]], " ")))
[seq(2,p[ii]*2,2)],
ind=ii, facto=im[ii])
newd=rbind(newd, kk1)
}
newd=newd[-1,]
# plot(Y~X, newd[newd$N==1,], typ="o", asp=1) #Check plotting
xyplot(Y~X, newd, groups=facto, typ="l", aspect="iso") #Check
plotting
newd
# if(scale){
# for(ii in 1:length(table(newd$ind))){
# kk=scale(newd[newd$ind==ii,1:2], center=T, scale=F)
# kk=kk/max(kk)
# newd[newd$ind==ii,1:2]= kk
# }}
}



vharmonics=function(dat){
# vharmonics goes for visualize harmonics
# function to plot original contour (on grey) vs. NEF and
efourier fitted contour, with varying number of harmonics (5-15).
# dat= data file of the original contour coordinates, X,Y
# note that your 'dat' should be scaled and centered (this bit,
I have some doubts
names(dat)=c("X","Y")
par(mar=c(2,2,2,2), mfrow=c(1,1))
plot(dat, typ="l", main="Harmonics reconstruction of outline:
EFA",
asp=1, axes=F, ylim=c(min(dat$Y), max(dat
$Y+3)), xlim=range(dat$X))
box(col="grey90")
for(jj in 1:11){
ii=seq(0,3,.3)[jj]; kk=c(5:16)[jj]
# Using efourier
f1=efourier(dat, n=kk)
f2=iefourier(f1$an, f1$bn, f1$cn, f1$dn, k=kk, n=100)
polygon(dat$X, dat$Y+ii, col="grey80", border=0)
lines(f2$x, f2$y+ii, col=2)
text(dat$X[1], c(dat$Y+ii)[1], paste(kk,"harm."), col=2,
pos=4, cex=.7)
# Using NEF instead
f1=NEF(dat, n=kk)
f2=iefourier(f1$A, f1$B, f1$C, f1$D, k=kk, n=100)
lines(f2$x, f2$y+ii, col=4)
legend("topleft", c("efourier","NEF"), text.col=c(2,4), cex=.
7, box.col=0)
}}



nharmonics=function(dat, P){
# nharmonics goes for Number of harmonics. Function to plot and
estimate the power and variance explained according to the number of
harmonics
# dat is a data frame with all outline's coordinates: i.e. dat=newd[,
1:2] # X and Y coordinates of all outlines
# P=100 # number of points in each outline.

# Store first 20 harmonics in a file
harm=matrix(NA,N,20*4)
for(ii in 1:N){
N= dim(dat)[1]/P #total number of individuals
kk=rep(1:N,e=P)
Ne=efourier(dat[kk==ii,], n=20) #taut[,,ii] are the xy coordinates
of the points.
# When we do the NEF, we obtain 50 harmonic
coefficients of each type (A, B, C and D).
# Thus, we get 200 variables * Number of
individuals
harm[ii,]=c(Ne$an, Ne$bn, Ne$cn, Ne$dn)
# OR
# Ne=NEF(newd1[newd1$ind==ii,c("X","Y")], n=20)
# harm[ii,]=c(Ne$A, Ne$B, Ne$C, Ne$D)
}
rm(Ne)

#1st approach
co=harm^2
tharn=apply(co, 2, sum)
Power=apply(matrix(tharn, 20, 4),1,sum)
po1=round(cumsum(Power[-1])/sum(Power[-1]),3)[2:19]

# Add a curve with Cumulative variance explained by the
coeffcients
vharn=apply(harm, 2, var)
variation=apply(matrix(vharn,20,4),1,sum)
var1=round(cumsum(variation[-1])/sum(variation[-1]),3)[2:19]

# and plot it
par(las=1, mar=c(1,4,1,1), mfrow=c(1,1))
plot(po1*100, typ="o", col=2, cex=.7, pch=16, ylab=c("%
explained"), axes=F,
ylim=c(min(var1*100),101), main="How many harmonics?")
box()
axis(2)
abline(h=100, col="lightblue", lwd=2)
abline(v=c(1:19), col=8, lty=3)
abline(h=seq(0,100,5), col="lightblue", lty=3)
text(1:19, rep(101,19), c(2:19), cex=.7)
points(2:19, var1*100, typ="o", col=4, cex=.7, pch=17, lty=2)
legend("bottomright", c("Total power","Variance explained"),
text.col=c(2,4), cex=.7, box.col=0)

kk=list("Total comulative power"=po1*100, "Variance
explained"=var1*100)
kk
}