Subject: morphometric functions for R users
From: marta rufino <marta.m.ruf...@gmail.com>
Date: Thu, 15 Jul 2010 08:56:37 -0700 (PDT)
To: morphmet_test_000 <morphmet_test_...@googlegroups.com>

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/)
========================================================

Reply via email to