R/gifti_writer.R

In freesurferformats: Read and Write 'FreeSurfer' Neuroimaging File Formats

# This file contains basic functions for generating valid GIFTI XML documents and writing them to files.
#
# Note that GIFTI is a very versatile and complex format, and one can write a lot of different things in different formats.
# While this is great, it also means that it is not guaranteed that a software x will be able to interprete the data stored in a GIFTI file by software y in the
# way intended by y (even though the file is technically a valid GIFTI file).
# The functions in here leave the user the ability to edit the generated XML trees (e.g., with the xml2 package), so you should be able to use these functions
# to write data exporters for any software.
#
# Some writers for how a specific software stores data of certain typpes (surface, measurements/morphometry, labels) in GIFTI format files are implemented in separate
# source files. E.g., writers for FreeSurfer can be found in gifti_writer_freesurfer.R.

#' @title Get GIFTI XML representation of data.
#'
#' @description Creates a GIFTI XML tree from your datasets (vectors and matrices). The tree can be further modified to add additional data, or written to a file as is to produce a valid GIFTI file (see \code{\link{gifti_xml_write}}).
#'
#' @param data_array list of data vectors and/or data matrices.
#'
#' @param intent vector of NIFTI intent strings for the data vectors in 'data_array' parameter, see \code{\link[gifti]{convert_intent}}. Example: 'NIFTI_INTENT_SHAPE'. See https://nifti.nimh.nih.gov/nifti-1/documentation/nifti1fields/nifti1fields_pages/group__NIFTI1__INTENT__CODES.html.
#'
#' @param datatype vector of NIFTI datatype strings. Example: 'NIFTI_TYPE_FLOAT32'. Should be suitable for your data.
#'
#' @param encoding vector of encoding definition strings. One of 'ASCII', 'Base64Binary', 'GZipBase64Binary'.
#'
#' @param endian vector of endian definition strings. One of 'LittleEndian' or 'BigEndian'. See \code{\link[gifti]{convert_endian}}.
#'
#' @param transform_matrix optional, a list of transformation matrices, one for each data_array. If one of the data arrays has none, pass `NA`. Each transformation matrix in the outer list has to be a 4x4 matrix or given as a named list with entries 'transform_matrix', 'data_space', and 'transformed_space'. Here is an example: \code{list('transform_matrix'=diag(4), 'data_space'='NIFTI_XFORM_UNKNOWN', 'transformed_space'='NIFTI_XFORM_UNKNOWN')}.
#'
#' @param force logical, whether to force writing the data, even if issues like a mismatch of datatype and data values are detected.
#'
#' @return xml tree, see xml2 package. One could modify this tree as needed using xml2 functions, e.g., add metadata.
#'
#' @references See https://www.nitrc.org/frs/download.php/2871/GIFTI_Surface_Format.pdf
#'
#' @note Unless you want to modify the returned tree manually, you should not need to call this function. Use \code{\link{gifti_writer}} instead.
#'
#' @seealso The example for \code{\link{gifti_xml_write}} shows how to modify the tree.
#'
#' @examples
#' \dontrun{
#'   my_data_sets = list(rep(3.1, 3L), matrix(seq(6)+0.1, nrow=2L));
#'   transforms = list(NA, list('transform_matrix'=diag(4), 'data_space'='NIFTI_XFORM_UNKNOWN',
#'    'transformed_space'='NIFTI_XFORM_UNKNOWN'));
#'   xmltree = gifti_xml(my_data_sets, datatype='NIFTI_TYPE_FLOAT32', transform_matrix=transforms);
#'   # Verify that the tree is a valid GIFTI file:
#'   gifti_xsd = "https://www.nitrc.org/frs/download.php/158/gifti.xsd";
#'   xml2::xml_validate(xmltree, xml2::read_xml(gifti_xsd));
#' }
#'
#' @importFrom xml2 xml_new_root xml_set_attr xml_add_child read_xml
#' @export
gifti_xml <- function(data_array, intent='NIFTI_INTENT_SHAPE', datatype='NIFTI_TYPE_FLOAT32', encoding='GZipBase64Binary', endian='LittleEndian', transform_matrix=NULL, force=FALSE) {

  if (requireNamespace("gifti", quietly = TRUE)) {
    if( ! is.list(data_array)) {
      stop("Parameter 'data_array' must be a list.");
    }

    supported_encodings = c('ASCII', 'Base64Binary', 'GZipBase64Binary');

    num_data_arrays = length(data_array);
    dataarray_contains_matrices = any(lapply((lapply(data_array, dim)), length) > 0L);

    num_transform_matrices = 0L;
    if(! is.null(transform_matrix)) {
      if( ! is.list(transform_matrix)) {
        stop("Parameter 'transform_matrix' must be NULL or a list.");
      }
      if( length(names(transform_matrix)) != 0L) {
        stop("Parameter 'transform_matrix' must not be a named list. Hint: When passing a single matrix, you need to enclose it in an outer list.");
      }
      num_transform_matrices = length(transform_matrix);
      if(num_transform_matrices != num_data_arrays) {
        stop(sprintf("Found %d data arrays, but %d transform matrices: mismatch. Pass NA if you have no matrix for a data array.\n", num_data_arrays, num_transform_matrices));
      }
    }

    if(num_data_arrays > 1L) {
      if(length(intent) == 1L) {
        intent = rep(intent, num_data_arrays);
      }
      if(length(datatype) == 1L) {
        datatype = rep(datatype, num_data_arrays);
      }
      if(length(encoding) == 1L) {
        encoding = rep(encoding, num_data_arrays);
      }
      if(length(endian) == 1L) {
        endian = rep(endian, num_data_arrays);
      }
    }

    dim0 = rep(1L, num_data_arrays);                               # gets filled later
    dim1 = rep(1L, num_data_arrays);                               # gets filled later
    dimensionality = rep(1L, num_data_arrays);                     # gets filled later
    array_indexing_order = rep("RowMajorOrder", num_data_arrays);  # currently fixed

    root = xml2::xml_new_root("GIFTI", 'Version' = "1.0", 'NumberOfDataArrays'=num_data_arrays);
    metadata = xml2::xml_add_child(root, read_xml("<MetaData></MetaData>"));
    xml2::xml_add_child(metadata, read_xml("<MD><Name>Generator</Name><Value>fsbrain</Value></MD>"));


    da_index = 1L;
    data_is_matrix = FALSE;
    for(da in data_array) {
      if(is.vector(da)) {
        dim0[da_index] = length(da);
        dim1[da_index] = 1L;
        dimensionality[da_index] = 1L;
      } else if(is.matrix(da)) {
        data_is_matrix = TRUE;
        dimensionality[da_index] = 2L;
        dim0[da_index] = dim(da)[1];
        dim1[da_index] = dim(da)[2];
        if(array_indexing_order[da_index] == "RowMajorOrder") {
          da = as.vector(t(da));
        } else if(array_indexing_order[da_index] == "ColumnMajorOrder") {
          da = as.vector((da));
        } else {
          stop(sprintf("Dataarray # %d: invalid array_indexing_order, must be 'RowMajorOrder' or 'ColumnMajorOrder'.\n", da_index));
        }
      } else {
        stop("The data_arrays must be of type vector or matrix.");
      }

      if(! encoding[da_index] %in% supported_encodings) {
        stop(sprintf("Dataarray # %d: invalid encoding '%s'.\n", da_index, encoding[da_index]));
      }
      if(! endian[da_index] %in% c('LittleEndian', 'BigEndian')) {
        stop(sprintf("Dataarray # %d: invalid endian '%s'.\n", da_index, endian[da_index]));
      }

      check_data_and_settings_consistency(da_index, da, datatype[da_index], intent[da_index], force=force);

      data_array_node = xml2::read_xml("<DataArray/>");
      xml2::xml_set_attr(data_array_node, 'Dimensionality', dimensionality[da_index]);
      xml2::xml_set_attr(data_array_node, 'Dim0', dim0[da_index]);

      if(dataarray_contains_matrices) {   # We only need this attribute if any dataarray has more than 1 dimension. But if so, we need it for all datasets.
        xml2::xml_set_attr(data_array_node, 'Dim1', dim1[da_index]);
      }

      xml2::xml_set_attr(data_array_node, 'Encoding', encoding[da_index]);
      xml2::xml_set_attr(data_array_node, 'DataType', datatype[da_index]);
      xml2::xml_set_attr(data_array_node, 'Intent', intent[da_index]);
      xml2::xml_set_attr(data_array_node, 'Endian', endian[da_index]);

      xml2::xml_set_attr(data_array_node, 'ExternalFileName', '');                   # not supported atm
      xml2::xml_set_attr(data_array_node, 'ExternalFileOffset', '');                 # not supported atm
      xml2::xml_set_attr(data_array_node, 'ArrayIndexingOrder', array_indexing_order[da_index]);

      data_array_node_added = xml2::xml_add_child(root, data_array_node);
      data_array_metadata = xml2::xml_add_child(data_array_node_added, xml2::read_xml("<MetaData></MetaData>"));
      encoded_data = gifti::data_encoder(da, encoding = encoding[da_index], datatype = datatype[da_index], endian = endian[da_index]);
      data_node = xml2::read_xml(sprintf("<Data>%s</Data>", encoded_data));

      if(num_transform_matrices > 0L) {
        tf = transform_matrix[[da_index]];
        if(is.list(tf)) {
          tf_node = xml_node_gifti_coordtransform(tf$transform_matrix, data_space=tf$data_space, transformed_space=tf$transformed_space);
          xml2::xml_add_child(data_array_node_added, tf_node);
        } else if(is.matrix(tf)) {
          tf_node = xml_node_gifti_coordtransform(tf);
          xml2::xml_add_child(data_array_node_added, tf_node);
        }
        else {
          if(! is.na(tf)) {
            stop(sprintf("Invalid transformation matrix at index %d: neither a named list, nor NA.\n", da_index));
          }
        }
      }

      xml2::xml_add_child(data_array_node_added, data_node);

      da_index = da_index + 1L;
    }
    return(root);
  } else {
    stop("Writing files in GIFTI format requires the 'gifti' package to be installed.");
  }
}


#' @title Warn about common errors in combining data and datatype.
#'
#' @param index positive integer, the dataarray index to report. Makes it easier for the user to find the broken one.
#'
#' @param data vector or matrix, the data to write to the GIFTI file. Checked against the datatype.
#'
#' @param datatype NIFTI datatype string, the datatype to use when writing to the GIFTI file. Checked against the data.
#'
#' @param intent NIFTI intent string, checked independently. In no way do we check whether it makes sense for the data.
#'
#' @note The checks in here are in no way exhaustive.
#'
#' @keywords internal
check_data_and_settings_consistency <- function(index, data, datatype, intent, force=FALSE) {
  msg = NULL;
  if(is.integer(data) & startsWith(datatype, "NIFTI_TYPE_FLOAT")) {
    msg = sprintf("Dataset # %d in file will be corrupted: integer R data passed, and written as '%s'.\n", index, datatype);
  }

  if(gifti::convert_intent(intent) == "unknown") {
    msg = sprintf("Dataset # %d: Invalid NIFTI intent '%s'.\n", index, intent);
  }

  # Should check more here.

  if(! is.null(msg)) {
    if(force) {
      warning(msg);
    } else {
      stop(msg);
    }
  }
}


#' @title Write XML tree to a gifti file.
#'
#' @param filepath path to the output gifti file
#'
#' @param xmltree XML tree from xml2
#'
#' @param options output options passed to \code{\link[xml2]{write_xml}}.
#'
#' @references \url{https://www.nitrc.org/frs/download.php/2871/GIFTI_Surface_Format.pdf}
#'
#' @examples
#' \dontrun{
#'   outfile = tempfile(fileext = '.gii');
#'   my_data_sets = list(rep(3.1, 3L), matrix(seq(6)+0.1, nrow=2L));
#'   xmltree = gifti_xml(my_data_sets, datatype='NIFTI_TYPE_FLOAT32');
#'   # Here we add global metadata:
#'   xmltree = gifti_xml_add_global_metadata(xmltree, list("User"="Me", "Day"="Monday"));
#'   # Validating your XML never hurts
#'   gifti_xsd = "https://www.nitrc.org/frs/download.php/158/gifti.xsd";
#'   xml2::xml_validate(xmltree, xml2::read_xml(gifti_xsd));
#'   gifti_xml_write(outfile, xmltree);  # Write your custom tree to a file.
#' }
#'
#' @export
#' @importFrom xml2 write_xml
gifti_xml_write <- function(filepath, xmltree, options=c('as_xml', 'format')) {
  return(invisible(xml2::write_xml(xmltree, file=filepath, options=options)));
}


#' @title Write data to a gifti file.
#'
#' @param filepath path to the output gifti file
#'
#' @param ... parameters passed to \code{\link{gifti_xml}}.
#'
#' @references \url{https://www.nitrc.org/frs/download.php/2871/GIFTI_Surface_Format.pdf}
#'
#' @examples
#' \dontrun{
#'   outfile = tempfile(fileext = '.gii');
#'   dataarrays = list(rep(3.1, 3L), matrix(seq(6), nrow=2L));
#'   gifti_writer(outfile, dataarrays, datatype=c('NIFTI_TYPE_FLOAT32', 'NIFTI_TYPE_INT32'));
#' }
#'
#' @export
gifti_writer <- function(filepath, ...) {
  if(! is.character(filepath)) {
    stop("Parameter 'filepath' must be a character string.");
  }
  xmltree = gifti_xml(...);
  return(invisible(gifti_xml_write(filepath, xmltree)));
}