R/1_repair_soma.R

In ComputationalIntelligenceGroup/3DSomaMS: 3DSomaMS: A univocal definition of the neuronal soma morphology using Gaussian mixture models

adjust_faces <- function(list_vertices, actual_faces) {
    # Vertices of the faces.
    dendritic_faces <- sort(unique(c(actual_faces)))
    # The list of vertices that were removed and still are defined in the faces matrix are obtained.
    vertices_deleted <- sort(setdiff(dendritic_faces, list_vertices))
    
    # Those faces that have a vertex that were deleted are removed.
    new_faces <- actual_faces
    new_faces <- new_faces[which(!(new_faces[, 1] %in% vertices_deleted)), ]
    new_faces <- new_faces[which(!(new_faces[, 2] %in% vertices_deleted)), ]
    new_faces <- c(new_faces[which(!(new_faces[, 3] %in% vertices_deleted)), ])
    
    # Isolated vertices, that is, those that do not belong to any face, are removed.
    verticesIsolated <- which(!(list_vertices %in% unique(c(new_faces))))
    if (length(verticesIsolated) > 0) {
      list_vertices <- list_vertices[-c(verticesIsolated)]
    }
    
    
    # Gaps in the numeration are removed reordering the index
    new_faces <- as.numeric(factor(new_faces))
    new_faces <- matrix(new_faces, ncol = 3)
    
    # Return vertices and faces of the mesh.
    return(list(vertices = list_vertices, faces = new_faces))
}

#' Repair soma morphology 
#' 
#' Repair the surface of the neuron to remove holes and vertices inside the mesh
#'
#' @param directory path to the folder where PLY files representing damaged soma mesh
#' @param output_ambient_occlusion path to the folder where somas with ambient occlusion are saved   
#' @param broken_mesh path to the folder where somas were saved after the GMM threshold was imposed and vertices inside the soma or forming cavities were removed
#' @param output_poisson_reconstruction path to the folder where somas were saved after the mesh closing process
#' 
#' @return None
#' 
#' @examples 
#' repair_soma(directory = file.path(tempdir(),"PLYs"), output_ambient_occlusion = file.path(tempdir(),"ambient_occlusion"), broken_mesh = file.path(tempdir(),"broken_mesh_ao"), output_poisson_reconstruction = file.path(tempdir(),"poisson_reconstruction_ao"))
#' 
#' @export
repair_soma <- function(directory, output_ambient_occlusion, broken_mesh, output_poisson_reconstruction) {
    if (!file.exists(file.path(directory))){
      stop("The reading directory does not exist")
    }
    
    if (!file.exists(file.path(output_ambient_occlusion))){
      if (output_ambient_occlusion == "") {
        stop("The ambient occlusion directory is not defined")
      }else{
        warning("The ambient occlusion directory does not exist, it was created automatically")
        dir.create(output_ambient_occlusion,showWarning = TRUE, recursive = TRUE)
      }
    }
    
    if (!file.exists(file.path(broken_mesh))){
      if (broken_mesh == "") {
        stop("The broken mesh directory is not defined")
      }else{
        warning("The broken mesh directory does not exist, it was created automatically")
        dir.create(broken_mesh,showWarning=TRUE, recursive = TRUE)
      }
    }
    
    if (!file.exists(file.path(output_poisson_reconstruction))){
      if (output_poisson_reconstruction == "") {
        stop("The poisson reconstruction directory is not defined")
      }else{
        warning("The poisson reconstruction directory does not exist, it was created automatically")
        dir.create(output_poisson_reconstruction, showWarning = TRUE, recursive = TRUE)
      }
    }
    
    PLY_files <- list.files(directory)[which(file_ext(list.files(directory)) == "ply")]
    
    for (file in PLY_files) {
        # Extract soma name from the file.
        soma_name <- file_path_sans_ext(file)
    
        # Execute ambient occlusion script.
        output <- execute_meshlab_script(input = file.path(directory, file), output = file.path(output_ambient_occlusion,soma_name), meshlab_script="ambient_occlusion")
      
        # Read the .ply file generated by the script and extract vertices, faces and ambient occlusion values.
        soma <- read_binary_PLY(output)

        # Clustering to differentiate surface points from interior points.
        quality_classification <- Mclust(soma$quality, G = 2)
        surface_vertices <- which(quality_classification$classification == 2)
        reparation <- adjust_faces(surface_vertices, as.matrix(soma$faces))
        adjusted_vertex <- soma$vertices[reparation$vertices, ]
        
  
        # Export results to a ply file
        vcgPlyWrite(tmesh3d(rbind(t(soma$vertices[reparation$vertices, ]), 1), t(reparation$faces)), file.path(broken_mesh,soma_name), addNormals=T)
        
        remove(adjusted_vertex)
        remove(reparation)
        
        output <- execute_meshlab_script(input = paste0(file.path(broken_mesh,soma_name),".ply"), output = paste0(file.path(output_poisson_reconstruction,soma_name),".ply"), meshlab_script="poisson_reconstruction")
    }
    return(0)
}

#' Repair single soma from data
#' 
#' Repair the surface of the soma to remove holes and vertices inside the mesh
#'
#' @param soma A filepath or the structure returned by read_VRML or read_PLY representing damaged soma mesh @seealso \code{read_VRML}
#' @param ao.path folder where intermediate files with ambient occlusion data are saved
#' @param broken.path Folder where intermediate files with broken somas after GMM threshold are saved
#' @param repaired.path Folder where final file is stored
#' @param name Name of the final file
#' @param returnSoma
#' 
#' @return Repaired soma (if returnSoma is True)
#' 
#' @useDynLib SomaMS
#' 
#' @export
repair.soma3d <- function( soma, ao.path = tempdir(), broken.path = tempdir(), repaired.path = tempdir(), name="soma", returnSoma = F ){
  
  ###
  # Validate inputs
  ###
  
  # Validate repaired
  stopifnot( length(name)>0 )
  
  # Validate tempdirs
  stopifnot( utils.validateDir(ao.path, pre="Ambient Occlusion", createDir = T) )
  stopifnot( utils.validateDir(broken.path, pre="Broken Mesh", createDir = T) )
  stopifnot( utils.validateDir(repaired.path, pre="Repaired Meshes", createDir = T) )
  
  ######

  # If soma is already a ply path keeps it, if its a vrml read and validates  and if is a list validates.
  soma.path <- utils.somaToPLY(soma, paste("original",name,sep="_"), moveToLocation = F)
  stopifnot(!is.null(soma.path))
  
  # Call meshlab
  soma.ao.path <- file.path(ao.path,paste0("ao_",name,".ply"))
  ret <- execute_meshlab_script(input = soma.path, output = soma.ao.path, meshlab_script="ambient_occlusion")
  if( ret != 0 )
    stop(sprintf("Error executing poisson reconstruction. Error code: %d - out: %s - in %s",ret,soma.ao.path,soma.path))
  
  # Read the .ply file generated by the script and extract vertices, faces and ambient occlusion values.
  soma.ao <- read_binary_PLY(soma.ao.path)

  # Clustering to differentiate surface points from interior points.
  quality_classification <- Mclust(soma.ao$quality, G = 2)
  surface_vertices <- which(quality_classification$classification == 2)
  reparation <- adjust_faces(surface_vertices, as.matrix(soma.ao$faces))
  adjusted_vertex <- soma.ao$vertices[reparation$vertices, ]
  
  # Export results to a ply file
  soma.broken.path <-file.path(broken.path,paste0("broken_",name))
  if( vcgPlyWrite(tmesh3d(rbind(t(soma.ao$vertices[reparation$vertices, ]), 1), t(reparation$faces)), soma.broken.path, addNormals=T) != 0)
    stop("Error while writing PLY broken soma file")
  ## update name name (.ply is added by vcgPlyWrite)
  soma.broken.path <- paste0(soma.broken.path,".ply")
  
  
  # Clean up a little
  remove(adjusted_vertex)
  remove(reparation)
  
  # repaired file 
  soma.repaired.path <- file.path(repaired.path,paste0("repaired_",name,".ply"))
  ret <- execute_meshlab_script(input = soma.broken.path, output = soma.repaired.path, meshlab_script="poisson_reconstruction")
  if( ret != 0)
    stop(sprintf("Error executing poisson reconstruction. Error code: %d",ret))
  
  # If returnSoma read and return it
  return( ifelse(returnSoma, read_binary_PLY(soma.repaired.path), soma.repaired.path) )
}