R/talairach.R

In dfsp-spirit/brainloc: Find Region Names and Coordinates for a Point on the Human Brain

# Functions for querying Talairach labels for coordinates in Talairach space.
# These functions require the Talairach lookup data (volume with label indices, and mapping of indices to label strings) from talairach.org.
# All

#' @title Retrieve Talairach labels for Talairach coordinates from \code{talairach.org} files.
#'
#' @param tal_coords \code{nx3} numeric matrix, the \code{n} query Talairach coordinates for which you want to retrieve the labels. (Three columns, each row represents the x,y,z coordinates of a point in Talairach space.)
#'
#' @param talairach_vol_file optional character string, the path to the NIFTI Talairach label volume file from \code{talairach.org}. It is named \code{talairach.nii} on the website. You can also gzip it and use a much smaller \code{.nii.gz} version of the file. Leave at \code{NULL} to auto-download from \code{talairach.org} if needed.
#'
#' @param lookup_table_file optional character string, the path to the Talairach label index text file from \code{talairach.org}. It is named \code{labels.txt} on the website. Leave at \code{NULL} to auto-download from \code{talairach.org} if needed.
#'
#' @param check_oob logical, whether to check for out-of-bounds voxels, which do not fall into the NIFTI volume during the computation (because the talairach coordinate is a bit off). Checking is slower, but if such coordinates exist in your query and you do not check, the function will stop with an error. Leave this alone if in doubt.
#'
#' @return data.frame describing labels for the coordinates. The following columns are included: \code{cx,cy,cz}: the query coordinates, as given in parameter 'tal_coords'. \code{v1,v2,v3}: The voxel indices in the talairach.nii file that map to the query coordinates. \code{label_lvl1,...,label_lvl5}: the label strings for the location, in a hierarchy with 5 levels (parsed from \code{label_full} for you as a convenience). \code{label_full}: The raw, full label string for the location. A star \code{*} means that no label name is available for the given coordinate at this level.
#'
#' @note If you leave the file path arguments \code{talairach_vol_file} and {lookup_table_file} at NULL, the files will be downloaded from \code{talairach.org} automatically if that is possible (e.g., you have a working internet connection and the server is up). If it fails, the function will stop in the next step.
#'
#' @note When using this function or the Talairach.org data in your research, please cite the following two publications: \code{Lancaster JL, Woldorff MG, Parsons LM, Liotti M, Freitas CS, Rainey L, Kochunov PV, Nickerson D, Mikiten SA, Fox PT, "Automated Talairach Atlas labels for functional brain mapping". Human Brain Mapping 10:120-131, 2000.} and \code{Lancaster JL, Rainey LH, Summerlin JL, Freitas CS, Fox PT, Evans AC, Toga AW, Mazziotta JC. Automated labeling of the human brain: A preliminary report on the development and evaluation of a forward-transform method. Hum Brain Mapp 5, 238-242, 1997.}
#'
#' @note Currently it is not possible to search 'in the vicinity' to compensate for a slightly inaccurate query coordinate that may be outside of the brain by just a tiny bit. Please open an issue if you feel this is required.
#'
#' @examples
#' \dontrun{
#' query_tal_coords = matrix(seq.int(15), nrow=5, ncol=3);
#' get_talairach_label(query_tal_coords);
#' }
#'
#' @importFrom utils read.table
#' @importFrom freesurferformats mghheader.ras2vox read.fs.volume doapply.transform.mtx
#'
#' @export
get_talairach_label <- function(tal_coords, talairach_vol_file=NULL, lookup_table_file=NULL, check_oob=TRUE) {

    if(is.null(talairach_vol_file) | is.null(lookup_table_file)) {
        download_talairach(accept_talairach_usage =  TRUE);
        talairach_vol_file = get_optional_data_filepath('talairach/talairach.nii');
        lookup_table_file = get_optional_data_filepath('talairach/labels.txt');
    }

    if(! file.exists(talairach_vol_file)) {
        stop(sprintf("Please download the 'talairach.nii' file from talairach.org and specify the correct path. Expected file not found at '%s'.\n", talairach_vol_file));
    }
    if(! file.exists(lookup_table_file)) {
        stop(sprintf("Please download the 'labels.txt' file from talairach.org and specify the correct path. Expected file not found at '%s'.\n", lookup_table_file));
    }

    if(is.vector(tal_coords)) {
        tal_coords = matrix(tal_coords, ncol = 3L);
    }

    if(is.matrix(tal_coords)) {
        if(ncol(tal_coords) != 3L) {
            stop("Matrix 'tal_coords' must have exactly 3 columns.");
        }
    }


    tal = freesurferformats::read.fs.volume(talairach_vol_file, with_header = TRUE);
    lab = read.table(lookup_table_file, sep = "\t", col.names = c("index", "region"));

    taldata = drop(tal$data);

    r2v = freesurferformats::mghheader.ras2vox(tal);
    voxels = floor(freesurferformats::doapply.transform.mtx(tal_coords, r2v, as_mat = TRUE));

    oob = NULL;
    if(check_oob) {
        # Find out-of-bounds voxels, which do not fall into the NIFTI volume. These need to be ignored.
        oob = which(voxels[,1] > dim(taldata)[1]);
        oob = c(oob, which(voxels[,2] > dim(taldata)[2]));
        oob = c(oob, which(voxels[,3] > dim(taldata)[3]));
        # We should also ignore voxels < 1 here.
        oob = c(oob, which(voxels[,1] < 1L));
        oob = c(oob, which(voxels[,2] < 1L));
        oob = c(oob, which(voxels[,3] < 1L));
    }

    if(length(oob) > 0L) {
        voxels[oob,] = c(1L,1L,1L); # replace
        message(sprintf("Ignoring %d of %d query coordinates that mapped to voxel indices outside of NIFTI volume.\n", length(oob), nrow(tal_coords)));
    }

    label_indices = taldata[voxels];

    if(length(oob) > 0L) {
        label_indices[label_indices < 1L] = 1L; # fix label_indices for out-of-bounds voxels from 0 to 1. 1 is an arbitrary but working index that ensures the computation finishes. The label values will be replaced with 'out-of-bounds' later.
    }

    voxel_labels = lab$region[label_indices]; # each label is single string that looks like: ''
    voxel_labels_split = strsplit(voxel_labels, split = ".", fixed = TRUE);

    res = data.frame("cx"=tal_coords[,1], "cy"=tal_coords[,2], "cz"=tal_coords[,3], "v1"=voxels[,1], "v2"=voxels[,2], "v3"=voxels[,3], stringsAsFactors = FALSE);

    ncoords = length(voxel_labels_split); # Could also use the number of query voxels.
    nlevels = 5L;
    for(level_idx in seq(nlevels)) {
        level_names = rep("*", ncoords);
        for(coord_idx in seq(ncoords)) {
            level_names[coord_idx] = voxel_labels_split[[coord_idx]][level_idx];
        }
        if(length(oob) > 0L) {
            level_names[oob] = "out-of-bounds";
        }
        key = paste("label_lvl", level_idx, sep = "");
        res[[key]] = level_names;
    }
    if(length(oob) > 0L) {
        voxel_labels[oob] = "out-of-bounds";
    }
    res$label_full = voxel_labels;
    return(res);
}