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]>wrote: > > Subject: morphometric functions for R users > From: marta rufino <[email protected]> > Date: Thu, 15 Jul 2010 08:56:37 -0700 (PDT) > To: morphmet_test_000 <[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/) > ======================================================== > >
