tests/testthat/test-write_fs_mgh.R

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

test_that("An one-dimensional MGH file of double values can be written", {
  output_file <- tempfile()

  # generate data
  data_length <- 149244
  data <- rep(1.25, data_length)
  data_array <- array(data)

  # write data to file
  write.fs.mgh(output_file, data_array)

  # load data again and check it
  mgh <- read.fs.mgh(output_file, with_header = TRUE)
  read_data <- mgh$data
  header <- mgh$header
  expect_equal(header$mr_params, c(0, 0, 0, 0, 0))
  expect_equal(header$dtype, 3) # MRI_FLOAT
  expect_equal(header$ras_good_flag, -1L)

  # Check data dimensions: we should get 4 dimensions back. Since we wrote a vector of data only, the last 3 should be 1.
  expect_equal(typeof(read_data), "double")
  expect_equal(length(dim(read_data)), 4L)
  expect_equal((dim(read_data)[1]), data_length)
  expect_equal((dim(read_data)[2]), 1L)
  expect_equal((dim(read_data)[3]), 1L)
  expect_equal((dim(read_data)[4]), 1L)

  # Check the data values
  expect_equal(length(data), length(as.vector(read_data)))
  expect_equal(data_length, length(as.vector(read_data)))
  expect_equal(data, as.vector(read_data))
})


test_that("A three-dimensional MGH file with realistic size can be written", {
  output_file <- tempfile()

  # generate data
  data_array <- array(data = rep(1, 256 * 256 * 256), dim = c(256, 256, 256))

  # write data to file
  write.fs.mgh(output_file, data_array)

  # load data again and check it
  read_data <- read.fs.mgh(output_file)

  # Check data dimensions: we should get 4 dimensions back. Since we wrote a vector of data only, the last 3 should be 1.
  expect_equal(length(dim(read_data)), 4)
  expect_equal((dim(read_data)[1]), 256)
  expect_equal((dim(read_data)[2]), 256)
  expect_equal((dim(read_data)[3]), 256)
  expect_equal((dim(read_data)[4]), 1)

  # Check the data values
  expect_equal(length(data_array), length(read_data))
  expect_equal(data_array, read_data[, , , 1])
})


test_that("An one-dimensional uncompressed MGH file of integer values can be written", {
  output_file <- tempfile(fileext = ".mgh")

  # generate data
  data_length <- 111111L
  data <- rep.int(1L, data_length)
  data[10000:20000] <- 0L # Set some of the values to zero.
  expect_true(is.integer(data))
  expect_equal(range(data), c(0L, 1L))

  # write data to file
  write.fs.mgh(output_file, data)

  # load data again and check it
  mgh <- read.fs.mgh(output_file, with_header = TRUE)
  read_data <- mgh$data
  header <- mgh$header
  expect_equal(header$mr_params, c(0, 0, 0, 0, 0))
  expect_equal(header$dtype, 1) # MRI_INT
  expect_equal(header$ras_good_flag, -1L)

  # Check data dimensions: we should get 4 dimensions back. Since we wrote a vector of data only, the last 3 should be 1.
  expect_equal(length(dim(read_data)), 4L)
  expect_equal(typeof(read_data), "integer")
  expect_equal((dim(read_data)[1]), data_length)
  expect_equal((dim(read_data)[2]), 1L)
  expect_equal((dim(read_data)[3]), 1L)
  expect_equal((dim(read_data)[4]), 1L)

  # Check the data values
  expect_equal(length(data), length(as.vector(read_data)))
  expect_equal(data_length, length(as.vector(read_data)))
  expect_equal(data, as.vector(read_data))
  expect_equal(range(read_data), c(0L, 1L))
})


test_that("An one-dimensional compressed MGZ file of integer values can be written", {
  # This test is very similar to the last one, but it uses MGZ instead of MGH.
  output_file <- tempfile(fileext = ".mgz")

  # generate data
  data_length <- 111111L
  data <- rep.int(1L, data_length)
  data[10000:20000] <- 0L # Set some of the values to zero.
  expect_true(is.integer(data))
  expect_equal(range(data), c(0L, 1L))

  # write data to file
  write.fs.mgh(output_file, data)

  # load data again and check it
  mgh <- read.fs.mgh(output_file, with_header = TRUE)
  read_data <- mgh$data
  header <- mgh$header
  expect_equal(header$mr_params, c(0, 0, 0, 0, 0))
  expect_equal(header$dtype, 1) # MRI_INT
  expect_equal(header$ras_good_flag, -1L)

  # Check data dimensions: we should get 4 dimensions back. Since we wrote a vector of data only, the last 3 should be 1.
  expect_equal(length(dim(read_data)), 4L)
  expect_equal(typeof(read_data), "integer")
  expect_equal((dim(read_data)[1]), data_length)
  expect_equal((dim(read_data)[2]), 1L)
  expect_equal((dim(read_data)[3]), 1L)
  expect_equal((dim(read_data)[4]), 1L)

  # Check the data values
  expect_equal(length(data), length(as.vector(read_data)))
  expect_equal(data_length, length(as.vector(read_data)))
  expect_equal(data, as.vector(read_data))
  expect_equal(range(read_data), c(0L, 1L))
})


test_that("Logical data can be written and re-read from MGH and MGZ files", {
  # This test is very similar to the last one, but it uses MGZ instead of MGH.
  output_file_mgh <- tempfile(fileext = ".mgh")
  output_file_mgz <- tempfile(fileext = ".mgz")

  # generate data
  data_int <- sample(0:1, 2000, replace = TRUE) # only zeroes and ones
  data_logical <- as.logical(data_int)

  # write data to files, re-read and check
  write.fs.mgh(output_file_mgh, data_logical)
  mgh <- read.fs.mgh(output_file_mgh, with_header = TRUE, drop_empty_dims = TRUE)
  expect_equal(mgh$data, data_int)
  expect_equal(mgh$header$dtype, translate.mri.dtype("MRI_UCHAR"))

  write.fs.mgh(output_file_mgz, data_logical, mr_params = c(0., 0., 0., 0.))
  mgz <- read.fs.mgh(output_file_mgz, with_header = TRUE, drop_empty_dims = TRUE)
  expect_equal(mgz$data, data_int)
  expect_equal(mgz$header$dtype, translate.mri.dtype("MRI_UCHAR"))

  # Also write integer data as MRI_SHORT and _UCHAR for completeness
  write.fs.mgh(tempfile(fileext = ".mgz"), data_int, mri_dtype = "MRI_SHORT")
  write.fs.mgh(tempfile(fileext = ".mgz"), data_int, mri_dtype = "MRI_UCHAR")

  # There should be warning for unsiotable data range / dtype combinations
  data_int_large <- c(12L, 288L, 343456L) # integer data with values outside of range for MRI_SHORT and _UCHAR
  expect_warning(write.fs.mgh(tempfile(fileext = ".mgz"), data_int_large, mri_dtype = "MRI_SHORT"))
  expect_warning(write.fs.mgh(tempfile(fileext = ".mgz"), data_int_large, mri_dtype = "MRI_UCHAR"))
  # ... and for invalid input data type.
  data_char <- rep("hi", 42L)
  expect_warning(write.fs.mgh(tempfile(fileext = ".mgz"), data_char, mri_dtype = "MRI_SHORT"))
  expect_warning(write.fs.mgh(tempfile(fileext = ".mgz"), data_char, mri_dtype = "MRI_UCHAR"))
  expect_warning(write.fs.mgh(tempfile(fileext = ".mgz"), data_char, mri_dtype = "MRI_INT"))
  expect_warning(write.fs.mgh(tempfile(fileext = ".mgz"), data_char, mri_dtype = "MRI_FLOAT"))
})


test_that("Improper use of write.fs.mgh leads to errors", {
  data_int <- sample(0:1, 2000, replace = TRUE)
  filepath <- tempfile(fileext = ".mgh")

  expect_error(write.fs.mgh(123, data_int)) # 123 is not a valid file path (not a string)
  expect_error(write.fs.mgh(filepath, data_int, mr_params = c(0., 0))) # invalid length of mr_params
  expect_error(write.fs.mgh(filepath, data_int, mr_params = "hi")) # invalid mr_params, must be double
  expect_error(write.fs.mgh(filepath, data_int, mr_params = "what")) # mr_params must be double
  expect_error(write.fs.mgh(filepath, data_int, vox2ras_matrix = "what")) # vox2ras_matrix must be double matrix
  expect_error(write.fs.mgh(filepath, data_int, vox2ras_matrix = matrix(seq(6), nrow = 2))) # vox2ras_matrix must be a 4x4 matrix
  expect_error(write.fs.mgh(filepath, data_int, mri_dtype = 3)) # mri_dtype must be a character string
  expect_error(write.fs.mgh(filepath, data_int, mri_dtype = "no such MRI dtype I guess")) # mri_dtype must be valid
  expect_error(write.fs.mgh(filepath, data = NULL)) # invalid data
  expect_error(write.fs.mgh(filepath, data = array(seq.int(6), dim = rep(1L, 6L)))) # invalid data: too many dimensions (>5)
  expect_error(write.fs.mgh(filepath, data = rep("hi", 100L))) # invalid data: character strings not supported, no suitable MRI_DTYPE
})


test_that("An MGH volume written with an anisotropic, rotated vox2ras matrix reads back with the same geometry", {
  # Regression test: the direction cosines (Mdc) written to the MGH header must be the columns of
  # the vox2ras linear part normalized by their own voxel size. A previous version normalized them
  # with the voxel sizes misaligned (row-wise instead of column-wise), which corrupted the stored
  # orientation for anisotropic volumes whose axes are not aligned with the coordinate axes.
  data <- array(as.integer(sample(0:255, 4 * 5 * 6, replace = TRUE)), dim = c(4, 5, 6, 1))

  # A rotation about the z axis combined with anisotropic voxel sizes (1, 2, 3) mm. The columns of
  # the resulting linear part have norms 1, 2 and 3 and contain off-axis entries, which exposes
  # misaligned column normalization.
  rot_z <- matrix(c(0, -1, 0, 1, 0, 0, 0, 0, 1), nrow = 3, byrow = TRUE)
  vox2ras <- diag(4)
  vox2ras[1:3, 1:3] <- rot_z %*% diag(c(1, 2, 3))
  vox2ras[1:3, 4] <- c(10, -20, 30)

  output_file <- tempfile(fileext = ".mgz")
  write.fs.mgh(output_file, data, vox2ras_matrix = vox2ras)

  mgh <- read.fs.mgh(output_file, with_header = TRUE)
  expect_equal(mgh$data, data)
  expect_equal(mghheader.vox2ras(mgh), vox2ras, tolerance = 1e-4)

  # And a second, differently signed/rotated anisotropic case for good measure.
  rot_z2 <- matrix(c(0, 1, 0, -1, 0, 0, 0, 0, -1), nrow = 3, byrow = TRUE)
  vox2ras2 <- diag(4)
  vox2ras2[1:3, 1:3] <- rot_z2 %*% diag(c(1.5, 2.0, 3.5))
  vox2ras2[1:3, 4] <- c(5, 6, 7)
  output_file2 <- tempfile(fileext = ".mgz")
  write.fs.mgh(output_file2, data, vox2ras_matrix = vox2ras2)
  mgh2 <- read.fs.mgh(output_file2, with_header = TRUE)
  expect_equal(mghheader.vox2ras(mgh2), vox2ras2, tolerance = 1e-4)
})