---
title: "Read imaging data"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Read imaging data}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include=FALSE, echo = FALSE, results='hide'}
if (identical(Sys.getenv("IEEGIO_PKGDOWN", unset = ""), "")) {
  knitr::opts_chunk$set(
    collapse = TRUE,
    comment = "#>",
    eval = FALSE
  )
} else {
  options("rgl.useNULL" = TRUE)
  library(ieegio)
  library(rgl)
  knitr::opts_chunk$set(
    collapse = TRUE,
    comment = "#>"
  )
  rgl::setupKnitr(autoprint = FALSE)
  cache_dir <- tools::R_user_dir("ieegio", "cache")
  dir.create(cache_dir, showWarnings = FALSE, recursive = TRUE)
  cache_dir <- normalizePath(cache_dir, mustWork = TRUE)
  options("ieegio.extract_path", cache_dir)
}

```

`ieegio` supports reading from and writing to multiple imaging formats:

* Volume: `NIfTI` & `FreeSurfer MGH/MGZ`
* Surface: `GIfTI` & `FreeSurfer` geometry, annotation, curvature/measurement, `w` format, `AFNI`/`SUMA` `NIML`, `VTK` polygon meshes
* Streamlines: `TRK`, `TCK`, `TT` (read-only), `VTK`, `VTP`, ...

Beyond reading and writing, `ieegio` can also map data between volumes and surfaces, chain coordinate transforms, and compare regions of interest; those are covered at the end of this article.

To start, please load `ieegio`. This vignette uses sample data which requires extra download.


```{r sample}
library(ieegio)

# volume file
nifti_file <- ieegio_sample_data("brain.demosubject.nii.gz")

# geometry
geom_file <- ieegio_sample_data(
  "gifti/icosahedron3d/geometry.gii")

# measurements
shape_file <- ieegio_sample_data(
  "gifti/icosahedron3d/rand.gii"
)

# time series
ts_file <- ieegio_sample_data(
  "gifti/icosahedron3d/ts.gii")

# streamlines
trk_file <- ieegio_sample_data(
  "streamlines/CNVII_R.trk")

tck_file <- ieegio_sample_data(
  "streamlines/CNVII_R.tck")

tt_file <- ieegio_sample_data(
  "streamlines/CNVII_R.tt")

# AFNI/SUMA std.141 geometry and matching annotation
std141_geom_file <- ieegio_sample_data(
  "gifti/std.141.lh.inf_200.gii")

niml_file <- ieegio_sample_data(
  "niml/std.141.lh.aparc.a2009s.annot.niml.dset")

# volumetric atlas
atlas_file <- ieegio_sample_data(
  "atlases/YBA/YBA690.nii.gz")
```


## Volume files

`ieegio::read_volume` and `ieegio::write_volume` provides high-level interfaces for reading and writing volume data such as `MRI`, `CT`. `fMRI`, etc.

Each volume data (`NIfTI`, `MGH`, ...) contains a `header`, a `data`, and a `transforms` list.

```{r read_volume}
volume <- read_volume(nifti_file)
volume
```

The transforms contain transforms from volume (column, row, slice) index to other coordinate systems. The most commonly used one is `vox2ras`, which is a `4x4` matrix mapping the voxels to scanner (usually `T1-weighted`) `RAS` (right-anterior-superior) system.

Accessing the image values via `[` operator. For example,

``` r
volume[128, , ]
```

Plotting the anatomical slices:

```{r plot_volume, fig.width=8, fig.height=3, out.width="95%"}
par(mfrow = c(1, 3), mar = c(0, 0, 3.1, 0))

ras_position <- c(-50, -10, 15)

ras_str <- paste(sprintf("%.0f", ras_position), collapse = ",")

for (which in c("coronal", "axial", "sagittal")) {
  plot(x = volume, position = ras_position, crosshair_gap = 10,
       crosshair_lty = 2, zoom = 3, which = which,
       main = sprintf("%s T1RAS=[%s]", which, ras_str))
}

```

## Surface files

Reading surface file using `read_surface` supports multiple data types


```{r read_surface}
library(ieegio)
# geometry
geometry <- read_surface(geom_file)

# measurements
measurement <- read_surface(shape_file)

# time series
time_series <- read_surface(ts_file)
```

You can merge them to a single object, making an object with multiple embedding data-sets:

```{r merge}
merged <- merge(geometry, measurement, time_series)
print(merged)
```

Plot the surfaces in `3D` viewer, colored by shape measurement

```{r plot_surface, webgl = TRUE, out.width="70%", fig.width = 7}
# plot the first column in measurements section
plot(merged, name = list("measurements", 1))
```

Plot the normalized time-series data

```{r time_series, webgl = TRUE, out.width="100%", fig.width = 7}
ts_demean <- apply(
  merged$time_series$value,
  MARGIN = 1L,
  FUN = function(x) {
    x - mean(x)
  }
)
merged$time_series$value <- t(ts_demean)
plot(
  merged, name = "time_series",
  col = c(
    "#053061", "#2166ac", "#4393c3",
    "#92c5de", "#d1e5f0", "#ffffff",
    "#fddbc7", "#f4a582", "#d6604d",
    "#b2182b", "#67001f"
  )
)
```

# Streamline files

Reading streamlines via universal entry function `read_streamlines`

```{r read_streamlines, results='hide'}
trk <- read_streamlines(trk_file, half_voxel_offset = TRUE)
tck <- read_streamlines(tck_file)
tt <- read_streamlines(tt_file)
```

To obtain the streamline data

```{r streamline_subset}
message("Total number of streamlines: ", length(trk))

head(trk[[1]]$coords)
```

To preview the streamline data

```{r streamline_plot, out.width="100%", fig.width = 9, fig.height=3}
pal <- colorRampPalette(c("navy", "grey", "red"))
plot(trk, col = pal(length(trk)))
```


To write the streamlines, for example, write `tck` file to `trk` file:

```{r streamline_write}
# Create a temporary file
tfile <- tempfile(fileext = ".trk")
write_streamlines(x = tck, con = tfile)
```

```{r streamline_cleanup, echo = FALSE, results='hide'}
if (file.exists(tfile)) {
  unlink(tfile)
}
```

## Surface annotations from `AFNI`/`SUMA`

`NIML` datasets (file names ending with `.niml.dset`) are read by the same `read_surface` entry point. The data type is resolved from the dataset itself: datasets carrying a label table are read as `annotations`, the rest as `measurements`.

```{r read_niml}
std141_geometry <- read_surface(std141_geom_file)

annotation <- read_surface(niml_file)
annotation
```

The sample geometry and annotation live on the same `std.141` mesh, so they merge directly:

```{r merge_niml}
labeled <- merge(std141_geometry, annotation)
labeled
```

```{r plot_niml, fig.width=7, fig.height=4, out.width="95%"}
plot(labeled, name = "annotations")
```

If you need the raw `NIML` element tree rather than a surface object, use the low-level `io_read_niml` together with `niml_find`.

## Regions of interest

A region of interest is described first and computed later. `as_ieegio_roi` records the criteria without applying them, and `resolve_roi_as` carries them out, returning geometry in world (`RAS`) coordinates.

```{r roi_describe}
atlas <- read_volume(atlas_file)

# describe: which voxels count as the region
roi <- as_ieegio_roi(atlas, threshold_lb = 1, threshold_ub = 5)
roi
```

```{r roi_resolve}
# compute: turn that description into geometry
resolve_roi_as(roi, "pointcloud")
```

Because both sides are resolved to world coordinates first, regions stored in different ways can be compared directly. Here the whole atlas is tested against the facial-nerve tracts read earlier:

```{r roi_overlap}
overlap <- detect_roi_overlap(
  as_ieegio_roi(atlas, threshold_lb = 1),
  trk,
  radius = 2
)
overlap
```

The result carries the annotated streamlines back in `overlap$annotated`, so each tract knows whether it reached the region, and `overlap$hit_ratio` reports the proportion that did.

## Volume to surface

`volume_to_surface` turns a mask or a set of atlas labels into a smoothed mesh:

```{r volume_to_surface}
volume_to_surface(atlas, threshold_lb = 1, threshold_ub = 5)
```

## Coordinate spaces and transforms

`new_space` names a coordinate space, and `surface_to_surface` moves a surface into it, recording the target on the surface transform list.

```{r spaces}
mni <- new_space("MNI152", orientation = "RAS")
mni

surface_to_surface(
  geometry,
  space_from = "scanner",
  space_to = mni,
  transform = diag(1, 4)
)
```

Transforms themselves can come from other tools: `io_read_ants_transform` reads `ANTs` affine and displacement field transforms, `io_read_flirt_transform` reads `FSL FLIRT` matrices, and `transform_flirt2ras` converts a `FLIRT` transform into world (`RAS`) coordinates.

```{r teardown, echo=FALSE, results='hide'}
cache_dir <- tools::R_user_dir("ieegio", "cache")
if (file.exists(cache_dir)) {
  unlink(cache_dir, recursive = TRUE, force = TRUE)
}
```