R/sketchPhylo3D.R

In treeDbalance: Computation of 3D Tree Imbalance

#' Sketch a phylo3D object
#'
#' \code{sketchPhylo3D} - Sketches a phylo3D object using functions of the
#' package 'rgl'. This is a much faster but rougher version of
#' \code{plotPhylo3D} as the resulting sketch is using lines instead of
#' cylinders and the user can specify in various ways that only a portion of
#' the edges will be depicted. This is intended to get a quick overview of a
#' phylo3D object that consists of thousands of edges and would therefore take
#' very long to plot regularly.
#'
#' @author Sophie Kersting
#'
#' @param tree A rooted tree in phylo3D format (no special node enumeration
#' required, except that nodes are numbered from 1 to |V| = the total number of
#' nodes). There must be at least 2 nodes, i.e. one edge. The attributes
#' 'node.coord' and 'edge.weight' are strictly required. \cr
#' Optional: Add attribute edge.color' (character vector, e.g.
#' c("red","blue", "red", "gray30",...)) to the pylo3D object to sketch
#' the edges differently. If it is not given, every edge is
#' depicted as "gray30" by default. \cr
#' The attribute edge.diam' can be added optionally to  set the edge diameters,
#' otherwise the edge diameters will be calculated based on the edge lengths
#' and the edge weights (treated as volume).
#' @param show_node_enum A boolean value (default FALSE). If true each node
#' of the visualized phylo3D object is marked with its number. This helps to
#' identify specific nodes and edges.
#' @param ratio_of_edges A numeric value (>0 and <=1, default 1) that
#' indicates the ratio of edges that should be visualized (beginning with
#' the widest edges). This can be ignored, i.e. set to NULL, if
#' 'number_of_edges' or 'min_edge_diam' is set. If the parameters are
#' conflicting: 'ratio_of_edges' trumps 'number_of_edges' trumps
#' 'min_edge_diam'.
#' @param number_of_edges An integer value (>0, default NULL) indicating how
#' many of the edges should be visualized (beginning with the widest edges).
#' This can be ignored, i.e. left at NULL, if 'ratio_of_edges' or
#' 'min_edge_diam' is set. If the parameters are conflicting: 'ratio_of_edges'
#' trumps 'number_of_edges' trumps 'min_edge_diam'.
#' @param min_edge_diam A numeric value (default NULL) that
#' indicates that only edges with equal or larger diameter should be visualized.
#' This can be ignored, i.e. left at NULL, if
#' 'number_of_edges' is set. If the parameters are conflicting: 'ratio_of_edges'
#' trumps 'number_of_edges' trumps 'min_edge_diam'.
#' @param lwd_factor A numeric vlaue (>0, default 100) to increase or decrease
#' the line width.
#'
#'
#' @export
#' @rdname sketchPhylo3D
#' @examples
#' tree <- treeDbalance::example3Dtrees$bean09
#' sketchPhylo3D(tree, show_node_enum = TRUE, ratio_of_edges = 0.8)
sketchPhylo3D <- function(tree, show_node_enum = FALSE,
                          ratio_of_edges = 1, number_of_edges = NULL,
                          min_edge_diam = NULL, lwd_factor = 100) {
  rgl::par3d(windowRect = c(50, 50, 800, 800))
  # Measure the tree to create the box for the sketch.
  square_width <- max(
    max(tree$node.coord[, 1]) - min(tree$node.coord[, 1]),
    max(tree$node.coord[, 2]) - min(tree$node.coord[, 2]),
    max(tree$node.coord[, 3]) - min(tree$node.coord[, 3])
  )
  my_epsilon <- square_width / 20
  square_mids <- c(
    (max(tree$node.coord[, 1]) + min(tree$node.coord[, 1])) / 2,
    (max(tree$node.coord[, 2]) + min(tree$node.coord[, 2])) / 2,
    (max(tree$node.coord[, 3]) + min(tree$node.coord[, 3])) / 2
  )
  # Construct basic sketch.
  rgl::plot3d(1, 1, 1,
    type = "n",
    xlab = "x", ylab = "y", zlab = "z",
    xlim = c(
      square_mids[1] - square_width / 2 - my_epsilon,
      square_mids[1] + square_width / 2 + my_epsilon
    ),
    ylim = c(
      square_mids[2] - square_width / 2 - my_epsilon,
      square_mids[2] + square_width / 2 + my_epsilon
    ),
    zlim = c(
      square_mids[3] - square_width / 2 - my_epsilon,
      square_mids[3] + square_width / 2 + my_epsilon
    ),
    surface = FALSE, axis.scales = FALSE, aspect = TRUE
  )
  addSketchPhylo3D(tree,
    offset = c(0, 0, 0), show_node_enum,
    ratio_of_edges, number_of_edges, min_edge_diam,
    lwd_factor
  )
}
#' Sketch a phylo3D object
#'
#' \code{addSketchPhylo3D} - This function sketches a phylo3D object without any
#' coordinate axis or adds the tree to an existing sketch
#' (e.g. for sketchPhylo3D). This is a much faster but rougher version of
#' \code{addPhylo3D} as the resulting sketch is using lines instead of
#' cylinders and the user can specify in various ways that only a portion of
#' the edges will be depicted. This is intended to get a quick overview of a
#' phylo3D object that consists of thousands of edges and would therefore take
#' very long to plot regularly.
#'
#' @param offset Numeric vector of length 3, contains 3D coordinates by which
#' the phylo object should be shifted (default = c(0,0,0), i.e. no shift).
#' @export
#' @rdname sketchPhylo3D
#' @examples
#' addSketchPhylo3D(tree,
#'   offset = c(1, 1, 0), ratio_of_edges = NULL,
#'   number_of_edges = 200
#' )
addSketchPhylo3D <- function(tree, offset = c(0, 0, 0), show_node_enum = FALSE,
                             ratio_of_edges = 1, number_of_edges = NULL,
                             min_edge_diam = NULL, lwd_factor = 100) {
  if (sum(c("node.coord") %in% attributes(tree)$names) != 1) {
    stop(paste(
      "Cannot calculate subtree centroids. The tree is missing",
      "this attribute: 'node.coord'."
    ))
  }
  rgl::par3d(windowRect = c(50, 50, 800, 800))
  tree$node.coord <- tree$node.coord + # shift all coordinates by offset
    outer(rep.int(1L, nrow(tree$node.coord)), offset)

  if (!"edge.diam" %in% attributes(tree)$names) {
    edge.len <- sapply(1:nrow(tree$edge), function(x) {
      sqrt(sum((tree$node.coord[tree$edge[x, 1], ] -
        tree$node.coord[tree$edge[x, 2], ])^2))
    })
    tree$edge.diam <- 2 * sqrt(tree$edge.weight / pi / edge.len)
    message("Edge diameters were computed based on edge length and weight.\n")
  }
  #------------------- Choose only the widest edges for visualization.
  edge_selection <- NULL
  if (is.null(ratio_of_edges)) {
    if (is.null(number_of_edges)) {
      if (!is.null(min_edge_diam)) {
        lower_bound <- min(min_edge_diam, max(tree$edge.diam))
        edge_selection <- which(tree$edge.diam >= lower_bound)
      }
    } else {
      if (number_of_edges > 0) {
        how_many <- min(round(number_of_edges), nrow(tree$edge))
        edge_selection <- order(tree$edge.diam, decreasing = T)[1:how_many]
      }
    }
  } else {
    if (ratio_of_edges > 0 && ratio_of_edges <= 1) {
      how_many <- max(1, round(ratio_of_edges * nrow(tree$edge)), na.rm = TRUE)
      edge_selection <- order(tree$edge.diam, decreasing = T)[1:how_many]
    }
  }
  if (is.null(edge_selection)) {
    stop(paste(
      "Set (at least one of) the parameters 'ratio_of_edges'",
      "and 'number of edges' appropriately."
    ))
  }
  #-------------------
  if ("edge.color" %in% attributes(tree)$names) {
    ecol <- tree$edge.color
  } else {
    ecol <- rep("gray30", nrow(tree$edge))
  }
  edges_to_draw <- which(tree$edge.diam > 0)
  edges_to_draw <- intersect(edges_to_draw, edge_selection)
  rgl::segments3d(tree$node.coord[as.vector(t(tree$edge[edges_to_draw, 1:2])), ],
    col = rep(ecol[edges_to_draw], each = 2),
    lwd = max(tree$edge.diam[edges_to_draw]) / 2 * lwd_factor
  )
  if (show_node_enum) {
    nodes_to_show <- sort(unique(as.vector(tree$edge[edges_to_draw, ])))
    for (v in nodes_to_show) {
      rgl::text3d(
        x = tree$node.coord[v, 1],
        y = tree$node.coord[v, 2],
        z = tree$node.coord[v, 3], texts = paste(v), pos = 5, offset = 1
      )
    }
  }
}