tests/testthat/test-read_analyze.R

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

# Tests for reading ANALYZE 7.5 files and NIFTI v1 pair files (see R/read_analyze.R).
#
# The fixtures in inst/extdata/analyze were generated by dev_tools/generate_analyze_test_data.py (with nibabel)
# and verified against nibabel and FreeSurfer by dev_tools/check_analyze_conversion.R. The larger and more
# exotic fixtures live in extra_test_data/analyze, which is not part of the package, and the tests that need
# them skip when the data is not available.

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

analyze_extra_file <- function(filename) {
  return(find_extra_test_data_file(file.path("analyze", filename)))
}

# The data of the 'tiny_*' fixtures: 4x3x2 values, 1 to 24 in file order.
analyze_reference_values <- as.numeric(1:24)


test_that("The ANALYZE header of a file is parsed correctly", {
  hdrfile <- analyze_test_file("tiny_u8.hdr")
  header <- read.analyze.header(hdrfile)

  expect_equal(header$sizeof_hdr, 348L)
  expect_equal(header$datatype, 2L) # unsigned 8 bit
  expect_equal(header$bitpix, 8L)
  expect_equal(header$dim, c(3L, 4L, 3L, 2L, 1L, 1L, 1L, 1L))
  expect_equal(header$pix_dim[2:4], c(1., 2., 3.))
  expect_equal(header$orient, 1L)
  expect_equal(header$descrip, "freesurferformats")
  expect_equal(header$vox_offset, 0.)
  expect_equal(header$funused1, 0.)
  expect_null(header$spm_origin)
  expect_equal(header$magic, "")
  expect_equal(header$header_format, "analyze")
  expect_equal(header$filepath_header, hdrfile)
  expect_equal(header$filepath_image, analyze_test_file("tiny_u8.img"))
})


test_that("The ANALYZE volume reader returns the data with the correct dimensions and values", {
  vol <- read.fs.volume.analyze(analyze_test_file("tiny_u8.hdr"), with_header = TRUE)

  expect_equal(class(vol), "fs.volume")
  expect_equal(names(vol), c("header", "data"))
  expect_equal(dim(vol$data), c(4L, 3L, 2L))
  expect_equal(as.vector(vol$data), analyze_reference_values)
  # ANALYZE does not store the orientation, so no matrix is reported by default.
  expect_null(vol$header$vox2ras_matrix)
})


test_that("Signed 16 bit ANALYZE data is read correctly", {
  # The values of this fixture include negative ones and values above 127, so a reader that uses the wrong
  # signedness returns different data.
  vol <- read.fs.volume.analyze(analyze_test_file("tiny_i16.hdr"), with_header = TRUE)
  expect_equal(vol$header$datatype, 4L)
  expect_equal(vol$header$bitpix, 16L)
  expect_equal(as.vector(vol$data), c(-1, 130, 3, 32767, -32768, 6, 7, 8, 9, 1000, -1000, 12, 13, 14, -250, 16, 17, 200, 19, 20, 21, 22, 23, 24))
})


test_that("A 4D ANALYZE file is read correctly", {
  vol <- read.fs.volume.analyze(analyze_test_file("tiny_4d.hdr"), with_header = TRUE)
  expect_equal(vol$header$dim[1], 4L)
  expect_equal(dim(vol$data), c(4L, 3L, 2L, 4L))
  # Frame f contains the values 1..24 plus 100*f, in file order.
  for (frame in 1L:4L) {
    expect_equal(as.vector(vol$data)[(1L:24L) + 24L * (frame - 1L)], analyze_reference_values + 100 * (frame - 1L))
  }
})


test_that("The parameters flatten, with_header and drop_empty_dims work for ANALYZE files", {
  hdrfile <- analyze_test_file("tiny_u8.hdr")

  flat <- read.fs.volume.analyze(hdrfile, flatten = TRUE)
  expect_true(is.vector(flat))
  expect_equal(flat, analyze_reference_values)

  with_header <- read.fs.volume.analyze(hdrfile, with_header = TRUE, flatten = TRUE)
  expect_equal(with_header$header$voldim, 24)
  expect_equal(with_header$header$filepath_image, analyze_test_file("tiny_u8.img"))

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


test_that("Compressed ANALYZE and NIFTI pair files are read", {
  vol <- read.fs.volume.analyze(analyze_test_file("tiny_u8_gz.hdr.gz"), with_header = TRUE)
  expect_equal(as.vector(vol$data), analyze_reference_values)
  expect_true(grepl("\\.hdr\\.gz$", vol$header$filepath_header))
  expect_true(grepl("\\.img\\.gz$", vol$header$filepath_image))

  pair <- read.fs.volume.analyze(analyze_test_file("pair_u8_gz.hdr.gz"), with_header = TRUE)
  expect_equal(dim(pair$data), c(4L, 3L, 2L))
  expect_equal(as.vector(pair$data), analyze_reference_values)
  # The affine of this fixture has non-unit voxel sizes, so the matrix is not a sign matrix.
  expect_true(!is.null(pair$header$vox2ras_matrix))
})


test_that("The two files of a pair can be given in several ways", {
  hdrfile <- analyze_test_file("tiny_u8.hdr")
  imgfile <- analyze_test_file("tiny_u8.img")
  basefile <- file.path(dirname(hdrfile), "tiny_u8")

  expect_equal(as.vector(read.fs.volume.analyze(hdrfile)), analyze_reference_values)
  expect_equal(as.vector(read.fs.volume.analyze(imgfile)), analyze_reference_values)
  expect_equal(as.vector(read.fs.volume.analyze(basefile)), analyze_reference_values)

  pair <- analyze.pair.files(basefile)
  expect_equal(pair$header, hdrfile)
  expect_equal(pair$image, imgfile)
  expect_true(pair$header_exists)
  expect_true(pair$image_exists)
})


test_that("The ANALYZE reader reports an error for a missing data file and for NIFTI files", {
  # A header file whose data file does not exist.
  hdrfile <- analyze_test_file("tiny_u8.hdr")
  tmpdir <- tempfile()
  dir.create(tmpdir)
  broken <- file.path(tmpdir, "broken.hdr")
  file.copy(hdrfile, broken)
  expect_error(read.fs.volume.analyze(broken), "does not exist")
  expect_error(read.analyze.data(broken), "does not exist")

  # The NIFTI v1 pair fixtures are NIFTI files, not ANALYZE files.
  expect_error(read.analyze.header(analyze_test_file("pair_u8.hdr")), "not an ANALYZE 7.5 file")
  # ...but the volume reader handles both variants.
  expect_equal(as.vector(read.fs.volume.analyze(analyze_test_file("pair_u8.hdr"))), analyze_reference_values)

  # A single file NIFTI image has no data file and cannot be read as a pair.
  single_file <- system.file("extdata", "lh.area.gz", package = "freesurferformats", mustWork = TRUE)
  expect_error(read.fs.volume.analyze(single_file), "does not exist|not in ANALYZE")
})


test_that("is.analyze.file distinguishes ANALYZE files from NIFTI v1 pairs", {
  expect_true(is.analyze.file(analyze_test_file("tiny_u8.hdr")))
  expect_true(is.analyze.file(analyze_test_file("tiny_spm.hdr")))
  expect_false(is.analyze.file(analyze_test_file("pair_u8.hdr")))
  # A single file NIFTI image is not a pair file either, it has the 'n+1' magic.
  single_file <- tempfile(fileext = ".nii")
  write.nifti1(single_file, array(1:24, dim = c(4L, 3L, 2L)))
  expect_false(is.analyze.file(single_file))
})


test_that("NIFTI v1 pair files are read with their geometry", {
  # This fixture stores an sform.
  vol <- read.fs.volume.analyze(analyze_test_file("pair_u8.hdr"), with_header = TRUE)
  expect_equal(as.vector(vol$data), analyze_reference_values)
  expect_equal(vol$header$vox2ras_source, "sform")
  expect_equal(vol$header$sform_code, 2L)
  expect_equal(unname(vol$header$vox2ras_matrix),
               matrix(c(-1., 0., 0., 10., 0., 0., 1., -20., 0., -1., 0., 30., 0., 0., 0., 1.), nrow = 4L, byrow = TRUE))

  # This one stores only a qform, so the quaternion code path is used.
  qvol <- read.fs.volume.analyze(analyze_test_file("pair_qform_i16.hdr"), with_header = TRUE)
  expect_equal(qvol$header$sform_code, 0L)
  expect_equal(qvol$header$qform_code, 1L)
  expect_equal(qvol$header$vox2ras_source, "qform")
  # The affine is a rotation with voxel sizes 1.5, 1.5 and 2.5 and a flip (qfac), so its determinant is negative.
  expect_true(det(qvol$header$vox2ras_matrix[1:3, 1:3]) < 0)
  expect_equal(abs(det(qvol$header$vox2ras_matrix[1:3, 1:3])), 1.5 * 1.5 * 2.5, tolerance = 1e-5)
})


test_that("The SPM interpretation of an ANALYZE header is available and off by default", {
  hdrfile <- analyze_test_file("tiny_spm.hdr")
  header <- read.analyze.header(hdrfile)
  expect_equal(header$funused1, 0.5) # the SPM data scale factor
  expect_equal(header$spm_origin, c(2L, 2L, 2L)) # the SPM image origin
  expect_equal(as.integer(charToRaw(header$originator)), c(2L, 2L, 2L))

  # Ignoring the SPM fields is the default, and it warns about them, since the data would be wrong otherwise.
  expect_warning(vol <- read.fs.volume.analyze(hdrfile, with_header = TRUE), "SPM")
  expect_equal(sum(vol$data), sum(analyze_reference_values)) # unscaled
  expect_null(vol$header$vox2ras_matrix)

  # With spm = TRUE, the scale factor is applied and a matrix is derived from the voxel sizes and the origin.
  spm_vol <- read.fs.volume.analyze(hdrfile, with_header = TRUE, spm = TRUE)
  expect_equal(sum(spm_vol$data), sum(analyze_reference_values) * 0.5)
  expect_equal(spm_vol$header$data_scaled_by, 0.5)
  expect_equal(spm_vol$header$vox2ras_source, "spm origin")
  # The origin is at voxel (1, 1, 1) (0-based), so that voxel maps to the world origin.
  expect_equal(as.numeric(spm_vol$header$vox2ras_matrix[1:3, 4]), c(1., -1., -1.))

  # A file without SPM fields does not warn, and the ANALYZE convention is used for the matrix.
  expect_silent(simple_vol <- read.fs.volume.analyze(analyze_test_file("tiny_u8.hdr"), with_header = TRUE, spm = TRUE))
  expect_equal(simple_vol$header$vox2ras_source, "analyze convention")
  expect_equal(as.numeric(simple_vol$header$vox2ras_matrix[1:3, 4]), c(1.5, -2., -1.5))
})


test_that("Reading a pair header with the NIFTI reader and vice versa reports helpful errors", {
  # This used to return the header bytes of the file as voxel data, without any error.
  expect_error(read.nifti1.data(analyze_test_file("tiny_u8.hdr")), "ANALYZE")
  # The NIFTI reader can read the header of a NIFTI pair, that is the same 348 byte layout.
  nii_header <- read.nifti1.header(analyze_test_file("pair_u8.hdr"))
  expect_equal(nii_header$magic, "ni1")
  expect_equal(nii_header$dim, c(3L, 4L, 3L, 2L, 1L, 1L, 1L, 1L))
})


test_that("The ANALYZE reader is reachable through the read.fs.volume dispatcher", {
  expect_equal(as.vector(read.fs.volume(analyze_test_file("tiny_u8.hdr"))), analyze_reference_values)
  expect_equal(as.vector(read.fs.volume(analyze_test_file("tiny_u8.img"))), analyze_reference_values)
  expect_equal(as.vector(read.fs.volume(analyze_test_file("pair_u8.hdr"))), analyze_reference_values)
  expect_equal(as.vector(read.fs.volume(analyze_test_file("tiny_u8.hdr"), format = "analyze")), analyze_reference_values)
  # The base name of a pair is accepted as well.
  expect_equal(as.vector(read.fs.volume(file.path(dirname(analyze_test_file("tiny_u8.hdr")), "tiny_u8"))), analyze_reference_values)

  expect_error(read.fs.volume(analyze_test_file("tiny_u8.hdr"), format = "bogus"), "Format must be one of")
})


test_that("The MATLAB sidecar file of an ANALYZE image is read as the transformation matrix", {
  # SPM and FreeSurfer write the transformation matrix into a MATLAB file next to the image. The fixture was
  # written by nibabel's Spm99AnalyzeImage class, which stores the affine that was passed in, and the check script
  # dev_tools/check_analyze_conversion.R verifies that we compute the same affine as nibabel does.
  hdrfile <- analyze_test_file("spm_mat_u8.hdr")
  vol <- read.fs.volume.analyze(hdrfile, with_header = TRUE)
  expect_equal(vol$header$vox2ras_source, "mat sidecar")
  expect_equal(unname(vol$header$vox2ras_matrix),
               matrix(c(-1., 0., 0., 10., 0., 0., 1., -20., 0., -1., 0., 30., 0., 0., 0., 1.), nrow = 4L, byrow = TRUE))
  # The data is unaffected by the geometry, and no warning is needed for this file.
  expect_equal(as.vector(vol$data), analyze_reference_values)

  # The sidecar file is found next to the data file.
  pair <- analyze.pair.files(hdrfile)
  expect_true(pair$mat_exists)
  expect_equal(pair$mat, analyze_test_file("spm_mat_u8.mat"))
})


test_that("The MATLAB v4 parser reads the matrices of a sidecar file", {
  variables <- read.matlab.v4.matrix(analyze_test_file("spm_mat_u8.mat"))
  expect_equal(sort(names(variables)), c("M", "mat"))
  expect_equal(dim(variables$M), c(4L, 4L))
  expect_equal(dim(variables$mat), c(4L, 4L))

  # nibabel writes both variables: 'mat' includes the flip of the first axis and the 1-based voxel indices of
  # MATLAB, 'M' is the same without the flip (nibabel's docs: "the 'M' matrix does not include flips").
  mat <- variables$mat
  m <- variables$M
  expect_equal(mat, diag(c(-1., 1., 1., 1.)) %*% m, tolerance = 1e-12)

  # Our matrix is 'mat' adjusted for 0-based voxel indices, which adds the sum of the rows of the rotation part.
  expected <- mat
  expected[1:3, 4L] <- expected[1:3, 4L] + rowSums(expected[1:3, 1:3])
  expect_equal(unname(analyze.mat.sidecar.to.vox2ras(analyze_test_file("spm_mat_u8.mat"))$vox2ras), expected)

  # A file that is not a MATLAB v4 file is reported as unreadable, with a reason.
  not_matlab <- tempfile(fileext = ".mat")
  writeLines("this is not a MATLAB file at all", not_matlab)
  expect_null(read.matlab.v4.matrix(not_matlab))
  expect_null(analyze.mat.sidecar.to.vox2ras(not_matlab)$vox2ras)
  expect_true(nchar(analyze.mat.sidecar.to.vox2ras(not_matlab)$reason) > 0L)

  # So is a valid v4 file that contains different variables.
  other_variables <- tempfile(fileext = ".mat")
  con <- file(other_variables, "wb")
  writeBin(as.integer(c(0L, 2L, 2L, 0L, 6L)), con, size = 4L) # type double, 2x2, no imaginary part, 6 bytes of name
  writeChar("other", con, eos = NULL)
  writeBin(as.raw(0), con)
  writeBin(as.double(1:4), con)
  close(con)
  expect_equal(sort(names(read.matlab.v4.matrix(other_variables))), "other")
  result <- analyze.mat.sidecar.to.vox2ras(other_variables)
  expect_null(result$vox2ras)
  expect_true(grepl("other", result$reason))
})


test_that("A sidecar file that cannot be used does not prevent reading the image", {
  tmpdir <- tempfile()
  dir.create(tmpdir)
  base <- file.path(tmpdir, "vol")
  file.copy(analyze_test_file("tiny_u8.hdr"), paste0(base, ".hdr"))
  file.copy(analyze_test_file("tiny_u8.img"), paste0(base, ".img"))
  writeLines("not a MATLAB file", paste0(base, ".mat"))

  expect_warning(vol <- read.fs.volume.analyze(base, with_header = TRUE), "sidecar")
  expect_equal(as.vector(vol$data), analyze_reference_values)
  expect_null(vol$header$vox2ras_matrix)
})


test_that("ANALYZE files with unusual headers are read", {
  badbitpix <- analyze_extra_file("vol_badbitpix.hdr")
  if (is.null(badbitpix)) {
    testthat::skip("Test data missing.")
  }
  # The file is an int16 file whose 'bitpix' field claims 32 bits per value. The data type is the reliable
  # field, so the reader uses it and warns.
  expect_warning(vol <- read.fs.volume.analyze(badbitpix), "bits per value")
  expect_equal(as.vector(vol), c(-1, 130, 3, 32767, -32768, 6, 7, 8, 9, 1000, -1000, 12, 13, 14, -250, 16, 17, 200, 19, 20, 21, 22, 23, 24))

  # A float32 file and a uint16 pair file.
  f32 <- analyze_extra_file("vol_f32.hdr")
  u16 <- analyze_extra_file("vol_u16_pair.hdr")
  if (!is.null(f32)) {
    expect_equal(as.vector(read.fs.volume.analyze(f32)), analyze_reference_values / 4., tolerance = 1e-6)
  }
  if (!is.null(u16)) {
    # Values above 32767 have to survive, which requires reading them as unsigned 16 bit values.
    expect_equal(as.vector(read.fs.volume.analyze(u16)), c(40000, 65535, 1, 50000, 32768, 6, 7, 8, 9, 1000, 60000, 12, 13, 14, 33000, 16, 17, 200, 65534, 20, 21, 22, 23, 24))
  }

  # The file that FreeSurfer wrote from the 'tiny.mgh' volume that comes with the package: the values have to
  # match, and the SPM origin it stored has to be found.
  fs_tiny <- analyze_extra_file("fs_tiny.hdr")
  if (!is.null(fs_tiny)) {
    tiny_mgh <- system.file("extdata", "tiny.mgh", package = "freesurferformats", mustWork = TRUE)
    fs_vol <- read.fs.volume.analyze(fs_tiny, with_header = TRUE, spm = TRUE)
    expect_equal(as.numeric(fs_vol$data), as.numeric(read.fs.mgh(tiny_mgh)))
    expect_equal(fs_vol$header$spm_origin, c(2L, 2L, 2L))
    expect_equal(fs_vol$header$orient, 4L)
  }
})