R/QSM2Phylo.R

In treeDbalance: Computation of 3D Tree Imbalance

#' Convert QSM format to a phylo3D object
#'
#' \code{qsm2phylo3D} - Wrapper function that converts the output .mat-file of
#' the software treeQSM (version 2.0 or 2.4.x) to a rooted tree in phylo3D
#' format.
#'
#' @author Sophie Kersting
#'
#' @param file File name (with path if not located in working directory),
#' e.g.:\cr
#' "C:/...path.../myMatlabQSMfile.mat"
#' @param version Specifies the version of the treeQSM .mat-file. Available
#' are \cr
#' - "2.0" for version 2.0 and
#' - "2.4.x" for version 2.4.1 and 2.4.2.
#' @param setConnection2zero Logical value (default TRUE) indicating if the
#' width of the edges, which connect two consecutive cylinders, should be
#' set to zero. If FALSE, then the radius is set to the radius of the
#' "child" cylinder.
#'
#' @return \code{qsm2phylo3D(file)} A rooted tree in phylo3D format.
#' @export
#' @rdname QSM2Phylo
#' @examples
#' mat_file <- system.file("extdata", "exampleQSM.mat",
#'   package = "treeDbalance"
#' )
#' new_phylo <- qsm2phylo3D(file = mat_file)
qsm2phylo3D <- function(file, version = "2.4.x", setConnection2zero = TRUE) {
  if (is.null(file)) {
    comment("Input file empty.")
    return(NULL)
  }
  if (!grepl("mat", file)) stop("File is not a .mat-file")
  matfile <- R.matlab::readMat(file)
  if (version == "2.0") {
    if (sum(c("Rad", "Len", "Sta", "Axe", "CPar", "CExt", "CiB", "BOrd") %in%
      names(matfile)) != 8) {
      warning("Input file not formatted correctly.")
      return(NULL)
    }
    cyl_branch_order <- rep(NA, nrow(matfile$Rad))
    for (i in seq(nrow(matfile$BOrd))) {
      cyl_branch_order[unlist(matfile$CiB[i, 1])] <- matfile$BOrd[i, 1]
    }
    qsm_list <- list(
      radius = matfile$Rad, length = matfile$Len,
      start = matfile$Sta, axis = matfile$Axe,
      parent = matfile$CPar, extension = matfile$CExt,
      cyl_branch_order = cyl_branch_order
    )
  } else if (version == "2.4.x") {
    ncyl <- length(matfile$QSM[[1]][[1]])
    if (ncyl < 1 || length(matfile$QSM[[1]]) < 6 ||
      sum(c(
        "radius", "length", "start", "axis",
        "parent", "extension", "BranchOrder"
      ) %in%
        rownames(matfile$QSM[[1]])) != 7) {
      warning("Input file not formatted correctly.")
      return(NULL)
    }
    qsm_list <- matfile$QSM[, , ]$cylinder[c(
      "radius", "length", "start",
      "axis", "parent", "extension", "BranchOrder"
    ), , ]
    names(qsm_list) <- gsub("BranchOrder", "cyl_branch_order", names(qsm_list))
  } else {
    stop("Unsupported version.")
  }
  tree <- qsmList2phylo3D(
    qsm_list = qsm_list,
    setConnection2zero = setConnection2zero
  )
  return(tree)
}
#' Convert QSM format to a phylo3D object
#'
#' \code{qsmList2phylo3D} - Converts a list containing QSM data to a rooted
#' tree in phylo3D format.
#'
#' @author Sophie Kersting
#'
#' @param qsm_list List containing the QSM-data, with the attributes
#' "radius", "length", "start", "axis", "parent" and "extension".
#'
#' @return \code{qsmList2phylo3D(qsmList)} A rooted tree in phylo3D format.
#' @export
#' @rdname QSM2Phylo
#' @examples
#' new_phylo <- qsmList2phylo3D(NULL)
qsmList2phylo3D <- function(qsm_list, setConnection2zero = TRUE) {
  if (is.null(qsm_list)) {
    comment("Input list empty.")
    return(NULL)
  }
  ncyl <- length(qsm_list$radius)
  if (ncyl < 1 || length(qsm_list$length) != ncyl ||
    !identical(dim(qsm_list$start), as.integer(c(ncyl, 3))) ||
    !identical(dim(qsm_list$axis), as.integer(c(ncyl, 3))) ||
    length(qsm_list$parent) != ncyl || length(qsm_list$extension) != ncyl ||
    sum(c(
      "radius", "length", "start", "axis", "parent", "extension",
      "cyl_branch_order"
    ) %in%
      names(qsm_list)) != 7) {
    warning("Input QSM list not formatted correctly.")
    return(NULL)
  }
  # Rename (extension/child == 0) to NA and compute end points of cylinders from
  # the axis and the edge length.
  qsm_list$extension[qsm_list$extension == 0] <- NA
  qsm_list$parent[qsm_list$parent == 0] <- NA
  qsm_list$end <- qsm_list$start + qsm_list$axis * matrix(qsm_list$length,
    byrow = FALSE, nrow = ncyl,
    ncol = 3
  )
  # Each cylinder can contribute two nodes since they are not precisely
  # connected. For the case that some do in fact meet in the same coordinates,
  # the duplicate nodes will later be identified and sorted out.
  # The node coordinates of the source nodes are already given, the end nodes
  # have been computed from the axis and the edge length. Here, the nodes
  # 1,...,ncyl are the starting points of the cylinders, and ncyl+1,...,2ncyl
  # the respective end points. The latter will be later replaced by the
  # starting points of the respective next cylinder if they coincide.
  ncoords <- rbind(qsm_list$start, qsm_list$end)
  # Reserve storage for twice as many edges than cylinders since there might
  # be connection edges needed between the cylinder edges. Here the edges
  # 1,...,ncyl are the given cylinders, and the edges ncyl+1,...,2ncyl
  # are the respective incoming edges -- again if needed.
  edgemat <- rbind(
    cbind(1:ncyl, (ncyl + 1):(2 * ncyl)),
    matrix(NA, nrow = ncyl, ncol = 2)
  )
  eradius <- c(qsm_list$radius, rep(NA, ncyl))
  elength <- c(qsm_list$length, rep(NA, ncyl))
  ebranch_order <- c(qsm_list$cyl_branch_order, rep(NA, ncyl))
  eoriginal_cyl <- c(rep(TRUE, ncyl), rep(FALSE, ncyl))
  edges_of_zero_weight <- which(eradius == 0 | elength == 0)
  # Now, create the phylo structure and track which nodes are leaves and which
  # nodes are actually needed (end coincides with start of next cylinder).
  is_leaf <- rep(FALSE, 2 * ncyl)
  nodes_needed <- rep(TRUE, 2 * ncyl)
  replaced_by <- 1:(2 * ncyl)
  edges_needed <- rep(TRUE, 2 * ncyl)
  for (curr_cyl in 1:ncyl) {
    # If the next node (extension=child) is NA, this cylinder is the last of
    # a chain of cylinders, i.e. ends with a leaf node.
    child <- qsm_list$extension[curr_cyl]
    if (is.na(child)) { # leaf cylinder
      is_leaf[curr_cyl + ncyl] <- TRUE # end node is leaf
    } else { # internal cylinder
      # Check if end coincides with child's start.
      if (identical(qsm_list$end[curr_cyl, ], qsm_list$start[child, ])) {
        # Is identical: Node curr_cyl+ncyl not needed and replaced by child,
        # no connection edge needed.
        edgemat[curr_cyl, ] <- c(curr_cyl, child)
        nodes_needed[curr_cyl + ncyl] <- FALSE
        replaced_by[curr_cyl + ncyl] <- child
        edges_needed[child + ncyl] <- FALSE
      } else {
        # Not identical: Node curr_cyl+ncyl needed and a connection edge from
        # that node to the child has to be added - note that this is the
        # incoming edge of the child edge.
        edgemat[child + ncyl, ] <- c(curr_cyl + ncyl, child)
        elength[child + ncyl] <- stats::dist(rbind(
          qsm_list$end[curr_cyl, ],
          qsm_list$start[child, ]
        ))
        ebranch_order[child + ncyl] <- ebranch_order[child]
        eradius[child + ncyl] <- eradius[child]
      }
    }
  }
  # The incoming connection edge of the first cylinder of each chain has to
  # be set differently:
  root <- NA
  # Check every still missing and still needed incoming connection edge.
  for (connect_edge in which(is.na(elength) & edges_needed)) {
    # Note curr_cyl+ncyl is the incoming connection edge of curr_cyl, so:
    curr_cyl <- connect_edge - ncyl
    parent_cyl <- qsm_list$parent[curr_cyl]
    if (is.na(parent_cyl)) {
      # No parent, then this is the root.
      root <- curr_cyl # since the cylinder number here equals the starting node
      edges_needed[root + ncyl] <- FALSE
    } else {
      # If there is a parent cylinder connect the current cylinder to the
      # parental cylinder's start point.
      # Check if that coincides with current cylinder's start.
      if (identical(ncoords[parent_cyl, ], qsm_list$start[curr_cyl, ])) {
        # Is identical: Node curr_cyl not needed and replaced by parental
        # connection node; no connection edge needed.
        edgemat[curr_cyl, 1] <- replaced_by[parent_cyl]
        nodes_needed[curr_cyl] <- FALSE
        replaced_by[curr_cyl] <- replaced_by[parent_cyl]
        edges_needed[connect_edge] <- FALSE
      } else {
        # Not identical: Node curr_cyl needed and a connection edge from
        # the parental cylinder's start node to curr_cyl has to be added - note
        # that this is the incoming edge of the current cylinder.
        edgemat[connect_edge, ] <- c(replaced_by[parent_cyl], curr_cyl)
        elength[connect_edge] <- stats::dist(rbind(
          qsm_list$start[curr_cyl, ],
          ncoords[parent_cyl, ]
        ))
        ebranch_order[connect_edge] <- ebranch_order[curr_cyl]
        eradius[connect_edge] <- eradius[curr_cyl]
      }
    }
  }
  # Finally, remove all edges which are not needed.
  edgemat <- edgemat[edges_needed, ]
  eradius <- eradius[edges_needed]
  elength <- elength[edges_needed]
  ebranch_order <- ebranch_order[edges_needed]
  eoriginal_cyl <- eoriginal_cyl[edges_needed]
  eradius_ConnSet <- eradius
  if (setConnection2zero) {
    eradius_ConnSet[!eoriginal_cyl] <- 0
  }
  # ------------------
  # Create the final object.
  phyl <- list(
    edge = edgemat, tip.label = rep("", sum(is_leaf)),
    node.coord = ncoords,
    edge.diam = 2 * eradius_ConnSet, edge.length = elength,
    edge.weight = eradius_ConnSet^2 * pi * elength,
    Nnode = sum(nodes_needed) - sum(is_leaf),
    edge.is_original_cyl = eoriginal_cyl,
    edge.branch_order = ebranch_order
  )
  attr(phyl, "edges_of_zero_weight") <- edges_of_zero_weight
  # ------------------
  # Compute the DBH -- Search for the cylinder around 1.3m of main stem.
  DBH <- NA
  if (max(ncoords[, 3]) - min(ncoords[, 3]) > 1.3) { # If the tree is large enough.
    DBH_candidates <- which((phyl$node.coord[phyl$edge[, 1], 3] -
      min(phyl$node.coord[, 3]) <= 1.3) &
      (phyl$node.coord[phyl$edge[, 2], 3] -
        min(phyl$node.coord[, 3]) >= 1.3) &
      phyl$edge.branch_order == 0)
    if (length(DBH_candidates) >= 1) {
      DBH <- 2 * mean(eradius[DBH_candidates])
    }
  }
  attr(phyl, "DBH") <- DBH
  # ------------------
  class(phyl) <- "phylo3D"
  phyl <- enum2_1toV(phyl) # Change the enumeration if necessary.
  return(phyl)
}