SpectralBrain
Spectral Shape Analysis for Brain Structures

---

SpectralBrain computes, analyzes, and visualizes spectral shape descriptors of brain structures — cortical surfaces, subcortical meshes, hippocampal subfields, white-matter tracts, and point clouds from volumetric segmentations. It connects spectral geometry (the Laplace–Beltrami operator) to clinical neuroimaging, with one pipeline from FreeSurfer / HippUnfold output through statistically rigorous analysis to publication-ready figures.

The core idea: any input geometry — cortical surface, subcortical/hippocampal mesh, tract, or point cloud — is passed through the Laplace–Beltrami operator to obtain its eigenpairs {λ, φ}, from which pose- and mesh-free spectral descriptors (ShapeDNA, HKS, SI-HKS, WKS, GPS, BKS, functional maps, wavelets) are read out.

Statement of need

Volumetric and thickness measures collapse a structure's shape to a few scalars and are sensitive to registration and voxel size. Intrinsic spectral descriptors derived from the Laplace–Beltrami operator (LBO) — ShapeDNA, the Heat/Wave Kernel Signatures, and relatives — characterize shape independently of pose and parameterization, capturing geometry that volume alone misses. They are well established in geometry processing but scattered across research code, rarely packaged with the I/O, multi-site harmonization, correct multiple-comparison statistics, and rendering that a neuroimaging study needs end to end. SpectralBrain fills that gap as a single, tested library, with a primary focus on the hippocampus in mesial temporal lobe epilepsy, while remaining general to any brain surface or point cloud.

Key capabilities

• Spectral descriptors — ShapeDNA, Heat Kernel Signature (HKS), Scale-Invariant HKS, Wave Kernel Signature (WKS), Global Point Signature (GPS), Bates–Kornfeld Signature (BKS) and its inverse, functional maps, and more — all from the LBO eigenpairs of a mesh or point cloud.
• Input-agnostic I/O — FreeSurfer surfaces and morphometry, GIfTI (.surf.gii / .func.gii / .shape.gii), NIfTI / MGZ volumes and labels, HippUnfold v1 & v2 outputs, .ply / .obj / .stl / .vtk, HDF5, and point clouds, with automatic format detection.
• Cohort loading — BIDS / derivatives, FreeSurfer SUBJECTS_DIR, or an explicit list, loaded in parallel and stacked for group analysis; FreeSurfer measures can be resampled onto a common template; TractSeg bundle masks import directly as point clouds or isosurface meshes.
• Statistics done right — vertex-wise tests with genuine family-wise error control (max-statistic permutation), FDR, partial correlations with correct degrees of freedom, TFCE, the analytic DeLong AUC test, BCa bootstrap, ComBat / ComBat-GAM harmonization, and six PyMC Bayesian models.
• Contiguous clustering & atlas comparison — a distance-dependent Chinese Restaurant Process (ddCRP) with a Normal-Inverse-Wishart collapsed marginal likelihood (spatial and fPCA-functional variants), consensus clustering, data-driven hyperparameter autotuning, and a cluster-vs-atlas suite (ARI/AMI/ NMI/VI, Dice/Jaccard, within-parcel homogeneity) reported against size-matched random parcellations, plus spatially-aware ARI (spARI) and eigenstrapping/ BrainSMASH spatial nulls for bounded surfaces.
• Publication figures — a template-free six-view 3D renderer (vedo), unfolded flat-maps, Bayesian-posterior plots, advanced 3D tractography (direction- encoded/scalar streamlines, bundle surfaces with HKS/WKS overlays, multi-POV montages), and parcellation-vs-clustering grids (views × labelings) for hippocampi, brains, and bundles.

Installation

pip install spectralbrain

Optional feature sets (extras):

pip install "spectralbrain[bayesian]"   # PyMC, nutpie, NumPyro, BlackJAX, ArviZ
pip install "spectralbrain[viz]"        # vedo, fury, trimesh, cmcrameri, …
pip install "spectralbrain[gpu]"        # torch, CuPy, JAX (CUDA)
pip install "spectralbrain[neuro]"      # nilearn, dipy, pybids, templateflow, …
pip install "spectralbrain[tuning]"     # optuna (ddCRP autotuning; random fallback)
pip install "spectralbrain[full]"       # everything above

Requires Python 3.11–3.12.

API at a glance

The core API is on the top-level package; heavier statistics and visualization live in submodules you import explicitly (mirroring scipy.stats):

import spectralbrain as sb               # meshes, descriptors, I/O
import spectralbrain.statistics as sbstats   # frequentist + Bayesian
import spectralbrain.viz as sbviz             # 3D / 2D figures

Quick start

1 — Mesh → eigenpairs → descriptors

import spectralbrain as sb

# A BrainMesh from vertices (N, 3) and faces (M, 3).
vertices, faces = sb.io.load_gifti_surface("path/to/surf/gii")
mesh   = sb.BrainMesh(vertices, faces)
decomp = mesh.decompose(k=100)                       # 100 LBO eigenpairs

hks = sb.compute_hks(decomp, t_values=[1.0, 10.0, 100.0])   # (N, 3)
wks = sb.compute_wks(decomp, n_energies=50)                 # (N, 50)
dna = sb.compute_shapedna(decomp)                           # (k-1,) global

Point clouds work identically — sb.BrainPointCloud(points).decompose(k=...).

2 — Compare two shapes

d = sb.shapedna_distance(dna_a, dna_b)   # pose-invariant spectral distance

3 — Vertex-wise group statistics with FWER control

import spectralbrain.statistics as sbstats

# controls, patients : (n_subjects, n_vertices) descriptor fields
res = sbstats.vertexwise_permutation(
    controls, patients,
    n_permutations=5000,
    correction="max",      # family-wise error via the max-statistic null
    seed=0,
)
significant = res.significant          # boolean mask, FWER-controlled

correction="fdr" and "none" are also available; vertexwise_ttest defaults to Welch's t-test.

4 — Compare two classifiers (analytic DeLong)

auc_new, auc_ref, p = sbstats.auc_comparison_delong(y_true, scores_new, scores_ref)

5 — Six-view 3D render (template-free)

import spectralbrain.viz as sbviz

fig = sbviz.plot_hippocampus_sixview(
    mesh, scalars=hks[:, 1],
    cmap="plasma", scalar_bar_title="HKS(t=10)",
    save="hipp_sixview.png",
)
# Pick any subset/order of the six canonical views:
fig = sbviz.plot_hippocampus_sixview(mesh, scalars=hks[:, 1],
                                     views=("superior", "left_lateral"))

Views: anterior, posterior, inferior, superior, left_lateral, right_lateral. It renders any surface — HippUnfold v2 den-8k, an aseg ROI mesh, or a whole cortical hemisphere — with no bundled template, so scalar↔vertex correspondence is guaranteed.

6 — Bayesian sparse regression (extra: [bayesian])

from spectralbrain.statistics import HorseshoeRegression

model = HorseshoeRegression(tau_prior=0.5).fit(X, y, sampler="nuts")
importance = model.feature_importance()   # sparse posterior shrinkage

7 — Contiguous ddCRP clustering + atlas comparison

The sampler maintains incremental per-component sufficient statistics (the cluster scatter is never recomputed inside the candidate loop) and caches the NIW marginal per cluster, so it scales to dense surfaces (~1-3 s/draw at ~7k vertices). Pass n_components to PCA-whiten the descriptors for a further speedup and better conditioning.

import spectralbrain.statistics as sbstats

# H : (V, d) per-vertex descriptors (e.g. fused HKS/WKS); vertices/faces from the mesh.
res = sbstats.cluster_ddcrp(H, faces=faces, vertices=vertices,
                            decay_kind="exponential")   # decay uses real edge distances
print(res)                                              # ClusterResult(method='ddcrp', n_clusters=…)

# Tune hyperparameters from the data (optuna if installed, else random search):
tuned = sbstats.autotune_ddcrp(H, faces=faces, vertices=vertices, n_trials=40)
final = tuned.cluster_result                            # refit ClusterResult

# Is the partition more atlas-like than a size-matched random parcellation?
report = sbstats.cluster_atlas_concordance(final.labels, atlas_labels,
                                           faces=faces, coords=vertices, n_null=1000)
print(report["metrics"]["ari"], report["ari_null"]["z"], report["spARI"]["spARI"])

8 — Tracts in 3D and parcellation-vs-clustering grids (extra: [viz])

import spectralbrain.viz as sbviz

# (a) advanced 3D tractography from several POVs (FURY; DEC orientation colours)
sl, _ = sbviz.load_streamlines("CST_left.trk", to_space="world")   # needs dipy
fig, _ = sbviz.streamlines_multiview(sl, views=("left", "anterior", "superior", "oblique"),
                                     out_path="cst_multiview.pdf")

# bundle surface from a TractSeg mask, with an HKS overlay
V, F = sbviz.mask_to_mesh(mask, affine=affine)
hks = sbviz.spectral_overlay(V, F, kind="hks")
sbviz.render_bundle_surface(V, F, scalars=hks, view="oblique", out_path="cst_hks.png")

# (b) 3D grid: columns = views, rows = a reference parcellation then each clustering
fig, meta = sbviz.plot_parcellation_vs_clusters(
    vertices, faces, parcellation=atlas_labels,
    clusterings={"ddCRP": res, "Leiden": leiden_res},   # ClusterResult or label arrays
    views=["left_lateral", "anterior", "superior"],
    save="parcellation_vs_clusters.png",
)

The grid works identically for hippocampi, whole brains, and bundle surfaces — it operates on any (vertices, faces) mesh plus a dict of per-vertex labelings.

Loading a cohort

import spectralbrain as sb

# BIDS / derivatives (one file per subject):
files = sb.discover_bids("/data/derivatives/hippunfold",
                         "sub-{sub}/surf/sub-{sub}_hemi-L_*thickness.shape.gii")
group = sb.load_group(files, mode="maps", n_jobs=8)
res   = sb.group_comparison(group, group.covariate("group"), test="ttest")

# FreeSurfer SUBJECTS_DIR, resampled to a common template:
group = sb.load_group_freesurfer("/data/fs", measure="thickness",
                                 template="fsaverage", n_jobs=8)

# TractSeg bundle masks → meshes ready for .decompose():
bundles = sb.load_tractseg("/data/sub-01/tractseg_output", output="mesh")
decomp  = bundles["CST_left"].decompose(k=80)

mode="pipeline" runs load → decompose → descriptor per subject (with an optional GPU backend=); mode="maps" stacks vertex-corresponded fields.

Compute backends

Eigen-decomposition and Bayesian sampling run on pluggable backends:

from spectralbrain.backends import TorchBackend       # or CupyBackend, JaxBackend
decomp = mesh.decompose(k=200, backend=TorchBackend())  # GPU eigsolve

Bayesian models accept sampler="auto" | "nuts" | "nutpie" | "numpyro" | "blackjax".

Structure & workflow

From inputs to inference, SpectralBrain is one coherent pipeline — load geometry → build the LBO & decompose → compute descriptors → group statistics → assess & visualize — assembled from a small set of focused, independently usable subpackages.

(a) inputs — geometry, operator, cohort, atlas, spectral descriptors, and inference; (b) the five-step main workflow, ending in a lateralization effect-size read-out; (c) the modular structure — core · spectral · io · statistics · viz · backends.

The pipeline in practice

The five stages above are not just a diagram — they run end to end on real study designs. The figure below carries two analyses through the whole pipeline side by side and adds a parallel Bayesian lane, so the same data are assessed with both frequentist tools (max-statistic permutation, TFCE, DeLong AUC) and Bayesian ones (hierarchical models, horseshoe priors, HDI + ROPE, LOO).

An end-to-end walkthrough. Columns are worked analyses — hippocampal lateralization in MTLE-HS (HippUnfold den-8k surfaces, L vs R), cortical morphometry (MTLE-HS vs. controls, FreeSurfer / Schaefer-200), and a Bayesian lane (hierarchical model → horseshoe priors → variable selection → posterior effect size with HDI + ROPE → posterior-predictive check and LOO model comparison). Rows are the five pipeline stages — load geometry → decompose → descriptors → group statistics → assess — tied on the left to the subpackages that implement them (io · core · spectral · statistics · viz). Panels are illustrative, generated from synthetic example data to show the shape and flow of an analysis, not empirical results.

Documentation map

Subpackage	What it provides
spectralbrain (top level)	BrainMesh, BrainPointCloud, decompose, all compute_* descriptors, distances, I/O, cohort loading
spectralbrain.io	loaders/savers, BIDS & FreeSurfer discovery, load_group, template resampling, TractSeg import, parcellation
spectralbrain.statistics	vertex-wise tests, TFCE, effect sizes, RSA, classification, ComBat(-GAM), normative models, bootstrap & null models, six Bayesian models
spectralbrain.backends	CPU / Torch / CuPy / JAX eigensolvers; PyMC / nutpie / NumPyro / BlackJAX samplers
spectralbrain.viz	six-view 3D renderer, unfolded flat-maps, cluster overlays, Bayesian-posterior and general scientific plots

Validation

examples/example_clustering_tracts.py is a runnable end-to-end pipeline (HippUnfold hippocampus -> HKS/WKS -> ddCRP -> atlas stats -> 3D figures; synthetic fallback if no data). validation/validate_spari.py checks the spatially-aware Rand index: analytical properties (identity, label invariance, symmetry, reduction to the exact ARI as the spatial scale vanishes, chance level, locality monotonicity) run anywhere, plus an R cell-for-cell bridge (--with-r) against the published spARI package that justifies flipping validated_against_R=True.

Development

git clone https://github.com/rdneuro/spectralbrain
cd spectralbrain
uv sync --group dev          # or: pip install -e ".[full]" + dev tools
uv run pytest                # run the test suite
uv run ruff check src/ tests/

Citing

If SpectralBrain contributes to your work, please cite it. Someday, maybe, if we feel lucky, a JOSS paper will be submited; until then, cite the archived release on Zenodo:

Debona, R. SpectralBrain: Spectral Shape Analysis for Brain Structures. Zenodo. https://doi.org/10.5281/zenodo.21090748

The DOI 10.5281/zenodo.21090748 always resolves to the latest release. See CITATION.cff for a machine-readable citation — GitHub's "Cite this repository" button reads it.

License

MIT — see LICENSE.