From 3e5e1937934ac2eddc92fca1b4e86d283f169e6f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Tim=20Sch=C3=A4fer?= Date: Mon, 21 Dec 2020 22:30:40 +0100 Subject: [PATCH] Add support for reading MRtrix3 scalar DTI track data in TSF format. --- NAMESPACE | 1 + R/read_dti_tck.R | 111 ++++++++++++++++++++++++++++++++++++++++++++ man/read.dti.tsf.Rd | 30 ++++++++++++ 3 files changed, 142 insertions(+) create mode 100644 man/read.dti.tsf.Rd diff --git a/NAMESPACE b/NAMESPACE index dc76226..42ccd53 100644 --- a/NAMESPACE +++ b/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) diff --git a/R/read_dti_tck.R b/R/read_dti_tck.R index 84751f9..258c657 100644 --- a/R/read_dti_tck.R +++ b/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); +} diff --git a/man/read.dti.tsf.Rd b/man/read.dti.tsf.Rd new file mode 100644 index 0000000..8d641f6 --- /dev/null +++ b/man/read.dti.tsf.Rd @@ -0,0 +1,30 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/read_dti_tck.R +\name{read.dti.tsf} +\alias{read.dti.tsf} +\title{Read DTI tracking per-coord data from file in MRtrix 'TSF' format.} +\usage{ +read.dti.tsf(filepath) +} +\arguments{ +\item{filepath}{character string, path to the \code{TSF} file to read.} +} +\value{ +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. +} +\description{ +Read DTI tracking per-coord data from file in MRtrix 'TSF' format. +} +\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. +} +\examples{ +\dontrun{ + tsff = "~/simple.tsf"; + tsf = read.dti.tsf(tsff); +} + +} +\seealso{ +\code{read.dti.tck} +}