Actually, thanks should go to Marta Rufino. I am but a humble medium.
-the morphmet mod (dslice)
-------- Original Message --------
Subject: Re: morphometric functions for R users
Resent-Date: Thu, 15 Jul 2010 15:31:30 -0400
Resent-From: [email protected]
Date: Thu, 15 Jul 2010 13:31:22 -0600
From: Jimena BoHe <[email protected]>
Reply-To: [email protected]
To: [email protected]
Dr. Slice
Thank you so much. They are definitely useful.
Jimena.
--
Jimena Bohórquez Herrera, M. Sc.
PhD Student on Marine Sciences.
Centro Interdisciplinario de Ciencias Marinas - Instituto Politécnico
Nacional.
La Paz, Baja California Sur, México.
Researcher
Fundación Colombiana para la Investigación y Conservación de Tiburones y
Rayas - SQUALUS
Cali - Colombia
On Thu, Jul 15, 2010 at 12:29 PM, Dennis E. Slice
<[email protected] <mailto:[email protected]>> wrote:
Subject: morphometric functions for R users
From: marta rufino <[email protected]
<mailto:[email protected]>>
Date: Thu, 15 Jul 2010 08:56:37 -0700 (PDT)
To: morphmet_test_000 <[email protected]
<mailto:[email protected]>>
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=harm2
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
}
--
Dennis E. Slice
Associate Professor
Dept. of Scientific Computing
Florida State University
Dirac Science Library
Tallahassee, FL 32306-4120
-
Guest Professor
Department of Anthropology
University of Vienna
-
Software worth having/learning/using...
Linux (Operating System: Ubuntu, CentOS, openSUSE, etc.)
OpenOffice (Office Suite: http://www.openoffice.org/)
R package (Stats/Graphics environment: http://www.r-project.org/)
Eclipse (Java/C++/etc IDE: http://www.eclipse.org/)
Netbeans (Java/C++/etc IDE: http://netbeans.org/)
Zotero (FireFox bibliographic extension: http://www.zotero.org/)
========================================================
