rm(list=ls())


library(geomorph)
source('ShapeRotator.R')


#Skulls
crania_data <- read.table("XXX.txt", header=TRUE, row.names=1, stringsAsFactors = FALSE)
rownames(crania_data)
head(crania_data)
is.numeric(crania_data)

coords_crania <- as.matrix(crania_data)
is.numeric(coords_crania)
dim(coords_crania)

skull <- arrayspecs(coords_crania, p=95, k=3)

str(skull)
skull[,,1]

plot.3D(skull, 10, colour = 'black')


#Mandibles

mandible_data <- read.table("XXX.txt", header=TRUE, row.names=1, stringsAsFactors = FALSE)
rownames(mandible_data)
head(mandible_data)
is.numeric(mandible_data)

coords_mandible <- as.matrix(mandible_data)
is.numeric(coords_mandible)
dim(coords_mandible)

mandible <- arrayspecs(coords_mandible, p=44, k=3)

str(mandible)
mandible[,,1]
plot.3D(mandible, 10, colour = 'red')



# Select the landmarks

land.a=XX
land.b=XX
land.c=XX 
land.d=XX

land.e=XX 
land.f=XX 
land.g=XX
land.h=XX


# Translate the data
data.1_t = translate(skull, land.a)
data.2_t = translate(mandible, land.e)

str(data.1_t)
str(data.2_t)

# We need to match the dim.names 
matched = match.datasets(data.1_t, data.2_t)

# The real datasets (BUT NOT JOINED YET) are now as follows, whre matched$matched1[i] is matched with mathched$matched2[i].
data.1 = matched$matched1
data.2 = matched$matched2

# The rotations return a list
dr0l = double.rotation(data.1, data.2, land.a, land.b, land.c, land.d, land.e, land.f, land.g, land.h, 0)
dr15l = double.rotation(data.1, data.2, land.a, land.b, land.c, land.d, land.e, land.f, land.g, land.h, 15)
dr45l = double.rotation(data.1, data.2, land.a, land.b, land.c, land.d, land.e, land.f, land.g, land.h, 45)
dr90l = double.rotation(data.1, data.2, land.a, land.b, land.c, land.d, land.e, land.f, land.g, land.h, 90)
dr120l = double.rotation(data.1, data.2, land.a, land.b, land.c, land.d, land.e, land.f, land.g, land.h, 120)


# Join.

dr0 = join.arrays (dr0l$rotated1, dr0l$rotated2)
dr15 = join.arrays (dr15l$rotated1, dr15l$rotated2)
dr45 = join.arrays (dr45l$rotated1, dr45l$rotated2)
dr90 = join.arrays (dr90l$rotated1, dr90l$rotated2)


# Examples of rotations for several specimens

# This set is with the skull translated to vector skull_translate, which will determine the position of skull relative to mandible
#Might have to change the y depending on the angle
skull_translate = c(2.2,-0.1, 0) #You can use the function norm3D to calculate distances, in the case you have measured the distance between
#structures in the skull and mandible

vector_1 <- c(1, 2, 3)
vector_1_length <- norm_3D(vector_1)
# 3.741657
# This way you can make sure to check it is at a distance you like.

dr0_st = join.arrays (dr0l$rotated1, translate(dr0l$rotated2,land.e, skull_translate))
plot.rotation.3D(dr0_st, data.1, data.2, 1) 

skull_translate = c(1.7,0.1, 0)
dr15_st = join.arrays (dr15l$rotated1, translate(dr15l$rotated2,land.e, skull_translate))
plot.rotation.3D(dr15_st, data.1, data.2, 1) 

skull_translate = c(1,1, 0) # You need to choose which position would be best. I just guessed what I thought that looks good
dr45_st = join.arrays (dr45l$rotated1, translate(dr45l$rotated2,land.e, skull_translate))

plot.rotation.3D(dr45_st, data.1, data.2, 1)
plot.rotation.3D(dr45_st, data.1, data.2, 2)
plot.rotation.3D(dr45_st, data.1, data.2, 4)
plot.rotation.3D(dr45_st, data.1, data.2, 5)
plot.rotation.3D(dr45_st, data.1, data.2, 6)
plot.rotation.3D(dr45_st, data.1, data.2, 7)

skull_translate = c(0,2, 0) # Again, just re-position it to where you prefer it.
dr90_st = join.arrays (dr90l$rotated1, translate(dr90l$rotated2,land.e, skull_translate))
plot.rotation.3D(dr90_st, data.1, data.2, 1)    


### FINAL DATASETS TO EXPORT

#Example with dr45_st (translated mandible at 45 degrees from skull)

writeland.tps(dr45_st, file="Heads_45degrees", scale = NULL, specID = TRUE)