tests/testthat/test-read_nrrd.R

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

# Tests for the read-only NRRD volume reader (see R/read_nrrd.R).
#
# The test data in inst/extdata/nrrd was generated by dev_tools/generate_nrrd_test_data.py, which also
# writes reference dumps of the values (read back with pynrrd) and of the geometry (computed with
# SimpleITK, i.e. ITK itself) into extra_test_data/nrrd/expected. The agreement of our reader with these
# two independent implementations is verified by dev_tools/check_nrrd_conversion.R. The tests below pin
# the behavior of the reader itself, and use the small fixtures that can be shipped with the package.

nrrd_test_file <- function(filename) {
  return(system.file("extdata", "nrrd", filename, package = "freesurferformats", mustWork = TRUE))
}

nrrd_type_test_file <- function(filename) {
  return(find_extra_test_data_file(file.path("nrrd", filename)))
}

# The values of all 'vol_u8_*' fixtures in inst/extdata/nrrd, in file order (NRRD stores the first axis
# fastest, which is the same order in which R fills an array).
nrrd_reference_values <- as.numeric(1:24)


test_that("The header of an uncompressed NRRD file is parsed correctly", {
  nrrd_file <- nrrd_test_file("vol_u8_raw.nrrd")
  header <- read.nrrd.header(nrrd_file)

  expect_equal(header$magic, "NRRD0005")
  expect_equal(header$type, "uint8")
  expect_equal(header$dimension, 3L)
  expect_equal(header$sizes, c(4L, 3L, 2L))
  expect_equal(header$encoding, "raw")
  expect_equal(header$endian, "little")
  expect_null(header$data_file)
  expect_null(header$data_files)
  expect_false(header$gzipped_file)
  expect_equal(header$data_offset, header$header_size)
  expect_equal(header$filepath, nrrd_file)
})


test_that("The space information of an NRRD file is parsed and turned into a voxel-to-RAS matrix", {
  header <- read.nrrd.header(nrrd_test_file("vol_u8_raw.nrrd"))

  expect_equal(header$space, "right-anterior-superior")
  expect_equal(as.vector(header$space_origin), c(0, 0, 0))
  expect_equal(unname(header$space_directions), diag(3L))
  expect_equal(header$vox2ras_source, "space directions")
  expect_equal(unname(header$vox2ras_matrix), diag(4L))
})


test_that("The voxel data of an NRRD file is read with the correct dimensions and values", {
  nrrd_file <- nrrd_test_file("vol_u8_raw.nrrd")
  data <- read.fs.volume.nrrd(nrrd_file)

  expect_true(is.array(data))
  expect_equal(dim(data), c(4L, 3L, 2L))
  expect_equal(as.vector(data), nrrd_reference_values)
  expect_equal(data[1, 1, 1], 1)
  expect_equal(data[4, 3, 2], 24) # The last voxel of the volume, stored last in the file.
})


test_that("The parameters flatten, with_header and drop_empty_dims work for NRRD files", {
  nrrd_file <- nrrd_test_file("vol_u8_raw.nrrd")

  flat <- read.fs.volume.nrrd(nrrd_file, flatten = TRUE)
  expect_true(is.vector(flat))
  expect_null(dim(flat))
  expect_equal(flat, nrrd_reference_values)

  with_header <- read.fs.volume.nrrd(nrrd_file, with_header = TRUE, flatten = TRUE)
  expect_equal(class(with_header), "fs.volume")
  expect_equal(names(with_header), c("header", "data"))
  expect_equal(with_header$data, nrrd_reference_values)
  expect_equal(with_header$header$magic, "NRRD0005")
  expect_equal(with_header$header$voldim, 24)

  # This volume has no dimension of length 1, so dropping empty dimensions must not change anything.
  expect_equal(dim(read.fs.volume.nrrd(nrrd_file, drop_empty_dims = TRUE)), c(4L, 3L, 2L))
})


test_that("All NRRD encodings and data locations give the same volume", {
  # The same 4x3x2 uint8 volume (values 1 to 24), stored in every way the format allows: compressed and
  # uncompressed, with the data attached to the header or in a separate file, and with skipped bytes.
  variants <- c("vol_u8_raw.nrrd", "vol_u8_ascii.nrrd", "vol_u8_gzip.nrrd", "vol_u8_bzip2.nrrd",
                "vol_u8_raw.nrrd.gz", "vol_u8_detached.nhdr", "vol_u8_detached_gz.nhdr",
                "vol_byteskip.nrrd", "vol_gzip_byteskip_minus1.nrrd")

  for (variant in variants) {
    data <- read.fs.volume.nrrd(nrrd_test_file(variant))
    expect_equal(dim(data), c(4L, 3L, 2L), info = variant)
    expect_equal(as.vector(data), nrrd_reference_values, info = variant)
  }

  # The header has to report how the data is stored, not only what the values are.
  expect_equal(read.nrrd.header(nrrd_test_file("vol_u8_ascii.nrrd"))$encoding, "ascii")
  expect_equal(read.nrrd.header(nrrd_test_file("vol_u8_gzip.nrrd"))$encoding, "gzip")
  expect_equal(read.nrrd.header(nrrd_test_file("vol_u8_bzip2.nrrd"))$encoding, "bzip2")
})


test_that("A gzipped whole NRRD file and a detached gzipped data file are read correctly", {
  # A NRRD file that is itself gzipped: the header is inside the compressed stream.
  gz_header <- read.nrrd.header(nrrd_test_file("vol_u8_raw.nrrd.gz"))
  expect_true(gz_header$gzipped_file)
  expect_equal(gz_header$type, "uint8")
  expect_equal(gz_header$sizes, c(4L, 3L, 2L))
  expect_equal(as.vector(gz_header$space_origin), c(0, 0, 0))

  # A detached header, i.e. the data lives in a separate file next to it.
  detached_header <- read.nrrd.header(nrrd_test_file("vol_u8_detached.nhdr"))
  expect_equal(detached_header$data_file, "vol_u8_detached.raw")
  expect_false(detached_header$gzipped_file)

  gz_detached_header <- read.nrrd.header(nrrd_test_file("vol_u8_detached_gz.nhdr"))
  expect_equal(gz_detached_header$data_file, "vol_u8_detached_gz.raw.gz")
  expect_false(gz_detached_header$gzipped_file) # Only the data file is gzipped, not the header file.
})


test_that("NRRD files with all supported data types are read correctly", {
  # A big endian int16 file is part of the shipped test data, the remaining types are only available in
  # the development test data (see the doc of the reader for the supported types).
  big_endian <- read.fs.volume.nrrd(nrrd_test_file("vol_i16_be.nrrd"))
  expect_equal(read.nrrd.header(nrrd_test_file("vol_i16_be.nrrd"))$endian, "big")
  expect_equal(as.vector(big_endian), nrrd_reference_values)

  float_file <- nrrd_test_file("vol_f32_oblique.nrrd")
  float_data <- read.fs.volume.nrrd(float_file)
  expect_equal(read.nrrd.header(float_file)$type, "float")
  expect_equal(as.vector(float_data), seq(0.25, 6.0, by = 0.25), tolerance = 1e-6)

  if (is.null(nrrd_type_test_file("vol_type_i8.nrrd"))) {
    testthat::skip("Test data missing.")
  }
  # For all integer and floating point types, the reader has to return the values 1 to 24, see the doc
  # of the 'type' header field: the reader returns doubles for all of them.
  for (typename in c("i8", "u16", "i32", "u32", "i64", "u64", "f32", "f64")) {
    data <- read.fs.volume.nrrd(nrrd_type_test_file(sprintf("vol_type_%s.nrrd", typename)))
    expect_equal(as.vector(data), nrrd_reference_values, info = typename)
  }
})


test_that("Integer NRRD types that R cannot represent natively are converted exactly", {
  if (is.null(nrrd_type_test_file("vol_type_u64.nrrd"))) {
    testthat::skip("Test data missing.")
  }
  # Unsigned 32 bit and both 64 bit integer types have to be read as double without losing precision for
  # the value range covered by the test data (the reader refuses to warn only for values above 2^53).
  for (typename in c("u32", "i64", "u64")) {
    nrrd_file <- nrrd_type_test_file(sprintf("vol_type_%s.nrrd", typename))
    expect_silent(data <- read.fs.volume.nrrd(nrrd_file))
    expect_equal(as.vector(data), nrrd_reference_values, info = typename)
  }
})


test_that("The geometry of oblique and LPS NRRD files is converted correctly", {
  # A rotated (oblique) basis, stored in RAS space.
  oblique_header <- read.nrrd.header(nrrd_test_file("vol_f32_oblique.nrrd"))
  expect_equal(as.vector(oblique_header$space_origin), c(10, -20, 30))
  expect_equal(unname(oblique_header$vox2ras_matrix),
               matrix(c(0, -1, 0, 10,
                        1, 0, 0, -20,
                        0, 0, 1, 30,
                        0, 0, 0, 1), nrow = 4L, byrow = TRUE))

  # The same geometry, but the space directions are given in LPS coordinates (left-posterior-superior).
  # The first two axes of the space directions and of the origin are flipped to get RAS coordinates.
  lps_header <- read.nrrd.header(nrrd_test_file("vol_u8_lps.nrrd"))
  expect_equal(lps_header$space, "left-posterior-superior")
  expect_equal(unname(lps_header$space_directions),
               matrix(c(0, 1, 0, -1, 0, 0, 0, 0, 1), nrow = 3L, byrow = TRUE))
  expect_equal(as.vector(lps_header$space_origin), c(5, 6, 7))
  expect_equal(unname(lps_header$vox2ras_matrix),
               matrix(c(0, 1, 0, -5,
                        -1, 0, 0, -6,
                        0, 0, 1, 7,
                        0, 0, 0, 1), nrow = 4L, byrow = TRUE))
})


test_that("A 4D NRRD file is read correctly", {
  nrrd_file <- nrrd_test_file("vol_u8_4d.nrrd")
  header <- read.nrrd.header(nrrd_file)

  expect_equal(header$dimension, 4L)
  expect_equal(header$sizes, c(4L, 3L, 2L, 5L))
  # The 4th axis has no space direction in this file, which is stored as a row of NA values.
  expect_equal(dim(header$space_directions), c(4L, 3L))
  expect_true(all(is.na(header$space_directions[4L, ])))

  data <- read.fs.volume.nrrd(nrrd_file)
  expect_equal(dim(data), c(4L, 3L, 2L, 5L))
  expect_equal(as.vector(data), as.numeric(1:120))
  expect_equal(data[1, 1, 1, 1], 1)
  expect_equal(data[4, 3, 2, 5], 120)
})


test_that("DWI metadata in NRRD files is parsed", {
  nrrd_file <- nrrd_test_file("vol_dwi.nrrd")
  header <- read.nrrd.header(nrrd_file)

  expect_equal(header$type, "uint16")
  expect_equal(header$dimension, 4L)
  expect_equal(header$sizes, c(4L, 3L, 2L, 6L))
  expect_equal(as.vector(read.fs.volume.nrrd(nrrd_file)), as.numeric(1:144))

  expect_equal(header$dwi$b_value, 1000)
  expect_equal(header$dwi$num_gradients, 6L)
  expect_equal(dim(header$dwi$bvec), c(6L, 3L))

  # The first gradient is the b=0 (or b-value scaling) entry, the other ones are the diffusion directions.
  expect_equal(as.vector(header$dwi$bvec[1L, ]), c(0, 0, 0))
  expect_equal(as.vector(header$dwi$bvec[2L, ]), c(1, 0, 0))
  expect_equal(as.vector(header$dwi$bvec[5L, ]), rep(1 / sqrt(3), 3L), tolerance = 1e-6)
  # The gradient field names from the file are used as row names, normalized to lower case like all
  # other header field keys.
  expect_equal(rownames(header$dwi$bvec)[5L], "dwmrigradient0004")

  # The measurement frame is the identity for this file.
  expect_equal(unname(header$dwi$measurement_frame), diag(3L))
})


test_that("NRRD files without space information are handled", {
  if (is.null(nrrd_type_test_file("vol_u8_no_geometry.nrrd"))) {
    testthat::skip("Test data missing.")
  }
  nrrd_file <- nrrd_type_test_file("vol_u8_no_geometry.nrrd")
  header <- read.nrrd.header(nrrd_file)

  expect_null(header$space)
  expect_null(header$space_directions)
  expect_null(header$space_origin)
  expect_null(header$vox2ras_matrix)
  # The data itself is still read fine.
  expect_equal(as.vector(read.fs.volume.nrrd(nrrd_file)), nrrd_reference_values)
})


test_that("NRRD files with a sheared basis, an ASCII dtype unit and line skips are read", {
  if (is.null(nrrd_type_test_file("vol_f64_sheared.nrrd"))) {
    testthat::skip("Test data missing.")
  }
  # A non-orthogonal basis cannot be stored by ITK, but the file format allows it.
  sheared <- read.fs.volume.nrrd(nrrd_type_test_file("vol_f64_sheared.nrrd"))
  expect_equal(dim(sheared), c(4L, 3L, 2L))
  expect_equal(as.vector(sheared), seq(1.5, 36.0, by = 1.5), tolerance = 1e-9)

  line_skip <- read.fs.volume.nrrd(nrrd_type_test_file("vol_lineskip.nrrd"))
  expect_equal(as.vector(line_skip), nrrd_reference_values)
})


test_that("A list of data files in the NRRD header is supported", {
  if (is.null(nrrd_type_test_file("vol_datafile_list.nhdr"))) {
    testthat::skip("Test data missing.")
  }
  # This file uses the 'data file: LIST' mode, i.e. the data is stored in several files whose contents
  # are concatenated. It contains the same volume as the fixtures above.
  nrrd_file <- nrrd_type_test_file("vol_datafile_list.nhdr")
  header <- read.nrrd.header(nrrd_file)

  expect_equal(length(header$data_files), 2L)
  expect_equal(basename(header$data_files[[1L]]), "vol_datafile_list_0.raw")
  expect_equal(basename(header$data_files[[2L]]), "vol_datafile_list_1.raw")
  expect_equal(as.vector(read.fs.volume.nrrd(nrrd_file)), nrrd_reference_values)
})


test_that("The NRRD reader reports errors for invalid input files", {
  expect_error(read.nrrd.header("file_that_does_not_exist.nrrd"), "does not exist")
  expect_error(read.fs.volume.nrrd("file_that_does_not_exist.nrrd"), "does not exist")

  # A file that is not an NRRD file at all.
  not_nrrd <- tempfile(fileext = ".nrrd")
  writeLines(c("hello", "", "world"), not_nrrd)
  expect_error(read.nrrd.header(not_nrrd), "not in NRRD format")

  # A header without the mandatory 'sizes' field.
  no_sizes <- tempfile(fileext = ".nrrd")
  writeLines(c("NRRD0005", "type: uint8", "encoding: raw", "endian: little", "", "abc"), no_sizes)
  expect_error(read.nrrd.header(no_sizes), "has no valid 'sizes' header field")

  # An unknown data type.
  bad_type <- tempfile(fileext = ".nrrd")
  writeLines(c("NRRD0005", "type: bogus", "dimension: 3", "sizes: 4 3 2", "encoding: raw",
               "endian: little", "", "abc"), bad_type)
  expect_error(read.nrrd.header(bad_type), "Unsupported NRRD data type 'bogus'")

  # An encoding that the format knows, but that this reader does not implement.
  bad_encoding <- tempfile(fileext = ".nrrd")
  writeLines(c("NRRD0005", "type: uint8", "dimension: 3", "sizes: 4 3 2", "encoding: hex",
               "endian: little", "", "abcd"), bad_encoding)
  expect_error(read.fs.volume.nrrd(bad_encoding), "Unsupported NRRD encoding 'hex'")

  # A file in which the data is shorter than the header claims.
  truncated <- tempfile(fileext = ".nrrd")
  writeLines(c("NRRD0005", "type: uint8", "dimension: 3", "sizes: 4 3 2", "encoding: raw",
               "endian: little", "", "abcd"), truncated)
  expect_error(read.fs.volume.nrrd(truncated), "is truncated")
})


test_that("The NRRD reader is reachable through the read.fs.volume dispatcher", {
  # By file extension, for both file name endings of the format.
  expect_equal(as.vector(read.fs.volume(nrrd_test_file("vol_u8_raw.nrrd"))), nrrd_reference_values)
  expect_equal(as.vector(read.fs.volume(nrrd_test_file("vol_u8_detached.nhdr"))), nrrd_reference_values)

  # And by explicitly requesting the format.
  expect_equal(as.vector(read.fs.volume(nrrd_test_file("vol_u8_raw.nrrd"), format = "nrrd")),
               nrrd_reference_values)

  expect_error(read.fs.volume(nrrd_test_file("vol_u8_raw.nrrd"), format = "bogus"),
               "Format must be one of")
})