Skip to content

Atlas

The atlas module provides unified atlas discovery, region listing, and overlap analysis for both surface (mesh) and volumetric (voxel) atlases. It is used internally by the analyzer, optimization, and statistics modules.

graph LR
    SEG([Segmentation Dir]) --> MESH[MeshAtlasManager]
    FS([FastSurfer / FreeSurfer mri/]) --> VOXEL[VoxelAtlasManager]
    MESH --> REGIONS([Region Lists])
    VOXEL --> REGIONS
    SIG([Significant Mask]) --> OVERLAP[atlas_overlap_analysis]
    VOXEL --> OVERLAP
    OVERLAP --> RESULTS([Overlap Results])
    style SEG fill:#1a3a5c,stroke:#48a,color:#fff
    style FS fill:#1a3a5c,stroke:#48a,color:#fff
    style SIG fill:#1a3a5c,stroke:#48a,color:#fff
    style MESH fill:#2d5a27,stroke:#4a8,color:#fff
    style VOXEL fill:#2d5a27,stroke:#4a8,color:#fff
    style OVERLAP fill:#2d5a27,stroke:#4a8,color:#fff
    style REGIONS fill:#1a5c4a,stroke:#4a8,color:#fff
    style RESULTS fill:#1a5c4a,stroke:#4a8,color:#fff

Surface (Mesh) Atlases

MeshAtlasManager discovers and queries FreeSurfer .annot parcellation files in the m2m_{subject}/segmentation/ directory. Built-in atlases (DK40, a2009s, HCP_MMP1) are always included in discovery results.

from tit.atlas import MeshAtlasManager

manager = MeshAtlasManager(seg_dir="/data/my_project/derivatives/SimNIBS/sub-001/m2m_001/segmentation")

# List all available mesh atlases
atlases = manager.list_atlases()
# ['DK40', 'HCP_MMP1', 'a2009s']

# List regions for a specific atlas
regions = manager.list_regions("DK40")
# ['bankssts-lh', 'bankssts-rh', 'caudalanteriorcingulate-lh', ...]

# Find the .annot file for a specific atlas and hemisphere
annot_path = manager.find_atlas_file("DK40", "lh")

# Get all atlas files for a hemisphere
all_atlases = manager.find_all_atlases("lh")
# {'DK40': '/path/to/lh.aparc_DK40.annot', ...}

Volumetric (Voxel) Atlases

VoxelAtlasManager discovers current FastSurfer outputs, optional or existing FreeSurfer outputs, SimNIBS segmentation, and custom label masks. Region labels are read with nibabel and NumPy and cached; no FreeSurfer executable is needed.

from tit.atlas import VoxelAtlasManager

manager = VoxelAtlasManager(
    fastsurfer_mri_dir="/data/my_project/derivatives/fastsurfer/sub-001/mri",
    freesurfer_mri_dir="/data/my_project/derivatives/freesurfer/sub-001/mri",  # optional FreeSurfer data
    masks_dir="/data/my_project/derivatives/SimNIBS/sub-001/m2m_001/masks",
    seg_dir="/data/my_project/derivatives/SimNIBS/sub-001/m2m_001/segmentation",
)

# Discover available voxel atlases
atlases = manager.list_atlases()
# [('aparc.DKTatlas+aseg.deep.mgz', '/path/to/aparc.DKTatlas+aseg.deep.mgz'), ...]

# List regions in a specific atlas
regions = manager.list_regions(atlases[0][1])
# ['Left-Cerebellum-Cortex (ID: 8)', 'Right-Hippocampus (ID: 53)', ...]

# Detect MNI-space atlases from the resources directory
mni_atlases = VoxelAtlasManager.detect_mni_atlases("/ti-toolbox/resources/atlas")

# Find the SimNIBS labeling LUT
lut_path = manager.find_labeling_lut()

Atlas Overlap Analysis

The atlas_overlap_analysis function identifies which atlas regions overlap with a binary mask of significant voxels. This is used by the statistics module to map cluster results back to anatomical regions. The companion check_and_resample_atlas function handles dimension mismatches by resampling the atlas to the reference image space.

import nibabel as nib
import numpy as np
from tit.atlas import atlas_overlap_analysis, check_and_resample_atlas

# Load a significance mask and reference image
sig_mask = np.zeros((182, 218, 182), dtype=bool)
sig_mask[80:100, 100:120, 80:100] = True

reference_img = nib.load("/path/to/subject_field.nii.gz")

# Run overlap analysis against multiple atlases
results = atlas_overlap_analysis(
    sig_mask=sig_mask,
    atlas_files=["aparc.DKTatlas+aseg.deep.mgz"],
    data_dir="/data/my_project/derivatives/fastsurfer/sub-001/mri",
    reference_img=reference_img,
)

# results is a dict: atlas_name -> list of overlap dicts
for atlas_name, overlaps in results.items():
    for region in overlaps[:5]:
        pct = 100 * region["overlap_voxels"] / region["region_size"]
        print(f"  Region {region['region_id']}: {region['overlap_voxels']} voxels ({pct:.1f}%)")

Built-in Atlases

Mesh (Surface) Atlases

These are always available after running subject_atlas during preprocessing (CHARM):

Atlas Description
DK40 Desikan-Killiany atlas -- 68 cortical regions
a2009s Destrieux atlas -- 148 cortical regions
HCP_MMP1 Human Connectome Project Multi-Modal Parcellation -- 360 cortical regions

Voxel (Volumetric) Atlases

Current FastSurfer segmentation is discovered from derivatives/fastsurfer/sub-<id>/mri/:

File Description
aparc.DKTatlas+aseg.deep.mgz or .nii.gz DKT cortical and subcortical deep segmentation

The following FreeSurfer outputs are discovered under derivatives/freesurfer/sub-<id>/mri/. Optional recon-all produces the cortical parcellations; subregions require an additional operation after reconstruction. Existing legacy outputs remain readable:

File Hemisphere Description
aparc.DKTatlas+aseg.mgz both DKT cortical + subcortical segmentation
aparc.a2009s+aseg.mgz both Destrieux cortical + subcortical segmentation
lh.hippoAmygLabels-T1.v22.mgz lh Left hippocampal/amygdala subfields
rh.hippoAmygLabels-T1.v22.mgz rh Right hippocampal/amygdala subfields
ThalamicNuclei.v13.T1.mgz both Thalamic nuclei segmentation

Custom Masks

Any integer label volume placed in the subject's mask directory is discovered automatically and offered alongside the curated atlases — no whitelist entry and no conversion required:

derivatives/SimNIBS/sub-<id>/m2m_<id>/masks/*.nii.gz | *.nii | *.mgz

Custom entries appear in the ROI pickers prefixed with masks/. Requirements:

  • The volume must be an integer label map in the subject space m2m_<id> was built from (same geometry as T1.nii.gz). A binary mask is label 1.
  • Region names come from the volume's own label values. Drop a <mask-stem>_LUT.txt colour table next to the mask to name them; otherwise labels are listed by id (using the standard FreeSurfer names where the id matches).

MNI-Space Atlases

Available from the container resources directory (/ti-toolbox/resources/atlas):

File Description
MNI_Glasser_HCP_v1.0.nii.gz Glasser HCP parcellation in MNI space (default)
massp2021-parcellation_decade-18to40.nii.gz MASSP subcortical parcellation in MNI space

API Reference

tit.atlas.mesh.MeshAtlasManager

MeshAtlasManager(seg_dir: str)

Discovers and queries FreeSurfer .annot mesh atlases.

Parameters:

Name Type Description Default
seg_dir str

Path to m2m_{subject}/segmentation/ directory.

required
Source code in tit/atlas/mesh.py
def __init__(self, seg_dir: str) -> None:
    self.seg_dir = seg_dir

list_atlases

list_atlases() -> list[str]

List available mesh atlas names.

Returns:

Type Description
list[str]

Sorted list of atlas names (always includes builtins).

Source code in tit/atlas/mesh.py
def list_atlases(self) -> list[str]:
    """List available mesh atlas names.

    Returns:
        Sorted list of atlas names (always includes builtins).
    """
    discovered: list[str] = []
    if os.path.isdir(self.seg_dir):
        for fname in os.listdir(self.seg_dir):
            if fname.startswith("lh.") and fname.endswith(".annot"):
                stem = fname[3:-6]
                atlas_name = stem.split("_", 1)[-1] if "_" in stem else stem
                if (
                    atlas_name not in discovered
                    and atlas_name not in BUILTIN_ATLASES
                ):
                    discovered.append(atlas_name)
    return sorted(set(BUILTIN_ATLASES + discovered))

list_regions

list_regions(atlas_name: str) -> list[str]

List regions for a mesh atlas from .annot files.

Returns:

Type Description
list[str]

Sorted list of region names in SimNIBS format (e.g. "lh.precentral").

Source code in tit/atlas/mesh.py
def list_regions(self, atlas_name: str) -> list[str]:
    """List regions for a mesh atlas from .annot files.

    Returns:
        Sorted list of region names in SimNIBS format (e.g. "lh.precentral").
    """
    if os.path.isdir(self.seg_dir):
        import nibabel.freesurfer as nfs

        regions: list[str] = []
        for hemi in ("lh", "rh"):
            matches = glob.glob(
                os.path.join(self.seg_dir, f"{hemi}.*{atlas_name}.annot")
            )
            if not matches:
                continue
            _, _, names = nfs.read_annot(matches[0])
            for n in names:
                name = n.decode("utf-8") if isinstance(n, bytes) else str(n)
                if name != "unknown":
                    regions.append(f"{hemi}.{name}")
        if regions:
            return sorted(regions)

    return []

find_atlas_file

find_atlas_file(atlas_name: str, hemisphere: str) -> str | None

Find the .annot file path for a given atlas and hemisphere.

Returns:

Type Description
str | None

Path to the .annot file, or None if not found.

Source code in tit/atlas/mesh.py
def find_atlas_file(self, atlas_name: str, hemisphere: str) -> str | None:
    """Find the .annot file path for a given atlas and hemisphere.

    Returns:
        Path to the .annot file, or None if not found.
    """
    if not os.path.isdir(self.seg_dir):
        return None
    matches = glob.glob(
        os.path.join(self.seg_dir, f"{hemisphere}.*{atlas_name}.annot")
    )
    return matches[0] if matches else None

find_all_atlases

find_all_atlases(hemisphere: str) -> dict[str, str]

Find all available atlas files for a hemisphere.

Returns:

Type Description
dict[str, str]

Dict mapping atlas display name to file path.

Source code in tit/atlas/mesh.py
def find_all_atlases(self, hemisphere: str) -> dict[str, str]:
    """Find all available atlas files for a hemisphere.

    Returns:
        Dict mapping atlas display name to file path.
    """
    atlas_map: dict[str, str] = {}
    if not os.path.isdir(self.seg_dir):
        return atlas_map
    pattern = f"{hemisphere}.*.annot"
    for atlas_file in sorted(glob.glob(os.path.join(self.seg_dir, pattern))):
        fname = os.path.basename(atlas_file)
        parts = fname.split(".")
        if len(parts) == 3 and parts[2] == "annot":
            atlas_full = parts[1]
            atlas_display = (
                atlas_full.split("_", 1)[-1] if "_" in atlas_full else atlas_full
            )
            atlas_map[atlas_display] = atlas_file
    return atlas_map

list_annot_regions

list_annot_regions(annot_path: str) -> list[tuple[int, str]]

List all regions in a .annot file.

Returns:

Type Description
list[tuple[int, str]]

List of (region_index, region_name) tuples.

Source code in tit/atlas/mesh.py
def list_annot_regions(self, annot_path: str) -> list[tuple[int, str]]:
    """List all regions in a .annot file.

    Returns:
        List of (region_index, region_name) tuples.
    """
    import nibabel.freesurfer.io as fsio

    labels, ctab, names = fsio.read_annot(annot_path)
    regions = []
    for i, name in enumerate(names):
        region_name = name.decode("utf-8") if isinstance(name, bytes) else str(name)
        regions.append((i, region_name))
    return regions

tit.atlas.voxel.VoxelAtlasManager

VoxelAtlasManager(freesurfer_mri_dir: str = '', seg_dir: str = '', masks_dir: str = '', fastsurfer_mri_dir: str = '')

Discovers and queries volumetric atlas files.

All discovery methods use the same canonical atlas-name lists from :mod:tit.atlas.constants so that the analyzer, flex-search and NIfTI viewer show identical atlases.

Search order for a name present in both trees: fastsurfer_mri_dir first, then freesurfer_mri_dir. The two never collide in practice -- FastSurfer writes aparc.DKTatlas+aseg.deep.* and recon-all writes aparc.DKTatlas+aseg.mgz -- but the order is fixed so a project that holds both offers the current pipeline's output first.

Parameters:

Name Type Description Default
freesurfer_mri_dir str

Path to a FreeSurfer mri/ directory (derivatives/freesurfer/sub-<id>/mri).

''
fastsurfer_mri_dir str

Path to the FastSurfer mri/ directory (derivatives/fastsurfer/sub-<id>/mri).

''
seg_dir str

Path to m2m_{subject}/segmentation/ directory.

''
masks_dir str

Path to m2m_{subject}/masks/, holding user-supplied custom label volumes. Optional; most subjects have no such directory.

''
Source code in tit/atlas/voxel.py
def __init__(
    self,
    freesurfer_mri_dir: str = "",
    seg_dir: str = "",
    masks_dir: str = "",
    fastsurfer_mri_dir: str = "",
) -> None:
    self.freesurfer_mri_dir = freesurfer_mri_dir
    self.fastsurfer_mri_dir = fastsurfer_mri_dir
    self.seg_dir = seg_dir
    self.masks_dir = masks_dir

list_atlases

list_atlases() -> list[tuple[str, str]]

Discover available voxel atlas files for a subject.

Checks FastSurfer mri/ first, then a legacy FreeSurfer mri/, then segmentation/ for charm's own labeling.nii.gz, then masks/ for any user-supplied label volume. Used by the analyzer tab, the flex subcortical tab and the NIfTI viewer.

Returns:

Type Description
list[tuple[str, str]]

List of (display_name, full_path) tuples, first match per name.

Source code in tit/atlas/voxel.py
def list_atlases(self) -> list[tuple[str, str]]:
    """Discover available voxel atlas files for a subject.

    Checks FastSurfer ``mri/`` first, then a legacy FreeSurfer ``mri/``,
    then ``segmentation/`` for charm's own ``labeling.nii.gz``, then
    ``masks/`` for any user-supplied label volume. Used by the analyzer
    tab, the flex subcortical tab and the NIfTI viewer.

    Returns:
        List of (display_name, full_path) tuples, first match per name.
    """
    results: list[tuple[str, str]] = []
    seen: set[str] = set()

    for mri_dir, names in (
        (self.fastsurfer_mri_dir, FASTSURFER_ATLAS_FILES),
        (self.freesurfer_mri_dir, FREESURFER_ATLAS_FILES),
    ):
        if not mri_dir or not os.path.isdir(mri_dir):
            continue
        for name in names:
            if name in seen:
                continue
            path = os.path.join(mri_dir, name)
            if os.path.isfile(path):
                results.append((name, path))
                seen.add(name)

    if self.seg_dir:
        labeling = os.path.join(self.seg_dir, "labeling.nii.gz")
        if os.path.isfile(labeling):
            results.append(("labeling.nii.gz", labeling))

    results.extend(self.list_custom_masks())

    return results

list_custom_masks

list_custom_masks() -> list[tuple[str, str]]

Discover user-supplied label volumes in the subject's masks/ dir.

Any *.nii.gz/*.nii/*.mgz file dropped into m2m_{subject}/masks/ is offered for targeting. Display names are prefixed with masks/ so custom entries are distinguishable from the curated atlases. A missing directory yields no entries.

Returns:

Type Description
list[tuple[str, str]]

Sorted list of (display_name, full_path) tuples.

Source code in tit/atlas/voxel.py
def list_custom_masks(self) -> list[tuple[str, str]]:
    """Discover user-supplied label volumes in the subject's masks/ dir.

    Any ``*.nii.gz``/``*.nii``/``*.mgz`` file dropped into
    ``m2m_{subject}/masks/`` is offered for targeting.  Display names are
    prefixed with ``masks/`` so custom entries are distinguishable from
    the curated atlases.  A missing directory yields no entries.

    Returns:
        Sorted list of (display_name, full_path) tuples.
    """
    if not self.masks_dir or not os.path.isdir(self.masks_dir):
        return []

    found: list[tuple[str, str]] = []
    for name in sorted(os.listdir(self.masks_dir)):
        if not name.lower().endswith(MASK_EXTENSIONS):
            continue
        path = os.path.join(self.masks_dir, name)
        if os.path.isfile(path):
            found.append((f"masks/{name}", path))
    return found

list_regions

list_regions(atlas_path: str) -> list[str]

List regions in a voxel atlas.

Pure-Python replacement for shelling out to FreeSurfer's mri_segstats: labels present in the volume come from :func:~tit.atlas.segstats.compute_segstats (nibabel + numpy), named via :func:~tit.atlas.segstats.resolve_lut_for_atlas (a sidecar colour table next to the atlas, else the bundled standard FreeSurfer colour table). Validated voxel-for-voxel against live mri_segstats output on real recon-all data -- see docs/dev/DECISIONS.md § 2026-09-03 (One Docker image and a real development loop).

Caches the label file next to the atlas, in mri_segstats --sum's own text layout (other modules parse this exact sidecar filename), so subsequent calls skip recomputation.

Returns:

Type Description
list[str]

Sorted list of "RegionName (ID: N)" strings.

Source code in tit/atlas/voxel.py
def list_regions(self, atlas_path: str) -> list[str]:
    """List regions in a voxel atlas.

    Pure-Python replacement for shelling out to FreeSurfer's
    ``mri_segstats``: labels present in the volume come from
    :func:`~tit.atlas.segstats.compute_segstats` (nibabel + numpy),
    named via :func:`~tit.atlas.segstats.resolve_lut_for_atlas` (a
    sidecar colour table next to the atlas, else the bundled standard
    FreeSurfer colour table). Validated voxel-for-voxel against live
    ``mri_segstats`` output on real recon-all data -- see
    ``docs/dev/DECISIONS.md § 2026-09-03 (One Docker image and a real development loop)``.

    Caches the label file next to the atlas, in ``mri_segstats
    --sum``'s own text layout (other modules parse this exact sidecar
    filename), so subsequent calls skip recomputation.

    Returns:
        Sorted list of "RegionName (ID: N)" strings.
    """
    atlas_bname = os.path.splitext(os.path.basename(atlas_path))[0]
    if atlas_bname.endswith(".nii"):
        atlas_bname = os.path.splitext(atlas_bname)[0]
    labels_file = os.path.join(
        os.path.dirname(atlas_path), f"{atlas_bname}_labels.txt"
    )

    if os.path.isfile(labels_file):
        return self._parse_labels_file(labels_file)

    lut = resolve_lut_for_atlas(atlas_path)
    stats = compute_segstats(atlas_path, lut)
    write_segstats_sum(stats, labels_file)
    return sorted({f"{s.name} (ID: {s.seg_id})" for s in stats})

detect_mni_atlases staticmethod

detect_mni_atlases(atlas_dir: str) -> list[str]

Detect available MNI atlases in an assets directory.

The list, its order and every atlas's kind/licence come from resources/atlas/manifest.json (:mod:tit.atlas.manifest); a manifest entry whose file is absent is skipped.

Parameters:

Name Type Description Default
atlas_dir str

Path to the atlas resources directory.

required

Returns:

Type Description
list[str]

List of full paths to found MNI atlas files.

Source code in tit/atlas/voxel.py
@staticmethod
def detect_mni_atlases(atlas_dir: str) -> list[str]:
    """Detect available MNI atlases in an assets directory.

    The list, its order and every atlas's kind/licence come from
    ``resources/atlas/manifest.json`` (:mod:`tit.atlas.manifest`); a
    manifest entry whose file is absent is skipped.

    Args:
        atlas_dir: Path to the atlas resources directory.

    Returns:
        List of full paths to found MNI atlas files.
    """
    from tit.atlas.manifest import mni_atlas_entries

    if not os.path.isdir(atlas_dir):
        return []
    return [entry["path"] for entry in mni_atlas_entries(atlas_dir)]

find_labeling_lut

find_labeling_lut() -> str | None

Find the LUT file for the SimNIBS labeling atlas.

Returns:

Type Description
str | None

Path to labeling_LUT.txt if it exists, else None.

Source code in tit/atlas/voxel.py
def find_labeling_lut(self) -> str | None:
    """Find the LUT file for the SimNIBS labeling atlas.

    Returns:
        Path to labeling_LUT.txt if it exists, else None.
    """
    lut_path = os.path.join(self.seg_dir, "labeling_LUT.txt")
    return lut_path if os.path.isfile(lut_path) else None

tit.atlas.overlap.atlas_overlap_analysis

atlas_overlap_analysis(sig_mask, atlas_files: list[str], data_dir: str, reference_img=None) -> dict[str, list]

Analyze overlap between significant voxels and atlas regions.

Parameters:

Name Type Description Default
sig_mask

Binary ndarray (x, y, z) of significant voxels.

required
atlas_files list[str]

List of atlas file names.

required
data_dir str

Directory containing atlas files.

required
reference_img

nibabel image for resampling (optional).

None

Returns:

Type Description
dict[str, list]

Dict mapping atlas names to lists of region overlap dicts.

Source code in tit/atlas/overlap.py
def atlas_overlap_analysis(
    sig_mask,
    atlas_files: list[str],
    data_dir: str,
    reference_img=None,
) -> dict[str, list]:
    """Analyze overlap between significant voxels and atlas regions.

    Args:
        sig_mask: Binary ndarray (x, y, z) of significant voxels.
        atlas_files: List of atlas file names.
        data_dir: Directory containing atlas files.
        reference_img: nibabel image for resampling (optional).

    Returns:
        Dict mapping atlas names to lists of region overlap dicts.
    """
    import numpy as np
    import nibabel as nib

    logger.info("\n" + "=" * 60)
    logger.info("ATLAS OVERLAP ANALYSIS")
    logger.info("=" * 60)

    results: dict[str, list] = {}

    for atlas_file in atlas_files:
        atlas_path = os.path.join(data_dir, atlas_file)
        if not os.path.exists(atlas_path):
            logger.warning("Atlas file not found - %s", atlas_file)
            continue

        logger.info("\n--- %s ---", atlas_file)
        atlas_img = nib.load(atlas_path)

        if reference_img is not None:
            atlas_data = check_and_resample_atlas(atlas_img, reference_img, atlas_file)
        else:
            atlas_data = atlas_img.get_fdata().astype(int)

        regions = np.unique(atlas_data[atlas_data > 0])

        region_counts = []
        for region_id in regions:
            region_mask = atlas_data == region_id
            overlap = np.sum(sig_mask & region_mask)

            if overlap > 0:
                region_counts.append(
                    {
                        "region_id": int(region_id),
                        "overlap_voxels": int(overlap),
                        "region_size": int(np.sum(region_mask)),
                    }
                )

        region_counts = sorted(
            region_counts, key=lambda x: x["overlap_voxels"], reverse=True
        )

        logger.info("\nTop regions by significant voxel count:")
        for i, r in enumerate(region_counts[:15], 1):
            pct = 100 * r["overlap_voxels"] / r["region_size"]
            logger.info(
                "%2d. Region %3d: %4d sig. voxels (%.1f%% of region)",
                i,
                r["region_id"],
                r["overlap_voxels"],
                pct,
            )

        results[atlas_file] = region_counts

    return results

tit.atlas.overlap.check_and_resample_atlas

check_and_resample_atlas(atlas_img, reference_img, atlas_name: str)

Check if atlas dimensions match reference, resample if needed.

Parameters:

Name Type Description Default
atlas_img

nibabel image of the atlas.

required
reference_img

nibabel image of the reference (subject data).

required
atlas_name str

Name of atlas for logging.

required

Returns:

Type Description

Atlas data as integer ndarray in correct dimensions.

Source code in tit/atlas/overlap.py
def check_and_resample_atlas(atlas_img, reference_img, atlas_name: str):
    """Check if atlas dimensions match reference, resample if needed.

    Args:
        atlas_img: nibabel image of the atlas.
        reference_img: nibabel image of the reference (subject data).
        atlas_name: Name of atlas for logging.

    Returns:
        Atlas data as integer ndarray in correct dimensions.
    """
    import numpy as np
    from nibabel.processing import resample_from_to
    import nibabel as nib

    atlas_shape = atlas_img.shape
    ref_shape = reference_img.shape

    logger.info("  Atlas shape: %s", atlas_shape)
    logger.info("  Reference shape: %s", ref_shape[:3])

    if atlas_shape[:3] != ref_shape[:3]:
        logger.info("  Dimensions don't match. Resampling atlas...")

        if len(ref_shape) > 3:
            ref_data_3d = reference_img.get_fdata()[:, :, :, 0]
        else:
            ref_data_3d = reference_img.get_fdata()

        ref_img_3d = nib.Nifti1Image(
            ref_data_3d.astype(np.float32),
            reference_img.affine[:4, :4],
            None,
        )

        atlas_data_raw = atlas_img.get_fdata()
        if len(atlas_data_raw.shape) > 3:
            atlas_data_raw = atlas_data_raw[:, :, :, 0]

        atlas_img_3d = nib.Nifti1Image(
            atlas_data_raw.astype(np.float32),
            atlas_img.affine[:4, :4],
            None,
        )

        resampled_atlas = resample_from_to(atlas_img_3d, ref_img_3d, order=0)
        atlas_data = resampled_atlas.get_fdata().astype(int)
        logger.info("  Resampled to: %s", atlas_data.shape)
    else:
        logger.info("  Dimensions match.")
        atlas_data = atlas_img.get_fdata().astype(int)
        if len(atlas_data.shape) > 3:
            atlas_data = atlas_data[:, :, :, 0]

    return atlas_data