Basic operations with 3D models and geometric morphometric analysis in R (GPA, PCA, clustering, visualization)

3d-geom-morph.R

## 3D model geometric morphometry

# start with NextEngine 3D scanner, acquire model and save in PLY format
# open PLY file in IDAV's free Landmark software and place landmarks
# export landmark point co-ords as a PTS file, then continue here:

## read in pst files of landmark co-ord that were made by Landmark (http://www.idav.ucdavis.edu/research/EvoMorph)
# placing landmarks can be done in R with 'shapes' or 'geomorph', but is easier with Landmark

# get file names
setwd("E:\\My Documents\\My UW\\Research\\0812 Buddha Scanning\\Landmarks\\pts")
pts <- list.files(pattern="pts$")

# operate on filenames by reading in each text file as a single vector
lmks1 <- lapply(pts,  scan, what = "character")

## clean, tidy and arrange the data

# drop non-data parts...
lmks2 <- lapply(1:length(lmks1), function(i) lmks1[[i]][-(1:3)])

# arrange vector into a 4-column table (where each row is 
# set of 3d landmarks and a label)...
lmks3 <- lapply(1:length(lmks2), function(i) matrix(lmks2[[i]], ncol = 4, nrow = length(lmks2[[i]])/4, byrow = TRUE)) 

# convert 3d co-ords to numeric...
lmks4 <- lapply(1:length(lmks3), function(i) apply(lmks3[[i]][,2:4], 2, as.numeric))
# just get 30 landmarks (since that's the min number for this particular sample)
lmks4 <-  lapply(1:length(lmks4), function(i) lmks4[[i]][ 1: 30  ,   ] )

# convert list of matrices to array (input format for gpa)
lmks5 <- simplify2array(lmks4)

## inspect the data

# view rotatable 3d plots of each specimen to check for outliers, etc
library(rgl)
for(i in 1:dim(lmks5)[3]) {
  open3d()
  points3d(lmks5[,,i])
}

## do calculations...
library(shapes)
library(geomorph)

# calculate Generalised Procrustes analysis (two methods)
out <- procGPA(lmks5)
out1 <- gpagen(lmks5) # makes a plot also

# visualize PC relationships
pch = c(rep("L", 3), rep("S", 4)) # specific sample names for this dataset
par(mfrow=c(2,2))
plot(out$rawscores[,1],out$rawscores[,2],xlab="PC1",ylab="PC2", pch = pch)
title("PC scores")
plot(out$rawscores[,2],out$rawscores[,3],xlab="PC2",ylab="PC3", pch = pch)
plot(out$rawscores[,1],out$rawscores[,3],xlab="PC1",ylab="PC3", pch = pch)
plot(out$size,out$rho,xlab="size",ylab="rho", pch = pch)
title("Size versus shape distance")

# visualize Canonical variate analysis 
dev.off()
shapes.cva(out$rotated, groups = c(rep("VTL", 3), rep("VVS", 4)) )



shapepca(out,pcno=3)



          
# functions from Claude, J. (2008). Morphometrics with R. New York, Springer.          
          
angle2d <- function(v1,v2)
  {v1<-complex(1,v1[1],v1[2])
     v2<-complex(1,v2[1],v2[2])
     (pi+Arg(v1)-Arg(v2))%%(2*pi)-pi}

angle3<-function(v1, v2)
   {a<-angle(v1, v2)
      b<-sign( det(rbind(1, v1, v2)) )
      if (a == 0 & b == 1){jo<-pi/2}
      else if (a == 0 & b == -1){jo<- - pi/2}
      else {jo<- a * b}
      (pi+jo)%%(2*pi)-pi}


aligne<-function(A)
  {B<-A
      n<-dim(A)[3]; k<-dim(A)[2]
      for (i in 1:n)
        {sv<-svd(var(A[,,i]))
           M<-A[,,i]%*%sv$u
           v1<-A[2,,i]-A[1,,i]; v2<-A[3,,i]-A[1,,i]
           V1<-M[2,]-M[1,]; V2<-M[3,]-M[1,]
          
           if (k ==2)
             {if (round(angle2d(v1,v2),3)!=
                       round(angle2d(V1,V2),3))
               {M[,1]=-M[,1]}}
           if (k ==3)
             {if (round(angle3(v1,v2),3)!=
                       round(angle2d(V1,V2),3))
               {M[,1]=-M[,1]}}
           B[,,i]<-M}
      B}

trans1<-function(M){scale(M,scale=F)}

centsiz<-function(M)
   {p<-dim(M)[1]
      size<-sqrt(sum(apply(M, 2,var))*(p-1))
      list("centroid_size" = size,"scaled" = M/size)}
          
ild2<-function(M1, M2){sqrt(apply((M1-M2)^2, 1, sum))}

pPsup<-function(M1,M2)
   {k<-ncol(M1)
      Z1<-trans1(centsiz(M1)[[2]])
      Z2<-trans1(centsiz(M2)[[2]])
     sv<-svd(t(Z2)%*%Z1)
      U<-sv$v; V<-sv$u; Delt<-sv$d
      sig<-sign(det(t(Z1)%*%Z2))
      Delt[k]<-sig*abs(Delt[k]) ; V[,k]<-sig * V[,k]
      Gam<-U%*%t(V)
      beta<-sum(Delt)
      list(Mp1=Z1%*%Gam,Mp2=Z2, rotation=Gam,
       DP=sqrt(sum(ild2(Z1%*%Gam, Z2)^2)),rho=acos(beta))}

mshape<-function(A){
   apply(A, c(1,2), mean)}




pgpa<-function(A)
 # Extract the number of landmarks, coordinate dimensions, and number of configurations contained
 # in the array.
 {p<-dim(A)[1];k<-dim(A)[2];n<-dim(A)[3]
   # Create an empty array for storing scaled and rotated configurations.
   temp2<-temp1<-array(NA, dim=c(p,k,n)); Siz<-numeric(n)
   # Translate every configuration by aligning their centroid with the origin, and scale them to unit
   # size.
    for (i in 1:n)
      {Acs<-centsiz(A[,,i])
         Siz[i]<-Acs[[1]]
         temp1[,,i]<-trans1(Acs[[2]])}
   # Define the quantity Qm that should be minimized.
    Qm1<-dist(t(matrix(temp1,k*p,n)))
    Q<-sum(Qm1); iter<-0
   # Loop until differences between shape coordinates do not decrease anymore.
    while (abs(Q)>0.00001)
      {for(i in 1:n){
       # Define the mean shape ignoring the configuration that is going to be rotated.
        M<-mshape(temp1[,,-i])
       # Perform a partial Procrustes superimposition between the mean shape and each configuration.
        temp2[,,i]<-pPsup(temp1[,,i],M)[[1]]}
          Qm2<-dist(t(matrix(temp2,k*p,n)))
          Q<-sum(Qm1)-sum(Qm2)
          Qm1<-Qm2
          iter=iter+1
         temp1<-temp2}
    list("rotated"=temp2,"it.number"=iter,"Q"=Q,"intereucl.dist"=
              Qm2,"mshape"=centsiz(mshape(temp2))[[2]],"cent.size"=Siz)}


        orp<-function(A)
            {p<-dim(A)[1];k<-dim(A)[2];n<-dim(A)[3]
                Y1<-as.vector(centsiz(mshape(A))[[2]])
                oo<-as.matrix(rep(1,n))%*%Y1
                I<-diag(1,k*p)
                mat<-matrix(NA, n, k*p)
                for (i in 1:n){mat[i,]<-as.vector(A[,,i])}
                Xp<-mat%*%(I-(Y1%*%t(Y1)))
                Xp1<-Xp+oo
                array(t(Xp1), dim=c(p, k, n))}

## end Claude's functions

          library(MASS)
          library(shapes)

          # my data: cluster analysis
          AP<-orp(pgpa(lmks5)$rotated)
          m<-t(matrix(AP, (dim(lmks5)[1] * dim(lmks5)[2]), dim(lmks5)[3]))
          n<-dim(m)[1]
          fact <- c(rep("VTL", 3), rep("VVS", 4))
          par(mar=c(0.5,2,1,1))
          layout(matrix(c(1,2),2,1))
          plot(hclust(dist(m), method="average"),main="UPGMA"
               ,labels=fact,cex=0.7)
          plot(hclust(dist(m), method="complete"),main="COMPLETE"
               ,labels=fact,cex=0.7)
          
          # my data: PCA
          pcs<-prcomp(m)
          par(mar=c(4,4,1,1))
          layout(matrix(1:4,2,2))
          plot(pcs$x, pch=pch)
          barplot(pcs$sdev^2/sum(pcs$sdev^2),ylab="% of variance")
          title(sub="PC Rank",mgp=c(0,0,0))
          
          
          mesh<-as.vector(mshape(AP))
          max1<-matrix(mesh+max(pcs$x[,1])*pcs$rotation[,1], dim(lmks5)[1], dim(lmks5)[2])
          min1<-matrix(mesh+min(pcs$x[,1])*pcs$rotation[,1], dim(lmks5)[1], dim(lmks5)[2])
          # joinline<-c(1,6:8,2:5,1,NA,7,4)
          plot(min1,axes=F,frame=F,asp=1,xlab="",ylab="",pch=22)
          points(max1,pch=17)
          title(sub="PC1",mgp=c(0,0,0))
          # lines(max1[joinline,],lty=1)
          # lines(min1[joinline,],lty=2)
          max2<-matrix(mesh+max(pcs$x[,2])*pcs$rotation[,2], dim(lmks5)[1], dim(lmks5)[2])
          min2<-matrix(mesh+min(pcs$x[,2])*pcs$rotation[,2], dim(lmks5)[1], dim(lmks5)[2])
          plot(min2,axes=F,frame=F,asp=1,xlab="",ylab="",pch=22)
          points(max2,pch=17)
          title(sub="PC2",mgp=c(0,0,0))
          # lines(max2[joinline,],lty=1)
          # lines(min2[joinline,],lty=2)