tests/testthat/test-read_fs_mgh.R

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

test_that("The dimensions of our demo 3D MGZ file are read correctly", {
  brain_image <- system.file("extdata", "brain.mgz", package = "freesurferformats", mustWork = TRUE)
  vd <- read.fs.mgh(brain_image)

  expect_equal(class(vd), "array")
  expect_equal(length(dim(vd)), 4) # It has 4 dimensions
  expect_equal(dim(vd), c(256, 256, 256, 1))
})

test_that("The data values in the demo 3D MGZ file are read as in the reference implementation", {
  # Tests that the data read is identical to the values returned by the matlab function $FREESURFER_HOME/matlab/MRIread.m
  # You can also check with the `mri_info` command line tool, just be sure to shift the
  # vertex indices (zero-based for mri_info, one-based in R). See comment below for an example.

  brain_image <- system.file("extdata", "brain.mgz", package = "freesurferformats", mustWork = TRUE)
  vd <- read.fs.mgh(brain_image)

  expect_equal(class(vd), "array")
  expect_equal(length(dim(vd)), 4) # It has 4 dimensions
  expect_equal(dim(vd), c(256, 256, 256, 1))

  expect_equal(vd[100, 100, 100, 1], 77) # try on command line: mri_info --voxel 99 99 99 inst/extdata/brain.mgz
  expect_equal(vd[110, 110, 110, 1], 71)
  expect_equal(vd[1, 1, 1, 1], 0)

  expect_equal(sum(vd), 121035479)
})


test_that("The demo 1D MGZ file is read correctly", {
  morph_file <- system.file("extdata", "lh.curv.fwhm10.fsaverage.mgz", package = "freesurferformats", mustWork = TRUE)
  vd <- read.fs.mgh(morph_file)

  expect_equal(class(vd), "array")
  expect_equal(length(dim(vd)), 4) # It has 4 dimensions
  num_vers_fsaverage <- 163842
  expect_equal(dim(vd), c(num_vers_fsaverage, 1, 1, 1))

  vd2 <- read.fs.volume(morph_file)
  expect_equal(class(vd2), "array")
  expect_equal(length(dim(vd2)), 4) # It has 4 dimensions
  expect_equal(dim(vd2), c(num_vers_fsaverage, 1, 1, 1))
})


test_that("The header can be read", {
  brain_image <- system.file("extdata", "brain.mgz", package = "freesurferformats", mustWork = TRUE)
  ret <- read.fs.mgh(brain_image, with_header = TRUE)
  expect_true(is.fs.volume(ret))

  ret2 <- read.fs.volume(brain_image, with_header = TRUE)
  expect_true(is.fs.volume(ret2))

  header <- ret$header
  expect_true(is.mghheader(header))
  expect_equal(header$dtype, 0) # MRI_UCHAR
  expect_equal(header$dof, 0)
  expect_equal(header$ras_good_flag, 1)
  expect_equal(length(header$internal$delta), 3)
  expect_equal(header$internal$delta, c(1, 1, 1), tolerance = 1e-2)
  expect_equal(length(header$internal$Mdc), 9)
  expect_equal(header$internal$Mdc, matrix(c(-1, 0, 0, 0, 0, -1, 0, 1, 0), nrow = 3), tolerance = 1e-2)
  expect_equal(length(header$internal$Pxyz_c), 3)
  expect_equal(header$internal$Pxyz_c, c(-0.5, 29.4, -48.9), tolerance = 1e-2)
  expect_equal(header$voldim, c(256, 256, 256, 1))
  expect_equal(length(header$has_mr_params), 1)
  expect_equal(length(header$mr_params), 5)
  expect_equal(header$mr_params, c(2300.000000, 0.157080, 2.010000, 900.000000, 256.000), tolerance = 1e-2)
  expect_equal(length(header$internal$D), 9) # 3x3
  expect_equal(length(header$internal$Pcrs_c), 3) # 3x1
  expect_equal(length(header$internal$Pxyz_0), 3) # 3x1
  expect_equal(length(header$vox2ras_matrix), 16)
  expect_equal(header$vox2ras_matrix, matrix(c(-1, 0, 0, 0, 0, 0, -1, 0, 0, 1, 0, 0, 127.5, -98.6273, 79.0953, 1.000), nrow = 4, byrow = FALSE), tolerance = 1e-2)

  header2 <- ret2$header
  expect_true(is.mghheader(header2))
  expect_equal(header2$dtype, 0) # MRI_UCHAR
  expect_equal(header2$dof, 0)
  expect_equal(header2$ras_good_flag, 1)
  expect_equal(length(header2$internal$delta), 3)


  vd <- ret$data
  expect_equal(class(vd), "array")
  expect_equal(length(dim(vd)), 4) # It has 4 dimensions
  expect_equal(typeof(vd), "integer")
  expect_equal(dim(vd), c(256, 256, 256, 1))
  expect_equal(vd[80, 80, 80, 1], 110) # value of voxel 80,80,80 must be 110 (known from ref. implementation)
  expect_equal(vd[100, 100, 100, 1], 77)
  expect_equal(vd[100, 100, 80, 1], 105)

  vd2 <- ret2$data
  expect_equal(class(vd2), "array")
  expect_equal(length(dim(vd2)), 4) # It has 4 dimensions
  expect_equal(typeof(vd2), "integer")
  expect_equal(dim(vd2), c(256, 256, 256, 1))
  expect_equal(vd2[80, 80, 80, 1], 110) # value of voxel 80,80,80 must be 110 (known from ref. implementation)
  expect_equal(vd2[100, 100, 100, 1], 77)
  expect_equal(vd2[100, 100, 80, 1], 105)
})


test_that("Data is flattened when requested", {
  brain_image <- system.file("extdata", "brain.mgz", package = "freesurferformats", mustWork = TRUE)
  vd <- read.fs.mgh(brain_image, flatten = TRUE)

  expect_equal(class(vd), "integer")
  expect_equal(length((vd)), 16777216)
})


test_that("The gzip status is guessed as expected from a filename", {
  expect_equal(guess.filename.is.gzipped("noway"), FALSE)
  expect_equal(guess.filename.is.gzipped("/there/is/noway"), FALSE)
  expect_equal(guess.filename.is.gzipped("file.tar.gz"), TRUE)
  expect_equal(guess.filename.is.gzipped("relative/path/to/file.tar.gz"), TRUE)
  expect_equal(guess.filename.is.gzipped("file.tar"), FALSE)
  expect_equal(guess.filename.is.gzipped("brain.mgz"), TRUE)
  expect_equal(guess.filename.is.gzipped("brain.mgh"), FALSE)
})


test_that("A real MGH can be read, rewritten, read again, and the data and header are preserved as expected.", {
  brain_image <- system.file("extdata", "brain.mgz", package = "freesurferformats", mustWork = TRUE)
  ret_orig <- read.fs.mgh(brain_image, with_header = TRUE)
  expect_true(is.fs.volume(ret_orig))

  new_copy <- tempfile(fileext = "mgz")
  write.fs.mgh(new_copy, ret_orig$data, vox2ras_matrix = ret_orig$header$vox2ras_matrix, mr_params = ret_orig$header$mr_params)

  ret <- read.fs.mgh(new_copy, with_header = TRUE)
  header <- ret$header
  expect_true(is.mghheader(header))
  expect_equal(header$dtype, 1) # IMPORTANT: The data type will have changed from MRI_UCHAR to MRI_INTEGER. This is fine with us for now.
  expect_equal(header$dof, 0)
  expect_equal(header$ras_good_flag, 1)
  expect_equal(length(header$internal$delta), 3)
  expect_equal(header$internal$delta, c(1, 1, 1), tolerance = 1e-2)
  expect_equal(length(header$internal$Mdc), 9)
  expect_equal(header$internal$Mdc, matrix(c(-1, 0, 0, 0, 0, -1, 0, 1, 0), nrow = 3), tolerance = 1e-2)
  expect_equal(length(header$internal$D), 9) # 3x3
  expect_equal(length(header$internal$Pcrs_c), 3) # 3x1
  expect_equal(length(header$internal$Pxyz_0), 3) # 3x1
  expect_equal(header$internal$Pxyz_c, c(-0.5, 29.4, -48.9), tolerance = 1e-2)
  expect_equal(header$voldim, c(256, 256, 256, 1))
  expect_equal(length(header$has_mr_params), 1)
  expect_equal(length(header$mr_params), 5)
  expect_equal(header$mr_params, c(2300.000000, 0.157080, 2.010000, 900.000000, 256.000), tolerance = 1e-2)
  expect_equal(header$vox2ras_matrix, matrix(c(-1, 0, 0, 0, 0, 0, -1, 0, 0, 1, 0, 0, 127.5, -98.6273, 79.0953, 1.000), nrow = 4, byrow = FALSE), tolerance = 1e-2)


  vd <- ret$data
  expect_equal(class(vd), "array")
  expect_equal(length(dim(vd)), 4) # It has 4 dimensions
  expect_equal(dim(vd), c(256, 256, 256, 1))
  expect_equal(vd[80, 80, 80, 1], 110) # value of voxel 80,80,80 must be 110 (known from ref. implementation)
  expect_equal(vd[100, 100, 100, 1], 77)
  expect_equal(vd[100, 100, 80, 1], 105)
})


test_that("Invalid volume format leads to errors", {
  brain_image <- system.file("extdata", "brain.mgz", package = "freesurferformats", mustWork = TRUE)
  expect_error(read.fs.volume(brain_image, format = "invalid format")) # invalid format
})


test_that("Trying to read a non-existent file that does not match any of the file extension checks with read.fs.volume throws an error (issue #25)", {
  non_existent_file <- tempfile(, fileext = "")
  testthat::expect_false(file.exists(non_existent_file)) # ensure it does not exist.

  testthat::expect_error(read.fs.volume(non_existent_file))
})