Masking¶
Run this tutorial
This page is rendered from the marimo notebook docs/tutorials/data-operations/02_masking.py. Click the badge to run it in the cloud (free, no install), or locally: download 02_masking.py and run uvx marimo edit --sandbox 02_masking.py. The outputs below were produced when this page was built.
A mask restricts a BrainData object to the voxels you care about. This
tutorial builds masks three ways — a sphere around a coordinate, a
parcellation split into its regions, and a thresholded statistic map — then
uses them to summarize data and to paint per-region results back onto the
brain.
Load data¶
We use the pain dataset throughout. The joblib cache keeps the docs build
from reloading 84 images on every run; you can call fetch_pain() directly.
from joblib import Memory
from nltools.datasets import fetch_pain
memory = Memory(".tutorial-cache", verbose=0)
data = memory.cache(fetch_pain)()
data
nltools.data.braindata.BrainData(data=(84, 238955), resolution=2.0mm, space=mni, mask=2mm-MNI152-2009fsl-mask.nii.gz)
A spherical mask¶
create_sphere draws binary spheres. Centers are MNI millimeter coordinates
and radius is in millimeters, so the same request covers the same physical
volume whatever grid you are on. apply_mask then keeps only the voxels
inside it.
from nltools.mask import create_sphere
sphere = create_sphere([0, 0, 0], radius=30)
masked_data = data.apply_mask(sphere)
masked_data.mean().plot()
Masking drops voxels rather than zeroing them, so the masked object is narrower than the original:
whole brain: (84, 238955) 30 mm sphere: (84, 14147)
Average within a region¶
extract_roi collapses each region of a mask to one number per image. With a
single binary region you get one value per image — the mean signal in that
sphere for each of the 84 images. The mask is resampled onto the object's own
grid first, so it does not have to match resolutions.
import matplotlib.pyplot as plt
roi_mean = data.extract_roi(sphere)
_fig, _ax = plt.subplots(figsize=(8, 3))
_ax.plot(roi_mean)
_ax.set(xlabel="image", ylabel="mean signal", title="30 mm sphere at [0, 0, 0]")
[Text(0.5, 0, 'image'), Text(0, 0.5, 'mean signal'), Text(0.5, 1.0, '30 mm sphere at [0, 0, 0]')]
Expand and collapse a parcellation¶
A parcellation is one image whose voxel values are integer region IDs. nltools ships several; here is a 50-region whole-brain parcellation, fetched from the package's data repository.
from nltools.data import BrainData
from nltools.templates import fetch_resource
parcellation = BrainData(
fetch_resource("masks/default/2mm-MNI152-2009fsl-k50.nii.gz")
)
parcellation.plot()
expand_mask turns those IDs into a stack of 50 binary masks, one per region:
from nltools.mask import expand_mask
regions = expand_mask(parcellation)
print(regions.shape)
regions[:3].plot()
(50, 238955)
[<Figure size 950x350 with 6 Axes>, <Figure size 950x350 with 6 Axes>, <Figure size 950x350 with 6 Axes>]
collapse_mask is the inverse: it folds a stack of binary masks back into one
labeled image, numbering the regions in stack order and dropping any overlap.
Masks from a statistic map¶
threshold cuts a map at an absolute value or a percentile. Here we average
the high-pain images and keep the tails outside the middle 95%.
binarize=True turns the survivors into a mask of ones:
regions goes further and splits a thresholded map into its spatially
contiguous blobs, one image per blob. Pass limit to draw them all:
blobs = high.threshold(lower="2.5%", upper="97.5%").regions()
print(f"{len(blobs)} contiguous regions")
blobs.plot(limit=len(blobs))
4 contiguous regions
[<Figure size 950x350 with 6 Axes>, <Figure size 950x350 with 6 Axes>, <Figure size 950x350 with 6 Axes>, <Figure size 950x350 with 6 Axes>]
Analyze regions, then map the answer back¶
Masks are most useful as a round trip: summarize each region, run an analysis over regions, and paint the result back into a brain image.
Here we compute a linear pain contrast per subject, correlate the 50 regions across subjects, threshold that correlation structure into a graph, and map each region's degree back onto the brain.
import numpy as np
# High minus low pain, one contrast image per subject
contrast = BrainData(
[
data[data.X["SubjectID"] == subject] * np.array([1, -1, 0])
for subject in data.X["SubjectID"].unique(maintain_order=True)
]
)
contrast
nltools.data.braindata.BrainData(data=(28, 238955), resolution=2.0mm, space=mni, mask=2mm-MNI152-2009fsl-mask.nii.gz)
extract_roi on a labeled parcellation gives regions by images — the profile
of each region across the 28 subjects:
region_profiles = contrast.extract_roi(parcellation)
print(region_profiles.shape) # (regions, subjects)
(50, 28)
Correlation distance between those profiles is a 50-node Adjacency.
Thresholding it keeps only the region pairs that covary across subjects:
from sklearn.metrics import pairwise_distances
from nltools.data import Adjacency
distance = Adjacency(
pairwise_distances(region_profiles, metric="correlation"),
matrix_type="distance",
)
connected = distance.threshold(upper=0.4, binarize=True)
connected.plot()
to_graph hands the thresholded matrix to networkx, where any graph metric
is available. Degree counts how many regions each region is tied to:
graph = connected.to_graph()
degree = np.array([d for _, d in sorted(graph.degree())])
print(f"{graph.number_of_nodes()} nodes, {graph.number_of_edges()} edges")
print(f"degree range: {degree.min()}-{degree.max()}")
50 nodes, 337 edges degree range: 0-40
roi_to_brain writes one value per region back into the expanded masks, giving
a brain image of degree centrality: