Inputs: NIfTI files, folders, BIDS lists, nilearn images and arrays

htfa.load turns every supported input into a list of htfa.Subject objects, each with a (V, T) data matrix and (V, D) coordinates. htfa.fit calls it, so everything here also works as the first argument of htfa.fit. The reference table is on the Inputs page of the documentation; this notebook runs each case on small synthetic images and prints what comes back.

The examples use synthetic 12 × 12 × 12 images with 3 mm voxels: a bright ball on a constant dark background, with noise inside the ball and a few time points. These simple shapes make loading and masking behavior visible; they are not anatomical models or measurements from an experiment.

# Install htfa if it is missing (for example on Colab); other packages
# install themselves the first time they are needed.
try:
    import htfa  # noqa: F401
except ImportError:
    %pip install -q htfa

import tempfile
import time
import warnings
from pathlib import Path

import nibabel as nib
import numpy as np

import htfa

NiftiMasker = htfa.require("nilearn.maskers").NiftiMasker

rng = np.random.default_rng(0)
shape, T, voxel_mm = (12, 12, 12), 10, 3.0
affine = np.diag([voxel_mm, voxel_mm, voxel_mm, 1.0])
grid = np.stack(np.meshgrid(*[np.arange(n) for n in shape], indexing="ij"), axis=-1)
inside = np.linalg.norm(grid - (np.array(shape) - 1) / 2, axis=-1) < 4.5
print("voxels inside the ball:", int(inside.sum()), "of", inside.size)


def brain(T):
    """Ball of ~100 on a constant background of 5; noise inside the ball."""
    base = np.where(inside, 100.0, 5.0)[..., None]
    return base + rng.normal(0, 2, size=shape + (T,)) * inside[..., None]


tmp = tempfile.TemporaryDirectory()
root = Path(tmp.name)

voxels inside the ball: 360 of 1728

1. A single file and a nibabel image

A path to a .nii or .nii.gz file is one subject; a 4D file is its time series. A nibabel image (for example one returned by nilearn) is also one subject. Coordinates are world coordinates computed from the affine, so with 3 mm voxels they step by 3. By default a voxel is kept only if it is finite at every time point and non-zero at least once.

single = root / "single.nii.gz"
nib.save(nib.Nifti1Image(brain(T), affine), single)

(s,) = htfa.load(str(single))
print(
    "subject from a path:",
    s.data.shape,
    "| coords",
    s.coords.shape,
    "| image_shape",
    s.image_shape,
)
print("first coordinates (mm):", s.coords[:2].tolist())
print("source:", Path(str(s.source)).name)

(s_img,) = htfa.load(nib.load(single))
print(
    "subject from a nibabel image:",
    s_img.data.shape,
    "| same data:",
    np.array_equal(s.data, s_img.data),
)

subject from a path: (1728, 10) | coords (1728, 3) | image_shape (12, 12, 12)
first coordinates (mm): [[0.0, 0.0, 0.0], [0.0, 0.0, 3.0]]
source: single.nii.gz
subject from a nibabel image: (1728, 10) | same data: True

2. Folders and nested conditions

A flat folder contains one subject per NIfTI file. Files are taken in natural filename order (sub-2 before sub-10), the extension check ignores case, and hidden files and non-image files such as .json sidecars are skipped. Subjects can differ in number of time points.

Subfolders describe conditions or nested groups. Their names form each subject’s group_path, relative to the supplied root. htfa.fit uses these paths to construct the hierarchy by default. Organize the folders according to the grouping you intend to model; a directory structure alone does not establish that observations are independent.

folder = root / "folder"
folder.mkdir()
for name, t in [("sub-10", 8), ("sub-2", 12), ("sub-1", 10)]:
    nib.save(nib.Nifti1Image(brain(t), affine), folder / f"{name}_bold.nii.gz")
(folder / "sub-1_bold.json").write_text("{}")
(folder / ".hidden.nii.gz").write_text("not a NIfTI file")
nib.save(nib.Nifti1Image(brain(9), affine), folder / "sub-3_bold.NII.GZ")

subjects = htfa.load(str(folder))
print("files in the folder:", sorted(p.name for p in folder.iterdir()))
print("subjects loaded:", len(subjects))
for s in subjects:
    print(" ", Path(str(s.source)).name, s.data.shape)

files in the folder: ['.hidden.nii.gz', 'sub-10_bold.nii.gz', 'sub-1_bold.json', 'sub-1_bold.nii.gz', 'sub-2_bold.nii.gz', 'sub-3_bold.NII.GZ']
subjects loaded: 4
  sub-1_bold.nii.gz (1728, 10)
  sub-2_bold.nii.gz (1728, 12)
  sub-3_bold.NII.GZ (1728, 9)
  sub-10_bold.nii.gz (1728, 8)

The flat folder contains four visible NIfTI files, a .json sidecar, and a hidden file. The four images load in natural order: sub-1, sub-2, sub-3, sub-10. The uppercase .NII.GZ extension is accepted.

For an illustrative hierarchy, the next tree groups six independent simulated participants first by acquisition site (site-east, site-west) and then by an experimental condition (condition-a, condition-b). These are folder labels for synthetic data, not results from an experiment. Both sites contain a condition with the same name; its full path identifies the corresponding subgroup.

The two-iteration fit below demonstrates input routing and group construction, not a converged analysis. Its stopping warning is printed explicitly; increase the optimization budget for an actual analysis.

nested = root / "nested"
for site, condition, participant in [
    ("site-east", "condition-a", "01"),
    ("site-east", "condition-a", "02"),
    ("site-east", "condition-b", "03"),
    ("site-west", "condition-a", "04"),
    ("site-west", "condition-b", "05"),
    ("site-west", "condition-b", "06"),
]:
    directory = nested / site / condition
    directory.mkdir(parents=True, exist_ok=True)
    nib.save(nib.Nifti1Image(brain(T), affine), directory / f"sub-{participant}.nii.gz")

nested_subjects = htfa.load(nested)
for subject in nested_subjects:
    print(Path(subject.source).name, "->", subject.group_path)

# Short illustrative fit: group names and parent-child paths come from folders.
with warnings.catch_warnings(record=True) as fit_warnings:
    warnings.simplefilter("always")
    nested_result = htfa.fit(
        nested, K=2, max_global_iter=2, max_local_iter=3, random_state=0
    )
print("hierarchical model:", type(nested_result.model).__name__)
# Pass groups=False to htfa.fit to use a single shared template instead.
for warning in fit_warnings:
    print(f"{warning.category.__name__}: {warning.message}")

sub-01.nii.gz -> ('site-east', 'condition-a')
sub-02.nii.gz -> ('site-east', 'condition-a')
sub-03.nii.gz -> ('site-east', 'condition-b')
sub-04.nii.gz -> ('site-west', 'condition-a')
sub-05.nii.gz -> ('site-west', 'condition-b')
sub-06.nii.gz -> ('site-west', 'condition-b')
hierarchical model: HTFA
ConvergenceWarning: HTFA did not converge in max_global_iter=2 iterations; increase max_global_iter or threshold

3. BIDS-style datasets: pass a list of paths

The Brain Imaging Data Structure (BIDS) organizes neuroimaging files by participant and acquisition metadata such as task and run. Its participant/acquisition folders describe storage, which need not be the scientific hierarchy you want to fit. Select the intended runs explicitly to avoid treating these folders as group labels.

The synthetic directory below contains files named task-rest and task-other to illustrate selection by filename. Those names are placeholders, not recordings of rest or another task. The pattern selects one file per simulated participant. Collect the intended paths with glob and pass the list; each path becomes one Subject, in the order given. For a real study, selecting several runs from the same person does not make those runs independent participants.

bids = root / "bids"
for sub in ["sub-01", "sub-02", "sub-03"]:
    func = bids / sub / "func"
    func.mkdir(parents=True)
    nib.save(nib.Nifti1Image(brain(T), affine), func / f"{sub}_task-rest_bold.nii.gz")
    nib.save(nib.Nifti1Image(brain(T), affine), func / f"{sub}_task-other_bold.nii.gz")

runs = sorted(bids.glob("sub-*/func/*_task-rest_bold.nii.gz"))
print("collected:", [p.relative_to(bids).as_posix() for p in runs])
subjects = htfa.load([str(p) for p in runs])
print("subjects:", len(subjects), "| shapes:", [s.data.shape for s in subjects])

t0 = time.perf_counter()
with warnings.catch_warnings(record=True) as fit_warnings:
    warnings.simplefilter("always")
    res = htfa.fit(
        [str(p) for p in runs], K=2, max_global_iter=2, max_local_iter=3, random_state=0
    )
print(
    f"htfa.fit on the list: {time.perf_counter() - t0:.1f} s ->",
    type(res.model).__name__,
    "with",
    len(res.subjects),
    "subjects",
)
for warning in fit_warnings:
    print(f"{warning.category.__name__}: {warning.message}")

collected: ['sub-01/func/sub-01_task-rest_bold.nii.gz', 'sub-02/func/sub-02_task-rest_bold.nii.gz', 'sub-03/func/sub-03_task-rest_bold.nii.gz']
subjects: 3 | shapes: [(1728, 10), (1728, 10), (1728, 10)]
htfa.fit on the list: 0.2 s -> HTFA with 3 subjects
ConvergenceWarning: HTFA did not converge in max_global_iter=2 iterations; increase max_global_iter or threshold

4. Volumes as samples: stack=True

A folder of 3D volumes is, by default, one subject per volume with a single sample. stack=True treats all volumes in the folder as the samples of a single subject’s time series.

vols = root / "volumes"
vols.mkdir()
for i in range(6):
    nib.save(nib.Nifti1Image(brain(1)[..., 0], affine), vols / f"vol-{i:02d}.nii.gz")

per_volume = htfa.load(str(vols))
stacked = htfa.load(str(vols), stack=True)
print("default :", len(per_volume), "subjects, first data", per_volume[0].data.shape)
print("stack=True:", len(stacked), "subject, data", stacked[0].data.shape)

default : 6 subjects, first data (1728, 1)
stack=True: 1 subject, data (1728, 6)

5. Masks

With mask=None, voxels that are non-finite at any time point or zero at every time point are dropped (2D images keep every finite pixel). The ball image has no zeros, so every voxel is kept by default. mask= selects voxels explicitly: a NIfTI mask path or image (nonzero voxels kept), a boolean array with the shape of the image grid, a fitted nilearn NiftiMasker, or 'epi' / 'background' for nilearn’s automatic masks, computed for each subject (NIfTI input only). A NaN or infinite value inside an explicit mask raises ValueError.

mask_path = root / "mask.nii.gz"
nib.save(nib.Nifti1Image(inside.astype(np.uint8), affine), mask_path)


masker = NiftiMasker(mask_img=nib.load(mask_path)).fit()

options = {
    "mask=None": None,
    "mask path": str(mask_path),
    "mask nibabel image": nib.load(mask_path),
    "boolean array": inside,
    "fitted NiftiMasker": masker,
    "'epi'": "epi",
    "'background'": "background",
}
print(f"voxels in the ball: {int(inside.sum())} of {inside.size}")
for name, mask in options.items():
    (s,) = htfa.load(str(single), mask=mask)
    within = int((s.mask & inside).sum()) if s.mask is not None else None
    print(
        f"{name:20s}: V = {s.data.shape[0]:4d} voxels kept, inside the ball: {within}"
    )

voxels in the ball: 360 of 1728
mask=None           : V = 1728 voxels kept, inside the ball: 360
mask path           : V =  360 voxels kept, inside the ball: 360
mask nibabel image  : V =  360 voxels kept, inside the ball: 360
boolean array       : V =  360 voxels kept, inside the ball: 360
fitted NiftiMasker  : V =  360 voxels kept, inside the ball: 360
'epi'               : V =  304 voxels kept, inside the ball: 304
'background'        : V =  360 voxels kept, inside the ball: 360

The printed voxel counts compare the default mask, explicit ball masks, and automatic Nilearn masks. The default keeps all finite, nonzero voxels; the explicit masks select the ball exactly. Automatic masks depend on image intensity and spatial structure, so they can select a different boundary.

A zero-filled background is dropped by the default mask. Below, the background of the same image is set to exactly zero (as in skull-stripped data). The last cell puts a NaN at one voxel inside the ball and passes the ball as an explicit mask.

stripped = brain(T) * inside[..., None]
nib.save(nib.Nifti1Image(stripped, affine), root / "stripped.nii.gz")
(s,) = htfa.load(str(root / "stripped.nii.gz"))
print(
    "zero background, mask=None: V =",
    s.data.shape[0],
    "(ball has",
    int(inside.sum()),
    "voxels)",
)

nan_img = brain(T)
nan_img[6, 6, 6, 0] = np.nan
nib.save(nib.Nifti1Image(nan_img, affine), root / "nan.nii.gz")
try:
    htfa.load(str(root / "nan.nii.gz"), mask=inside)
except ValueError as err:
    print("NaN inside an explicit mask ->", type(err).__name__ + ":", err)

zero background, mask=None: V = 360 (ball has 360 voxels)
NaN inside an explicit mask -> ValueError: non-finite values inside the mask; adjust the mask or data

6. nilearn images

Any nibabel image works, including those returned by nilearn functions. Here nilearn smooths the image and concatenates volumes; the result is passed to htfa.fit (or htfa.load) as it is. The affine gives the coordinates, and HTFAResults.to_nifti returns images on the same grid.

image = htfa.require("nilearn.image")

smoothed = image.smooth_img(nib.load(single), fwhm=6)
concat = image.concat_imgs(
    [nib.Nifti1Image(brain(1)[..., 0], affine) for _ in range(5)]
)
for name, img in [("smooth_img", smoothed), ("concat_imgs", concat)]:
    (s,) = htfa.load(img)
    print(
        f"{name:12s}: image shape {img.shape} -> data {s.data.shape}, coordinate step "
        f"{np.diff(s.coords[:2], axis=0)[0].tolist()} mm"
    )

smooth_img  : image shape (12, 12, 12, 10) -> data (1728, 10), coordinate step [0.0, 0.0, 3.0] mm
concat_imgs : image shape (12, 12, 12, 5) -> data (1728, 5), coordinate step [0.0, 0.0, 3.0] mm

Both nilearn images load as one subject on the same 12 x 12 x 12 grid with 3 mm steps.

7. Arrays

• (X, Y, Z, T), (X, Y, Z) or (H, W) with no coords: an image grid. Coordinates are the grid indices of the kept voxels times spacing (a scalar or one value per axis).

• (V, T) plus coords (V, D): voxel-by-sample data at the given coordinates.

• A list of arrays (with a matching list of coordinate arrays): one subject per element.

arr4 = brain(T)
(s,) = htfa.load(arr4)
print(
    "(X,Y,Z,T) array     :",
    s.data.shape,
    "| first coords",
    s.coords[:2].tolist(),
    "| affine:",
    s.affine,
)
(s,) = htfa.load(arr4, spacing=(2.0, 2.0, 5.0))
print("with spacing=(2,2,5):", "first coords", s.coords[:2].tolist())
(s,) = htfa.load(arr4[..., 0])
print("(X,Y,Z) array       :", s.data.shape)
(s,) = htfa.load(rng.random((20, 30)))
print("(H,W) array         :", s.data.shape, "(one sample)")

V = 40
coords = rng.random((V, 2)) * 10
(s,) = htfa.load(rng.random((V, 15)), coords=coords)
print(
    "(V,T) + coords      :",
    s.data.shape,
    "| coords",
    s.coords.shape,
    "| image_shape:",
    s.image_shape,
)

Xs = [rng.random((V, t)) for t in (15, 12)]
Rs = [rng.random((V, 2)) * 10 for _ in Xs]
subjects = htfa.load(Xs, coords=Rs)
print("lists of arrays     :", [s.data.shape for s in subjects])
try:
    htfa.load(rng.random((V, 15)), coords=rng.random((V + 1, 2)))
except ValueError as err:
    print("mismatched coords ->", type(err).__name__ + ":", err)

(X,Y,Z,T) array     : (1728, 10) | first coords [[0.0, 0.0, 0.0], [0.0, 0.0, 1.0]] | affine: None
with spacing=(2,2,5): first coords [[0.0, 0.0, 0.0], [0.0, 0.0, 5.0]]
(X,Y,Z) array       : (1728, 1)
(H,W) array         : (600, 1) (one sample)
(V,T) + coords      : (40, 15) | coords (40, 2) | image_shape: None
lists of arrays     : [(40, 15), (40, 12)]
mismatched coords -> ValueError: with coords, data must be (V, T) matching coords (V, D); got data (40, 15), coords (41, 2)

htfa.fit accepts the same inputs; for arrays that are not on the generating scale, leave standardize at its default.

8. Inspect the loaded image and mask

These middle-plane views show the synthetic ball, the explicit mask, and the pixels recovered through Subject.to_image. Pixels outside the mask are NaN and appear as display background, not fitted observations. This visual check complements the shape and coordinate checks above; for real data, inspect anatomical alignment as well.

plt = htfa.require("matplotlib.pyplot")
(masked,) = htfa.load(single, mask=inside)
original = nib.load(single).get_fdata()[..., 0]
restored = masked.to_image(masked.data[:, 0]).get_fdata()
plane = shape[2] // 2
fig, axes = plt.subplots(1, 3, figsize=(10, 3.5))
for ax, values, title in zip(
    axes,
    [original[:, :, plane], inside[:, :, plane], restored[:, :, plane]],
    ["Synthetic input", "Explicit mask", "Loaded voxels on image grid"],
    strict=True,
):
    ax.imshow(
        values.T,
        origin="lower",
        cmap="gray",
        vmin=0,
        vmax=1 if values.dtype == bool else original.max(),
    )
    ax.set(title=title, xlabel="Voxel index x", ylabel="Voxel index y")
fig.tight_layout()
plt.show()

tmp.cleanup()