R/antsrSurf.R

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

#' Render a function onto a surface.
#'
#' Use a system call to \code{antsSurf} in order to render a function onto
#' a surface.  This requires \code{antsSurf} to be in the environment as well
#' as \code{ConvertScalarImageToRGB}, both available in Advanced Normalization
#' Tools software on github.
#'
#' @param x input antsImage defining the surface on which to render
#' @param y input antsImage list defining the function to render
#' on the surface. these image(s) should be in the same space as \code{x}.
#' @param z input antsImage list mask for each \code{y} function to render
#' on the surface. these image(s) should be in the same space as \code{y}.
#' @param quantlimits lower and upper quantile limits for overlay
#' @param colormap character, one of: grey, red, green, blue, copper, jet, hsv,
#' spring, summer, autumn, winter, hot, cool, overunder, custom
#' @param alpha transparency vector for underlay and each overlay, default zero
#' @param inflationFactor number of inflation iterations to run
#' @param smoothingSigma gaussian smooth the overlay by this sigma
#' @param rotationParams 3 Rotation angles expressed in degrees or a matrix of
#' rotation parameters that will be applied in sequence.
#' @param overlayLimits absolute lower and upper limits for functional overlay.
#' this parameter will override \code{quantlimits}.  Currently, this will set
#' levels above \code{overlayLimits[2]} to \code{overlayLimits[2]}. Can be a
#' list of length of y.
#' @param backgroundColor a three entry vector where each entry is in the range of 0 to 255 and corresponds to red green and blue color channels. Alternatively, pass a string white or black.
#' @param filename prefix filename for output pngs
#' @param antspath pass the ANTSPATH here otherwise we try to detect it from environment
#' @param verbose prints the command used to call \code{antsSurf}
#' @return no output
#' @author Avants BB, Tustison NJ
#' @examples
#' \dontrun{
#'
#' myantspath <- Sys.getenv("ANTSPATH")
#' surf_ex <- file.path(myantspath, "antsSurf")
#' if (file.exists(surf_ex)) {
#'   ch2i <- antsImageRead(getANTsRData("ch2"))
#'   ch2seg <- thresholdImage(ch2i, "Otsu", 3)
#'   wm <- thresholdImage(ch2seg, 3, 3)
#'   wm2 <- smoothImage(wm, 1) %>% thresholdImage(0.5, Inf)
#'   kimg <- weingartenImageCurvature(ch2i, 1.5) %>% smoothImage(1)
#'   rp1 <- matrix(c(90, 180, 90), ncol = 3)
#'   rp2 <- matrix(c(90, 180, 270), ncol = 3)
#'   rp3 <- matrix(c(90, 180, 180), ncol = 3)
#'   rp <- rbind(rp1, rp3, rp2)
#'   antsrSurf(
#'     x = wm2, y = list(kimg), z = list(wm2 %>% iMath("MD", 3)),
#'     inflationFactor = 255, overlayLimits = c(-0.3, 0.3), verbose = TRUE,
#'     rotationParams = rp, filename = tempfile()
#'   )
#'
#'   fn <- getANTsRData("surf")
#'   img <- antsImageRead(fn) # see antsSurf on github for data
#'   wm <- thresholdImage(img, 3, 4)
#'   wm <- thresholdImage(img, 3, 4) %>% iMath("FillHoles")
#'   wm <- wm %>% thresholdImage(1, 2)
#'   wms <- smoothImage(wm, 1.0)
#'   wmk <- weingartenImageCurvature(wms, 1.5, 0)
#'   # will display to screen
#'   antsrSurf(
#'     x = wm, y = list(wmk %>% smoothImage(1)), z = list(wm %>% iMath("MD", 1)),
#'     rotationParams = c(270, 0, 90), overlayLimits = c(-0.4, 0.4)
#'   )
#'   # surface and two overlays
#'   blob <- antsImageRead(getANTsRData("blob"))
#'   blob[1:266, 1:266, 1:100] <- 0
#'   z <- list(wm %>% iMath("MD", 1), blob %>% smoothImage(1))
#'   antsrSurf(
#'     x = wm, y = list(wmk %>% smoothImage(1), blob), z = z,
#'     colormap = c("jet", "blue"), alpha = c(1, 0.5, 1),
#'     rotationParams = c(270, 0, 90),
#'     overlayLimits = list(c(-0.4, 0.4), c(0.9, 1.001))
#'   )
#'   # separate pos and neg curvature
#'   y <- list(
#'     thresholdImage(wmk, 0.00, Inf),
#'     thresholdImage(wmk, -100, -0.00)
#'   )
#'   z <- list(
#'     y[[1]] %>% iMath("MD", 1),
#'     y[[2]] %>% iMath("MD", 1)
#'   )
#'   antsrSurf(
#'     x = wm, y = y, z = z, smoothingSigma = 0.5, alpha = c(1, 1, 1),
#'     colormap = c("red", "blue"),
#'     inflationFactor = 155, overlayLimits = list(c(0.5, 1.0001), c(0.5, 1.0001)),
#'     verbose = TRUE, rotationParams = rp[1, ]
#'   )
#' }
#' }
#'
#' @export antsrSurf
antsrSurf <- function(x, y, z,
                      quantlimits = c(0.1, 0.9),
                      colormap = "jet",
                      alpha = NA,
                      inflationFactor = 25,
                      smoothingSigma = 0.0,
                      rotationParams = c(270, 0, 270),
                      overlayLimits = NA,
                      backgroundColor,
                      filename = NA,
                      antspath = NA,
                      verbose = FALSE) {
  domainImageMap <- NULL
  if (any(is.na(alpha))) {
    alpha <- rep(1, length(x) + length(y))
  }
  if (length(z) != length(y)) stop("each y must have a mask in z")
  if (is.numeric(overlayLimits)) {
    overlayLimits <- list(overlayLimits)
  }
  # #' @param domainImageMap resamples surf and func to this domain FIXME
  # check for needed programs
  # first get antspath
  if (is.na(antspath)) {
    myantspath <- Sys.getenv(c("ANTSPATH"))
  } else {
    myantspath <- antspath
  }
  if (nchar(myantspath) == 0 & is.na(antspath)) {
    stop("Please set ANTSPATH in your environment")
  }
  if (!is.na(antspath)) {
    myantspath <- antspath
  }
  ass <- paste(myantspath, "antsSurf", sep = "/")
  cvt <- paste(myantspath, "ConvertScalarImageToRGB", sep = "/")
  if (!file.exists(ass)) {
    stop("This function needs antsSurf in the environment.")
  }
  if (!file.exists(cvt)) {
    stop("This function needs ConvertScalarImageToRGB in the environment.")
  }

  if (!is.null(domainImageMap)) {
    if (is.antsImage(domainImageMap)) {
      tx <- new("antsrTransform",
        precision = "float",
        type = "AffineTransform", dimension = x@dimension
      )
      x <- applyAntsrTransformToImage(tx, x, domainImageMap)
      if (!missing("y")) {
        if (is.antsImage(y)) y <- list(y)
        for (i in 1:length(y)) {
          y[[i]] <-
            applyAntsrTransformToImage(tx, y[[i]], domainImageMap,
              interpolation = "Linear"
            )
        }
      }
    }
  }
  xfn <- tempfile(fileext = ".nii.gz")
  antsImageWrite(x, xfn)
  if (is.na(filename)) filename <- tempfile()
  if (!is.matrix(rotationParams)) {
    rotationParams <- matrix(rotationParams, ncol = 3)
  }
  if (nrow(rotationParams) > 1) {
    if (!usePkg("abind") | !usePkg("png") | !usePkg("grid")) {
      print("Need abind, grid and png")
      return(NULL)
    }
  }
  pngs <- rep(NA, nrow(rotationParams))
  if (missing(backgroundColor)) {
    backgroundColor <- "white"
  }
  if (any(backgroundColor == "black")) {
    backgroundColor <- paste("0x0x0x", alpha[1], sep = "")
  } else if (any(backgroundColor == "white")) {
    backgroundColor <- paste("255x255x255x", alpha[1], sep = "")
  } else {
    backgroundColor <-
      paste(c(backgroundColor, alpha[1]), collapse = "x")
  }

  backgroundColorS <- "255x255x255x1"
  for (myrot in 1:nrow(rotationParams))
  {
    asscmd <- paste(ass, " -s [ ", xfn, ",", backgroundColorS, "] ")
    if (!missing(y)) {
      ct <- 0
      if (length(colormap) != length(y)) {
        colormap <- rep(colormap, length.out = length(y))
      }
      for (overlay in y)
      {
        ct <- ct + 1
        wms <- smoothImage(overlay, smoothingSigma)
        myquants <- quantile(wms[abs(wms) > 0], quantlimits)
        if (!all(is.na(overlayLimits))) {
          myquants <- overlayLimits[[ct]]
          if (all(myquants < 0)) {
            #      wms = wms * -1.0
            #      myquants = rev( myquants ) * ( -1.0 )
            #      wms = wms * thresholdImage( wms, myquants[1] , myquants[2] )
          }
          if (all(myquants > 0)) {
            #      wms = wms * thresholdImage( wms, myquants[1] , myquants[2] )
          }
        } else {
          #      wms[ wms <  quantlimits[1] ] = quantlimits[1]
          #      wms[ wms >  quantlimits[2] ] = quantlimits[2]
        }
        if (verbose) {
          print(paste("overlay quantiles for overlay", ct))
          print(myquants)
          print("range")
          print(range(wms))
        }
        kblobfn <- tempfile(fileext = ".nii.gz")
        antsImageWrite(z[[ct]], kblobfn)
        overlayfn <- tempfile(fileext = ".nii.gz")
        antsImageWrite(wms, overlayfn)
        csvlutfn <- tempfile(fileext = ".csv")
        overlayrgbfn <- tempfile(fileext = ".nii.gz")
        if (verbose) print(colormap[ct])
        cvtcmd <- paste(
          cvt, " 3 ", overlayfn, overlayrgbfn,
          kblobfn, colormap[ct], " none ", myquants[1], myquants[2], " 0 255", csvlutfn
        )
        sss <- system(cvtcmd)
        if (verbose) print(cvtcmd)
        if (verbose) cat("\n")
        alphaloc <- alpha[min(c(ct + 1, length(alpha)))]
        if (verbose) print(paste("alpha", alphaloc))
        asscmd <- paste(asscmd, "-f [ ", overlayrgbfn, ", ", kblobfn, ", ", alphaloc, " ] ")
      }
    }
    if (nrow(rotationParams) == 1) {
      asscmd <- paste(
        asscmd, " -i ", inflationFactor,
        " -d [", paste(rotationParams[myrot, ], collapse = "x"), ",", backgroundColor, "]"
      )
    } else {
      pngext <- myrot
      if (myrot < 10) pngext <- paste("0", pngext, sep = "")
      if (myrot < 100) pngext <- paste("0", pngext, sep = "")
      pngfnloc <- paste(filename, pngext, ".png", sep = "")
      system(paste("rm", pngfnloc))
      asscmd <- paste(asscmd, " -i ", inflationFactor, " -d ", pngfnloc,
        "[", paste(rotationParams[myrot, ], collapse = "x"), ",", backgroundColor, "] ",
        sep = ""
      )
    }
    if (verbose) print(asscmd)
    sss <- system(asscmd)
    Sys.sleep(3)
    if (nrow(rotationParams) > 1) pngs[myrot] <- pngfnloc
  }
  Sys.sleep(3)
  if (nrow(rotationParams) > 1) {
    mypng <- png::readPNG(pngs[1])
    for (i in 2:length(pngs))
    {
      mypng <- abind::abind(mypng, png::readPNG(pngs[i]), along = 2)
    }
    png(paste(filename, ".png", sep = ""), width = dim(mypng)[2], height = dim(mypng)[1])
    grid::grid.raster(mypng)
    dev.off()
  }
}









#' Render an image volume with or without overlay.
#'
#' Use a system call to \code{antsVol} in order to render an image.  This
#' requires \code{antsVol} to be in the environment as well as
#' \code{ConvertScalarImageToRGB}, both available in Advanced Normalization
#' Tools software on github.
#'
#' @param x input antsImage defining the image to render
#' @param y input antsImage list defining the function to render
#' on the image. these image(s) should be in the same space.
#' @param quantlimits lower and upper quantile limits for overlay
#' @param colormap character, one of: grey, red, green, blue, copper, jet, hsv,
#' spring, summer, autumn, winter, hot, cool, overunder, custom
#' @param rotationParams 3 Rotation angles expressed in degrees or a matrix of
#' rotation parameters that will be applied in sequence.
#' @param overlayLimits absolute lower and upper limits for functional overlay.
#' this parameter will override \code{quantlimits}.  Currently, this will set
#' levels above \code{overlayLimits[2]} to \code{overlayLimits[2]}.
#' @param magnificationFactor zooms in on image during rendering
#' @param intensityTruncation lower and upper quantiles at which to truncate intensity
#' @param filename prefix filename for output pngs
#' @param antspath pass the ANTSPATH here otherwise we try to detect it from environment
#' @param verbose prints the command used to call \code{antsVol}
#' @return no output
#' @author Avants BB, Tustison NJ
#' @examples
#' \dontrun{
#' ch2i <- antsImageRead(getANTsRData("mni"))
#' ch2seg <- thresholdImage(ch2i, "Otsu", 3)
#' wm <- thresholdImage(ch2seg, 3, 3)
#' kimg <- weingartenImageCurvature(ch2i, 1.5) %>% smoothImage(1)
#' ap <- path.expand("~/code/ants-src/bin/")
#' rp0 <- matrix(c(90, 180, 90), ncol = 3)
#' rp2 <- matrix(c(0, 0, 0), ncol = 3)
#' rp3 <- matrix(c(270, 90, 90), ncol = 3)
#' rp <- rbind(rp0, rp2, rp3) # pass these as rotationParams
#' antsrVol(wm, list(kimg), quantlimits = c(0.01, 0.99))
#' }
#'
#' @export antsrVol
antsrVol <- function(x, y,
                     quantlimits = c(0.1, 0.9),
                     colormap = "jet",
                     rotationParams = c(270, 0, 270),
                     overlayLimits = NA,
                     magnificationFactor = 1.0,
                     intensityTruncation = c(0.0, 1.0),
                     filename = NA,
                     antspath = NA,
                     verbose = FALSE) {
  # check for needed programs
  # first get antspath
  myantspath <- Sys.getenv(c("ANTSPATH"))
  if (nchar(myantspath) == 0 & is.na(antspath)) {
    stop("Please set ANTSPATH in your environment")
  }
  if (!is.na(antspath)) myantspath <- antspath
  ass <- file.exists(paste(myantspath, "antsVol", sep = "/"))
  cvt <- file.exists(paste(myantspath, "ConvertScalarImageToRGB", sep = "/"))
  if (length(ass) == 0) stop("This function needs antsSurf in the environment.")
  if (length(cvt) == 0) stop("This function needs ConvertScalarImageToRGB in the environment.")
  ass <- paste(myantspath, "antsVol", sep = "/")
  cvt <- paste(myantspath, "ConvertScalarImageToRGB", sep = "/")
  xfn <- tempfile(fileext = ".nii.gz")
  xmod <- antsImageClone(x)
  if (intensityTruncation[1] > 0 | intensityTruncation[2] < 1) {
    xmod <- iMath(x, "TruncateIntensity", intensityTruncation[1], intensityTruncation[2])
  }
  antsImageWrite(xmod, xfn)
  if (is.na(filename)) filename <- tempfile()
  if (!is.matrix(rotationParams)) {
    rotationParams <- matrix(rotationParams, ncol = 3)
  }
  if (nrow(rotationParams) > 1) {
    if (!usePkg("abind") | !usePkg("png") | !usePkg("grid")) {
      print("Need abind, grid and png")
      return(NULL)
    }
  }
  pngs <- rep(NA, nrow(rotationParams))
  for (myrot in 1:nrow(rotationParams))
  {
    assvcmd <- paste(ass, " -i  ", xfn, " ", sep = "")
    if (!missing(y)) {
      ct <- 0
      if (length(colormap) != length(y)) {
        colormap <- rep(colormap, length.out = length(y))
      }
      for (overlay in y)
      {
        ct <- ct + 1
        wms <- smoothImage(overlay, 1.0)
        myquants <- quantile(overlay[abs(overlay) > 0], quantlimits)
        if (!all(is.na(overlayLimits))) {
          myquants <- overlayLimits
          overlay[overlay < myquants[1]] <- 0
          overlay[overlay > myquants[2]] <- myquants[2]
          if (verbose) print(myquants)
        }
        kblob <- thresholdImage(wms, myquants[1], Inf)
        kblobfn <- tempfile(fileext = ".nii.gz")
        antsImageWrite(kblob, kblobfn)
        overlayfn <- tempfile(fileext = ".nii.gz")
        antsImageWrite(overlay, overlayfn)
        csvlutfn <- tempfile(fileext = ".csv")
        overlayrgbfn <- tempfile(fileext = ".nii.gz")
        if (verbose) print(colormap[ct])
        cvtcmd <- paste(
          cvt, " 3 ", overlayfn, overlayrgbfn,
          kblobfn, colormap[ct], " none ", myquants[1], myquants[2], " 0 255", csvlutfn
        )
        sss <- system(cvtcmd)
        if (verbose) print(cvtcmd)
        if (verbose) cat("\n")
        assvcmd <- paste(assvcmd, "-f [ ", overlayrgbfn, ", ", kblobfn, " ] ")
      }
    }
    if (is.na(filename)) {
      assvcmd <- paste(
        assvcmd,
        " -d [", magnificationFactor, ",", paste(rotationParams[myrot, ], collapse = "x"), ",255x255x255] "
      )
    } else {
      pngext <- myrot
      if (myrot < 10) pngext <- paste("0", pngext, sep = "")
      if (myrot < 100) pngext <- paste("0", pngext, sep = "")
      pngfnloc <- paste(filename, pngext, ".png", sep = "")
      system(paste("rm", pngfnloc))
      assvcmd <- paste(assvcmd, " -d ", pngfnloc,
        "[", magnificationFactor, ",",
        paste(rotationParams[myrot, ], collapse = "x"), ",255x255x255] ",
        sep = ""
      )
    }
    if (verbose) print(assvcmd)
    sss <- system(assvcmd)
    if (nrow(rotationParams) > 1) pngs[myrot] <- pngfnloc
  }
  if (nrow(rotationParams) > 1) {
    mypng <- png::readPNG(pngs[1])
    for (i in 2:length(pngs))
    {
      mypng <- abind::abind(mypng, png::readPNG(pngs[i]), along = 2)
    }
    png(paste(filename, ".png", sep = ""), width = dim(mypng)[2], height = dim(mypng)[1])
    grid::grid.raster(mypng)
    dev.off()
  }
}