natverse / nat

NeuroAnatomy Toolbox: An R package for the (3D) visualisation and analysis of biological image data, especially tracings of single neurons.
https://natverse.org/nat/
64 stars 29 forks source link

Consider adding geodesic_distances function #434

Open jefferis opened 4 years ago

jefferis commented 4 years ago

See https://groups.google.com/forum/#!topic/nat-user/3oozKMYg6zY

#' Find geodesic distance along neuronal tree from the soma to every synapse
#'
#' @description Find the geodesic distance (i.e. path length along) the neuron
#'   from the the root (usually soma) to all synapses or other selected nodes.
#' @param x A \code{\link{neuron}}
#' @param root The index of the node to use as the root for the calculation.
#'   Defaults to the soma.
#' @param nodes Optional argument specifying to which nodes you want to measure
#'   distances. By default this will be all connectors in the neuron, but you
#'   can also specify all nodes, or pass in a any data frame whose first columm
#'   corresponds to (a subset) of the node identifiers of the neuron.
#' @param rval What kind of object to return
#' @return A vector of distances, a data.frame or a \code{\link{neuron}} object
#'   with an extra column in the connectors data.frame entitled
#'   \code{rootdistance}.
#' @export
#'
#' @examples
#' # nb pick first neuron out of neuron list returned by neuprint_read_neurons
#' dl4=neuprint_read_neurons('DL4')[[1]]
#' hist(geodesic_distances(dl4), br=100)
#' df=geodesic_distances(dl4, rval='dataframe')
#' # the correction factor of 125 scales to microns
#' plot(dl4/125, col='black', WithNodes=F)
#' points(xyzmatrix(df)[,c(1,2)]/125, col=topo.colors(10)[cut(df$rootdistance, 10)], pch=19)
#' # plot the soma position
#' points(xyzmatrix(dl4)[rootpoints(dl4),,drop=FALSE]/125, col='red', pch=19)
geodesic_distances <- function(x, root=rootpoints(x),
                               nodes=c("synapses", "nodes"),
                               rval=c("distance", "dataframe", 'neuron')) {
  rval=match.arg(rval)
  if(is.character(nodes)) {
    nodes=match.arg(nodes)
    nodes <- if(nodes=="synapses") x$connectors else x$d
  }
  if(is.null(nodes))
    return(NULL)
  if(!is.data.frame(nodes))
    stop("I cannot understand the nodes argument. See docs for help!")
  g=as.ngraph(x, weights=TRUE)
  # Some nodes may make multiple connections, igraph::distances won't cope
  havedups=any(duplicated(nodes[[1]]))
  tns <- if(havedups) unique(nodes[[1]]) else nodes[[1]]
  # nb converting the to nodes to character handles both the cases where they are
  # 64 bit ids (typical of catmaid neurons) and the case where they are indices
  # into the point array (typical of neuprint neurons)
  dd=igraph::distances(g, v=root, to=as.character(tns), mode='all')
  if(havedups) {
    # re-expand to match the incoming data.frame of connectors
    dd=dd[match(nodes[[1]], tns)]
  }
  if(rval=='distance')
    return(dd)
  nodes$rootdistance=dd
  if(rval=='dataframe') nodes else {x$connectors=nodes; x}
}
jefferis commented 3 years ago

@dokato I forgot to mention this issue which includes a specimen geodesic_distances function that I have been meaning to add to nat. You might want to try this out.