R/ants_motion_estimation.R

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

#' @rdname motion_correction
#' @title Simple motion_correction function.
#'
#' @description motion corrects 4D time series imaging data
#'
#' @param img 4D antsImage
#' @param fixed target fixed image
#' @param moreaccurate 0, 1 or 2 with higher values being more accurte
#'  (use 2 for real applications, 0 for testing)
#' @param txtype Transformation type
#' @param verbose verbosity boolean
#' @param num_threads will execute
#' \code{Sys.setenv(ITK_GLOBAL_DEFAULT_NUMBER_OF_THREADS = num_threads)} before
#' running to attempt a more reproducible result.  See
#' \url{https://github.com/ANTsX/ANTs/wiki/antsRegistration-reproducibility-issues}
#' for discussion.  If \code{NULL}, will not set anything.
#' @param seed will execute
#' \code{Sys.setenv(ANTS_RANDOM_SEED = seed)} before
#' running to attempt a more reproducible result.  See
#' \url{https://github.com/ANTsX/ANTs/wiki/antsRegistration-reproducibility-issues}
#' for discussion.  If \code{NULL}, will not set anything.
#'
#' @note This function may give different results on multiple runs.
#'
#' @return list of outputs
#' @author Avants BB
#' @examples
#' set.seed(1000)
#' testimg <- makeImage(c(10, 10, 10, 5), rnorm(5000))
#' testimg <- iMath(testimg, "PadImage", 5)
#' mocorr <- .motion_correction(testimg, num_threads = 1, seed = 10)
#' aimg_to_array <- function(x) {
#'   if (is.antsImage(x)) {
#'     x <- as.array(x)
#'   }
#'   x
#' }
#' amocorr <- lapply(mocorr, aimg_to_array)
#' mocorr2 <- .motion_correction(testimg, num_threads = 1, seed = 10)
#' amocorr2 <- lapply(mocorr2, aimg_to_array)
#'
#' # This function may give different results on multiple runs
#' # without setting the seed
#' @export
.motion_correction <- function(img,
                               fixed = NULL,
                               moreaccurate = 1,
                               txtype = "Affine",
                               verbose = FALSE,
                               num_threads = 1,
                               seed = NULL) {
  ants_random_seed <- itk_threads <- NULL
  if (!is.null(seed)) {
    ants_random_seed <- Sys.getenv("ANTS_RANDOM_SEED")
    Sys.setenv(ANTS_RANDOM_SEED = seed)
  }
  if (!is.null(num_threads)) {
    itk_threads <- Sys.getenv("ITK_GLOBAL_DEFAULT_NUMBER_OF_THREADS")
    Sys.setenv(ITK_GLOBAL_DEFAULT_NUMBER_OF_THREADS = num_threads)
  }
  on.exit({
    if (!is.null(ants_random_seed)) {
      Sys.setenv(ANTS_RANDOM_SEED = ants_random_seed)
    }
    if (!is.null(itk_threads)) {
      Sys.setenv(ITK_GLOBAL_DEFAULT_NUMBER_OF_THREADS = itk_threads)
    }
  })
  if (is.character(img)) {
    if (length(img) != 1) {
      print("'img' should be only one filename")
      return(NULL)
    }
    img <- antsImageRead(img, dimension = 4, "float")
    inpixeltype <- img@pixeltype
  } else if (inherits(img, "antsImage")) {
    inpixeltype <- img@pixeltype
    if (img@pixeltype != "float") {
      print("'img' must have pixeltype  'float' ")
      img <- antsImageClone(img, "float")
    }
    if (img@dimension != 4) {
      print("'img' must have pixeltype 'float' and dimension '4'")
      return(NULL)
    }
  } else {
    print("'img' must be a filename or an 'antsImage'")
    return(NULL)
  }

  if (is.null(fixed)) {
    fixed <- getAverageOfTimeSeries(img)
  } else {
    if (is.character(fixed)) {
      if (length(fixed) != 1) {
        print("'fixed' should be only one filename")
        return(NULL)
      }
      fixed <- antsImageRead(fixed, dimension = 3, "float")
    } else if (is.antsImage(fixed)) {
      if (fixed@pixeltype != "float") {
        print("'fixed' must have pixeltype  'float' ")
        fixed <- antsImageClone(fixed, "float")
      }
      if (fixed@dimension != 3) {
        print("'fixed' must have pixeltype 'float' and dimension '3'")
        return(NULL)
      }
    } else {
      print("'fixed' must be a filename or an 'antsImage'")
      return(NULL)
    }
  }

  n <- dim(img)[4]
  if (n > 10) {
    n <- 10
  }
  avg_img <- getAverageOfTimeSeries(img)
  moco_img <- new("antsImage", "float", 4)
  moco_params <- new("antsMatrix", "double")
  mibins <- 20
  if (moreaccurate == 3) {
    antsMotionCorr(list(
      d = 3, o = list(moco_params, moco_img, avg_img),
      m = list(name = "MI", fixed, img, 1, mibins),
      t = paste(txtype, "[0.1,3,0]", sep = ""),
      i = "100x50x20",
      u = 1, e = 1, s = "2x1x0", f = "4x2x1",
      n = n, v = as.numeric(verbose)
    ))
  }
  if (moreaccurate == 2) {
    antsMotionCorr(list(
      d = 3, o = list(moco_params, moco_img, avg_img),
      m = list(name = "MI", fixed, img, 1, mibins, "regular", 0.25),
      t = paste(txtype, "[0.1]", sep = ""),
      i = "100x50x30",
      u = 1, e = 1, s = "2x1x0", f = "4x2x1",
      n = n, v = as.numeric(verbose)
    ))
  }
  if (moreaccurate == "intraSubjectBOLD") {
    antsMotionCorr(list(
      d = 3, o = list(moco_params, moco_img, avg_img),
      m = list(name = "MI", fixed, img, 1, mibins, "regular", 0.2),
      t = paste(txtype, "[0.25]", sep = ""),
      i = "50x20",
      u = 1, e = 1, s = "1x0", f = "2x1", n = n, v = as.numeric(verbose)
    ))
  }
  if (moreaccurate == 1) {
    antsMotionCorr(list(
      d = 3, o = list(moco_params, moco_img, avg_img),
      m = list(name = "MI", fixed, img, 1, mibins, "regular", 0.25),
      t = paste(txtype, "[0.1]", sep = ""),
      i = "100",
      u = 1, e = 1, s = 0, f = 1, n = n, v = as.numeric(verbose)
    ))
  }
  if (moreaccurate == 0) {
    antsMotionCorr(list(
      d = 3, o = list(moco_params, moco_img, avg_img),
      m = list(name = "MI", fixed, img, 1, mibins, "regular", 0.02),
      t = paste(txtype, "[0.1]", sep = ""),
      i = 3, u = 1,
      e = 1, s = 0, f = 1, n = n, v = as.numeric(verbose)
    ))
  }
  moco_params <- as.data.frame(moco_params)
  mynames <- c(
    "MetricPre", "MetricPost",
    paste("MOCOparam", 1:(ncol(moco_params) - 2), sep = "")
  )
  names(moco_params) <- mynames
  return(
    list
    (
      moco_img = moco_img,
      moco_params = moco_params,
      moco_avg_img = antsImageClone(avg_img, inpixeltype)
    )
  )
}