R/camino_modelfit.R

In muschellij2/rcamino: Port of the Camino Software

#' @title Wrapper for Camino \code{modelfit} function
#' @description Performs the Camino \code{modelfit} function
#' @param infile Input 4D image
#' @param outfile Output filename for diffusion tensor
#' @param scheme file of Camino scheme with b-values and b-vectors
#' @param mask Provides the name of a file containing a brain / background mask.
#' The file can be raw binary or a NIFTI image.
#' Raw binary files must be big endian; the default data type is 16-bit shorts,
#'  but can be changed using the -maskdatatype option.
#'  The program does not process background voxels, but outputs the same
#'  number of values in background voxels and foreground voxels.
#'  Each value is zero in background voxels apart from the exit code which is -1.
#' @param model Type of model used: 1-tensor: dt (linear diffusion tensor, default),
#' ldt - same as dt
#' nldt_pos (nonlinear optimization, constrained to be positive semi-definite),
#' nldt     (unconstrained nonlinear optimization),
#' ldt_wtd  (weighted linear)
#' @param inputdatatype Input data type
#' @param maskdatatype Specifies the type of the mask file;
#' must be big-endian byte ordering. Ignored if a NIFTI mask is used.
#' @param outputdatatype Output data type, ignored if \code{outfile} is specified.
#' @param startpoint  Specifies the starting values for the model.
#' @param outliermap Specifies the name of the file to contain the outlier map generated by the RESTORE algorithm. See restore(1).
#' @param noisemap Specifies the name of the file to contain the estimated noise variance on the diffusion-weighted signal, generated by a weighted tensor fit. The data type of this file is big-endian double.
#' @param residualmap Specifies the name of the file to contain the weighted residual
#' errors after computing a weighted linear tensor fit. One value is produced
#' per measurement, in voxel order. The data type of this file is big-endian double.
#'  Images of the residuals for each measurement can be extracted with shredder.
#' @param sigma Specifies the standard deviation of the noise in the data.
#' Required by the RESTORE algorithm. See restore(1). See datastats(1)
#' for help on how to compute sigma for specific data sets.
#' @param bgthresh Sets a threshold on the average q=0 measurement to
#' separate foreground and background. The program does not process
#' background voxels, but outputs the same number of values in background
#' voxels and foreground voxels. Each value is zero in background voxels
#' apart from the exit code which is -1.
#' @param csfthresh Sets a threshold on the average q=0 measurement to
#' determine which voxels are CSF. This program does not treat CSF voxels any
#' different to other voxels.
#' @param gradadj file of gradient adjustments (relevant for HCP).
#' @param verbose print diagnostic messages
#' @export
camino_modelfit = function(
  infile,
  outfile = NULL,
  scheme,
  mask,
  model = c("dt", "ldt", "nldt_pos", "nldt", "ldt_wtd", "restore"),
  inputdatatype = c("float", "char", "short",
                    "int", "long", "double"),
  maskdatatype = c("float", "char", "short",
                   "int", "long", "double"),
  outputdatatype = c("double", "float", "char", "short",
                     "int", "long"),
  startpoint = NULL,
  outliermap = NULL,
  noisemap = NULL,
  residualmap = NULL,
  sigma = NULL,
  bgthresh = NULL,
  csfthresh = NULL,
  gradadj = NULL,
  verbose = TRUE
) {
  # Checking Inputs
  infile = checkimg(infile)
  mask = checkimg(mask)
  checkimg_null = function(x){
    if (is.null(x)) {
      return(x)
    } else {
      return(checkimg(x))
    }
  }
  gradadj = checkimg_null(gradadj)
  outliermap = checkimg_null(outliermap)

  inputdatatype = match.arg(inputdatatype)
  maskdatatype = match.arg(maskdatatype)
  outputdatatype = match.arg(outputdatatype)
  model = match.arg(model)
  if (is.null(outfile)) {
    outfile = tempfile(fileext = paste0(".B", outputdatatype))
  }

  names(scheme) = NULL
  mask = unname(mask)
  maskdatatype = unname(maskdatatype)
  model = unname(model)
  startpoint = unname(startpoint)
  outliermap = unname(outliermap)
  noisemap = unname(noisemap)
  residualmap = unname(residualmap)
  sigma = unname(sigma)
  bgthresh = unname(bgthresh)
  csfthresh = unname(csfthresh)
  gradadj = unname(gradadj)

  opts = c("-inputfile" = shQuote(infile),
           "-outputfile" = shQuote(outfile),
           "-inputdatatype" = inputdatatype,
           "-schemefile" = scheme,
           "-bgmask" = mask,
           "-maskdatatype" = maskdatatype,
           "-model" = model,
           "-startpoint" = startpoint,
           "-outliermap" = outliermap,
           "-noisemap" = noisemap,
           "-residualmap" = residualmap,
           "-sigma" = sigma,
           "-bgthresh" = bgthresh,
           "-csfthresh" = csfthresh,
           "-gradadj" = gradadj
  )
  opts = paste(names(opts), opts, collapse = " ")

  cmd = camino_cmd("modelfit")
  cmd = paste(cmd, opts)
  if (verbose) {
    message(cmd)
  }
  res = system(cmd)
  if (res != 0) {
    warning("Result is non-zero, may not work")
  }
  return(outfile)
}
# modelfit -inputfile dwi.Bfloat -schemefile 4Ddwi_b1000_bvector.scheme -model ldt -bgmask brain_mask.nii.gz \
# -outputfile dt.Bdouble