---
title: "brainloc"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{brainloc}
  %\VignetteEngine{knitr::rmarkdown}
  \usepackage[utf8]{inputenc}
---

## About brainloc

The brainloc package serves to describe a location on a brain surface mesh that is in MNI305 space, like the surfaces of the fsaverage standard brain templates used by FreeSurfer. Typical use cases of this package in the field of computational neuroimaging include:

* Finding the mesh vertex closest to a given coordinate (using Euclidean distance).
* Finding the coordinates of a vertex on in various standard coordinate systems (MNI305, MNI152 and Talairach).
* Finding the brain region the vertex lies in (according to any brain atlas available for the respective surface, like the Desikan atlas).
* Finding the distances of all atlas regions to the given vertex (using Euclidean or geodesic distance).

Statistical results are often presented as significant clusters of activation on a brain surface mesh (sets of connected vertices). Typical use cases for working with clusters are:

* Finding the extreme value of a cluster (min, max, or absolute max) and the vertex at which it occurs.
* Finding the peaks within a cluster.
* Finding the brain regions that overlap with the cluster (including the degree of overlap).

This vignette describes how to perform these tasks with the brainloc package.

### Preparation: Getting the data required to run the examples

You will need a directory with the brain templates you want to use, typically these are the fsaverage subjects that come with the FreeSurfer neuroimaging software suite. There are 2 ways to get them:

#### Using the template subjects from a local FreeSurfer installation

When FreeSurfer is installed and configured properly on your system, there is nothing you need to do. In that case, the environment variable `SUBJECTS_DIR` points at the `subjects` directory under the installation directory (`FREESURFER_HOME`), and it will be used to find the templates by the `get_subjects_dir` function. You can also set `SUBJECTS_DIR` to your custom template directory, of course, and that will be used instead.

If you cannot or do not want to set environment variables, you can set the option `brainloc.subjects_dir` in R to any directory that contains the expected directory structure like this: `options("brainloc.subjects_dir"="~/my_fs_home/subjects/")`. This would mean that `~/my_fs_home/subjects/` exists on your system and holds the directories of the template subjects. We assume that you have the `fsaverage` subject in there. Setting this option takes precedence over any environment variables.

#### Have the brainloc package download them for you from the internet

This requires that you have an internet connection and accept the FreeSurfer software license. See the code in the first example of this vignette below, which downloads the templates. The templates will be downloaded only if they cannot be found, and downloaded ones will be found on the next call to the function, of course.


## Usage Examples Part I: Location information for a vertex or coordinate


### Creating a brainparc object

Let's say we have a coordinate in MNI305 surface space and want to find the closest vertex. We first load `brainloc` and create a `brainparc` object, which internally is a named list holding the surfaces (the left and right hemisphere meshes) and the surface parcellations (according to some brain atlas) for these meshes. The parcellations assign to each vertex of a mesh a single region.

One can construct a brainparc manually from arbitrary loaded meshes and parcellations using the `brainparc` function, but here we use the function `brainparc_fs`, which is more convenient if you want to read it from files in a `SUBJECTS_DIR` of `recon-all` output data. We create a brainparc for the `fsaverage` template:


```{r, eval = TRUE}
library("brainloc");
# Searches for FreeSurfer SUBJECTS_DIR on machine and downloads if needed.
sjd = get_subjects_dir(allow_download = TRUE, accept_freesurfer_license = TRUE); 
bp = brainparc_fs(sjd, "fsaverage", surface = "white", atlas = "aparc");
```

By default, many functions in the `brainloc` package print information to STDOUT, which is handy for creating command line clients from the library. We do not want that here, so we disable it globally here:
```{r, eval = TRUE}
options("brainloc.silent"=TRUE);
```

Note: One can also decide for each function that prints something, they all have a `silent` parameter.


### Finding the vertex closest to a query coordinate

Let's say we have some query coordinates in surface space and want to find the closest vertex in the `brainparc` for each of the coordinates. Here is how we can do that:

```{r, eval = TRUE}
query_coords = matrix(seq(9)+0.1, ncol = 3, byrow = TRUE);
query_coords

ccv = coord_closest_vertex(query_coords, get_surface(bp));
ccv
```

This assumes that the query coordinates are in the surfaces space and uses Euclidean distance. It works for all subjects and meshes.

### Finding the MNI152 and Talairach coordinates for a point or vertex in MNI305 space

This assumes that you query vertex (or coordinate) is in MNI305 space. This is true for all FreeSurfer standard template subjects, like `fsaverage`, `fsaverage6`, etc. Our `brainparc` is for `fsaverage`, so we can use the coordinates of its vertices to query.

```{r, eval = TRUE}
query_vertices = c(10L, 145029L);
query_hemis = c("lh", "rh"); # vertex 10 on left hemi and 145029L on right hemi.
query_coords = get_surface_coords(bp, query_vertices, query_hemis);
query_coords

coord_info = coord_MNI305_info(query_coords, surface = get_surface(bp));
coord_info
```

Note: If you have the `regfusionr` package installed, have a look at the `method` parameter of the `coord_MNI305_info` function: you can use the more accurate `regfusionr` method to transform from MNI305 to MNI152 then. Open the docs with `?coord_MNI305_info` to learn more.


### Finding the brain atlas region a vertex belongs to, and the distance to all other brain atlas regions from the vertex

We can easily find the regions a vertex is assigned to in all atlases of a brain parcellation:

```{r, eval = TRUE}
query_vertices = c(10L, 145029L);
query_hemis = c("lh", "rh"); # vertex 10 on left hemi and 145029L on right hemi.
regions = vertex_regions(bp, query_vertices, query_hemis);

regions
```


One can also compute the distance of a vertex to all atlas regions. Several definitions for the distance of a vertex to a region are available, here we go with the defaults:
```{r, eval = TRUE}
query_vertices = c(10L, 145029L);
query_hemis = c("lh", "rh"); # vertex 10 on left hemi and 145029L on right hemi.
region_dists = vertex_closest_regions(bp, query_vertices, query_hemis, num_regions_to_report=5L);

region_dists
```


One can change the distance definition with the parameters `distance` and `linkage` of the `vertex_closest_regions` function. Open the docs with `?vertex_closest_regions` to learn more.


## Usage Examples Part II: Location information for a cluster

A cluster is a set of vertices on a surface, and the statistical values that are assigned to the cluster's vertices. Clusters are used in vertex-wise analyses to define activated brain regions (where the exact meaning of *activated* depends on the context). A clusters has a freeform label or name, but typically only consecutive numbers are used.

The data structure used to store clusters in `brainloc` is called `clusterinfo`. The following two concepts are relevant to understanding how a `clusterinfo` instance is typically constructed:

* a `cluster overlay` is an integer vector that assigns to each surface vertex a cluster number (zero is assigned to background/non-cluster vertices).
* a `statistical map` contains arbitrary values, one for each surface vertex (a.k.a. per-vertex data). A *thresholded* map only contains values for the cluster vertices, the values for all non-cluster vertices are set to an arbitrary value (typically *0.0*).

A `clusterinfo` object can be constructed for a subject in two ways:

* if you have both a `cluster overlay` and the `statistical map` (thresholded or unthresholded), use the `clusterinfo` function.
* if you have only a thresholded `statistical map`, use the `clusterinfo_from_thresholded_overlay` function instead.

In the following example, we construct a `clusterinfo` instance from a cluster overlay and the statistical map:

```{r, eval = TRUE}
lh_tmap_file = system.file("extdata", "lh.tmap.mgh", package = "brainloc", mustWork = TRUE);
rh_tmap_file = system.file("extdata", "rh.tmap.mgh", package = "brainloc", mustWork = TRUE);
lh_overlay_file = system.file("extdata", "lh.cluster.overlayID.mgh", package = "brainloc", mustWork = TRUE);
rh_overlay_file = system.file("extdata", "rh.cluster.overlayID.mgh", package = "brainloc", mustWork = TRUE);
clinfo = clusterinfo(lh_overlay_file, rh_overlay_file, lh_tmap_file, rh_tmap_file, template_subject = "fsaverage", subjects_dir = sjd);
is.clusterinfo(clinfo);
```

One can extract the clusters from the `clusterinfo`` instance like this:

```{r, eval = TRUE}
clusters = get_clusters(clinfo);
```

This returns a named list, where the keys are cluster names and the values are integer vectors containing the vertices that belong the the respective cluster. See `?get_clusters` for more options, like returning only the clusters for a specific hemisphere.

The `clusterinfo` object can now be used to compute various information on the clusters. We will start by computing information on the location of the clusters:

```{r, eval = TRUE}
cl_details = cluster_location_details(clinfo);
cl_details;
```

Note that the `clusterinfo` object contains a brain parcellation (a `brainparc` object, see at the top). We can also compute the overlap of the cluster with its atlas regions:

```{r, eval = TRUE}
overlap = cluster_region_overlap(clinfo);
overlap;
```

Sometimes one wants to find all peaks of a cluster. A peak exists at vertices which have the most extreme value in their neighborhood. For a brain mesh, the neighborhood of a vertex is defined as all vertices which are reachable by one hop along the mesh edges.

```{r, eval = TRUE}
peaks = cluster_peaks(clinfo);
peaks
```

This concludes the package vignette, thanks for reading.