migration from cesvima

DESCRIPTION

@@ -0,0 +1,18 @@
+Package: SomaMS
+Type: Package
+Title: 3DSomaMS: A univocal definition of the neuronal soma morphology
+        using Gaussian mixture models
+Version: 1.1
+Date: 2015-03-30
+Author: Sergio Luengo-Sanchez <sluengo@fi.upm.es>
+Maintainer: Sergio Luengo_Sanchez <sluengo@fi.upm.es>
+Description: The definition of the soma is fuzzy, as there is no clear line demarcating the soma of the labeled neurons and the origin of the dendrites and axon. Thus, the morphometric analysis of the neuronal soma is highly subjective. In this paper, we provide a mathematical definition and an automatic segmentation method to delimit the neuronal soma. We applied this method to the characterization of pyramidal cells, which are the most abundant neurons in the cerebral cortex. Since there are no benchmarks with which to compare the proposed procedure, we validated the goodness of this automatic segmentation method against manual segmentation by experts in neuroanatomy to set up a framework for comparison. We concluded that there were no significant differences between automatically and manually segmented somata, i.e., the proposed procedure segments the neurons more or less as an expert does. It also provides univocal, justifiable and objective cutoffs. Thus, this study is a means of characterizing pyramidal neurons in order to objectively compare the morphometry of the somata of these neurons in different cortical areas and species.
+License: GPL (>= 2)
+Depends: R (>= 3.0.0), mclust, Morpho, Rcpp, rgl, Rvcg, tools,
+        matrixStats, pracma, geometry, misc3d, bnlearn, caret
+Suggests: doSNOW, foreach, gplots, parallel
+LinkingTo: Rcpp
+LazyLoad: yes
+Archs: x64
+NeedsCompilation: yes
+Packaged: 2015-10-26 17:18:50 UTC; luis

NAMESPACE

@@ -0,0 +1,20 @@
+# Generated by roxygen2 (4.1.1): do not edit by hand
+
+export(BN_clustering)
+export(BN_structure_learning)
+export(MAQ_S)
+export(RMSE_computation)
+export(characterize.soma3d)
+export(compute_meshes_volumes)
+export(convert_somas_to_PLY)
+export(read_VRML)
+export(read_binary_PLY)
+export(repair.soma3d)
+export(repair_soma)
+export(segment.soma3d)
+export(segment_somas)
+export(simulate_new_somas)
+export(simulation_3D_mesh)
+export(soma_caracterization)
+export(wilcoxon_RMSE)
+useDynLib(SomaMS)

R/0_convert_somas_to_PLY.R

@@ -0,0 +1,84 @@
+#' Read VRMLs of somas 
+#' 
+#' Read all the VRMLs of somas in a folder and export them to PLY format in other folder
+#'
+#' @param read_directory path to the folder where VRML files are placed
+#' @param write_directory path to the folder where PLY files will be saved    
+#' 
+#' @return None
+#' 
+#' @examples 
+#' convert_somas_to_PLY(read_directory = system.file("test/VRMLs",package="SomaMS"), write_directory = file.path(tempdir(),"PLYs"), parallel = TRUE)
+#' 
+#' @export
+convert_somas_to_PLY <- function(read_directory, write_directory, parallel = TRUE) {
+    if (!file.exists(file.path(read_directory))){
+      stop("The reading directory does not exist")
+    }
+
+    if (!file.exists(file.path(write_directory))){
+      if (write_directory == "") {
+        stop("The writing directory is not defined")
+      }else{
+        warning("The writing directory does not exist, it was created automatically")
+        dir.create(write_directory, showWarning = TRUE, recursive = TRUE)
+      }
+    }
+
+    VRML_files <- list.files(file.path(read_directory))[which(file_ext(list.files(read_directory)) == "vrml")]
+
+    if (length(VRML_files)==0){
+      warning('There are not VRMLs in the reading directory')
+    }
+
+    if(!parallel)
+    {
+      ncores <- 1
+      cl <- NULL
+    }else{
+      ncores <- detectCores()
+      if(is.na(ncores)){
+        ncores <- 1
+        cl <- NULL
+      }else{
+        ncores <- ncores - 1
+        cl <- makeCluster(ncores, type = "SOCK")
+        registerDoSNOW(cl)
+      }
+    }
+
+    if(!parallel){
+      foreach (file = VRML_files) %do% {
+        # Number 1 denotes that it must be choosen that mesh with more vertices. It is useful if there are more
+        # than one mesh save in the same file. p.e: 1 mesh with soma+dendrites and 1 mesh with the soma
+        # segmented by an expert.
+        soma <- read_VRML(file.path(read_directory, file))
+
+        soma_name <- file_path_sans_ext(file)
+        soma_name <- gsub(" ", "_", soma_name)
+
+        vcgPlyWrite(tmesh3d(rbind(t(soma$vertices), 1), t(soma$faces)), file.path(write_directory, soma_name))
+      }
+    }else{
+      foreach (file = VRML_files, .packages=c("Rvcg","Morpho","tools", "rgl"), .export = "read_VRML") %dopar% {
+        # Number 1 denotes that it must be choosen that mesh with more vertices. It is useful if there are more
+        # than one mesh save in the same file. p.e: 1 mesh with soma+dendrites and 1 mesh with the soma
+        # segmented by an expert.
+        soma <- read_VRML(file.path(read_directory, file))
+
+        soma_name <- file_path_sans_ext(file)
+        soma_name <- gsub(" ", "_", soma_name)
+
+        vcgPlyWrite(tmesh3d(rbind(t(soma$vertices), 1), t(soma$faces)), file.path(write_directory, soma_name))
+      }
+    }
+
+    if(parallel)
+    {
+      stopCluster(cl)
+    }
+
+    return(0)
+} 
+
+

R/1_repair_soma.R

@@ -0,0 +1,176 @@
+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) )
+}