# Tests for mghheader functions

testthat::test_that("The native vox2ras matrix can be computed from a conformed volume", {
  brain_image = system.file("extdata", "brain.mgz", package = "freesurferformats", mustWork = TRUE);
  mgh = read.fs.mgh(brain_image, with_header=TRUE);
  testthat::expect_true(is.fs.volume(mgh));

  # Run `mri_info --vox2ras <volume>` on the OS command line to find the reference value.
  known_vox2ras = matrix(c(-1, 0, 0, 0, 0, 0, -1, 0, 0, 1, 0, 0, 127.5005, -98.62726, 79.09527, 1), nrow=4, byrow = FALSE);

  vox2ras = mghheader.vox2ras(mgh$header);
  testthat::expect_equal(vox2ras, known_vox2ras, tolerance=1e-4);
})


testthat::test_that("The native ras2vox matrix can be computed from a conformed volume", {
  brain_image = system.file("extdata", "brain.mgz", package = "freesurferformats", mustWork = TRUE);
  mgh = read.fs.mgh(brain_image, with_header=TRUE);
  testthat::expect_true(is.fs.volume(mgh));

  # Run `mri_info --ras2vox <volume>` on the OS command line to find the reference value.
  known_ras2vox = matrix(c(-1, 0, 0, 0, 0, 0, 1, 0, 0, -1, 0, 0, 127.5005, 79.09527, 98.62726, 1), nrow=4, byrow = FALSE);

  ras2vox = mghheader.ras2vox(mgh$header);
  testthat::expect_equal(ras2vox, known_ras2vox, tolerance=1e-4);
})


testthat::test_that("The tkregister vox2ras matrix can be computed from a conformed volume", {
  brain_image = system.file("extdata", "brain.mgz", package = "freesurferformats", mustWork = TRUE);
  mgh = read.fs.mgh(brain_image, with_header=TRUE);
  testthat::expect_true(is.fs.volume(mgh));

  # Run `mri_info --vox2ras-tkr <volume>` on the OS command line to find the reference value.
  known_vox2ras_tkr = matrix(c(-1, 0, 0, 0, 0, 0, -1, 0, 0, 1, 0, 0, 128., -128., 128., 1), nrow=4, byrow = FALSE);

  vox2ras_tkr = mghheader.vox2ras.tkreg(mgh$header);
  testthat::expect_equal(vox2ras_tkr, known_vox2ras_tkr, tolerance=1e-4);
})


testthat::test_that("The tkregister ras2vox matrix can be computed from a conformed volume", {
  brain_image = system.file("extdata", "brain.mgz", package = "freesurferformats", mustWork = TRUE);
  mgh = read.fs.mgh(brain_image, with_header=TRUE);
  testthat::expect_true(is.fs.volume(mgh));

  # Run `mri_info --ras2vox-tkr <volume>` on the OS command line to find the reference value.
  known_ras2vox_tkr = matrix(c(-1, 0, 0, 0, 0, 0, 1, 0, 0, -1, 0, 0, 128., 128., 128., 1), nrow=4, byrow = FALSE);

  ras2vox_tkr = mghheader.ras2vox.tkreg(mgh$header);
  testthat::expect_equal(ras2vox_tkr, known_ras2vox_tkr, tolerance=1e-4);
})


testthat::test_that("The slice direction and orientation can be computed from full fs.volume", {
  brain_image = system.file("extdata", "brain.mgz", package = "freesurferformats", mustWork = TRUE);
  mgh = read.fs.mgh(brain_image, with_header=TRUE);
  testthat::expect_true(mghheader.is.conformed(mgh));
  testthat::expect_equal(mghheader.primary.slice.direction(mgh), 'coronal');
  testthat::expect_equal(mghheader.crs.orientation(mgh), 'LIA');

  header_noras = mgh$header;
  header_noras$ras_good_flag = 0L;
  testthat::expect_false(mghheader.is.conformed(header_noras));
})


testthat::test_that("The slice direction and orientation can be computed from the volume header", {
  brain_image = system.file("extdata", "brain.mgz", package = "freesurferformats", mustWork = TRUE);
  mgh = read.fs.mgh(brain_image, with_header=TRUE);
  header = mgh$header;
  testthat::expect_true(mghheader.is.conformed(header));
  testthat::expect_equal(mghheader.primary.slice.direction(header), 'coronal');
  testthat::expect_equal(mghheader.crs.orientation(header), 'LIA');
})


testthat::test_that("The mgh header fields can be initialized from a vox2ras matrix", {
  brain_image = system.file("extdata", "brain.mgz", package = "freesurferformats", mustWork = TRUE);
  mgh = read.fs.mgh(brain_image, with_header=TRUE);
  header = mgh$header;

  vox2ras = mghheader.vox2ras(header);
  new_header = mghheader.update.from.vox2ras(header, vox2ras);   # Update with the matrix it already has, i.e., set basic fields from matrix.
  new_vox2ras = mghheader.vox2ras(new_header);                   # Recompute matrix from set basic fields.

  testthat::expect_equal(vox2ras, new_vox2ras);                            # This should lead to the same matrix.
})


testthat::test_that("The RAS coordinate of the center voxel (CRAS) can be computed from the RAS coordinate for the first voxel (P0 RAS).", {
  brain_image = system.file("extdata", "brain.mgz", package = "freesurferformats", mustWork = TRUE);
  mgh = read.fs.mgh(brain_image, with_header=TRUE);
  header = mgh$header;

  vox2ras = mghheader.vox2ras(header);
  first_voxel_CRS = c(1L, 1L, 1L);
  first_voxel_RAS = vox2ras[,4][0:3];  # the P0 RAS
  known_first_voxel_RAS = c(127.5, -98.6273, 79.0953);   # known from: `mri_info --p0  path/to/brain.mgz` (in system shell)
  testthat::expect_equal(first_voxel_RAS, known_first_voxel_RAS, tolerance=1e-2);

  center_RAS = mghheader.centervoxelRAS.from.firstvoxelRAS(header, first_voxel_RAS);  # center RAS is also known as 'CRAS'.

  known_center_RAS = c(-0.499954, 29.3727, -48.9047);     # known from: `mri_info --cras  path/to/brain.mgz` (in system shell)

  testthat::expect_equal(center_RAS, known_center_RAS, tolerance=1e-2);
})


testthat::test_that("MGH header can be checked for valid RAS", {
  brain_image = system.file("extdata", "brain.mgz", package = "freesurferformats", mustWork = TRUE);
  mgh = read.fs.mgh(brain_image, with_header=TRUE);
  header = mgh$header;

  testthat::expect_true(mghheader.is.ras.valid(header));
  testthat::expect_true(mghheader.is.ras.valid(mgh));

  testthat::expect_false(mghheader.is.ras.valid(NULL));
  testthat::expect_false(mghheader.is.ras.valid("what"));

  header_noras = header;
  header_noras$ras_good_flag = 0L;
  testthat::expect_false(mghheader.is.ras.valid(header_noras));

  header_noras$ras_good_flag = NULL;
  testthat::expect_false(mghheader.is.ras.valid(header_noras));
})


testthat::test_that("MGH header can be used to compute tkreg2scanner", {
  brain_image = system.file("extdata", "brain.mgz", package = "freesurferformats", mustWork = TRUE);
  mgh = read.fs.mgh(brain_image, with_header=TRUE);
  header = mgh$header;

  tkm = mghheader.tkreg2scanner(header);
  tkm2 = mghheader.tkreg2scanner(mgh);

  header_noras = header;
  header_noras$ras_good_flag = 0L;
  testthat::expect_error(mghheader.tkreg2scanner(header_noras));
})


testthat::test_that("MGH header can be used to compute scanner2tkreg", {
  brain_image = system.file("extdata", "brain.mgz", package = "freesurferformats", mustWork = TRUE);
  mgh = read.fs.mgh(brain_image, with_header=TRUE);
  header = mgh$header;

  tkm = mghheader.scanner2tkreg(header);
  tkm2 = mghheader.scanner2tkreg(mgh);

  header_noras = header;
  header_noras$ras_good_flag = 0L;
  testthat::expect_error(mghheader.scanner2tkreg(header_noras));
})


testthat::test_that("Two MGH headers can be used to compute vox2vox", {
  brain_image = system.file("extdata", "brain.mgz", package = "freesurferformats", mustWork = TRUE);
  mgh = read.fs.mgh(brain_image, with_header=TRUE);
  header = mgh$header;
  header2 = header;

  v2v = mghheader.vox2vox(header, header2);
  v2v = mghheader.vox2vox(mgh, mgh);

  header_noras = header;
  header_noras$ras_good_flag = 0L;
  testthat::expect_error(mghheader.vox2vox(header_noras, header));
  testthat::expect_error(mghheader.vox2vox(header, header_noras));
  testthat::expect_error(mghheader.vox2vox(header_noras, header_noras));
})


testthat::test_that("Transformation matrices can be adapted to 0-based and 1-based indices", {
  brain_image = system.file("extdata", "brain.mgz", package = "freesurferformats", mustWork = TRUE);
  mgh = read.fs.mgh(brain_image, with_header=TRUE);
  vox2ras = mghheader.vox2ras(mgh$header);

  vox2ras_for_1based = sm0to1(vox2ras);
  vox2ras_for_0based = sm1to0(vox2ras_for_1based);
  testthat::expect_equal(vox2ras, vox2ras_for_0based, tolerance=1e-2);
})


testthat::test_that("Transformation from surface RAS to RAS and back works", {
  brain_image = system.file("extdata", "brain.mgz", package = "freesurferformats", mustWork = TRUE);
  mgh = read.fs.mgh(brain_image, with_header=TRUE);
  mghheader = mgh$header;

  # get surface RAS data: we use mesh coords of a surface.
  surface_file = system.file("extdata", "lh.tinysurface", package = "freesurferformats", mustWork = TRUE);
  surface = read.fs.surface(surface_file);

  ras = surfaceras.to.ras(mghheader, surface$vertices);
  # transform back
  sras = ras.to.surfaceras(mghheader, ras);
  testthat::expect_equal(surface$vertices, sras);
})


testthat::test_that("Transformation from Talairach RAS to RAS and back works", {
  brain_image = system.file("extdata", "brain.mgz", package = "freesurferformats", mustWork = TRUE);
  mgh = read.fs.mgh(brain_image, with_header=TRUE);
  mghheader = mgh$header;

  # get RAS data: we use mesh coords of a surface (in sras) and transform to RAS.
  surface_file = system.file("extdata", "lh.tinysurface", package = "freesurferformats", mustWork = TRUE);
  surface = read.fs.surface(surface_file);
  ras = surfaceras.to.ras(mghheader, surface$vertices);

  # transform to Talairach
  talairach = system.file("extdata", "talairach.xfm", package = "freesurferformats", mustWork = TRUE);
  talras = ras.to.talairachras(ras, talairach = talairach);

  # transform back
  orig_ras = talairachras.to.ras(talras, talairach = talairach);
  testthat::expect_equal(ras, orig_ras);
})


testthat::test_that("Transformation from surface vertex coord to Talairach space works", {
  brain_image = system.file("extdata", "brain.mgz", package = "freesurferformats", mustWork = TRUE);
  mgh = read.fs.mgh(brain_image, with_header=TRUE);
  mghheader = mgh$header;

  # get RAS data: we use mesh coords of a surface (in sras) and transform to RAS.
  surface_file = system.file("extdata", "lh.tinysurface", package = "freesurferformats", mustWork = TRUE);
  surface = read.fs.surface(surface_file);
  ras = surfaceras.to.ras(mghheader, surface$vertices);

  # transform to Talairach
  talairach = system.file("extdata", "talairach.xfm", package = "freesurferformats", mustWork = TRUE);
  talras = surfaceras.to.talairach(surface$vertices, header_cras = mghheader, talairach = talairach);
  testthat::expect_equal(nrow(ras), nrow(talras));
})