R/read_merge_info.R

In natverse/fafbseg: Support Functions for Analysis of FAFB EM Segmentation

#' Read CSV files containing raw agglomeration (merge) information
#'
#' @details Peter Li provided 10 CSV files with merge information in Aug 2018
#'   against which this function has been tested. This merge information has
#'   been processed to identify which groups segments should be merged together
#'   in the \code{\link[fafbsegdata]{fafbsegdata}} package, which can be
#'   accessed using the \code{\link{find_merged_segments}} function.
#'
#' @param x Input directory
#' @param pattern Optional pattern to identify CSV files
#' @param voxdims Optional voxel dimensions to convert raw pixel coordinates
#'   into physical (nm) coordinates
#' @param ... Additional arguments passed to \code{readr::\link{read_csv}}
#'
#' @return a \code{data.frame} with columns \itemize{
#'
#'   \item \code{id1,id2} Segment ids to be merged
#'
#'   \item \code{x,y,z} Location (in nm) of merge point
#'
#'   }
#' @export
#' @seealso \code{\link{find_merged_segments}},
#'   \code{\link[fafbsegdata]{fafbsegdata}} package
#'
#' @examples
#' \dontrun{
#' csvdir="~/projects/fafbseg/fafb14_v00c_split3xfill2x_merges/ffnreseg16nm_ms1000_md0.02_c0.6_iou0.7/"
#' mergeinfo <- read_mergeinfo(csvdir)
#'
#' # Show 3D position of merge locations for a large merge group containing a
#' # tangential cell (widefield visual interneuron)
#' tang_segs <- find_merged_segments(22139217567)
#' tang_mi <- subset(mergeinfo, id1%in% tang_segs | id2 %in%tang_segs)
#' points3d(xyzmatrix(tang_mi))
#'
#' # plot brain surface for context
#' library(elmr)
#' plot3d(FAFB)
#' }
#' @importFrom readr read_csv col_double col_integer cols
read_mergeinfo <- function(x, pattern = "10\\.csv$", voxdims = c(16, 16, 40), ...) {
  read_merge_csv <- function(x, voxdims) {
    res = readr::read_csv(
      x,
      col_names = c("id1", 'id2', 'x', 'y', 'z'),
      col_types = readr::cols(
        `id1` = readr::col_double(),
        `id2` = readr::col_double(),
        `x` = readr::col_integer(),
        `y` = readr::col_integer(),
        `z` = readr::col_integer()
      )
    )
    res$x = res$x * voxdims[1]
    res$y = res$y * voxdims[2]
    res$z = res$z * voxdims[3]
    res
  }

  ff = dir(x, pattern = pattern, full.names = TRUE)
  stopifnot(length(ff)>0)
  mergel = lapply(ff, read_merge_csv, voxdims = voxdims)
  dplyr::bind_rows(mergel)
}

#' Find all merged segment ids for a given starting google segment id
#'
#' This is a naive implementation that only depends only on the \code{mergeinfo}
#' object created by \code{\link{read_mergeinfo}}. It is retained only to verify
#' correctness. It is slow but does not require any pre-processing of the table
#' of merge information.
#'
#' @param x a segment id or any other input that can be interpreted by
#'   \code{\link{ngl_segments}}
#' @param mergeinfo The merge information data.frame created by
#'   \code{\link{read_mergeinfo}}
#' @return vector of segment ids
#' @seealso \code{\link{find_merged_segments}}
find_merged_segments_slow <- function(x, mergeinfo) {
  x=ngl_segments(x)
  while(TRUE) {
    cat('.')
    w1=which(mergeinfo$id1%in%x)
    w2=which(mergeinfo$id2%in%x)
    w=union(w1,w2)
    if(length(w) < 1) break
    segs=union(mergeinfo$id1[w], mergeinfo$id2[w])
    newsegs=setdiff(segs, x)
    if(length(newsegs)==0) break
    x=union(x, newsegs)
  }
  x
}


#' Graph representation of the agglomeration (merge) between raw segments
#'
#' @description \code{make_merge_graph} generates a graph representation of all
#'   agglomeration merges.
#'
#' @param x The \code{mergeinfo} \code{data.frame} generated by
#'   \code{\link{read_mergeinfo}}
#' @param g Either a graph generated by \code{make_merge_graph} or the
#'   \code{mergeinfo} \code{data.frame} generated by
#'   \code{\link{read_mergeinfo}}.
#' @param n For testing purposes, the number of rows of merge information to
#'   process - the default processes all rows.
#'
#' @return \code{make_merge_graph} returns an \code{igraph} object,
#'   \code{merge_graph_components} returns a list with elements describing the
#'   mapping between segment ids and merge groups.
#' @export
#'
#' @examples
#' \dontrun{
#' g=make_merge_graph(mergeinfo)
#' gc=merge_graph_components(g)
#' str(gc)
#' }
make_merge_graph <- function(x, n=NA) {
  el=x[,1:2]
  if(!is.na(n))
    el=el[seq_len(n),,drop=FALSE]
  el=data.matrix(el)
  vertexnames=sort(unique(c(el[,1], el[,2])))
  rawel=match(t(el), vertexnames)
  if(any(is.na(rawel))) stop("Bad vertices!")
  g=igraph::graph(rawel, n=length(vertexnames), directed=FALSE)
  igraph::V(g)$id=vertexnames
  g
}

#' @rdname make_merge_graph
#' @description \code{merge_graph_components} uses the full merge graph to find
#'   each connected subgraph corresponding to isolated merge groups.
merge_graph_components <- function(g, n=NA) {
  if(!igraph::is.igraph(g)) g=make_merge_graph(g, n=n)
  gc=igraph::components(g)
  gc$id=igraph::vertex_attr(g, "id")
  gc$membership=as.integer(gc$membership)
  gc
}