## 3D model geometric morphometry # This code is for use with the R programming language. The software # for running this code can be freely downloaded from http://www.r-project.org/ # The free integrated development environment, RStudio, can also # be used with this code and downloaded from http://www.rstudio.com/ # Data used in the Buddha scanning project are freely available # online here: https://digital.lib.washington.edu/dspace/handle/1773/22739 # To acquire data, 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 read in # pst files that were made by Landmark (http://www.idav.ucdavis.edu/research/EvoMorph) # get file names, assuming that all data files are in a folder called 'pts' # You will need to change the file location ("E:\\My Docs, etc.) # to suit your computer 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) 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]) } ## load libraries library(shapes) library(geomorph) # calculate Generalised Procrustes analysis (two methods) out <- procGPA(lmks5) shapepca(out,pcno=3) # plot out1 <- gpagen(lmks5) # makes a plot also # visualize PC relationships dev.off() # clear previous graphics pch = c(rep("L", 3), rep("S", 4)) 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() # clear previous graphics shapes.cva(out$rotated, groups = c(rep("VTL", 3), rep("VVS", 4)) ) # Tests to examine differences in mean shape # between two independent populations A <- lmks5[,,1:3] # VTL samples B <- lmks5[,,4:7] # VVS samples testmeanshapes(A, B, resamples = 1000) # 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)} # p. 163 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)} Hotellingsp<-function(SSef, SSer, dfef, dfer, exact=F) {library(MASS) # p corresponds to the number of shape space dimensions. p <- qr(SSef+SSer)$rank k<-dfef; w<-dfer s<-min(k,p) m<-(w-p-1)/2 t1<-(abs(p-k)-1)/2 Ht<-sum(diag(SSef%*%ginv(SSer))) Fapprox<-Ht*(2 * (s*m+1))/(s^2*(2*t1+s+1)) ddfnum<-s*(2*t1+s+1) ddfden<-2*(s*m+1) pval= 1-pf(Fapprox, ddfnum, ddfden) if (exact) {b<-(p+2*m)*(k+2*m)/((2*m+1)*(2*m-2)) c1<-(2+(p*k+2)/(b-1))/(2*m) Fapprox<-((4+(p*k+2)/(b-1))/(p*k))*(Ht/c1) ddfnum<-p*k ddfden<-4+(p*k+2)/(b-1)} unlist(list("dfeffect"=dfef,"dferror"=dfer,"T2"=Ht, "Approx_F"=Fapprox,"df1"=ddfnum,"df2"=ddfden,"p"=pval))} 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))} tps2d<-function(M, matr, matt) {p<-dim(matr)[1]; q<-dim(M)[1]; n1<-p+3 P<-matrix(NA, p, p) for (i in 1:p) {for (j in 1:p){ r2<-sum((matr[i,]-matr[j,])^2) P[i,j]<- r2*log(r2)}} P[which(is.na(P))]<-0 Q<-cbind(1, matr) L<-rbind(cbind(P,Q), cbind(t(Q),matrix(0,3,3))) m2<-rbind(matt, matrix(0, 3, 2)) coefx<-solve(L)%*%m2[,1] coefy<-solve(L)%*%m2[,2] fx<-function(matr, M, coef) {Xn<-numeric(q) for (i in 1:q) {Z<-apply((matr-matrix(M[i,], p, 2, byrow=T))^2, 1, sum) Xn[i]<-coef[p+1]+coef[p+2]*M[i,1]+coef[p+3]* M[i,2]+sum(coef[1:p]*(Z*log(Z)))} Xn} } tps<-function(matr, matt, n){ xm<-min(matt[,1]) ym<-min(matt[,2]) xM<-max(matt[,1]) yM<-max(matt[,2]) rX<-xM-xm; rY<-yM-ym a<-seq(xm-1/5*rX, xM+1/5*rX, length=n) b<-seq(ym-1/5*rX, yM+1/5*rX, by=(xM-xm)*7/(5*(n-1))) m<-round(0.5+(n-1)*(2/5*rX+ yM-ym)/ 11 (2/5*rX+ xM-xm)) M<-as.matrix(expand.grid(a,b)) ngrid<-tps2d(M,matr,matt) plot(ngrid, cex=0.2,asp=1,axes=F,xlab="",ylab="") for (i in 1:m){lines(ngrid[(1:n)+(i-1)*n,])} for (i in 1:n){lines(ngrid[(1:m)*n-i+1,])}} ## end Claude's functions library(MASS) library(shapes) # my data AP<-orp(pgpa(lmks5)$rotated) m<-t(matrix(AP, (dim(lmks5)[1] * dim(lmks5)[2]), dim(lmks5)[3])) n<-dim(m)[1] fact <- as.factor(c(rep("VTL", 3), rep("VVS", 4))) # calculate the plot hierarchical cluster analysis plot(hclust(dist(m), method="complete"), main="", labels=fact,cex=0.7,xlab=NA, sub=NA) # calculate the plot Principle Components Analysis 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) # test for differences between within-group and between-group # variance-covariance matrices with the adapted # Hotelling-Lawley multivariate test. mod1<-lm(m~as.factor(fact)) dfef<- length(levels(fact))-1 dfer<- n - length(levels(fact)) SSef<-(n-1)*var(mod1$fitted.values) SSer<-(n-1)*var(mod1$residuals) # adapted Hotelling-Lawley multivariate test. Hotellingsp(SSef, SSer, dfef, dfer)