R code: Reading TPS and PLY files

[morphmet]

-------- Original Message --------
Subject: R code: Reading TPS and PLY files
Date: Fri, 13 Jan 2012 17:24:59 -0500
From: Dean Adams <dcad...@iastate.edu>
To: morphmet@morphometrics.org <morphmet@morphometrics.org>

Hi all,

Below are a few R-routines I've written that I thought folks might find
useful.

The first is a function for reading PLY files from NextEngine and David
Scanners, and getting the data into R. For other scanners, the code may
need to be adjusted somewhat, as different scanners output slightly
different PLY files.  The second is a function for reading TPS landmark
data files.  The pixel coordinates are multiplied by the scale factor if
one is included (note: an earlier post by M. Rubio provided code for
reading outline data in TPS format).

I also have a general GPA routine for performing 2D and 3D GPA using
points, curves, surfaces, or any combination. For semilandmarks, the
routine implements both sliding approaches described in the literature
(i.e., via bending energy and via Procrustes distance). I will be
cleaning this code up soon and hope to distribute it in the near future.

All the best,

Dean

--
Dr. Dean C. Adams
Associate Professor
Department of Ecology, Evolution, and Organismal Biology
Department of Statistics
Iowa State University
Ames, Iowa
50011
www.public.iastate.edu/~dcadams/
phone: 515-294-3834


###############
read.ply<-function(file){            #Dean C. Adams
  plyfile<-scan(file=file,what="char",sep="\n",strip.white=T,quiet=T)
    #header section
  xline<-unlist(strsplit(grep(c("vertex "),plyfile, value=T)," "))
  npoints<-as.numeric(xline[grep(c("vertex"),xline)+1])
  yline<-unlist(strsplit(grep(c("element face"),plyfile, value=T)," "))
  npoly<-as.numeric(yline[grep(c("face"),yline)+1])
  headerend<-grep(c("end_header"),plyfile)
    #3D points
  ncolpts<-(length(grep(c("property"),plyfile))-1)
  cols<-grep(c("property"),plyfile, value=T)      #finds x,y,z cols
from points section
  x<-grep(c(" x"),cols);y<-grep(c(" y"),cols);z<-grep(c(" z"),cols)

points<-as.matrix(as.numeric(unlist(strsplit(plyfile[(headerend+1):(headerend+npoints)]," 

"))))
  dim(points)<-c(ncolpts,npoints)
  points<-t(points)
  xpts<-points[,x];ypts<-points[,y];zpts<-points[,z]
  points<-cbind(xpts,ypts,zpts)
    #polygons
  plyrest<-plyfile[(headerend+npoints+1):(headerend+npoints+npoly)]
  size1<-as.matrix(as.numeric(unlist(strsplit(plyrest[1]," "))))[1]
  sizecheck<-function(A){
size<-as.matrix(as.numeric(unlist(strsplit(A," "))))[1]}
  sizeall<-unlist(lapply(plyrest,sizecheck))
  if (min(sizeall)!=max(sizeall)) print("Polygons not saved: faces have
different number of vertices") else {

poly<-as.matrix(as.numeric(unlist(strsplit(plyfile[(headerend+npoints+1):(headerend+npoints+npoly)]," 

"))))
    dim(poly)<-c((poly[1]+1),npoly) #poly[1] = #pts per face
    poly<-poly[-1,]
    ifelse(min(poly)==0,poly<-poly+1,poly<-poly)  #b/c DavidScanner
begins count @ 0,NULL
  }
return(list(coords=points,polygons=poly))
}

##########

readland.tps<-function(file){            #Dean C. Adams: inspired by M.
Rufino's code on Morphmet
  tpsfile<-scan(file=file,what="char",sep="\n",strip.white=T,quiet=T)
  lmdata<-grep("LM=",tpsfile) # landmark points
  nspecs<-length(lmdata)
  scale<-as.numeric(sub("SCALE=","",tpsfile[grep("SCALE",tpsfile)]))
  if (is.null(scale)){scale=array(1,nspecs)}          #NOTE: assumes
scale is present for all, or none
  nland<-as.numeric(sub("LM=","",tpsfile[lmdata]))
    #landmark data, scaled by scale values
  landdata<-NULL                #extract landmark coordinates, multiply
by scale factor
  for(i in 1:nspecs){
    for(j in 1:nland[1]){
    tmp<-as.numeric(unlist(strsplit(tpsfile[lmdata[i]+j]," ")))*scale[i]
    landdata<-rbind(landdata,tmp)
     }
   }
  imageID<-(sub("IMAGE=","",tpsfile[grep("IMAGE",tpsfile)]))
    imageID<-sub(".jpg","",imageID)          #note: done for .jpg
  return(list(coords=landdata,specID=imageID,nspecs=nspecs,nland=nland[1]))
}