#############################################################
#              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