Load multi-pose recordings¶

This example shows how to assemble a multi-pose fUSI recording that was not acquired through an Iconeus SCAN file, by loading each pose as a separate NIfTI file and stacking them into a single VoxelData array with pose-dependent voxel-to-world geometry.

The Pepe, Mariani et al. (2026) dataset contains transcranial mouse resting-state recordings acquired with a linear probe stepped across several positions. In the fUSI-BIDS export, each probe position is one file, distinguished by the chunk- entity (chunk-0, chunk-1, ...); this is exactly the "Other Systems" case described in the Multi-Pose Imaging guide — data must be assembled manually rather than loaded as one file, unlike Iconeus SCAN's 3Dscan/4Dscan modes.

Fetch one recording's chunks¶

We select one subject, session, and acquisition so only the four chunks belonging to a single recording are downloaded, rather than the full multi-gigabyte dataset.

import re
from pathlib import Path

import matplotlib as mpl
import xarray as xr

import confusius as cf

# Adapt background color to the current Matplotlib style.
bg_color = mpl.colors.to_hex(mpl.rcParams["figure.facecolor"])

xr.set_options(display_expand_data=False)

subject = "cp230420a"
session = "1MEDISOses5"
acq = "3dfusi"

bids_root = cf.datasets.fetch_pepe_mariani_2026(
    datasets="rawdata",
    subjects=subject,
    sessions=session,
    acqs=acq,
    datatypes="fusi",
)

If you use this dataset in your work, please cite the following source:

Pepe, C., Mariani, J.-C., Urosevic, M., Gini, S., Stuefer, A., Ricci, F., Galbusera, A., Iurilli, G., & Gozzi, A. 
(2026). Structural and dynamic embedding of the mouse functional connectome revealed by functional ultrasound 
imaging (Fusi). https://doi.org/10.64898/2026.02.05.704055

Locate and load every chunk¶

Each chunk is an ordinary single-pose fUSI recording — (time, k, j, i) with a singleton k, since a linear probe images one elevation slice per pose — with its own voxel_to_world affine derived from that file's own NIfTI header. We sort by the chunk- index so poses end up in acquisition order, then load and compute each one (they're small enough to fit comfortably in memory).

fusi_dir = Path(bids_root) / f"sub-{subject}" / f"ses-{session}" / "fusi"
chunk_glob = (
    f"sub-{subject}_ses-{session}_task-rest_acq-{acq}_probe-linear_chunk-*_pwd.nii.gz"
)
chunk_paths = sorted(
    fusi_dir.glob(chunk_glob),
    key=lambda p: int(re.search(r"chunk-(\d+)", p.name).group(1)),
)
chunks = [cf.load(p).compute() for p in chunk_paths]

npose = len(chunks)
print(f"{npose} poses, each shaped {chunks[0].dims} = {chunks[0].shape}")

4 poses, each shaped ('time', 'k', 'j', 'i') = (1500, 1, 82, 128)

Each chunk's own origin confirms the poses sit at different physical positions along elevation (z), stepped by about 1 mm, and were acquired with slightly offset start times — this is the per-pose geometry we want to preserve, not collapse into one shared affine. y/x origins don't vary across chunks, so we only show z/time.

origins = [chunk.fusi.origin for chunk in chunks]
print("z origins:", tuple(origin["z"] for origin in origins))
print("time origins:", tuple(origin["time"] for origin in origins))

z origins: (0.00559999980032444, -0.9943992495536804, -1.9943994283676147, -2.9953997135162354)
time origins: (0.0, 0.6, 1.2, 1.8)

Stack the poses into one pose-dependent DataArray¶

stack_poses concatenates independently loaded single-pose grids into one DataArray with a pose dimension, one voxel-to-world affine per pose. Poses were also acquired sequentially rather than simultaneously, so each chunk's time coordinate is offset from the others by a fraction of the repetition time; when per-pose timestamps differ, stack_poses keeps them as a pose-dependent (time, pose)-shaped time coordinate rather than collapsing them into one shared time axis.

multipose = cf.multipose.stack_poses(chunks)

Pose-transparent vs. pose-specific geometry¶

Some geometry queries are well-defined without picking a pose: voxel spacing must be identical across poses (a stacked affine with differing scale is rejected at construction), so .fusi.spacing works directly on the multi-pose array.

Origin and direction, however, are inherently single-grid concepts—there is no one answer for "the origin" of a stack of differently-positioned grids—so they require selecting a scalar pose first.

pose0_origin = multipose.isel(pose=0).fusi.origin
pose3_origin = multipose.isel(pose=3).fusi.origin
print("spacing (pose-transparent):", multipose.fusi.spacing)
print("z origins (pose 0, 3):", (pose0_origin["z"], pose3_origin["z"]))
print("time origins (pose 0, 3):", (pose0_origin["time"], pose3_origin["time"]))

spacing (pose-transparent): {'time': 2.400000000000091, 'pose': 1.0, 'k': 0.4000000059604645, 'j': 0.10000000149011612, 'i': 0.10999999940395355}
z origins (pose 0, 3): (0.00559999980032444, -2.9953997135162354)
time origins (pose 0, 3): (0.0, 1.8)

Visualize each pose¶

We plot the temporal mean of each pose side by side with plot_volume's slice_mode="pose", which labels each panel by its pose value.

mean_db = multipose.mean("time").fusi.scale.db().isel(k=0)

plotter = cf.plotting.plot_volume(
    mean_db,
    slice_mode="pose",
    cbar_label="Power Doppler (dB)",
    bg_color=bg_color,
    yincrease=True,
)

Consolidate into a single volume¶

Since every pose here shares the same rotation and is offset by a pure translation along elevation, consolidate_poses can merge pose and the swept voxel dimension (k, the default) into one physically ordered axis — the same operation load_scan-produced 3Dscan/4Dscan data goes through, now working directly off the DataArray's own per-pose geometry rather than a separately stored affine.

consolidated = cf.multipose.consolidate_poses(multipose)
consolidated

The consolidated volume has a genuine k extent of npose elevation slices (each pose here contributes exactly one, since the probe is linear) and no more pose dimension, positioned by real physical location rather than acquisition order.

plotter = cf.plotting.plot_volume(
    consolidated.mean("time").fusi.scale.db(),
    slice_mode="z",
    cmap="gray",
    cbar_label="Power Doppler (dB)",
    bg_color=bg_color,
    yincrease=True,
)

Total running time: 12.9 s

Launch in Binder Download .py Download .ipynb