Working with BrainData¶

BrainData is the object almost every nltools analysis starts from: imaging data as an images-by-voxels matrix, one row per image and one column per voxel inside the mask. That shape is what makes it behave like a dataframe — index it, slice it, do arithmetic on it, iterate over it, and the metadata follows.

Load¶

fetch_pain() retrieves the pain dataset from Chang et al., 2015: 28 subjects with one contrast image each at low, medium and high pain intensity, cached locally on first use. The joblib cache below keeps the docs build from reloading 84 images on every run; call fetch_pain() directly in your own code.

Printing an object reports its shape, grid and mask, and len() is the number of images. .data is the numpy array underneath, .X a polars DataFrame of per-image metadata, and .Y outcomes or labels the same way:

from joblib import Memory

from nltools.data import BrainData
from nltools.datasets import fetch_pain

memory = Memory(".tutorial-cache", verbose=0)
data = memory.cache(fetch_pain)()

print(data)
print(f".data is a {type(data.data).__name__} of shape {data.data.shape}")

nltools.data.braindata.BrainData(data=(84, 238955), resolution=2.0mm, space=mni, mask=2mm-MNI152-2009fsl-mask.nii.gz)
.data is a ndarray of shape (84, 238955)

.X, one row per image:

data.X.head()

Your own data loads by path, and many files load together as a list. BrainData also takes a URL, a nibabel image, or a list of any of those.

one = BrainData("sub-01_pain-high.nii.gz")
many = BrainData(["sub-01_pain-high.nii.gz", "sub-02_pain-high.nii.gz"])

Pass mask to put an object on a specific grid — see Brain space and resolution for what decides the grid when you do not.

Index and slice¶

Indexing works the Python way, and every form carries .X and .Y along with the images it keeps: an integer gives one image, a slice a range of them on the same voxels and grid, a list of integers any images in any order, and a boolean array — usually from a column of .X — filters images. append concatenates along the image axis, and objects are iterable:

One image, by integer:

print(data[0])

nltools.data.braindata.BrainData(data=(238955,), resolution=2.0mm, space=mni, mask=2mm-MNI152-2009fsl-mask.nii.gz)

A range of images, by slice:

print(data[:5])

nltools.data.braindata.BrainData(data=(5, 238955), resolution=2.0mm, space=mni, mask=2mm-MNI152-2009fsl-mask.nii.gz)

Any images in any order, by list:

print(data[[0, 10, 20, 30]])

nltools.data.braindata.BrainData(data=(4, 238955), resolution=2.0mm, space=mni, mask=2mm-MNI152-2009fsl-mask.nii.gz)

The high-pain images, by boolean filter:

high_pain = data[data.X["PainLevel"] == 3]

print(f"{len(high_pain)} high-pain images")

28 high-pain images

Two objects concatenated:

print(data[:2].append(data[4]))

nltools.data.braindata.BrainData(data=(3, 238955), resolution=2.0mm, space=mni, mask=2mm-MNI152-2009fsl-mask.nii.gz)

Iterating over a stack, one image at a time:

[round(float(image.data.mean()), 2) for image in data[:5]]

[-0.32, -0.21, -0.11, -0.1, -0.07]

Arithmetic and statistics¶

Scalars broadcast over every voxel, and two objects add and subtract voxelwise, so subtracting two images gives one contrast map. Multiplying by a vector as long as the stack is a weighted sum across images instead, which is how you write a contrast over more than two conditions. mean and std reduce across images by default (axis=0), giving one value per voxel; axis=1 reduces across voxels. Methods chain, so temporal signal-to-noise ratio is one expression.

standardize z-scores or centers each voxel across images. threshold zeros everything below upper, everything above lower, or everything between the two when you pass both, resolving a percentile string over the nonzero voxels. binarize=True turns the survivors into ones, which is how you get a mask out of a map:

print((data + 10) * 2)
print(data[1] - data[0])

nltools.data.braindata.BrainData(data=(84, 238955), resolution=2.0mm, space=mni, mask=2mm-MNI152-2009fsl-mask.nii.gz)
nltools.data.braindata.BrainData(data=(238955,), resolution=2.0mm, space=mni, mask=2mm-MNI152-2009fsl-mask.nii.gz)

mean_map = data.mean()
tsnr = mean_map / data.std()

tsnr.plot(title="Temporal signal-to-noise ratio")

z_scored = data.standardize(method="zscore")

print(f"z-scored mean across voxels: {z_scored.data.mean():.6f}")

z-scored mean across voxels: 0.000000

top_voxels = mean_map.threshold(upper="95%", binarize=True)

print(f"top 5% of the mean map: {top_voxels.data.sum():.0f} voxels")

top 5% of the mean map: 11693 voxels

mean_map.threshold(upper="95%").plot(cmap="Blues", title="Top 5% of voxels")

Save¶

write saves to NIfTI, or to HDF5 when the file name ends in .h5 or .hdf5. Both load back with BrainData, but only HDF5 preserves .X, .Y and the mask: a NIfTI file holds the images and nothing else, so a round trip through it loses the metadata.

import tempfile
from pathlib import Path

tmpdir = Path(tempfile.mkdtemp())
nifti_path = tmpdir / "pain_subset.nii.gz"
h5_path = tmpdir / "pain_subset.h5"

data[:3].write(str(nifti_path))
data[:3].write(str(h5_path))

from_nifti = BrainData(str(nifti_path))
from_h5 = BrainData(str(h5_path))

print(f"NIfTI {nifti_path.stat().st_size / 1e6:.1f} MB")
print(f"  -> {from_nifti.shape}, X has {from_nifti.X.shape[1]} columns")
print(f"HDF5  {h5_path.stat().st_size / 1e6:.1f} MB")
print(f"  -> {from_h5.shape}, X has {from_h5.X.shape[1]} columns")

NIfTI 2.7 MB
  -> (3, 238955), X has 0 columns
HDF5  2.7 MB
  -> (3, 238955), X has 8 columns

Plot¶

plot draws a glass brain by default, and method='slices' draws cross-sections, one row per axis in view; the remaining keywords go to the matching nilearn plot. On a stack it draws one figure per image up to limit, which defaults to 3 so an 84-image object does not silently produce 84 figures. plot_surf projects the volume onto the fsaverage pial surface as a lateral and medial montage, plot_flatmap onto a flattened surface, and to_nifti hands back a nibabel image — 3D for one image, 4D for a stack — for anything these do not cover:

mean_map.plot(title="Mean activation")

data[:2].plot(limit=2)

[<Figure size 950x350 with 6 Axes>, <Figure size 950x350 with 6 Axes>]