R/minc_graphs.R

In Mouse-Imaging-Centre/RMINC: Statistical Tools for Medical Imaging NetCDF (MINC) Files

#' Object to Graph
#' 
#' Convert a bic_obj object to a graph for handy graph algorithms
#' 
#' @param bic_obj The \code{bic_obj} object, generated by \link{read_obj} to be converted
#' @return an \code{igraph} object representing the object's mesh
#' @export
obj_to_graph <-
  function(bic_obj){
    edges <-
      with(bic_obj, {
        cbind(triangle_matrix[1:2,]
              , triangle_matrix[2:3]
              , triangle_matrix[c(1,3),])})
    
    graph <- igraph::make_undirected_graph(edges)
    igraph::vertex_attr(graph, "x") <- bic_obj$vertex_matrix[1,]
    igraph::vertex_attr(graph, "y") <- bic_obj$vertex_matrix[2,]
    igraph::vertex_attr(graph, "z") <- bic_obj$vertex_matrix[3,]
    igraph::vertex_attr(graph, "name") <- 1:length(igraph::V(graph))
    
    graph
  }

#' Find connected components of a object graph
#' 
#' Given vertex data and a threshold, find the connected
#' component of a graph representing a \code{bic_obj}
#' 
#' @param obj_graph The \code{igraph} object generated by \link{obj_to_graph}
#' @param data_map A set of data corresponding to the vertices
#' @param data_range A two-component vector giving lower and upper bounds
#' for vertices to include when generating connected components
#' @return a data frame of vertices after thresholding indicating the
#' vertex id, the cluster it's a member of, and the size of its cluster
#' @export
connected_components <- 
  function(obj_graph, data_map, data_range){
    subgraph <- 
      igraph::induced_subgraph(obj_graph
                               , igraph::V(obj_graph)[between(data_map
                                                              , data_range[1]
                                                              , data_range[2])])
    
    con_coms <- 
      igraph::components(subgraph)
    
    con_com_frame <-
      with(con_coms, {
        data_frame(vertex = as.numeric(names(membership))
                   , clust_id = membership) %>%
          inner_join(data_frame(clust_size = csize
                                , clust_id = 1:length(csize))
                     , by = "clust_id")
      })
    
    con_com_frame
  }

#' Components to map
#' 
#' Convert a \link{connected_components} data.frame into a vertex map
#' for visualizing the results.
#' 
#' @param graph the parent \code{igraph} generated by \code{obj_to_graph}
#' and used with \code{connected_components}
#' @param components the data.frame results of \code{connected_components}
#' @param type whether to return the size of the component or its label
#' @return A vector corresponding to each vertex in \code{graph} with either
#' zero, or the size/label from the connected components calculation
#' @export 
components_to_map <- 
  function(graph, components, type = c("size", "label")){
    type <- match.arg(type)
    col <- switch(type, 
                  size = "clust_size",
                  label = "clust_id")
    
    data_map <- rep(0, length(igraph::V(graph)))
    data_map[components$vertex] <- unlist(components[,col])
    
    data_map
  }

#' Vertex find peaks
#' 
#' Find extrema of set of values given a connectivity graph
#' @param data_map A vector or text file of values to search for peaks
#' @param graph An igraph object describing the connectivity, in most cases this will be
#' a mesh graph produced by \link{obj_to_graph}
#' @param mindist the search radius (in number of neighbours) to consider when finding a peak. A
#' node is considered a peak if it is more extreme than any element in it's neighbourhood up to
#' \code{mindist} edges away
#' @param direction Either positive, negative, or both indicating whether to consider maxima,
#' minima or both.
#' @param threshold threshold above or below which to consider peaks. In the both direction case
#' a two element vector may be passed indicating positive and negative thresholds respectively. 
#' If a single element is passed in the both direction case, it is treated as \code{c(threshold, -threshold)}.
#' Defaults to 0.
#' @param output Either mask or indices, indicating whether a logical mask or a vector of indices is desired.
#' @details Note that if coordinates of peaks are desired, they can be accessed from the parent \code{bic_obj}
#' by \code{obj$vertex_matrix[,peaks]} 
#' @return Either a logical vector indicating peaks or a vector of peak indices determined by \code{output}
#' @export
vertexFindPeaks <- function(data_map, graph
                            , mindist = 1
                            , direction = c("both", "positive", "negative")
                            , threshold = 0
                            , output = c("mask", "indices")){

  direction <- match.arg(direction)
  output <- match.arg(output)
  
  if(length(threshold) == 1 && direction == "both"){
    threshold <- c(threshold, -threshold)
  }
  
  find_graph_extrema <- function(data_map, comparator, thresh){
    vapply(seq_along(data_map), function(i){
      cur_val <- data_map[i]
      
      if(!comparator(cur_val, thresh))
        return(FALSE)
      
      nebs <- unlist(igraph::ego(graph, mindist, i))
      nebs <- nebs[nebs != i]
      nebs_vals <- data_map[nebs]
      
      all(comparator(cur_val, nebs_vals), na.rm = FALSE)
    }, FUN.VALUE = logical(1))
  }
  
  res <-
    switch(direction
           , "positive" = find_graph_extrema(data_map, `>`, threshold)
           , "negative" = find_graph_extrema(data_map, `<`, threshold)
           , "both" = find_graph_extrema(data_map, `>`, threshold[1]) | find_graph_extrema(data_map, `<`, threshold[2]))
  
  if(output == "indices") res <- which(res)
  
  res
}