############################################################# # STAT 718A Computer Practical 1 ############################################################# library(shapes) # load the `shapes' package (may need installing # on your machine by downloading from CRAN) # the digit 3 data were described in Chapter 1 help(digit3.dat) # gives some information about the data digit3.dat # Will print out the coordinates of the digit3 data # which is an array of size 13 x 2 x 30, # and we have k=13 landmarks in m=2 dimensions for each # of n=30 observations digit3.dat[,,10] # these are the coordinates of the 10th observation plot(digit3.dat[,,10]) # plots the 10th observation, but note # the x and y axes are not to scale. plotshapes(digit3.dat[,,10]) # this plot is to scale, with a different # symbol for each landmark plotshapes(digit3.dat) # all 30 observations - what a mess! plotshapes(digit3.dat,joinline=1:13) # a bit less of a mess #make the plots square par(pty="s") #plot the first 12 digits par(mfrow= c(3,4) ) for (i in 1:12){ plot(digit3.dat[,,i], xlim=c(0,50), ylim= c(-50,0) ) } #plot the next 12 digits par(mfrow= c(3,4) ) for (i in 13:24){ plot(digit3.dat[,,i], xlim=c(0,50), ylim= c(-50,0) ) } #plot the last 6 digits par(mfrow= c(3,4) ) for (i in 25:30){ plot(digit3.dat[,,i], xlim=c(0,50), ylim= c(-50,0) ) } ###################################### # Centroid size ###################################### k<-13 m<-2 n<-30 size <- rep(0,times=n) for (i in 1:n){ size[i] <- centroid.size( digit3.dat[,,i] ) } size plot(size) hist(size) qqnorm(size) # check for normality qqnorm(log(size)) summary(size) ######################################## # Bookstein shape variables ####################################### bookst <- array( 0, c(k,m,n) ) for (i in 1:n){ bookst[,,i] <- bookstein.shpv( digit3.dat[,,i] ) } #plot Bookstein shape variables par(pty="s") par(mfrow=c(1,1)) plot(bookst[,,1], xlim=c(-6,6), ylim=c(-4,8) xlab="x", ylab="y") for (i in 1:n){ points(bookst[,,i]) lines(bookst[,,i]) } # alternatively plotshapes(bookst) # calculate a mean shape digit3.mean <- matrix( 0, k, m) for (i in 1:n){ digit3.mean <- digit3.mean + bookst[,,i] } digit3.mean <- digit3.mean/n plotshapes(digit3.mean) ## there is a command in R to do the above # bookstein2d() library(help=shapes) # in order to find the commands help(bookstein2d) # tells you about this specific command ans <- bookstein2d( digit3.dat) ans\$bshpv # contains the Bookstein shape variables ans\$mshape # contains the Bookstein mean shape #warning : only use bookstein2d() with the default baseline choice (l1=1,l2=2) # as I just found a bug in the program for other choices!! # some pairwise plots x <- cbind( size , t(ans\$bshpv[1:4,1,]) , t(ans\$bshpv[1:4,2,])) pairs( x , label = c("s", "u1", "u2", "u3", "u4", "v1", "v2", "v3", "v4") ) ################################### # Procrustes match ################################### plotshapes( digit3.dat[,,1], digit3.dat[,,2]) # match 1st on to 2nd y <- digit3.dat[,,2] w <- digit3.dat[,,1] help(procOPA) # match w onto y (using translation, scale and rotation) proc <- procOPA( y, w ) proc plotshapes( proc\$Ahat , proc\$Bhat) # match y onto w (using translation, scale and rotation) proc <- procOPA( w, y ) proc plotshapes( proc\$Ahat , proc\$Bhat) # NB compare fitted scales, fitted rotation matrices, using both methods # Finally : if you have time then you could try the same commands # out on another dataset, e.g. mouse vertebrae: qset2.dat