# Packages
library(Morpho)
library(rgl)
library(geomorph)
library(Rvcg)
library(mvnormtest)

# Template creation
# Here we detailed only the creation of the humerus template. For other bones, we only give the point rows and the specimens selected as template at the end of the script.

template.lm <- read.pts("Ceratotherium_simum_85_32_M_Humerus_LM_surface.pts") # Point data import
template.mesh <- ply2mesh(file="Ceratotherium_simum_85_32_M_Humerus.ply") # Mesh import

fixed <- as.integer(c(1:35)) # Anatomical landmark (LM) rows
patch <- as.matrix(template.lm[675:1233,]) # Matrix with sliding surface semi-LM
template.lm2 <- template.lm[1:674,] # Anatomical LM and sliding curve semi-LM
slidings <- as.integer(c(36:1233)) # Rows of all sliding semi-LM (curves and surfaces)
surface <- as.integer(c(675:1233)) # Rows of all sliding surface semi-LM

# Definition of the point rows of all the curves
SC1d <- as.integer(c(36:155)) # Bicipital groove
SC2d <- as.integer(c(156:181)) # Lesser tubercle crest
SC3d <- as.integer(c(182:294)) # Humeral head
SC4d <- as.integer(c(295:320)) # Greater tubercle crest
SC5d <- as.integer(c(321:372)) # M. infraspinatus insertion
SC6d <- as.integer(c(373:398)) # Deltoid tuberosity
SC7d <- as.integer(c(399:451)) # Tricipital line
SC8d <- as.integer(c(452:477)) # Epicondylar crest
SC9d <- as.integer(c(478:622)) # Distral trochlea ridge
SC10d <- as.integer(c(623:648)) # Trochlea / capitulum limit
SC11d <- as.integer(c(649:674)) # Medial development of the trochlea

curves <- list(SC1d,SC2d,SC3d,SC4d,SC5d,SC6d,SC7d,SC8d,SC9d,SC10d,SC11d) # Concatenation of all curves into a single object

atlas <- createAtlas(template.mesh, as.matrix(template.lm2), patch = patch, corrCurves = curves, patchCurves = NULL, keep.fix = fixed) # Atlas creation
plotAtlas(atlas, pt.size = 0.5, render = "s", point="s", legend=F) # Atlas visualization for point checking

# Data import
Ceratotherium_simum_85_32_M_Humerus <- read.pts("Ceratotherium_simum_85_32_M_Humerus.pts")
Ceratotherium_simum_2005_297_Humerus_miroir <- read.pts("Ceratotherium_simum_2005_297_Humerus_miroir.pts")
Ceratotherium_simum_2018_143_Humerus <- read.pts("Ceratotherium_simum_2018_143_Humerus.pts")
Ceratotherium_simum_3086_Humerus_miroir <- read.pts("Ceratotherium_simum_3086_Humerus_miroir.pts")
Ceratotherium_simum_19904_Humerus <- read.pts("Ceratotherium_simum_19904_Humerus.pts")
Ceratotherium_simum_35208_Humerus <- read.pts("Ceratotherium_simum_35208_Humerus.pts")
Ceratotherium_simum_cottoni_C20_Humerus <- read.pts("Ceratotherium_simum_cottoni_C20_Humerus.pts")
Ceratotherium_simum_cottoni_C32_Humerus <- read.pts("Ceratotherium_simum_cottoni_C32_Humerus.pts")
Ceratotherium_simum_cottoni_C37_Humerus_miroir <- read.pts("Ceratotherium_simum_cottoni_C37_Humerus_miroir.pts")
Ceratotherium_simum_cottoni_C40_Humerus <- read.pts("Ceratotherium_simum_cottoni_C40_Humerus.pts")
Ceratotherium_simum_cottoni_C110_Humerus <- read.pts("Ceratotherium_simum_cottoni_C110_Humerus.pts")
Ceratotherium_simum_cottoni_C112_Humerus <- read.pts("Ceratotherium_simum_cottoni_C112_Humerus.pts")
Ceratotherium_simum_NSM_Z_2010_44_Humerus <- read.pts("Ceratotherium_simum_NSM_Z_2010_44_Humerus.pts")
Ceratotherium_simum_RG35146_Humerus_miroir <- read.pts("Ceratotherium_simum_RG35146_Humerus_miroir.pts")

Dicerorhinus_sumatrensis_1204_Humerus_miroir <- read.pts("Dicerorhinus_sumatrensis_1204_Humerus_miroir.pts")
Dicerorhinus_sumatrensis_1879_6_14_2_Humerus_miroir <- read.pts("Dicerorhinus_sumatrensis_1879_6_14_2_Humerus_miroir.pts")
Dicerorhinus_sumatrensis_1894_9_24_1_Humerus <- read.pts("Dicerorhinus_sumatrensis_1894_9_24_1_Humerus.pts")
Dicerorhinus_sumatrensis_1903_300_Humerus <- read.pts("Dicerorhinus_sumatrensis_1903_300_Humerus.pts")
Dicerorhinus_sumatrensis_1908_571_Humerus_miroir <- read.pts("Dicerorhinus_sumatrensis_1908_571_Humerus_miroir.pts")
Dicerorhinus_sumatrensis_1931_5_28_1_Humerus <- read.pts("Dicerorhinus_sumatrensis_1931_5_28_1_Humerus.pts")
Dicerorhinus_sumatrensis_1948_12_20_1_Humerus <- read.pts("Dicerorhinus_sumatrensis_1948_12_20_1_Humerus.pts")
Dicerorhinus_sumatrensis_1949_1_11_1_Humerus <- read.pts("Dicerorhinus_sumatrensis_1949_1_11_1_Humerus.pts")
Dicerorhinus_sumatrensis_2004_23_Humerus <- read.pts("Dicerorhinus_sumatrensis_2004_23_Humerus.pts")
Dicerorhinus_sumatrensis_3082_Humerus <- read.pts("Dicerorhinus_sumatrensis_3082_Humerus.pts")
Dicerorhinus_sumatrensis_A_7967_Humerus_miroir <- read.pts("Dicerorhinus_sumatrensis_A_7967_Humerus_miroir.pts")
Dicerorhinus_sumatrensis_H_6392_Humerus_miroir <- read.pts("Dicerorhinus_sumatrensis_H_6392_Humerus_miroir.pts")

Diceros_bicornis_1936_644_Humerus <- read.pts("Diceros_bicornis_1936_644_Humerus.pts")
Diceros_bicornis_1944_278_Humerus_miroir <- read.pts("Diceros_bicornis_1944_278_Humerus_miroir.pts")
Diceros_bicornis_1961_186_Humerus <- read.pts("Diceros_bicornis_1961_186_Humerus.pts")
Diceros_bicornis_1961_187_Humerus <- read.pts("Diceros_bicornis_1961_187_Humerus.pts")
Diceros_bicornis_1962_166_Humerus <- read.pts("Diceros_bicornis_1962_166_Humerus.pts")
Diceros_bicornis_9714_Humerus <- read.pts("Diceros_bicornis_9714_Humerus.pts")
Diceros_bicornis_50002040_Humerus <- read.pts("Diceros_bicornis_50002040_Humerus.pts")
Diceros_bicornis_50002046_Humerus <- read.pts("Diceros_bicornis_50002046_Humerus.pts")
Diceros_bicornis_RG2133_Humerus <- read.pts("Diceros_bicornis_RG2133_Humerus.pts")

Rhinoceros_sondaicus_1205F_Humerus_miroir <- read.pts("Rhinoceros_sondaicus_1205F_Humerus_miroir.pts")
Rhinoceros_sondaicus_1861_3_11_1_Humerus_miroir <- read.pts("Rhinoceros_sondaicus_1861_3_11_1_Humerus_miroir.pts")
Rhinoceros_sondaicus_1871_12_29_7_Humerus <- read.pts("Rhinoceros_sondaicus_1871_12_29_7_Humerus.pts")
Rhinoceros_sondaicus_1921_5_15_1_Humerus <- read.pts("Rhinoceros_sondaicus_1921_5_15_1_Humerus.pts")
Rhinoceros_sondaicus_50002041_Humerus_miroir <- read.pts("Rhinoceros_sondaicus_50002041_Humerus_miroir.pts")
Rhinoceros_sondaicus_50002043_Humerus_miroir <- read.pts("Rhinoceros_sondaicus_50002043_Humerus_miroir.pts")
Rhinoceros_sondaicus_A_7970_Humerus <- read.pts("Rhinoceros_sondaicus_A_7970_Humerus.pts")
Rhinoceros_sondaicus_A_7971_Humerus <- read.pts("Rhinoceros_sondaicus_A_7971_Humerus.pts")

Rhinoceros_unicornis_1208_Humerus_miroir <- read.pts("Rhinoceros_unicornis_1208_Humerus_miroir.pts")
Rhinoceros_unicornis_1884_1_22_1_2_Humerus_miroir <- read.pts("Rhinoceros_unicornis_1884_1_22_1_2_Humerus_miroir.pts")
Rhinoceros_unicornis_1885_734_Humerus <- read.pts("Rhinoceros_unicornis_1885_734_Humerus.pts")
Rhinoceros_unicornis_1932_49_Humerus <- read.pts("Rhinoceros_unicornis_1932_49_Humerus.pts")
Rhinoceros_unicornis_1950_10_18_5_Humerus <- read.pts("Rhinoceros_unicornis_1950_10_18_5_Humerus.pts")
Rhinoceros_unicornis_1960_59_Humerus <- read.pts("Rhinoceros_unicornis_1960_59_Humerus.pts")
Rhinoceros_unicornis_1961_5_10_1_Humerus <- read.pts("Rhinoceros_unicornis_1961_5_10_1_Humerus.pts")
Rhinoceros_unicornis_1967_101_Humerus <- read.pts("Rhinoceros_unicornis_1967_101_Humerus.pts")
Rhinoceros_unicornis_1972_822_Humerus <- read.pts("Rhinoceros_unicornis_1972_822_Humerus.pts")
Rhinoceros_unicornis_33382_Humerus <- read.pts("Rhinoceros_unicornis_33382_Humerus.pts")

# Creation of an array with all specimens
array.lm <- bindArr(as.matrix(Ceratotherium_simum_85_32_M_Humerus),as.matrix(Ceratotherium_simum_2005_297_Humerus_miroir),as.matrix(Ceratotherium_simum_2018_143_Humerus),as.matrix(Ceratotherium_simum_3086_Humerus_miroir),as.matrix(Ceratotherium_simum_19904_Humerus),as.matrix(Ceratotherium_simum_35208_Humerus),as.matrix(Ceratotherium_simum_cottoni_C20_Humerus),as.matrix(Ceratotherium_simum_cottoni_C32_Humerus),as.matrix(Ceratotherium_simum_cottoni_C37_Humerus_miroir),as.matrix(Ceratotherium_simum_cottoni_C40_Humerus),as.matrix(Ceratotherium_simum_cottoni_C110_Humerus),as.matrix(Ceratotherium_simum_cottoni_C112_Humerus),as.matrix(Ceratotherium_simum_cottoni_M51854_Humerus),as.matrix(Ceratotherium_simum_cottoni_M51855_Humerus),as.matrix(Ceratotherium_simum_cottoni_M51857_Humerus),as.matrix(Ceratotherium_simum_cottoni_M51858_Humerus),as.matrix(Ceratotherium_simum_NSM_Z_2010_44_Humerus),as.matrix(Ceratotherium_simum_RG35146_Humerus_miroir),as.matrix(Ceratotherium_simum_simum_M81815_Humerus),as.matrix(Dicerorhinus_sumatrensis_1204_Humerus_miroir),as.matrix(Dicerorhinus_sumatrensis_1879_6_14_2_Humerus_miroir),as.matrix(Dicerorhinus_sumatrensis_1894_9_24_1_Humerus),as.matrix(Dicerorhinus_sumatrensis_1903_300_Humerus),as.matrix(Dicerorhinus_sumatrensis_1908_571_Humerus_miroir),as.matrix(Dicerorhinus_sumatrensis_1931_5_28_1_Humerus),as.matrix(Dicerorhinus_sumatrensis_1948_12_20_1_Humerus),as.matrix(Dicerorhinus_sumatrensis_1949_1_11_1_Humerus),as.matrix(Dicerorhinus_sumatrensis_2004_23_Humerus),as.matrix(Dicerorhinus_sumatrensis_3082_Humerus),as.matrix(Dicerorhinus_sumatrensis_A_7967_Humerus_miroir),as.matrix(Dicerorhinus_sumatrensis_H_6392_Humerus_miroir),as.matrix(Dicerorhinus_sumatrensis_M54763_Humerus),as.matrix(Dicerorhinus_sumatrensis_M81892_Humerus),as.matrix(Diceros_bicornis_1936_644_Humerus),as.matrix(Diceros_bicornis_1944_278_Humerus_miroir),as.matrix(Diceros_bicornis_1961_186_Humerus),as.matrix(Diceros_bicornis_1961_187_Humerus),as.matrix(Diceros_bicornis_1962_166_Humerus),as.matrix(Diceros_bicornis_9714_Humerus),as.matrix(Diceros_bicornis_50002040_Humerus),as.matrix(Diceros_bicornis_50002046_Humerus),as.matrix(Diceros_bicornis_bicornis_M81805_Humerus_miroir),as.matrix(Diceros_bicornis_M27757_Humerus),as.matrix(Diceros_bicornis_M113776_Humerus),as.matrix(Diceros_bicornis_M113777_Humerus_miroir),as.matrix(Diceros_bicornis_M113778_Humerus_miroir),as.matrix(Diceros_bicornis_RG2133_Humerus),as.matrix(Rhinoceros_sondaicus_1205F_Humerus_miroir),as.matrix(Rhinoceros_sondaicus_1861_3_11_1_Humerus_miroir),as.matrix(Rhinoceros_sondaicus_1871_12_29_7_Humerus),as.matrix(Rhinoceros_sondaicus_1885_734_Humerus),as.matrix(Rhinoceros_sondaicus_1921_5_15_1_Humerus),as.matrix(Rhinoceros_sondaicus_50002041_Humerus_miroir),as.matrix(Rhinoceros_sondaicus_50002043_Humerus_miroir),as.matrix(Rhinoceros_sondaicus_A_7970_Humerus),as.matrix(Rhinoceros_sondaicus_A_7971_Humerus),as.matrix(Rhinoceros_unicornis_1208_Humerus_miroir),as.matrix(Rhinoceros_unicornis_1884_1_22_1_2_Humerus_miroir),as.matrix(Rhinoceros_unicornis_1932_49_Humerus),as.matrix(Rhinoceros_unicornis_1950_10_18_5_Humerus),as.matrix(Rhinoceros_unicornis_1960_59_Humerus),as.matrix(Rhinoceros_unicornis_1961_5_10_1_Humerus),as.matrix(Rhinoceros_unicornis_1967_101_Humerus),as.matrix(Rhinoceros_unicornis_1972_822_Humerus),as.matrix(Rhinoceros_unicornis_33382_Humerus),as.matrix(Rhinoceros_unicornis_M35759_Humerus),as.matrix(Rhinoceros_unicornis_M54454_Humerus),as.matrix(Rhinoceros_unicornis_M54456_Humerus),along=3)
                                                                       
# Definition of the rawnames of the array
dimnames(array.lm)[[3]] <- c("Ceratotherium_simum_85_32_M_Humerus","Ceratotherium_simum_2005_297_Humerus_miroir","Ceratotherium_simum_2018_143_Humerus","Ceratotherium_simum_3086_Humerus_miroir","Ceratotherium_simum_19904_Humerus","Ceratotherium_simum_35208_Humerus","Ceratotherium_simum_cottoni_C20_Humerus","Ceratotherium_simum_cottoni_C32_Humerus","Ceratotherium_simum_cottoni_C37_Humerus_miroir","Ceratotherium_simum_cottoni_C40_Humerus","Ceratotherium_simum_cottoni_C110_Humerus","Ceratotherium_simum_cottoni_C112_Humerus","Ceratotherium_simum_cottoni_M51854_Humerus","Ceratotherium_simum_cottoni_M51855_Humerus","Ceratotherium_simum_cottoni_M51857_Humerus","Ceratotherium_simum_cottoni_M51858_Humerus","Ceratotherium_simum_NSM_Z_2010_44_Humerus","Ceratotherium_simum_RG35146_Humerus_miroir","Ceratotherium_simum_simum_M81815_Humerus","Dicerorhinus_sumatrensis_1204_Humerus_miroir","Dicerorhinus_sumatrensis_1879_6_14_2_Humerus_miroir","Dicerorhinus_sumatrensis_1894_9_24_1_Humerus","Dicerorhinus_sumatrensis_1903_300_Humerus","Dicerorhinus_sumatrensis_1908_571_Humerus_miroir","Dicerorhinus_sumatrensis_1931_5_28_1_Humerus","Dicerorhinus_sumatrensis_1948_12_20_1_Humerus","Dicerorhinus_sumatrensis_1949_1_11_1_Humerus","Dicerorhinus_sumatrensis_2004_23_Humerus","Dicerorhinus_sumatrensis_3082_Humerus","Dicerorhinus_sumatrensis_A_7967_Humerus_miroir","Dicerorhinus_sumatrensis_H_6392_Humerus_miroir","Dicerorhinus_sumatrensis_M54763_Humerus","Dicerorhinus_sumatrensis_M81892_Humerus","Diceros_bicornis_1936_644_Humerus","Diceros_bicornis_1944_278_Humerus_miroir","Diceros_bicornis_1961_186_Humerus","Diceros_bicornis_1961_187_Humerus","Diceros_bicornis_1962_166_Humerus","Diceros_bicornis_9714_Humerus","Diceros_bicornis_50002040_Humerus","Diceros_bicornis_50002046_Humerus","Diceros_bicornis_bicornis_M81805_Humerus_miroir","Diceros_bicornis_M27757_Humerus","Diceros_bicornis_M113776_Humerus","Diceros_bicornis_M113777_Humerus_miroir","Diceros_bicornis_M113778_Humerus_miroir","Diceros_bicornis_RG2133_Humerus","Rhinoceros_sondaicus_1205F_Humerus_miroir","Rhinoceros_sondaicus_1861_3_11_1_Humerus_miroir","Rhinoceros_sondaicus_1871_12_29_7_Humerus","Rhinoceros_sondaicus_1885_734_Humerus","Rhinoceros_sondaicus_1921_5_15_1_Humerus","Rhinoceros_sondaicus_50002041_Humerus_miroir","Rhinoceros_sondaicus_50002043_Humerus_miroir","Rhinoceros_sondaicus_A_7970_Humerus","Rhinoceros_sondaicus_A_7971_Humerus","Rhinoceros_unicornis_1208_Humerus_miroir","Rhinoceros_unicornis_1884_1_22_1_2_Humerus_miroir","Rhinoceros_unicornis_1932_49_Humerus","Rhinoceros_unicornis_1950_10_18_5_Humerus","Rhinoceros_unicornis_1960_59_Humerus","Rhinoceros_unicornis_1961_5_10_1_Humerus","Rhinoceros_unicornis_1967_101_Humerus","Rhinoceros_unicornis_1972_822_Humerus","Rhinoceros_unicornis_33382_Humerus","Rhinoceros_unicornis_M35759_Humerus","Rhinoceros_unicornis_M54454_Humerus","Rhinoceros_unicornis_M54456_Humerus")

# Projection step
patched <- placePatch(atlas, array.lm, path="./") # Projection of the template on all specimens
checkLM(patched, path="./", atlas = atlas, render = "s") # Point checking after the projection step

# Relaxation step
template.tot <- rbind(as.matrix(template.lm2),patch) # Matrix with all the points of the template
relaxed.template <- array(dim=dim(patched)) # Relaxation of the points on all the meshes
for (i in 1:dim(patched)[3])
{
  relaxed.template[,,i]<-relaxLM(patched[,,i], template.tot, 
                                 slidings, outlines = curves,
                                 surp = surface, sur.name= paste("./",dimnames(patched)[[3]][i],".ply",sep=""), 
                                 mesh = NULL, tol = 1e-09,
                                 deselect = FALSE, inc.check = TRUE, iterations = 5, fixRepro = TRUE,
                                 missing = NULL, bending = TRUE, stepsize = 1,
                                 use.lm = NULL)
}
dimnames(relaxed.template)[[3]]<-dimnames(patched)[[3]]
checkLM(relaxed.template, path="./", atlas=atlas, render="s",alpha=1) # Point checking after the relaxation step

# Sliding step
slided<-slider3d(relaxed.template, slidings, outlines = curves, surp = surface,
                 sur.path = ".", sur.name = NULL, meshlist = NULL, ignore = NULL,
                 sur.type = "ply", tol = 1e-07, deselect = FALSE, inc.check = T,
                 recursive = TRUE, iterations = 2, initproc = TRUE,
                 pairedLM = 0, mc.cores = 2,
                 fixRepro = TRUE) 
checkLM(slided$dataslide, path="./", atlas=atlas, render="s") # Point checking after the sliding step

couleurs <- c("grey","grey","grey","grey","grey","grey","grey","grey","grey","grey","grey","grey","grey","grey","red","red","red","red","red","red","red","red","red","red","red","red","black","black","black","black","black","black","black","black","black","blue","blue","blue","blue","blue","blue","blue","blue","green","green","green","green","green","green","green","green","green","green") # Manual colour definition
labels_humerus <- c("Cs_85_32_M","Cs_2005_297","Cs_2018_143","Cs_3086","Cs_19904","Cs_35208","Csc_C20","Csc_C32","Csc_C37","Csc_C40","Csc_C110","Csc_C112","Cs_Z_2010_44","Cs_RG35146","Ds_1204","Ds_1879_6_14_2","Ds_1894_9_24_1","Ds_1903_300","Ds_1908_571","Ds_1931_5_28_1","Ds_1948_12_20_1","Ds_1949_1_11_1","Ds_2004_23","Ds_3082","Ds_A_7967","Ds_H_6392","Db_1936_644","Db_1944_278","Db_1961_186","Db_1961_187","Db_1962_166","Db_9714","Db_50002040","Db_50002046","Db_RG2133","Rs_1205F","Rs_1861_3_11_1","Rs_1871_12_29_7","Rs_1921_5_15_1","Rs_50002041","Rs_50002043","Rs_A_7970","Rs_A_7971","Ru_1208","Ru_1884_1_22_1_2","Ru_1885_734","Ru_1932_49","Ru_1950_10_18_5","Ru_1960_59","Ru_1961_5_10_1","Ru_1967_101","Ru_1972_822","Ru_33382") # Manual label definition

gpa <- gpagen(slided$dataslide) # Procrustes superimposition
acp <- plotTangentSpace(gpa$coords) # PCA
plot(acp$pc.scores,pch=19,col=couleurs, asp = 1) # PCA plot

# Correlation test between the five first PCs and the centroid size (only the two first PCs - and the third for the fibula - are given in the article)
PC1 <- as.vector(pca$PCscores[,1])
log_CS <- as.vector(log10(pca$size))
summary(lm(PC1~log_CS)) # Linear regression summary
correlation_PC1 <- cor.test(PC1,log_CS, method="pearson") # Correlation test

PC2 <- as.vector(pca$PCscores[,2])
summary(lm(PC2~log_CS))
correlation_PC2 <- cor.test(PC2,log_CS, method="pearson")

PC3 <- as.vector(pca$PCscores[,3])
summary(lm(PC3~log_CS))
correlation_PC3 <- cor.test(PC3,log_CS, method="pearson")

PC4 <- as.vector(pca$PCscores[,4])
summary(lm(PC4~log_CS))
correlation_PC4 <- cor.test(PC4,log_CS, method="pearson")

PC5 <- as.vector(pca$PCscores[,5])
summary(lm(PC5~log_CS))
correlation_PC5 <- cor.test(PC5,log_CS, method="pearson")

# Neighbour Joining Tree Computation
distances = dist(pca$PCscores)
arbre_nj = nj(distances)
plot.phylo(arbre_nj, type = 'u', tip.color = couleurs, show.tip.label = FALSE)
tiplabels(labels_humerus,  col = couleurs, frame = "none", bg = NULL)

# Theoretical shapes visualization
meanshape <- as.matrix(gpa$consensus) # Meanshape  coordinates
closest_lm <- as.matrix(slided$dataslide[,,49]) # Slided coordinates of the specimen closest to the meanshape (here Rhinoceros_unicornis_1960_59_Humerus)
closest_mesh <- ply2mesh(file="Rhinoceros_unicornis_1960_59_Humerus.ply")
humerus_meanshape <- tps3d(closest_mesh, closest_lm, meanshape) # Mapping of the specimen (reference) towards the meanshape (target)
shade3d(humerus_meanshape, col = "ivory") # Visualization
vcgPlyWrite(humerus_meanshape, filename = "humerus_meanshape") # Export of the meanshape mesh

# Extraction of the coordinates of the minimal and maximal theoretical shapes for the PC1
Min_PC1 <- as.matrix(acp$pc.shapes$PC1min)
Max_PC1 <- as.matrix(acp$pc.shapes$PC1max)

humerus_Min_PC1 <- tps3d(humerus_meanshape, meanshape, Min_PC1) # Mapping of the meanshape (reference) towards the Min PC1 (target)
shade3d(humerus_Min_PC1, col = "ivory") # Visualization
spheres3d(Min_PC1, radius = 0.001) # Landmark projection for checking
vcgPlyWrite(humerus_Min_PC1, filename = "humerus_Min_PC1") # Export of the Min PC1 mesh
humerus_Max_PC1 <- tps3d(humerus_meanshape, meanshape, Max_PC1) # Mapping of the meanshape (reference) towards the Max PC1 (target)
shade3d(humerus_Max_PC1, col = "ivory") # Visualization
spheres3d(Max_PC1, radius = 0.001) # Landmark projection for checking
vcgPlyWrite(humerus_Max_PC1, filename = "humerus_Max_PC1") # Export of the Max PC1 mesh

# Extraction of the coordinates of the minimal and maximal theoretical shapes for the PC2
Min_PC2 <- as.matrix(acp$pc.shapes$PC2max)
Max_PC2 <- as.matrix(acp$pc.shapes$PC2min)

humerus_Min_PC2 <- tps3d(closest_mesh, closest_lm, Min_PC2) # Mapping of the meanshape (reference) towards the Min PC2 (target)
shade3d(humerus_Min_PC2, col = "ivory") # Visualization
spheres3d(Min_PC2, radius = 0.001) # Landmark projection for checking
vcgPlyWrite(humerus_Min_PC2, filename = "humerus_Min_PC2") # Export of the Min PC2 mesh
humerus_Max_PC2 <- tps3d(closest_mesh, closest_lm, Max_PC2) # Mapping of the meanshape (reference) towards the Max PC2 (target)
shade3d(humerus_Max_PC2, col = "ivory") # Visualization
spheres3d(Max_PC2, radius = 0.001) # Landmark projection for checking
vcgPlyWrite(humerus_Max_PC2, filename = "humerus_Max_PC2") # Export of the Max PC2 mesh

# Allometry
species <- as.factor(c("simum","simum","simum","simum","simum","simum","simum","simum","simum","simum","simum","simum","simum","simum","sumatrensis","sumatrensis","sumatrensis","sumatrensis","sumatrensis","sumatrensis","sumatrensis","sumatrensis","sumatrensis","sumatrensis","sumatrensis","sumatrensis","bicornis","bicornis","bicornis","bicornis","bicornis","bicornis","bicornis","bicornis","bicornis","sondaicus","sondaicus","sondaicus","unicornis","sondaicus","sondaicus","sondaicus","sondaicus","sondaicus","unicornis","unicornis","unicornis","unicornis","unicornis","unicornis","unicornis","unicornis","unicornis")) # Manual definition of the species groups
gdf <- geomorph.data.frame(gpa, species = species) # Creation of a geomorph dataframe
allometry_fem <- procD.allometry(coords~Csize, ~species, logsz = TRUE, data = gdf) # Procrustes ANOVA of shape scores (computed from the coordinates) against the centroid size, taking into account the group affiliation (species)
plot(allometry_fem, method = "RegScore", warpgrids = FALSE, label = labels_humerus, pch = 19, cex = 1)
summary(allometry_fem)

# Extraction of the coordinates of the theoretical shapes for minimal and maximal centroid size values
Min_Csize <- as.matrix(allometry_fem$Ahat.at.min)
Max_Csize <- as.matrix(allometry_fem$Ahat.at.max)

humerus_Min_Csize <- tps3d(humerus_meanshape, meanshape, Min_Csize) # Mapping of the meanshape (reference) towards the Min centroid size (target)
shade3d(humerus_Min_Csize, col = "ivory")
spheres3d(allometry_fem$Ahat.at.min, radius = 0.001)
vcgPlyWrite(humerus_Min_Csize, filename = "humerus_Min_Csize")
humerus_Max_Csize <- tps3d(humerus_meanshape, meanshape, Max_Csize) # Mapping of the meanshape (reference) towards the Max centroid size (target)
shade3d(humerus_Max_Csize, col = "ivory")
spheres3d(allometry_fem$Ahat.at.max, radius = 0.001)
vcgPlyWrite(humerus_Max_Csize, filename = "humerus_Max_Csize")

mass <- as.matrix(c(2300,2300,2300,2300,2300,2300,2300,2300,2300,2300,2300,2300,2300,775,775,775,775,775,775,775,775,775,775,775,775,1050,1050,1050,1050,1050,1050,1050,1050,2300,1350,1350,1350,2000,1350,1350,1350,1350,1350,2000,2000,2000,2000,2000,2000,2000,2000,2000)) # Manual definition of the mean body mass
mass <- mass^(1/3) # Cube root of the mean body mass
gdf_mass <- geomorph.data.frame(gpa, mass = mass) # Creation of a geomorph dataframe
allometry_mass <- procD.allometry(coords~mass, logsz = TRUE, data = gdf_mass)
plot(allometry_mass)
summary(allometry_mass)

# Extraction of the coordinates of the theoretical shapes for minimal and maximal mean body mass values
Min_mass <- as.matrix(allometry_mass$Ahat.at.min)
Max_mass <- as.matrix(allometry_mass$Ahat.at.max)

humerus_Min_mass <- tps3d(humerus_meanshape, meanshape, Min_mass) # Mapping of the meanshape (reference) towards the Min mean body mass (target)
shade3d(humerus_Min_mass, col = "ivory")
spheres3d(allometry_mass$Ahat.at.min, radius = 0.001)
vcgPlyWrite(humerus_Min_mass, filename = "humerus_Min_mass")
humerus_Max_mass <- tps3d(humerus_meanshape, meanshape, Max_mass)
shade3d(humerus_Max_mass, col = "ivory")
spheres3d(allometry_mass$Ahat.at.max, radius = 0.001)
vcgPlyWrite(humerus_Max_mass, filename = "humerus_Max_mass")

# Visualization of the differences between theoretical shapes obtained for centroid size and mean body mass (used in Figure 9)
deformGrid3d(Min_mass, Min_Csize, col1 = 2, col2 = 4, ngrid = 0, lwd = 0.5, ask = FALSE, size = 0.0002)
deformGrid3d(Max_mass, Max_Csize, col1 = 2, col2 = 4, ngrid = 0, lwd = 0.5, ask = FALSE, size = 0.0002)

# Computation of the meanshape for each species
# Ceratotherium simum
liste_simum <- c(1:14)
gpa_simum <- gpagen(slided$dataslide[,,liste_simum])
meanshape_simum <- as.matrix(gpa_simum$consensus)
closest_lm_simum <- as.matrix(slided$dataslide[,,11])
closest_mesh_simum <- ply2mesh(file="Ceratotherium_simum_cottoni_C110_Humerus.ply")
mesh_meanshape_simum <- tps3d(closest_mesh_simum, closest_lm_simum, meanshape_simum)
shade3d(mesh_meanshape_simum, col = "ivory")
spheres3d(meanshape_simum, radius = 0.001)
vcgPlyWrite(mesh_meanshape_simum, filename = "humerus_meanshape_simum")

# Dicerorhinus sumatrensis
gpa_sumatrensis <- gpagen(slided$dataslide[,,(15:26)])
meanshape_sumatrensis <- as.matrix(gpa_sumatrensis$consensus)
closest_lm_sumatrensis <- as.matrix(slided$dataslide[,,21])
closest_mesh_sumatrensis <- ply2mesh(file="Dicerorhinus_sumatrensis_1948_12_20_1_Humerus.ply")
mesh_meanshape_sumatrensis <- tps3d(closest_mesh_sumatrensis, closest_lm_sumatrensis, meanshape_sumatrensis)
shade3d(mesh_meanshape_sumatrensis, col = "ivory")
spheres3d(meanshape_sumatrensis, radius = 0.001)
vcgPlyWrite(mesh_meanshape_sumatrensis, filename = "humerus_meanshape_sumatrensis")

# Diceros bicornis
gpa_bicornis <- gpagen(slided$dataslide[,,(27:35)])
meanshape_bicornis <- as.matrix(gpa_bicornis$consensus)
closest_lm_bicornis <- as.matrix(slided$dataslide[,,33])
closest_mesh_bicornis <- ply2mesh(file="Diceros_bicornis_50002040_Humerus.ply")
mesh_meanshape_bicornis <- tps3d(closest_mesh_bicornis, closest_lm_bicornis, meanshape_bicornis)
shade3d(mesh_meanshape_bicornis, col = "ivory")
spheres3d(meanshape_bicornis, radius = 0.001)
vcgPlyWrite(mesh_meanshape_bicornis, filename = "humerus_meanshape_bicornis")

# Rhinoceros sondaicus
gpa_sondaicus <- gpagen(slided$dataslide[,,(36:43)])
meanshape_sondaicus <- as.matrix(gpa_sondaicus$consensus)
closest_lm_sondaicus <- as.matrix(slided$dataslide[,,43])
closest_mesh_sondaicus <- ply2mesh(file="Rhinoceros_sondaicus_A_7971_Humerus.ply")
mesh_meanshape_sondaicus <- tps3d(closest_mesh_sondaicus, closest_lm_sondaicus, meanshape_sondaicus)
shade3d(mesh_meanshape_sondaicus, col = "ivory")
spheres3d(meanshape_sondaicus, radius = 0.001)
vcgPlyWrite(mesh_meanshape_sondaicus, filename = "humerus_meanshape_sondaicus")

# Rhinoceros unicornis
gpa_unicornis <- gpagen(slided$dataslide[,,(44:53)])
meanshape_unicornis <- as.matrix(gpa_unicornis$consensus)
closest_lm_unicornis <- as.matrix(slided$dataslide[,,53])
closest_mesh_unicornis <- ply2mesh(file="Rhinoceros_unicornis_33382_Humerus.ply")
mesh_meanshape_unicornis <- tps3d(closest_mesh_unicornis, closest_lm_unicornis, meanshape_unicornis)
shade3d(mesh_meanshape_unicornis, col = "ivory")
spheres3d(meanshape_unicornis, radius = 0.001)
vcgPlyWrite(mesh_meanshape_unicornis, filename = "humerus_meanshape_unicornis")


### Landmark definition for all other bones

# Radius
fixed <- as.integer(c(1:23))
patch <- as.matrix(template.lm[417:1336,])
template.lm2 <- template.lm[1:416,]
slidings <- as.integer(c(24:1336))
surface <- as.integer(c(417:1336))
SC1d <- as.integer(c(24:125)) # Proximal glenoid cavity
SC2d <- as.integer(c(126:177)) # Lateral synovial articular surface
SC3d <- as.integer(c(178:202)) # Medial synovial articular surface
SC4d <- as.integer(c(203:246)) # Interosseous crest
SC5d <- as.integer(c(247:281)) # Transverse crest
SC6d <- as.integer(c(282:366)) # Articular surface for the scaphoid
SC7d <- as.integer(c(367:399)) # Articular surface for the semilunar
SC8d <- as.integer(c(400:416)) # Distal articular surface for the ulna

# Ulna
fixed <- as.integer(c(1:21))
patch <- as.matrix(template.lm[365:1186,])
template.lm2 <- template.lm[1:364,]
slidings <- as.integer(c(22:1186))
surface <- as.integer(c(365:1186))
SC1d <- as.integer(c(22:116)) # Trochlear notch
SC2d <- as.integer(c(117:167)) # Articular surface for the radius
SC3d <- as.integer(c(168:184)) # Medial ridge of the proximal articular surface for the radius
SC4d <- as.integer(c(185:208)) # Articular surface for the semilunar
SC5d <- as.integer(c(209:241)) # Articular surface for the triquetrum
SC6d <- as.integer(c(242:258)) # Articular surface for the pisiform
SC7d <- as.integer(c(259:320)) # Caudal border of the ulna
SC8d <- as.integer(c(321:364)) # Interosseous crest

# Femur
patch <- as.matrix(template.lm[640:1670,])
template.lm2 <- template.lm[1:639,]
slidings <- as.integer(c(28:1670))
surface <- as.integer(c(640:1670))
SC1d <- as.integer(c(28:105)) # Greater trochanter
SC2d <- as.integer(c(106:185)) # Femoral head
SC3d <- as.integer(c(186:211)) # Lesser trochanter
SC4d <- as.integer(c(212:255)) # Third trochanter
SC5d <- as.integer(c(256:317)) # Lateral line
SC6d <- as.integer(c(318:379)) # Medial line
SC7d <- as.integer(c(380:481)) # Trochlea ridge
SC8d <- as.integer(c(482:560)) # Medial condyle
SC9d <- as.integer(c(561:639)) # Lateral condyle

# Tibia
fixed <- as.integer(c(1:24))
patch <- as.matrix(template.lm[409:1262,])
template.lm2 <- template.lm[1:408,]
slidings <- as.integer(c(25:1262))
surface <- as.integer(c(409:1262))
SC1d <- as.integer(c(25:100)) # Articular surface of the lateral condyle
SC2d <- as.integer(c(101:159)) # Articular surface of the medial condyle
SC3d <- as.integer(c(160:203)) # Lateral part of the tibial tuberosity
SC4d <- as.integer(c(204:237)) # Proximal articular surface for the fibula
SC5d <- as.integer(c(238:306)) # Distal articular surface for the fibula
SC6d <- as.integer(c(307:382)) # Distal articular surface for the talus
SC7d <- as.integer(c(383:408)) # Interosseous crest

# Fibula
fixed <- as.integer(c(1:12))
patch <- as.matrix(template.lm[282:735,])
template.lm2 <- template.lm[1:281,]
slidings <- as.integer(c(13:735))
surface <- as.integer(c(282:735))
SC1d <- as.integer(c(13:46)) # Proximal articular surface for the tibia
SC2d <- as.integer(c(47:98)) # Distal articular surface for the tibia
SC3d <- as.integer(c(99:148)) # Distal articular surface for the talus
SC4d <- as.integer(c(149:219)) # Laterocaudal crest
SC5d <- as.integer(c(220:281)) # Laterocranial crest