Add support for reading MRtrix3 scalar DTI track data in TSF format.

NAMESPACE

@@ -58,6 +58,7 @@ export(ras.to.surfaceras)
 export(ras.to.talairachras)
 export(read.dti.tck)
 export(read.dti.trk)
+export(read.dti.tsf)
 export(read.fs.annot)
 export(read.fs.annot.gii)
 export(read.fs.colortable)

R/read_dti_tck.R

@@ -104,3 +104,114 @@ read.dti.tck <- function(filepath) {
   return(tck);
 }
 
+
+#' @title Read DTI tracking per-coord data from file in MRtrix 'TSF' format.
+#'
+#' @param filepath character string, path to the \code{TSF} file to read.
+#'
+#' @note The data in such a file is one value per track point, the tracks are not part of the file but come in the matching TCK file.
+#'
+#' @seealso \code{read.dti.tck}
+#'
+#' @examples
+#' \dontrun{
+#'  tsff = "~/simple.tsf";
+#'  tsf = read.dti.tsf(tsff);
+#' }
+#'
+#' @return named list with entries 'header' and 'scalars'. The scala data are available in 2 representations: 'merged': a vector of all values (requires external knowledge on track borders), and 'scalar_list': organized into a list of vectors. Each vector represents the values for the points of one track.
+#'
+#' @export
+read.dti.tsf <- function(filepath) {
+
+  tsf = list('header' = list('derived' = list()), 'scalars' = list());
+
+  all_lines = suppressWarnings(readLines(filepath));
+  if(length(all_lines) < 4L) {
+    stop("File not in TSF format: too few lines.");
+  }
+
+  tsf$header$id = all_lines[1];
+  if(tsf$header$id != "mrtrix track scalars") {
+    stop("File not in TSF format: Invalid first line.");
+  }
+
+  for(line_idx in 2L:length(all_lines)) {
+    current_line = all_lines[line_idx];
+    if(current_line == "END") {
+      break;
+    } else {
+      line_parts = unlist(strsplit(current_line, ':'));
+      lkey = trimws(line_parts[1]);
+      lvalue = trimws(line_parts[2]);
+      tsf$header[[lkey]] = lvalue;
+      if(lkey == "file") {
+        file_parts = unlist(strsplit(lvalue, ' '));
+        tsf$header$derived$filename_part = trimws(file_parts[1]); # for future multi-file support, currently always '.'
+        if(tsf$header$derived$filename_part != ".") {
+          stop("Multi-file TSF files not supported.");
+        }
+        tsf$header$derived$data_offset = as.integer(trimws(file_parts[2]));
+      }
+    }
+  }
+
+
+  all_lines = NULL; # free, no longer needed.
+  valid_datatypes = c('Float32BE', 'Float32LE', 'Float64BE', 'Float64LE');
+  if(! tsf$header$datatype %in% valid_datatypes) {
+    stop("Invalid datatype in TSF file header");
+  }
+
+  if(is.null(tsf$header$derived$data_offset)) {
+    stop("Invalid TSF file, missing file offset header entry.");
+  }
+
+  # Determine endianness of following binary data.
+  tsf$header$derived$endian = "little";
+  if(endsWith(tsf$header$datatype, 'BE')) {
+    tsf$header$derived$endian = "big";
+  }
+
+  # Determine size of entries in bytes.
+  tsf$header$derived$dsize = 4L; # default to 32bit
+  if(startsWith(tsf$header$datatype, 'Float64')) {
+    tsf$header$derived$dsize = 8L;
+  }
+
+  # Read binary scalar data.
+  fs = file.size(filepath);
+  num_to_read = (fs - tsf$header$derived$data_offset) / tsf$header$derived$dsize;
+
+  fh = file(filepath, "rb");
+  on.exit({ close(fh) }, add=TRUE);
+
+  seek(fh, where = tsf$header$derived$data_offset, origin = "start");
+  scalar_rawdata = readBin(fh, numeric(), n = num_to_read, size = tsf$header$derived$dsize, endian = tsf$header$derived$endian);
+  # NaN and Inf values are track separators.
+
+  # Generate single vector representation. The user will have to split by tracks based on
+  # knowledge about track borders from the TCK file.
+  data_indices = which(!(is.nan(scalar_rawdata) | is.infinite(scalar_rawdata)));
+  tsf$scalars$merged = scalar_rawdata[data_indices];
+
+  # Generate the alternative list representation of the scalar data:
+  # Filter separators and end marker, organize into tracks list (of matrices).
+  tsf$scalars$scalar_list = list();
+  current_track_idx = 1L;
+  for(value_idx in 1L:length(scalar_rawdata)) {
+    current_value = scalar_rawdata[value_idx];
+    if(is.nan(current_value) | is.infinite(current_value)) {
+      current_track_idx = current_track_idx + 1L;
+      next;
+    } else {
+      if(length(tsf$scalars$scalar_list) < current_track_idx) {
+        tsf$scalars$scalar_list[[current_track_idx]] = current_value;
+      } else {
+        tsf$scalars$scalar_list[[current_track_idx]] = c(tsf$scalars$scalar_list[[current_track_idx]], current_value);
+      }
+    }
+  }
+
+  return(tsf);
+}