inst/doc/freesurferformats.R

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

## -----------------------------------------------------------------------------
    library("freesurferformats")
    mgh_file = system.file("extdata", "brain.mgz", package = "freesurferformats", mustWork = TRUE)
    brain = read.fs.mgh(mgh_file)
    cat(sprintf("Read voxel data with dimensions %s. Values: min=%d, mean=%f, max=%d.\n", paste(dim(brain), collapse = 'x'), min(brain), mean(brain), max(brain)));

## -----------------------------------------------------------------------------
    brain_with_hdr = read.fs.mgh(mgh_file, with_header = TRUE);
    brain = brain_with_hdr$data;   # as seen before, this is what we got in the last example (the data).
    header = brain_with_hdr$header;   # the header

## ---- eval = FALSE------------------------------------------------------------
#      header$dtype           # int, one of: 0=MRI_UCHAR; 1=MRI_INT; 3=MRI_FLOAT; 4=MRI_SHORT
#      header$ras_good_flag   # int, 0 or 1. Whether the file contains a valid vox2ras matrix and ras_xform (see header$vox2ras_matrix below)
#      header$has_mr_params   # int, 0 or 1. Whether the file contains mr_params (see header$mr_params below)
#      header$voldim          # integer vector or length 4. The volume (=data) dimensions. E.g., c(256, 256, 256, 1) for 3D data.
#  

## -----------------------------------------------------------------------------
    if(header$has_mr_params) {
        mr_params = header$mr_params;  
        cat(sprintf("MR acquisition parameters: TR [ms]=%f, filp angle [radians]=%f, TE [ms]=%f, TI [ms]=%f\n", mr_params[1], mr_params[2], mr_params[3], mr_params[4]));
    }

## -----------------------------------------------------------------------------
    if(header$ras_good_flag) {
        print(header$vox2ras_matrix);
    }

## -----------------------------------------------------------------------------
    if(header$ras_good_flag) {
        print(header$ras_xform);
    }

## ---- eval = FALSE------------------------------------------------------------
#      mgh_file = system.file("mystudy", "subject1", "surf", "lh.thickness.fwhm25.fsaverage.mgh")
#      cortical_thickness_standard = read.fs.mgh(mgh_file)

## ---- eval = FALSE------------------------------------------------------------
#      oro.nifti::orthographic(oro.nifti::nifti(brain))

## -----------------------------------------------------------------------------
    library("freesurferformats")
    curvfile = system.file("extdata", "lh.thickness", package = "freesurferformats", mustWork = TRUE)
    ct = read.fs.curv(curvfile)

## -----------------------------------------------------------------------------
    cat(sprintf("Read data for %d vertices. Values: min=%f, mean=%f, max=%f.\n",  length(ct), min(ct), mean(ct), max(ct)))
    hist(ct, main="lh cortical thickness", xlab="Cortical thickness [mm]", ylab="Vertex count")

## -----------------------------------------------------------------------------
    morphfile1 = system.file("extdata", "lh.thickness", package = "freesurferformats", mustWork = TRUE)
    thickness_native = read.fs.morph(morphfile1)

## -----------------------------------------------------------------------------
    morphfile2 = system.file("extdata", "lh.curv.fwhm10.fsaverage.mgz", package = "freesurferformats", mustWork = TRUE)
    curv_standard = read.fs.morph(morphfile2)
    curv_standard[curv_standard < -1] = 0; # remove extreme outliers
    curv_standard[curv_standard > 1] = 0;
    hist(curv_standard, main="lh std curvature", xlab="Mean Curvature [mm^-1], fwhm10", ylab="Vertex count")

## -----------------------------------------------------------------------------
    annotfile = system.file("extdata", "lh.aparc.annot.gz", package = "freesurferformats", mustWork = TRUE);
    annot = read.fs.annot(annotfile);

## -----------------------------------------------------------------------------
    num_vertices_total = length(annot$vertices);
    for (vert_idx in c(1, 5000, 123456)) {
        cat(sprintf("Vertex #%d with zero-based index %d has label code '%d' which stands for atlas region '%s'\n", vert_idx, annot$vertices[vert_idx], annot$label_codes[vert_idx], annot$label_names[vert_idx]));
    }

## -----------------------------------------------------------------------------
    ctable = annot$colortable$table;
    regions = annot$colortable$struct_names;
    for (region_idx in seq_len(annot$colortable$num_entries)) {
        cat(sprintf("Region #%d called '%s' has RGBA color (%d %d %d %d) and code '%d'.\n", region_idx, regions[region_idx], ctable[region_idx,1], ctable[region_idx,2], ctable[region_idx,3], ctable[region_idx,4], ctable[region_idx,5]));
    }

## -----------------------------------------------------------------------------
    r_index = 50;                       # one-based index as used by R and Matlab
    fs_index = annot$vertices[r_index];  # zero-based index as used in C, Java, Python and many modern languages
    cat(sprintf("Vertex at R index %d has FreeSurfer index %d and lies in region '%s'.\n", r_index, fs_index, annot$label_names[r_index]));

## -----------------------------------------------------------------------------
    region = "bankssts"
    thickness_in_region = thickness_native[annot$label_names == region]
    cat(sprintf("Region '%s' has %d vertices and a mean cortical thickness of %f mm.\n", region, length(thickness_in_region), mean(thickness_in_region)));

## -----------------------------------------------------------------------------
    surface_file = system.file("extdata", "lh.tinysurface", package = "freesurferformats", mustWork = TRUE);
    surf = read.fs.surface(surface_file);
    cat(sprintf("Loaded surface consisting of %d vertices and %d faces.\n", nrow(surf$vertices), nrow(surf$faces)));

## -----------------------------------------------------------------------------
    vertex_index = 5;
    v5 = surf$vertices[vertex_index,];
    cat(sprintf("Vertex %d has coordinates (%f, %f, %f).\n", vertex_index, v5[1], v5[2], v5[3]));

## -----------------------------------------------------------------------------
    face_index = 2;
    f2 = surf$faces[face_index,];
    cat(sprintf("Face %d consistes of the vertices %d, %d, and %d.\n", face_index, f2[1], f2[2], f2[3]));

## ---- eval = FALSE------------------------------------------------------------
#  pmesh = read.fs.surface("~/path/to/mesh.ply");
#  my_mesh3d = rgl::tmesh3d(c(t(pmesh$vertices)), c(t(pmesh$faces)), homogeneous=FALSE);

## ---- eval = FALSE------------------------------------------------------------
#  curv = Rvcg::vcgCurve(my_mesh3d);

## -----------------------------------------------------------------------------
    labelfile = system.file("extdata", "lh.entorhinal_exvivo.label", package = "freesurferformats", mustWork = TRUE);
    label = read.fs.label(labelfile);
    cat(sprintf("The label consists of %d vertices, the first one is vertex %d.\n", length(label), label[1]));

## -----------------------------------------------------------------------------
lutfile = system.file("extdata", "colorlut.txt", package = "freesurferformats", mustWork = TRUE);
colortable = read.fs.colortable(lutfile, compute_colorcode=TRUE);
head(colortable);

## -----------------------------------------------------------------------------
annotfile = system.file("extdata", "lh.aparc.annot.gz", package = "freesurferformats", mustWork = TRUE);
annot = read.fs.annot(annotfile);
colortable = colortable.from.annot(annot);
head(colortable);