R/alphashape.R

In alexanderbates/catnat: Analyse neuronal morphology and connectivity (extends catmaid and nat)

#' Combine alphashapes generated by the alphashapes package
#'
#' @param ashapelist list of alphashape objects
#' @param ... Additional arguments (currently ignored)
#'
#' @references Lafarge T, Pateiro-López B, Possolo A, Dunkers J (2014) R
#'   Implementation of a Polyhedral Approximation to a 3D Set of Points Using
#'   the alpha-Shape. J Stat Softw 56
#'
#' @return An alpha shape object
#' @export
#' @seealso \code{\link{transform3dalphashape}},
#'   \code{\link{WriteVTKalphashape}}
combine.alphashape <- function (ashapelist, ...){
  initial = ashapelist[[1]]
  for (a in 2:length(ashapelist))
  {
    ashape = ashapelist[[a]]
    count = nrow(initial$x)
    initial$x <- rbind(initial$x, ashape$x)
    ashape$tetra[,1:4] <- ashape$tetra[,1:4] + count
    initial$tetra <- rbind(initial$tetra, ashape$tetra)
    ashape$triang[,1:3] <- ashape$triang[,1:3] + count
    initial$triang <- rbind(initial$triang, ashape$triang)
    ashape$edge[,1:2] <- ashape$edge[,1:2] + count
    initial$edge <- rbind(initial$edge, ashape$edge)
    ashape$vertex[,1] <- ashape$vertex[,1] + count
    initial$vertex <- rbind(initial$vertex, ashape$vertex)
  }
  return (initial)
}

#' Transform the 3D vertices of an alphashape
#'
#' @param ashape list of alphashape objects
#' @param transformations A transformation that can be accepted by nat::xform(). I.e. A registration defined by a matrix, a function, a cmtkreg object, or a character vector specifying a path to one or more registrations on disk
#' @param ... additional arguments passed to methods
#'
#' @references Lafarge T, Pateiro-López B, Possolo A, Dunkers J (2014) R Implementation of a Polyhedral Approximation to a 3D Set of Points Using the alpha-Shape. J Stat Softw 56
#'
#' @return An alpha shape object
#' @export
#' @seealso \code{\link{combine.alphashape}} \code{\link{WriteVTKalphashape}}
transform3dalphashape = function (ashape, transformations, ...){
  positions = ashape$x
  if (is.list(transformations) == F){
    cat("Single transformation")
    positions <- nat::xform(positions, transformations)
  }
  if (is.list(transformations) == T){
    for (transformation in transformations){
      positions <- nat::xform(positions, transformation)
    }
  }
  shape = alphashape3d::ashape3d(positions, alpha = ashape$alpha)
  return (shape)
}

#' Write alphashape as a .vtk file
#'
#' @param ashape List of alphashape objects
#' @param filename Path to and name of desired file. Should end in .vtk
#' @param title Title of the .vtk file
#' @param datatype Either float or double
#' @param ... additional arguments passed to methods
#'
#' @return A .vtk file saved to a given location
#' @export
#' @seealso \code{\link{combine.alphashape}} \code{\link{transform3dalphashape}}
#' @importFrom utils write.table
WriteVTKalphashape <-function(ashape, filename, title = filename, datatype=c("float","double")){
  d = ashape$x
  if(ncol(d)!=3) stop("Expect N rows x 3 cols of 3d points")
  nummarkers=nrow(d)
  datatype=match.arg(datatype)
  if(missing(title)) title=paste("Data written from R by WriteVTKLandmarks at",Sys.time())

  cat("# vtk DataFile Version 2.0",
      title,
      "ASCII",
      "DATASET POLYDATA",
      paste("POINTS",nummarkers,datatype),sep="\n",file=filename)

  write.table(d,col.names=F,row.names=F,file=filename,append=TRUE)

  data = ashape$triang
  keeps <- apply(data, 1, function(x) {( any(as.numeric(x[9]) > 1))} ) # Removes rows for triangles not included in the alphashape, for the chosen alpha. Includes other simplexes: interior, regular and singular.
  mx = data.matrix(data[keeps,][,1:3]-1) # VTK files are 0 indexed
  numpoints = rep(3, nrow(mx))
  mx = cbind(numpoints, mx)
  cat(paste("POLYGONS",nrow(mx),nrow(mx)*4),sep="\n",file=filename, append = TRUE)
  write.table(mx,col.names=F,row.names=F,file=filename,append=TRUE)
}

#' Find neurites inside of an alphashape
#'
#' @param shape An alpha shape or mesh3d object
#' @param db Neuronlist object that serves as the search database
#' @param min_nodes Number of nodes a neuron has to have within the alpha shape
#'   to be selected
#' @param ... additional arguments passed to methods
#'
#' @return A neuronlist
#' @export
#' @seealso \code{\link{combine.alphashape}},
#'   \code{\link{transform3dalphashape}}, \code{\link[nat]{find.neuron}}
neurons.inside <- function(shape, db, min_nodes = 1) {
  if (is.null(shape$triang)) {
    selection = db[unlist(nat::nlapply(db, function(x)
      sum(pointsinsidemesh(
        shape, indexAlpha = 1, nat::xyzmatrix(x)
      )) > min_nodes))]
  } else{
    selection = db[unlist(nat::nlapply(db, function(x)
      sum(
        pointsinsidemesh(
          surf = shape,
          indexAlpha = 1,
          x = nat::xyzmatrix(x)
        )
      ) > min_nodes))]
  }
  selection
}

#' Convert an alpha shape to a mesh3D object
#'
#' @param ashape List of alphashape objects
#' @param remove.interior.points whether or not to remove points that are not part of a bounding triangle for the mesh
#' @param ... additional arguments passed to methods
#'
#' @return A mesh3D object
#' @export
#' @seealso \code{\link{combine.alphashape}} \code{\link{transform3dalphashape}}
ashape2mesh3d <- function(ashape, remove.interior.points = TRUE){
  triangles = ashape$triang[apply(ashape$triang, 1, function(x) {( any(as.numeric(x[9]) > 1))} ),][,1:3]
  if(remove.interior.points){
    if(!requireNamespace('pbapply', quietly = TRUE))
      stop("Please install suggested pbapply package")
    vertices = unique(as.vector(unique(triangles)))
    kept = 1:length(vertices)
    names(kept) = vertices
    vert = t(ashape$x)[,vertices]
    tri <- pbapply::pbapply(triangles, 1, function(x) kept[as.character(x)])
    mesh3d = rgl::tmesh3d(vertices=vert, indices = tri, homogeneous = F)
  }else{
    mesh3d = rgl::tmesh3d(t(ashape$x), t(triangles), homogeneous = F)
  }
  mesh3d
}