R/get_meshes.R

In amgfernandes/fishatlas_R: R client utilities for interacting with the MPI Fish Brain Atlas project

#' @title Read 3D larval zebrafish template brain from Kunst et al. 2019
#'
#' @description Read the template brain for a 6 day old larval zebrafish, generated by Kunst et al. 2019,
#' straight from https://fishatlas.neuro.mpg.de.
#' 3D triangular  meshes are returned as a \code{nat} package \code{\link[nat]{hxsurf}}
#' object, which mimics the Amira surface format. These can be be plotted in 3D using \code{rgl} and analysed with tools from the \code{nat} ecosystem.
#' This incldue subseting by neuropil, i.e.. if you only want to visualise or analyse the forebrain.
#' @param type which version of the template brain you want to retreive. See details. Defaults to retreiving an
#' expert segmnetation of the template brain into 51 neuropils
#' @param fishatlas_url the web location of the fish brain atlas project
#' @param surf Use if you want to add the downloaded brain to a hxsruf object already in your environment. Defaults to null, a new object is made from scratch
#' @details A single 3D brain object orits neuropils are read from .txt files hosted on https://fishatlas.neuro.mpg.de, and
#' turned into a \code{\link[nat]{hxsurf}} object in R, which mimics the Amira surface format.
#' As of May 2019, data available is:
#'\itemize{
#'  \item 	\emph{outline}	an outline from the zebrafish brain
#'  \item 	\emph{fiber}	 a 3D mesh based on a tract segmentation of the brain
#'  \item 	\emph{cell_bodies}	 a 3D mesh based on segmenting cell bodies of the brain
#'  \item 	\emph{brain_regions}	 51 expert annotated sub-regions of the brain. These include the major
#'  subdivisions beneth the forebrain, midbrain and hindbrain.
#'  Named in $Regionlist as 'major.neuropil_minor.neuropil_sub.region', e.g Forebrain_Diencephalon_Habenula
#'  }
#'  Note that since neurons are reconstructed from many different neural species,
#'  there is no 'standard' orientation between species, but within a species these neurons are registered to a
#'  template brain, usually using elastix.
#' @return a \code{nat} package \code{\link[nat]{hxsurf}} object, which mimics the Amira surface format, replete with metadata that can be
#' accessed using \code{$}
#' @seealso \code{\link{fishatlas_download_neurons}}, \code{\link{fishatlas_cell_types}}, \code{\link{fishatlas_read_saved_neurons}}
#' @examples
#' \dontrun{
#' ## Lets have a look at these neuropils!
#' zfbrain = fishatlas_read_brain(type = "brain_regions")
#' plot3d(zfbrain)
#'
#' ## Maybe we just want to plot for forebrain
#' clear3d()
#' plot3d(subset(zfbrain, "Forebrain"), alpha = 0.5, col = "red")
#'
#' ## Maybe within the forebrain, we are interested in plotting the habenula
#' plot3d(subset(zfbrain, "Habenula"), alpha = 0.3, col = "red")
#'
#' }
#' @references Kunst, Michael, Eva Laurell, Nouwar Mokayes, Anna Kramer, Fumi Kubo, António M. Fernandes, Dominique Förster, Marco Dal Maschio, and Herwig Baier. 2019. “A Cellular-Resolution Atlas of the Larval Zebrafish Brain.” Neuron, May. https://doi.org/10.1016/j.neuron.2019.04.034.
#' @export
#' @rdname fishatlas_read_brain
fishatlas_read_brain <- function( type = c("brain_regions", "outline","fiber", "cell_bodies"),
                                  fishatlas_url = "https://fishatlas.neuro.mpg.de",
                                  surf = NULL){
  type = match.arg(type)
  if(type %in% c("outline","fiber", "cell_bodies")){
    brain_info = fishatlas_fetch(path = "/neurons/get_brain",
                                 body = NULL,
                                 parse.json = TRUE,
                                 simplifyVector=FALSE,
                                 include_headers = FALSE,
                                 fishatlas_url = "https://fishatlas.neuro.mpg.de")
    brain.points = read_model_vertices(point.file = paste0(fishatlas_url,brain_info$brain[[type]]$file))
    brain.mesh = build_mesh(brain.points, surf = surf, region = type)
  }else{
    brain.mesh = fishatlas_read_brain_regions(fishatlas_url=fishatlas_url, surf = surf)
  }
  brain.mesh
}

# hidden
fishatlas_read_brain_regions <- function( fishatlas_url = "https://fishatlas.neuro.mpg.de",
                                          surf = NULL){
  brain_regions = fishatlas_fetch(path = "neurons/get_brain_regions",
                                  body = NULL,
                                  parse.json = TRUE,
                                  simplifyVector=FALSE,
                                  include_headers = FALSE,
                                  fishatlas_url = "https://fishatlas.neuro.mpg.de")
  roi.mesh = surf
  roi.files = nams = c()
  for(i in 1:length(brain_regions$brain_regions)){
    major.region = brain_regions$brain_regions[[i]]$name
    for(j in 1:length(brain_regions$brain_regions[[i]]$sub_regions)){
      middle.region = brain_regions$brain_regions[[i]]$sub_regions[[j]]$name
      for(k in 1:length(brain_regions$brain_regions[[i]]$sub_regions[[j]]$sub_regions) ){
        sub.region = brain_regions$brain_regions[[i]]$sub_regions[[j]]$sub_regions[[k]]$name
        roi.file = brain_regions$brain_regions[[i]]$sub_regions[[j]]$sub_regions[[k]]$file
        roi.files = c(roi.files, roi.file)
        nams = c(nams, paste(major.region, middle.region, sub.region, sep = "_"))
      }
    }
  }
  max = length(roi.files)
  cols = grDevices::rainbow(length(roi.files))
  for(i in 1:max){
    roi = roi.files[i]
    col = cols[i]
    region = nams[i]
    roi.points = read_model_vertices(point.file = paste0(fishatlas_url,roi))
    roi.mesh = build_mesh(vertices = roi.points,
                          surf = roi.mesh,
                          region = region,
                          color = col)
    fishatlas_progress(x = i, max=max, message = paste0("downloading ", region))
  }
  names(roi.mesh$RegionList) = NULL
  roi.mesh
}

# hidden
read_model_vertices <- function(point.file){
  txt = utils::read.table(point.file)
  vertices = data.frame()
  vertices = matrix(unlist(txt), byrow = TRUE, ncol=3)
  colnames(vertices) = c("X","Y","Z")
  vertices = as.data.frame(vertices)
  vertices
}

# hidden
build_mesh <- function(vertices,
                       surf = NULL,
                       region = "Interior",
                       color = "#FFCC66"){
  # are we adding to an extant surf object?
  if(is.null(surf)){
    surf = list()
    offset = 0
  }else{
    offset = max(surf$Vertices$PointNo)
  }
  # get indices for triangular mesh
  if(is.null(vertices$PointNo)){
    vertices$PointNo = 1:nrow(vertices)
  }
  V = matrix(vertices$PointNo, byrow = TRUE, ncol = 3)
  colnames(V) = c("V1","V2","V3")
  rownames(V) = 1:nrow(V)
  # there are many duplicated points, let's remove
  mm = match(data.frame(t(vertices)), data.frame(t(vertices)))
  vertices$PointNo = mm + offset
  V = matrix(mm[V] + offset, byrow = FALSE, ncol = 3)
  vertices = vertices[!duplicated(vertices$PointNo),]
  # build mesh
  surf$Vertices = rbind(surf$Vertices, vertices)
  surf$Regions[[region]] = V
  surf$RegionList = c(surf$RegionList, region)
  surf$RegionColourList = c(surf$RegionColourList, color)
  class(surf) = "hxsurf"
  surf
}