#############################################################
#              STAT 718A  Computer Practical 2
#############################################################

library(shapes)  

# we'll do some more analysis of the digit 3 data

plotshapes(digit3.dat, joinline=c(1:13) )

# Carry out Procrustes analysis (calculated using complex eigenvectors)

out <- procGPA( digit3.dat , eigen2d = TRUE ) 


out$rotated     # contains the Procrustes registered data

plotshapes ( out$rotated , joinline=c(1:13) )   # procrustes registered data

out$mshape    # contains the Procrustes mean shape

lines( out$mshape, col=2 )  # adds the mean shape to the plot (in red)

plot( out$size , out$rho )  # NOTE `rho' is called the Riemannian distance
                        # and it is given by   sin rho = d_F , 
                        # where d_F is the full Procrustes distance 
                        # given in our notes. 
                        # here out$rho is a vector of shape distances to the mean

text( out$size , out$rho , c(1:30) )

# is there a possible outlier ? 

plotshapes( digit3.dat , digit3.dat[,,1] , joinline= 1:13)


# To see all the calculations given by the command procGPA() 
# type 

summary( out )

#and look at

help( procGPA )


# Shape variability

print( out$rmsd1 ) #full Procrustes distances
print( out$rho )   #Riemannian distances


# We'll now look at the first few principal components

out$percent

# plot first two PC scores

plot(out$rawscores[,1],out$rawscores[,2], xlab="pc1",ylab="pc2")

# pairwise plots - do you see any interesting structure ?

pairs( cbind( out$size, out$rho, out$rawscores[,1:3] ) , 
           label=c("size","rho","pc1","pc2","pc3") )


# Let us interpret the structure of shape variability in the 
# first 3 PCs

shapepca ( out , 1:3 , type="r" , joinline=1:13 )  # row plot
  
shapepca ( out , 1:3 , type="v" , joinline=1:13 )  #vector plot

shapepca ( out , 1:3 , type="s" , joinline=1:13 )  #superposition plot

shapepca ( out , 1:3 , type="m" , joinline=1:13 )  # movie plot


##################################################################

# Now let us remove the 1st digit as it looks a bit strange. 

newdigit3.dat <- digit3.dat[,,-1]

# check dimension 

dim(newdigit3.dat)

#### carry out the above analysis on the digit3 data but with the 
#### outlier removed. 

#Is there a big difference in mean shape and shape variability  
#with or without the outlier ?

#################################################################### 

#Let's now compare the Procrustes mean shape and Bookstein mean shapes.

#Consider the female Chimpanzee data

help(panf.dat)

dim(panf.dat) 

proc <- procGPA( panf.dat )
plotshapes( panf.dat , proc$mshape )

book <- bookstein2d( panf.dat ) 
plotshapes( book$bshpv , book$mshape ) 


final <- procOPA( proc$mshape , book$mshape )

plotshapes( proc$mshape , final$Bhat) 

print( riemdist( proc$mshape, book$mshape ) )

# is there much difference ?? 

#####################################################################