Create atlas from whole-brain volumetric parcellation

Build a brain atlas from a single volumetric parcellation (NIfTI/MGZ) that contains both cortical and subcortical regions. Cortical regions are projected onto the fsaverage5 surface via FreeSurfer's mri_vol2surf and rendered as surface views, while subcortical regions go through the mesh-based subcortical pipeline.

Requires FreeSurfer.

Usage

create_wholebrain_from_volume(
  input_volume,
  verbose = get_verbose(),
  ...,
  input_lut = NULL,
  atlas_name = NULL,
  output_dir = NULL,
  labels = list(),
  projection_opts = list(),
  cortical_opts = list(),
  subcortical_opts = list(),
  cerebellar_opts = list(),
  steps = NULL,
  min_coverage = NULL,
  cleanup = NULL,
  skip_existing = NULL
)

Arguments

input_volume

Path to volumetric parcellation in MNI152 space (.mgz, .nii, .nii.gz).

verbose

How much to print: 0 silent, 1 progress (the default), 2 adds FreeSurfer's own output. TRUE and FALSE mean 1 and 0. Falls back to options("ggseg.extra.verbose"), then GGSEG_EXTRA_VERBOSE.

...

Not used. Present so a mistyped or misplaced argument is reported against this function, naming the list that should hold it, rather than silently ignored; anything passed here is an error.

input_lut

Path to FreeSurfer-style colour lookup table, or a data.frame with columns idx, label, R, G, B, A. An optional type column with values "cortical", "subcortical" or "cerebellar" controls label classification (see Label classification); any other value is an error, and NA leaves the label to the other classification steps. An optional hemi column with values "left" or "right" says which hemisphere a parcel belongs to (see Declared hemisphere). Voxel IDs not listed in the LUT are automatically zeroed out before surface projection (see Volume pre-processing). If NULL, generic names and no palette.

atlas_name

Name for the atlas. If NULL, derived from the input filename. It also names the folder of intermediate files inside output_dir, so it must be a single name without path separators.

output_dir

Where to put the intermediate files. Defaults to tempdir(), from options("ggseg.extra.output_dir") or GGSEG_EXTRA_OUTPUT_DIR.

labels

Named list routing labels to a sub-pipeline, overriding the lookup table's type column and the vertex-count heuristic. Entries:

cortical: label names to force as cortical.
subcortical: label names to force as subcortical.
cerebellar: label names to send through the cerebellar SUIT flatmap pipeline, using the bundled surfaces from suit_flatmap_path() and suit_3d_path().

Unnamed or unknown entries error.

projection_opts

Named list of volume-to-surface projection settings, with these entries and defaults:

subject ("fsaverage5") and registration ("header"): the target surface subject and how the volume is aligned to it. Read Registration below before relying on the default.
projfrac (0.5) and projfrac_range (c(0, 1, 0.1)): how deep through the cortical ribbon to sample, passed on to mri_vol2surf. Set projfrac_range to NULL to sample the single depth projfrac instead.
min_vertices (50): how much surface a label needs before the vertex-count heuristic calls it cortical. Only reached for labels that labels and a type column leave unclassified; see Label classification.

Unknown entries error.

cortical_opts

Named list of extra arguments forwarded to the cortical sub-pipeline. Allowed entry: views. Unknown entries error. Leave empty to use defaults.

subcortical_opts

Named list of extra arguments forwarded to create_subcortical_from_volume(). Any argument of that function may be set here except those managed by the wholebrain pipeline (input_volume, input_lut, atlas_name, output_dir, verbose, cleanup, skip_existing). Use this to tune vertex_size_limits, decimate, slabs.

cerebellar_opts

Named list of extra arguments forwarded to create_cerebellar_from_volume(). Any argument of that function may be set here except those managed by the wholebrain pipeline (input_volume, input_lut, atlas_name, output_dir, verbose, cleanup, skip_existing).

One entry is not forwarded: cerebellar_space is consumed here, and names the space input_volume's cerebellar labels are in. The flatmap they are drawn on is in SUIT space, so anything else is transformed with the matching suit_deformation_field() first. Nothing in a NIfTI header records a volume's space, so "suit" is an assumption rather than a detection: set it if your volume is in an MNI space. See suit_deformation_field() for the spaces available.

steps

Which pipeline steps to run. Default NULL runs all steps. Steps are:

1: Project volume onto surface
2: Split labels into cortical/subcortical/cerebellar
3: Run cortical pipeline
4: Run subcortical pipeline
5: Run cerebellar pipeline

Use steps = 1:2 to run projection and split only.

min_coverage

Smallest share of core labels that must end up with 2D geometry, between 0 and 1. A build that falls below it warns and names the labels that have no shape; the warning fires whatever verbose is set to, because a region that cannot be drawn is data loss rather than progress chatter. The build itself still returns, so a long pipeline is never thrown away over a threshold. Use 0 to turn the check off. Defaults to 0.8, or to options("ggseg.extra.min_coverage"), then GGSEG_EXTRA_MIN_COVERAGE.

cleanup

Remove the intermediate files afterwards. Default TRUE, from options("ggseg.extra.cleanup") or GGSEG_EXTRA_CLEANUP.

The intermediate files live in a folder named after the atlas inside output_dir, and the whole folder is removed. A build therefore refuses to start if that folder already holds files it did not write, since removing it would take them too. Choose another output_dir or atlas_name, or pass FALSE to build there and keep everything.

skip_existing

Reuse a step's intermediate files when they already exist, instead of running the step. Default FALSE, from options("ggseg.extra.skip_existing") or GGSEG_EXTRA_SKIP_EXISTING.

A cache records which ggseg.extra wrote it, not what it was built from, so reuse cannot tell that the volume or lookup table has changed. Reuse is therefore asked for rather than assumed, and is reported when it happens. To continue an interrupted build, either pass TRUE or leave the finished steps out of steps.

Value

For a full run, a named list with elements cortical, subcortical and cerebellar, each a ggseg_atlas object (or NULL if no regions of that type exist). A run that stops at or before step 2 instead returns, invisibly, the label split to inspect: cortical_labels, subcortical_labels, cerebellar_labels and vertex_counts.

Label classification

The pipeline must know which labels are cortical (rendered on the surface) and which are subcortical (rendered as 3D meshes / 2D slices). Three mechanisms are available, applied in priority order:

The labels argument (highest priority): labels = list(cortical = ..., subcortical = ..., cerebellar = ...) overrides everything for the labels it names.
LUT type column: If the colour lookup table has a type column with values "cortical", "subcortical" or "cerebellar", that classification is used for any labels labels does not name. This is the recommended approach for reproducible atlas creation. Any other value is an error, raised before the projection starts; NA leaves a label to the next step.
Vertex-count heuristic (fallback): Labels with at least projection_opts$min_vertices vertices on the surface projection are classified as cortical; the rest as subcortical. It measures how much surface a label covers rather than where the label sits, so a small cortical parcel and a deep structure look the same to it. It warns whenever it runs; treat that warning as a request to declare the labels instead.

lut_classify_anatomy() writes the type column for a lookup table that has none, by reading each label's position in FreeSurfer's aparc+aseg. Run it once while authoring the atlas and commit the column it returns: a declared classification is reviewable in a diff, where an inferred one is not.

Volume pre-processing

Before surface projection, the volume is filtered so that only voxel IDs listed in the LUT are kept; all other non-zero voxels are zeroed out. This prevents unlisted structures (e.g. white matter, ventricle masks) from bleeding onto the cortical surface during label dilation.

The subcortical pipeline also gets a brain-outline reference to draw as grey context behind its structures, under FreeSurfer's cortex labels 3 (left) and 42 (right). By default that is the atlas's own cortical voxels, split by hemisphere at the volume's midline: left-hemisphere voxels map to label 3, right-hemisphere to label 42.

A parcellation that covers both banks of every sulcus makes a solid mantle that way, and no amount of polishing can put sulci into a silhouette that never had any. When the cortical labels hold more than 1.5 times the voxels of a cortical ribbon, the context is taken from FreeSurfer's aseg instead, where sulcal CSF is unlabelled: its ribbon is resampled onto this volume's own grid through the two headers (mri_vol2vol --regheader --nearest) and written wherever no structure claims the voxel. An atlas whose cortical labels are already a ribbon, such as one derived from a surface, keeps its own.

A ribbon is 2.5-3 mm thick, so a volume sampled more coarsely than that cannot hold one: the resampled ribbon breaks into islands and the silhouette with it. The substitution is declined when the resampled ribbon is not at least a voxel thick through most of itself, and the atlas's own cortical labels are used, as they were before.

Either way, the aseg cerebellar cortex and brain stem are added to the context, so the posterior fossa is not drawn as empty space behind an atlas that reaches below the tentorium or one whose cerebellum lives in a separate atlas. Cerebellar white matter is left out: without it the cerebellum stays a foliated shell rather than a solid lump that merges with the occipital lobe.

When no usable aseg is available - no FreeSurfer, no aseg.mgz for the subject, a failed resampling, or a ribbon that does not land inside this volume, which is what a volume in some other space looks like - the midline split is used and the pipeline warns.

Human oversight

This is the most complex pipeline in ggseg.extra and the one most likely to need manual correction. Recommended workflow:

Run steps = 1:2 first to project the volume and classify labels.
Inspect result$cortical_labels and result$subcortical_labels. If the automatic split is wrong, either add a type column to the LUT or name the labels in labels.
Run the full pipeline once you are satisfied with the split.
Visually inspect the resulting atlas with plot() for a quick overview, or ggplot() + geom_brain() and ggseg3d() for a figure.

The cortical surface projection uses FreeSurfer's cortex label ({hemi}.cortex.label) to prevent label dilation into the medial wall. This file ships with fsaverage5 and is always required. Cortex vertices left without a listed label take the most common label of their neighbours. Everything outside the cortex label becomes the unknown medial wall, which the cortical atlas keeps as grey context geometry rather than as a region.

Declared hemisphere

A cortical label carries no hemisphere of its own. The lh_/rh_ prefix on a finished region comes from whichever surface its vertices landed on, so a parcel whose voxels cross the midline is sampled onto both surfaces and becomes two regions the parcellation never had – the larger one real, the smaller one built from the spill.

Give the lookup table a hemi column ("left" or "right", or lh/rh) and the parcel is kept only on the surface it belongs to. Where there is no column the label's own name is read, which covers Left-Thalamus, lh.something, region_L and FreeSurfer's ctx-lh-superiorfrontal. A value the column cannot mean – a misspelling, say – is an error rather than a row quietly treated as undeclared.

A label that declares nothing is left alone, on purpose: a lookup table that names each structure once for both hemispheres is meant to produce lh_ and rh_. So this is opt-in, and a table without the column behaves exactly as before.

The column is passed on to the subcortical and cerebellar pipelines, where it sets the region's hemisphere instead of the label's name doing so. There it may also say "midline". "vermis" and "middle" are accepted there and recorded as "midline", since hemi carries only left/right/midline.

The column can arrive in a data.frame or in a lookup table file: write_lut() stores it and read_lut() brings it back.

Hemisphere is never inferred from the voxels. The volume does know which side a parcel sits on, but a majority vote is a guess made silently on every build, whereas a column is a fact you can see and correct – the same reasoning as the type column and lut_classify_anatomy().

Registration

fsaverage lives in MNI305, so an MNI152 volume projected with "header" samples roughly 2 mm off, uniformly, anteriorly and inferiorly. Every atlas in the ggsegverse was built that way. For reference maps meant for visualisation that error is usually acceptable, which is why it remains the default.

"mni152" removes it, but only for volumes on the grid mni152.register.dat was built for: 1 mm, left-handed (LAS). tkregister coordinates are derived from the volume's own voxel order, size and field of view, so the same matrix applied to a right-handed volume mirrors left and right, and applied to an LAS volume of another resolution mislocates regions (7.7% of vertices at 1.5 mm, about 28% at 4 mm, measured against the same volume resampled to 1 mm first). Both mistakes produce an atlas that still looks anatomically plausible, so an unsuitable grid is refused rather than warned about. To use "mni152" with such a volume, resample it onto the 1 mm LAS MNI152 grid first.

No header identifies the space an arbitrary volume is in, so the registration you name is applied to whatever you pass. A volume sitting on the target subject's exact voxel grid warns, that being a claim a header does support, but a volume in some other non-MNI152 space cannot be detected.

mni152.register.dat targets the FSL/SPM MNI152 (NLin6) 1 mm template. Volumes in other MNI152 variants, such as the NLin2009cAsym template ggsegJulich uses, keep a sub-millimetre residual rather than the roughly 2 mm the transform corrects. It is registered against fsaverage, so subject must share fsaverage's conformed geometry, which the fsaverageN subjects do; for any other subject, supply your own register.dat or LTA file.

FreeSurfer's own INFO and WARNING lines about the registration are only visible with verbose = TRUE.

Examples

if (FALSE) { # \dontrun{
# --- Recommended: LUT with type column ---
lut <- data.frame(
  idx   = c(10, 11, 49, 50, 101:148),
  label = c("Left-Thalamus", "Left-Caudate",
            "Right-Thalamus", "Right-Caudate",
            paste0("cortical_region_", 1:48)),
  type  = c(rep("subcortical", 4), rep("cortical", 48)),
  R = sample(50:220, 52, TRUE),
  G = sample(50:220, 52, TRUE),
  B = sample(50:220, 52, TRUE),
  A = 255L
)

result <- create_wholebrain_from_volume(
  input_volume = "my_atlas.nii.gz",
  input_lut = lut,
  atlas_name = "my_atlas",
  subcortical_opts = list(vertex_size_limits = c(3e6, 3e7))
)
result$cortical   # surface-based cortical atlas
result$subcortical # mesh-based subcortical atlas

# --- Without type column: automatic classification ---
result <- create_wholebrain_from_volume(
  input_volume = "atlas.nii.gz",
  input_lut = "atlas_LUT.txt",
  atlas_name = "auto_atlas"
)

# --- Inspect classification before full run ---
result <- create_wholebrain_from_volume(
  input_volume = "atlas.nii.gz",
  input_lut = "atlas_LUT.txt",
  steps = 1:2
)
result$cortical_labels
result$subcortical_labels

# --- Overriding classification and projection ---
# `labels` forces a label down a particular pipeline; `projection_opts`
# holds everything about getting the volume onto the surface.
result <- create_wholebrain_from_volume(
  input_volume = "atlas.nii.gz",
  input_lut = "atlas_LUT.txt",
  labels = list(
    cerebellar = c("Left-Cerebellum-Cortex", "Right-Cerebellum-Cortex")
  ),
  projection_opts = list(registration = "mni152", subject = "fsaverage6"),
  cerebellar_opts = list(cerebellar_space = "MNI152NLin6AsymC")
)
} # }