R/antsApplyTransformsToPoints.R

In stnava/ANTsR: ANTs in R: Quantification Tools for Biomedical Images

#' Apply transforms to points.
#'
#' Apply a transform list to map a pointset from one domain to another. In
#' registration, one computes mappings between pairs of domains.
#' These transforms are often a sequence of increasingly complex maps, e.g.
#' from translation, to rigid, to affine to deformation.  The list of such
#' transforms is passed to this function to interpolate one image domain
#' into the next image domain, as below.  The order matters strongly and the
#' user is advised to familiarize with the standards established in examples.
#' Importantly, point mapping goes the opposite direction of image mapping,
#' for both reasons of convention and engineering.
#'
#' @param dim dimensionality of the transformation.
#' @param points moving point set with n-points in rows of at least dim
#' columns - we maintain extra information in additional columns. this
#' may be either a dataframe or a 2D antsImage - the latter may be better
#' for large pointsets.
#' @param transformlist character vector of transforms generated by
#' antsRegistration where each transform is a filename.
#' @param whichtoinvert vector of booleans, same length as transforms
#' @return same type as input. 1 -- Failure
#' @author Avants BB
#' @examples
#'
#' fixed <- antsImageRead(getANTsRData("r16"), 2)
#' moving <- antsImageRead(getANTsRData("r64"), 2)
#' fixed <- resampleImage(fixed, c(64, 64), 1, 0)
#' moving <- resampleImage(moving, c(64, 64), 1, 0)
#' mytx <- antsRegistration(
#'   fixed = fixed, moving = moving,
#'   typeofTransform = c("SyN"), verbose = FALSE
#' )
#' pts <- data.frame(
#'   x = c(110.5, 120, 130), y = c(108.1, 121.0, 130),
#'   label = c(1, 2, 3)
#' )
#' wpts <- antsApplyTransformsToPoints(
#'   dim = 2, points = pts,
#'   transformlist = mytx$fwdtransforms
#' )
#' wptsi <- antsApplyTransformsToPoints(
#'   dim = 2, points = wpts,
#'   transformlist = mytx$invtransforms
#' ) # close to pts
#'
#' \dontrun{
#' fixed <- antsImageRead(getANTsRData("r16"), 2)
#' moving <- antsImageRead(getANTsRData("r64"), 2)
#' fpts <- kmeansSegmentation(fixed, 3)$segmentation %>%
#'   thresholdImage(1, 1) %>%
#'   labelClusters(5) %>%
#'   getCentroids(5)
#' wpts <- antsApplyTransformsToPoints(
#'   dim = 2, points = fpts,
#'   transformlist = mytx$fwdtransforms
#' )
#' labimgf <- fixed * 0
#' labimgm <- moving * 0
#' for (p in 1:nrow(wpts))
#' {
#'   pt <- as.numeric(wpts[p, 1:2])
#'   idx <- round(antsTransformPhysicalPointToIndex(moving, pt))
#'   labimgm[idx[1], idx[2]] <- p
#'   pt <- as.numeric(fpts[p, 1:2])
#'   idx <- round(antsTransformPhysicalPointToIndex(fixed, pt))
#'   labimgf[idx[1], idx[2]] <- p
#' }
#' plot(fixed, labimgf %>% iMath("GD", 2))
#' plot(moving, labimgm %>% iMath("GD", 2))
#' }
#'
#' @seealso \code{\link{antsRegistration}}
#' @export antsApplyTransformsToPoints
antsApplyTransformsToPoints <- function(
    dim,
    points,
    transformlist = "",
    whichtoinvert = NA) {
  ttexists <- TRUE
  for (i in 1:length(transformlist))
  {
    if (!file.exists(transformlist[i])) ttexists <- FALSE
  }
  if (ttexists) {
    mytx <- list()
    if (all(is.na(whichtoinvert))) {
      whichtoinvert <- rep(FALSE, length(transformlist))
    }
    for (i in c(1:length(transformlist)))
    {
      ismat <- FALSE
      if ((i == 1 & length(transformlist) > 1) |
        whichtoinvert[i] == TRUE) {
        if (length(grep(".mat", transformlist[i])) == 1) {
          ismat <- TRUE
        }
      }
      if (!ismat) {
        mytx <- list(mytx, "-t", transformlist[i])
      } else if (ismat) {
        mytx <- list(
          mytx, "-t",
          paste("[", transformlist[i], ",1]", sep = "")
        )
      }
    }
    if (!inherits(points, "antsImage")) {
      usepts <- as.antsImage(data.matrix(points))
    } else {
      usepts <- antsImageClone(points)
    }
    if (usepts@dimension != 2) stop("must be 2d antsImage")
    pointsout <- antsImageClone(usepts)
    args <- list(d = dim, i = usepts, o = pointsout, unlist(mytx))
    myargs <- .int_antsProcessArguments(c(args))
    for (jj in c(1:length(myargs))) {
      if (!is.na(myargs[jj])) {
        if (myargs[jj] == "-") {
          myargs2 <- rep(NA, (length(myargs) - 1))
          myargs2[1:(jj - 1)] <- myargs[1:(jj - 1)]
          myargs2[jj:(length(myargs) - 1)] <- myargs[(jj + 1):(length(myargs))]
          myargs <- myargs2
        }
      }
    }
    ANTsRCore::antsApplyTransformsToPoints(c(myargs, "-f", 1, "--precision", 0))

    if (inherits(points, "antsImage")) {
      return(pointsout)
    }
    pointsout <- data.frame(as.matrix(pointsout))
    colnames(pointsout) <- colnames(points)
    if (ncol(pointsout) > dim) {
      pointsout[, (dim + 1):ncol(points)] <- points[, (dim + 1):ncol(points)]
    }
    return(pointsout)
  }
  if (!ttexists) {
    stop("transforms may not exist")
  }
}