Getting example data¶

Download notebook

dmri predict reads the four-file FSL layout: data.nii.gz, nodif_brain_mask.nii.gz, bvals and bvecs. Public diffusion datasets rarely ship exactly that, so this example takes an open dataset from download to a finished prediction.

Two datasets are covered. Prefer CFIN for a realistic run: it is a real brain and its shells match the b = 1000 / 2000 the pretrained checkpoints were trained around. The ISBI 2013 phantom is useful for a different reason -- it has known ground-truth crossing fibres -- but see the rounding trap below.

	`cfin_2shell`	`isbi_phantom`
content	real brain	software phantom
matrix	96 x 96 x 19	50 x 50 x 50
volumes	67	64
shells	b = 0, 1000, 2000	b = 0, 1500, 2500
voxels in mask	55,314	125,000
download	174 MB	22 MB

For comparison, a full HCP-style subject is around 285 MB with 105 volumes. Network work scales with the number of measurements, so fewer volumes is a faster run.

This notebook is not executed during the documentation build -- it downloads 174 MB and runs a full prediction. Run it yourself to reproduce the outputs.

Download¶

CFIN is published as one 496-volume acquisition: 33 gradient directions repeated across 15 shells from b = 200 to b = 3000. Download it first, then keep only the two shells you need.

mkdir -p ~/dmri_data_small/cfin_raw && cd ~/dmri_data_small/cfin_raw

BASE=https://digital.lib.washington.edu/researchworks/bitstream/handle//1773/38488
N=__DTI_AX_ep2d_2_5_iso_33d_20141015095334_4

wget -c $BASE/$N.nii
wget -c $BASE/$N.bval
wget -c $BASE/$N.bvec

The double slash in handle//1773/38488 is not a typo -- it is what the host serves, and what dipy's own fetcher uses. A single slash returns 404.

The phantom, if you want it too:

mkdir -p ~/dmri_data_small/isbi_raw && cd ~/dmri_data_small/isbi_raw

BASE=https://digital.lib.washington.edu/researchworks/bitstream/handle/1773/38465

wget -c $BASE/phantom64.nii.gz
wget -c $BASE/phantom64.bval
wget -c $BASE/phantom64.bvec

Convert to the FSL layout¶

Four differences have to be reconciled, and each one is a silent failure if missed:

b-value files are named .bval / .bvec, and b-vectors may be N-by-3 rather than the 3-by-N that FSL uses;
b0 volumes are sometimes labelled b = 5 rather than 0, but the normalisation step tests bvals == 0 exactly;
there is no brain mask;
multi-shell acquisitions carry far more volumes than a prediction needs.

from pathlib import Path

import nibabel as nb
import numpy as np


def brain_mask(b0):
    """A binary mask from the mean b0 volume.

    Uses dipy's median_otsu, which is what FSL's `bet` approximates for this
    purpose; falls back to a plain Otsu threshold so the function does not hard
    depend on dipy.
    """
    try:
        from dipy.segment.mask import median_otsu

        _, mask = median_otsu(b0, median_radius=2, numpass=1)
        return mask.astype(np.uint8)
    except ImportError:
        from scipy.ndimage import binary_fill_holes, gaussian_filter
        from skimage.filters import threshold_otsu

        smooth = gaussian_filter(b0.astype(np.float32), 1.0)
        return binary_fill_holes(smooth > threshold_otsu(smooth)).astype(np.uint8)


def to_fsl_layout(data, bval, bvec, out, shells=None, tolerance=60.0):
    """Write data.nii.gz, bvals, bvecs and nodif_brain_mask.nii.gz into `out`.

    Args:
        data: 4-D NIfTI of diffusion-weighted volumes.
        bval, bvec: the b-value and b-vector text files beside it.
        out: destination folder, created if needed.
        shells: keep only these nominal b-values; b0 is always kept. None keeps
            everything.
        tolerance: how far a b-value may sit from a nominal shell and still
            count as part of it.
    """
    img = nb.load(data)
    bvals = np.loadtxt(bval).ravel()
    bvecs = np.loadtxt(bvec)
    if bvecs.shape[0] != 3:
        bvecs = bvecs.T
    if bvecs.shape[1] != bvals.size:
        raise ValueError(f"bvecs {bvecs.shape} do not match {bvals.size} bvals")

    volumes = np.asanyarray(img.dataobj, dtype=np.float32)
    if volumes.shape[-1] != bvals.size:
        raise ValueError(f"{volumes.shape[-1]} volumes but {bvals.size} bvals")

    keep = bvals <= tolerance  # b0 always survives
    if shells:
        for shell in shells:
            keep |= np.abs(bvals - shell) <= tolerance
        print(f"dropping {int((~keep).sum())} of {bvals.size} volumes")
    else:
        keep[:] = True

    volumes, bvals, bvecs = volumes[..., keep], bvals[keep], bvecs[:, keep]
    # The loader tests `bvals == 0`, so snap near-zero labels down.
    bvals[bvals <= tolerance] = 0.0
    # And zero the direction on unweighted volumes, as FSL does.
    bvecs[:, bvals == 0] = 0.0

    mask = brain_mask(volumes[..., bvals == 0].mean(axis=-1))

    out = Path(out).expanduser()
    out.mkdir(parents=True, exist_ok=True)
    nb.Nifti1Image(volumes, img.affine, img.header).to_filename(out / "data.nii.gz")
    nb.Nifti1Image(mask, img.affine, img.header).to_filename(
        out / "nodif_brain_mask.nii.gz"
    )
    np.savetxt(out / "bvals", bvals[None, :], fmt="%g")
    np.savetxt(out / "bvecs", bvecs, fmt="%.6f")

    shell_names, counts = np.unique(np.round(bvals / 100) * 100, return_counts=True)
    print(f"wrote {out}")
    print(f"  data   {volumes.shape}")
    print(
        f"  shells {dict(zip(shell_names.astype(int).tolist(), counts.tolist(), strict=True))}"
    )
    print(f"  mask   {int(mask.sum()):,} of {mask.size:,} voxels in brain")
    return out

Run it on the CFIN download, keeping only b = 1000 and b = 2000:

raw = Path("~/dmri_data_small/cfin_raw").expanduser()
stem = "__DTI_AX_ep2d_2_5_iso_33d_20141015095334_4"

folder = to_fsl_layout(
    data=raw / f"{stem}.nii",
    bval=raw / f"{stem}.bval",
    bvec=raw / f"{stem}.bvec",
    shells=[1000, 2000],
    out="~/dmri_data_small/cfin_2shell",
)

dropping 429 of 496 volumes


wrote /home/manuel/dmri_data_small/cfin_2shell
  data   (96, 96, 19, 67)
  shells {0: 1, 1000: 33, 2000: 33}
  mask   55,314 of 175,104 voxels in brain

Predict¶

dmri predict ~/dmri_data_small/cfin_2shell \
  --output-subdir out_fast --quality fast --non-interactive

The run prints a clickable viewer link at the end; open out_fast/view_results.html to page through the maps.

The phantom needs `round_bvals=false`¶

The loader rounds b-values to the nearest 1000 by default. That is what lets a real acquisition labelled b = 5 be recognised as b0, but applied to the phantom's 1500 and 2500 it is destructive: 1500 rounds up, 2500 rounds down (numpy uses banker's rounding, so round(2.5) == 2), and both shells collapse onto a single nominal b = 2000 -- exactly the multi-shell contrast a ball-and-stick fit depends on.

from dmri.eval.load_data import process_bvals

for name, values in [
    ("cfin_2shell", np.array([0.0, 1000.0, 2000.0])),
    ("isbi_phantom", np.array([0.0, 1500.0, 2500.0])),
    ("HCP-style (b0 labelled 5)", np.array([5.0, 990.0, 2010.0])),
]:
    print(f"{name:26} {values.tolist()} -> {process_bvals(values, True).tolist()}")

cfin_2shell                [0.0, 1000.0, 2000.0] -> [0.0, 1000.0, 2000.0]
isbi_phantom               [0.0, 1500.0, 2500.0] -> [0.0, 2000.0, 2000.0]
HCP-style (b0 labelled 5)  [5.0, 990.0, 2010.0] -> [0.0, 1000.0, 2000.0]

So the phantom has to be run with rounding disabled:

dmri predict ~/dmri_data_small/isbi_phantom \
  --output-subdir out_fast --quality fast --non-interactive \
  --set evaluation.input.round_bvals=false

Its b0 is saturated in every voxel -- it is a synthetic block with signal throughout -- so the generated mask covers all 125,000 voxels. That is correct here, not a masking failure, but there is no background to skip.

Compatibility caveats¶

Both datasets sit inside the b-value range the checkpoints were trained on (0 to 4000, with typical shells at 1000, 2000 and 3000). That is necessary, not sufficient -- see the compatibility warning in the prediction guide.

Two specifics worth knowing:

Each has a single b0 volume, against five in a typical HCP-style subject, so the b0 normalisation is noisier.
The phantom's 1500 / 2500 shells are in range but are not among the typical values the training distribution favours, so results there are less representative than the CFIN run.