#' @title Script related to 3D meshes dataset of sagittal otoliths from red mullet in the Mediterranean Sea
#' @author Nicolas Andrialovanirina, Remi Laffont, Emilie Poisson Caillault
#' @date revised by Nicolas Andrialovanirina 2024.05.21 

# This script identifies meshes with holes,
# delete non otolith artefacts and
# put landmarks points on 3D otoliths

# packages and work directory ####

# how to install the necessary packages
# install.packages("processx")
# install.packages("rlang")
# install.packages("ps")
# install.packages("devtools")
# devtools::install_github("morphOptics/digit3DLand", local=FALSE)
# install.packages("digit3DLand")

# load packages
library("devtools")
library("digit3DLand")

# changing work directory (folder meshes location)
setwd("your_directory") # changes depending on the directory in which your meshes are stored

# delete non otolith artefact/count genus/save ply ####

# list of files with the .stl extension in work directory
stl_names <- list.files(pattern=".stl")

# loop to explore all meshes in order to delelte non-otolith artefact and identify mesh with hole

for (i in 1:length(stl_names)){ # browse all meshes by file names
  
  if (file.exists(stl_names[i])){ # test if mesh file names exist in the work directory
    
    print(stl_names[i])
    
    msh <- vcgImport(stl_names[i], silent=TRUE) # import 3D mesh into the msh object

    # remove non otolith artefacts
    msh <- vcgIsolated(msh, silent = TRUE)
    
    # index calculation in Euler formula
    
    #V: number of vertices
    V <- ncol(msh$vb)
    
    #E: number of edges
    E <- nrow(vcgGetEdge(msh))
    
    #F: number of faces
    Fa <- ncol(msh$it)
    
    
    
    #b: number of borders of the 3D mesh
    b <- sum(vcgBorder(msh)$vb)
    
    if (b==0){
      #g: genus, number of holes
      g <- -1 * (0.5 * (V - E + Fa + b) - 1)
      
      if (g==0){ # if no hole
        
        path <- paste0(getwd(), "/", "ply_meshes") 
        
        if (!dir.exists("ply_meshes")){ # test if "ply_meshes" exists in work directory
          dir.create(path) # create if "ply_meshes" does not exist
        }
        # save mesh as ply in "ply_meshes"
        vcgPlyWrite(msh, paste0(path, "/", gsub(".stl", ".ply", stl_names[i]))) 
        
      }else{
        print(paste0(stl_names[i], " has genus: ", g)) # show number of holes
        
      }
      
      
    }else{
      print(paste0(stl_names[i], " has border elements: ", b)) # display number border
    }
  }
}



# put landmarks on meshes ####

# define number of landmarks to be taken
TeOpt <- setTemplOptions(fixed=6, template=FALSE)
# manually place 6 landmarks as shown in Figure 3 (follow the instruction in the console)
FiOpt <- setFileOptions("ply_meshes", sdir="your_directory") # changes depending on the directory in which your meshes are stored
r <- digitMesh("tps_file_name", fixed=6, FiOpt=FiOpt, TeOpt=TeOpt)