parcellate

lacuna.prepare.parcellate

Reduce a whole-brain connectome to a parcellated NxN connectivity matrix.

Functional branch only for now: load a voxelwise fMRI HDF5 connectome (same format consumed by lacuna run fnm), project timeseries onto parcel means using the atlas resampled to the connectome grid, compute a per-subject Pearson correlation matrix, average across subjects via Fisher z, and write the resulting group-level :class:~lacuna.core.data_types.ConnectivityMatrix as a BIDS-style TSV + JSON sidecar (mirrors the layout produced by lacuna run snm for its disconnectivity matrix).

Exposed via the lacuna prepare afnm CLI target.

ResolvedParcellation dataclass

An atlas ready for use by the parcellate pipeline.

Source code in src/lacuna/prepare/parcellate.py

@dataclass
class ResolvedParcellation:
    """An atlas ready for use by the parcellate pipeline."""

    short_name: str
    image: nib.Nifti1Image
    labels: dict[int, str]
    space: str | None
    source: str  # "registry" or "custom"

parcellate_functional(connectome_path, parcellations, output_dir, *, overwrite=False, verbose=False)

Compute group-average parcellated FC matrices from a voxelwise HDF5 connectome.

One output TSV+JSON pair is written per parcellation under output_dir. Returns the list of TSV paths.

Source code in src/lacuna/prepare/parcellate.py

def parcellate_functional(
    connectome_path: Path,
    parcellations: list[ResolvedParcellation],
    output_dir: Path,
    *,
    overwrite: bool = False,
    verbose: bool = False,
) -> list[Path]:
    """Compute group-average parcellated FC matrices from a voxelwise HDF5 connectome.

    One output TSV+JSON pair is written per parcellation under ``output_dir``.
    Returns the list of TSV paths.
    """
    connectome_path = Path(connectome_path)
    output_dir = Path(output_dir)
    output_dir.mkdir(parents=True, exist_ok=True)

    batch_files = list_connectome_batch_files(connectome_path)
    mask_info = read_mask_info(batch_files[0])
    mi = mask_info["mask_indices"]
    mask_affine = mask_info["mask_affine"]
    mask_shape = mask_info["mask_shape"]

    # Pre-resample every atlas once and build a weight matrix per parcellation.
    prepared: list[dict] = []
    for parc in parcellations:
        resampled = _resample_atlas_to_mask_grid(parc.image, mask_affine, mask_shape)
        atlas_values = resampled[mi[0], mi[1], mi[2]]
        W, region_ids = _build_parcel_weight_matrix(atlas_values)
        region_labels = [parc.labels.get(rid, f"Region{rid}") for rid in region_ids]
        N = len(region_ids)
        prepared.append(
            {
                "parc": parc,
                "W": W,
                "region_ids": region_ids,
                "region_labels": region_labels,
                "sum_z": np.zeros((N, N), dtype=np.float64),
                "n_subjects": 0,
            }
        )

    hasher = hashlib.sha256()
    for bf in batch_files:
        hasher.update(str(bf.name).encode())
        hasher.update(str(bf.stat().st_size).encode())

        if verbose:
            print(f"[parcellate] batch {bf.name}")
        with h5py.File(bf, "r") as hf:
            ts = hf["timeseries"][:]  # (n_sub, n_tp, n_vox)
        n_sub = ts.shape[0]

        for item in prepared:
            parcel_ts = ts @ item["W"]  # (n_sub, n_tp, N)
            corr = _per_subject_correlation(parcel_ts)
            corr = np.clip(corr, -0.999999, 0.999999)
            item["sum_z"] += np.arctanh(corr).sum(axis=0)
            item["n_subjects"] += n_sub
        del ts

    source_hash = hasher.hexdigest()[:16]
    written: list[Path] = []
    timestamp = time.strftime("%Y-%m-%dT%H:%M:%S")

    for item in prepared:
        parc: ResolvedParcellation = item["parc"]
        mean_z = item["sum_z"] / max(item["n_subjects"], 1)
        np.fill_diagonal(mean_z, 0.0)
        mean_z_f32 = mean_z.astype(np.float32)

        matrix_r = np.tanh(mean_z).astype(np.float32)
        np.fill_diagonal(matrix_r, 1.0)

        shared_meta = {
            "connectome_path": str(connectome_path),
            "parcellation_name": parc.short_name,
            "parcellation_source": parc.source,
            "parcellation_space": parc.space,
            "n_subjects": item["n_subjects"],
            "n_batches": len(batch_files),
            "modality": "functional",
            "source_hash": source_hash,
            "timestamp": timestamp,
            "mask_shape": list(mask_shape),
        }

        outputs = [
            (
                matrix_r,
                f"groupFC_{parc.short_name}",
                "pearson_r",
                f"method-parcellate_atlas-{parc.short_name}" f"_desc-fcgroupr_connmatrix",
            ),
            (
                mean_z_f32,
                f"groupFCz_{parc.short_name}",
                "fisher_z",
                f"method-parcellate_atlas-{parc.short_name}" f"_desc-fcgroupz_connmatrix",
            ),
        ]

        for matrix, name, value_type, base in outputs:
            connmat = ConnectivityMatrix(
                name=name,
                matrix=matrix,
                region_labels=item["region_labels"],
                matrix_type="functional",
                metadata=dict(shared_meta, value_type=value_type),
            )

            tsv_path = output_dir / f"{base}.tsv"
            json_path = output_dir / f"{base}.json"

            if tsv_path.exists() and not overwrite:
                raise FileExistsError(f"Output exists: {tsv_path}. Pass --overwrite to replace.")

            df = pd.DataFrame(
                connmat.matrix,
                index=connmat.region_labels,
                columns=connmat.region_labels,
            )
            df.to_csv(tsv_path, sep="\t")

            sidecar = {
                "Description": connmat.name,
                "MatrixType": connmat.matrix_type,
                "ValueType": value_type,
                "RegionLabels": connmat.region_labels,
                "Shape": list(connmat.matrix.shape),
                "Metadata": connmat.metadata,
            }
            with open(json_path, "w") as f:
                json.dump(sidecar, f, indent=2)

            written.append(tsv_path)
            if verbose:
                print(f"[parcellate] wrote {tsv_path}")

    return written

resolve_parcellations(parcel_atlases, custom_parcellation)

Turn CLI --parcel-atlases / --custom-parcellation into objects.

Source code in src/lacuna/prepare/parcellate.py

def resolve_parcellations(
    parcel_atlases: list[str] | None,
    custom_parcellation: list[list[str]] | None,
) -> list[ResolvedParcellation]:
    """Turn CLI ``--parcel-atlases`` / ``--custom-parcellation`` into objects."""
    resolved: list[ResolvedParcellation] = []

    for name in parcel_atlases or []:
        parc = load_parcellation(name)
        resolved.append(
            ResolvedParcellation(
                short_name=name,
                image=parc.image,
                labels=parc.labels,
                space=getattr(parc.metadata, "space", None),
                source="registry",
            )
        )

    for entry in custom_parcellation or []:
        short_name, nifti_path, labels_path, space = entry
        image = nib.load(nifti_path)
        labels = _load_labels_file(Path(labels_path))
        resolved.append(
            ResolvedParcellation(
                short_name=short_name,
                image=image,
                labels=labels,
                space=space,
                source="custom",
            )
        )

    if not resolved:
        raise ValueError("At least one of --parcel-atlases or --custom-parcellation must be given.")
    return resolved

run_parcellate_functional_cli(args)

Entry point for lacuna prepare afnm.

Source code in src/lacuna/prepare/parcellate.py

def run_parcellate_functional_cli(args) -> int:
    """Entry point for ``lacuna prepare afnm``."""
    try:
        parcellations = resolve_parcellations(
            getattr(args, "parcel_atlases", None),
            getattr(args, "custom_parcellation", None),
        )
    except (ValueError, KeyError, FileNotFoundError) as e:
        print(f"Error: {e}", file=__import__("sys").stderr)
        return 1

    try:
        parcellate_functional(
            connectome_path=args.connectome_path,
            parcellations=parcellations,
            output_dir=args.output,
            overwrite=getattr(args, "overwrite", False),
            verbose=getattr(args, "verbose_count", 0) > 0,
        )
    except (FileNotFoundError, FileExistsError, ValueError) as e:
        print(f"Error: {e}", file=__import__("sys").stderr)
        return 1
    return 0