tests/testthat/test-cifti-connectome-official.R

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

# Tests for the real connectome files in extra_test_data (item I.2g).
#
# These files are written by Connectome Workbench from the official Conte69 example files of
# the HCP, so unlike the tiny fixtures in inst/extdata they have *real* geometry: a reduced
# grayordinate mapping with non-contiguous vertex indices (30,424 of 32,492 left indices in
# the source file), thousands of vertices per parcel, and a matrix that was computed by the
# reference implementation. They are generated by
# dev_tools/generate_cifti_connectome_test_data.py and are only available in a git checkout
# of the repository, so the tests skip when they are missing.

cifti.test.official.file <- function(filename) {
  return(find_extra_test_data_file(file.path("cifti", filename)))
}

test_that("The reduced dense mapping of a real grayordinates file is read", {
  testthat::skip_on_cran()
  restricted_file <- cifti.test.official.file("conte69.restricted300.dtseries.nii")
  testthat::skip_if(is.null(restricted_file),
                    message = "extra_test_data not available in this environment (a git checkout of the repository is required).")

  cii <- read.cifti.header(restricted_file)
  expect_equal(cii$matrix$dim_sizes, c(2L, 306L))
  structures <- cifti.structures(cii, dim = 1L)
  expect_equal(structures$index_count, c(153L, 153L))
  expect_equal(structures$surface_number_of_vertices, c(32492L, 32492L))

  # The vertex indices are a subset of the surface, i.e. neither contiguous nor starting at
  # 0 for the right hemisphere (which is the real-world case that a naive reader gets wrong).
  grayordinates <- cifti.grayordinates(cii, dim = 1L)
  expect_equal(nrow(grayordinates), 306L)
  left_indices <- grayordinates$vertex_index[grayordinates$structure_short == "CORTEX_LEFT"]
  right_indices <- grayordinates$vertex_index[grayordinates$structure_short == "CORTEX_RIGHT"]
  expect_equal(length(left_indices), 153L)
  expect_true(all(diff(left_indices) > 1L))
  expect_equal(left_indices[1L], 0L)
  expect_true(all(right_indices <= 32491L))
  expect_equal(sum(is.na(grayordinates$vertex_index)), 0L)

  # The per-vertex readers reconstruct the data for the complete surface, with NA for the
  # vertices the file does not contain: this is the NA pattern of a real grayordinates file.
  morph_lh <- read.fs.morph.cifti(restricted_file, "lh")
  expect_length(morph_lh, 32492L)
  expect_equal(sum(!is.na(morph_lh)), 153L)
  expect_equal(which(!is.na(morph_lh)) - 1L, left_indices)
  series_lh <- read.fs.series.cifti(restricted_file, "lh")
  expect_equal(dim(series_lh), c(32492L, 2L))
  expect_equal(colnames(series_lh), c("0", "1"))
})

test_that("A real dense connectome is read", {
  testthat::skip_on_cran()
  restricted_file <- cifti.test.official.file("conte69.restricted300.dtseries.nii")
  dconn_file <- cifti.test.official.file("conte69.restricted300.dconn.nii")
  testthat::skip_if(is.null(restricted_file) || is.null(dconn_file),
                    message = "extra_test_data not available in this environment.")

  conn <- read.fs.connectome.cifti(dconn_file)
  expect_equal(dim(conn$data), c(306L, 306L))
  expect_true(is.na(conn$parcels_dim))
  expect_equal(conn$grayordinates_dim, 0L)
  expect_equal(nrow(conn$grayordinates), 306L)
  expect_true(isSymmetric(conn$data, tol = 1e-5))
  expect_true(all(diag(conn$data) > 0.999))

  # The values are the correlations that Workbench computed from the restricted file, and the
  # file stores both halves of the matrix and the diagonal.
  expected <- cor(read.cifti(restricted_file)$data)
  expect_equal(unname(conn$data), unname(expected), tolerance = 1e-5)

  # The row reader returns the same values (the rows of a connectome are the same
  # brainordinates as its columns).
  part <- read.cifti.rows(dconn_file, rows = c(1L, 306L), columns = c(1L, 2L))
  expect_equal(unname(part$data), unname(conn$data[c(1L, 306L), c(1L, 2L), drop = FALSE]))
})

test_that("A real parcellated connectome is read", {
  testthat::skip_on_cran()
  ptseries_file <- cifti.test.official.file("Conte69.MyelinAndCorrThickness.32k_fs_LR.ptseries.nii")
  pconn_file <- cifti.test.official.file("conte69.ptseries.corr.pconn.nii")
  testthat::skip_if(is.null(ptseries_file) || is.null(pconn_file),
                    message = "extra_test_data not available in this environment.")

  conn <- read.fs.connectome.cifti(pconn_file)
  expect_equal(dim(conn$data), c(54L, 54L))
  expect_equal(conn$parcels_dim, 0L)
  expect_equal(length(conn$parcel_names), 54L)
  expect_equal(conn$parcel_names[1L], "MEDIAL.WALL")
  # The parcels are the real Conte69 parcels, with their vertex lists for both hemispheres.
  expect_equal(conn$parcels$num_vertices[2L], 977L)
  expect_equal(unname(conn$data), unname(cor(read.cifti(ptseries_file)$data)), tolerance = 1e-5)
})

test_that("The real connectome files round trip through the writer", {
  testthat::skip_on_cran()
  dconn_file <- cifti.test.official.file("conte69.restricted300.dconn.nii")
  pconn_file <- cifti.test.official.file("conte69.ptseries.corr.pconn.nii")
  testthat::skip_if(is.null(dconn_file) || is.null(pconn_file),
                    message = "extra_test_data not available in this environment.")

  for (source_file in c(dconn_file, pconn_file)) {
    conn <- read.fs.connectome.cifti(source_file)
    out_file <- file.path(tempdir(), paste0("roundtrip_", basename(source_file)))
    write.fs.connectome.cifti(out_file, conn$data, template = source_file)
    back <- read.fs.connectome.cifti(out_file)
    expect_equal(unname(back$data), unname(conn$data), info = basename(source_file))
    expect_equal(back$parcel_names, conn$parcel_names, info = basename(source_file))
    expect_equal(nrow(back$grayordinates), nrow(conn$grayordinates), info = basename(source_file))
    expect_equal(back$header$niiheader$intent_code, conn$header$niiheader$intent_code,
                 info = basename(source_file))
    unlink(out_file)
  }
})