A ground-truth testbed for tractography algorithms — with a faithful port of MRtrix3's probabilistic tracker, phantoms whose answer is known analytically, a plug-in interface for new algorithms, and an experimental quantum-optimisation track that reports where quantum methods lose as readily as where they win.
Built by Transconnectome (SNU Connectome Lab).
Tractography's central problem is that you usually cannot tell whether the answer is right. Whole-brain reconstructions look plausible and are substantially wrong; validation studies consistently find large false-positive rates that are invisible without an independent reference. This repository collects references of three kinds and scores against them with one shared metric implementation, so that a new algorithm is measured rather than admired.
| tier | ground truth | certainty | cost |
|---|---|---|---|
| A — synthetic | the generating curve itself, exact | highest | seconds, offline, no downloads |
| B — challenge data | simulated whole-brain and physical phantoms | high for simulated | ~250 MB download |
| C — biological | chemical tracing, histology, polarised light imaging | most realistic, least directly comparable | registration, often per-request |
Tier C carries a caveat no metric repairs: tracer connectivity is graded and directed; tractography is neither. A perfect score is unattainable in principle wherever forward and backward projections differ.
git clone https://github.com/Transconnectome/tractobench
cd tractobench
pip install -e ".[quantum,dev]"
pytest -q # 19 tests, ~2 sfrom tractobench import phantoms, signal, fod, metrics
from tractobench.tracker import FODField, TrackingConfig, track
from tractobench.algorithms import IFOD1
import numpy as np
ph = phantoms.crossing(angle_deg=90.0) # exact ground truth, generated here
scheme = signal.DiffusionScheme.single_shell(n_dirs=60, bval=3000.0)
data = ph.simulate(scheme, snr=30.0) # multi-tensor forward model + Rician noise
coeffs = fod.fit_csd(data, scheme, lmax=8, mask=ph.mask)
coeffs, _ = fod.normalise_fod_scale(coeffs, mask=ph.mask)
field = FODField(coeffs, ph.affine, mask=ph.mask)
cfg = TrackingConfig.for_voxel_size(ph.voxel_mm, "ifod1") # MRtrix3 defaults, scaled
sls = track(IFOD1(field, cfg), ph.all_seeds(60, np.random.default_rng(0)),
np.random.default_rng(1))
print(metrics.summarise(ph, [s.points for s in sls])["connections"])Subclass DirectionGetter and implement two methods. Propagation, bidirectional
tracking, length bookkeeping, mask exit and termination accounting are shared, so a
comparison isolates the direction-choosing step.
from tractobench.tracker import DirectionGetter, StepProposal, Termination
class MyTracker(DirectionGetter):
name = "my_tracker"
def initial_direction(self, pos_mm, rng):
... # a starting direction, or None if the seed is unusable
def propose(self, pos_mm, direction, rng):
... # StepProposal(direction, Termination.CONTINUE)StepProposal.positions lets an arc-based method (iFOD2-style) commit several vertices
at once. Positions are in world millimetres everywhere; only FODField knows the
affine.
crossing, kissing, fanning, bending, bottleneck — each isolating one geometric
ambiguity. Two of them are deliberately unsolvable by local methods:
kissing— two arcs tangent at the origin (verified: centreline separation 1×10⁻⁴ mm, tangent angle 0.24°). At the contact point the FOD contains a single population, so a tracker has no data-driven basis for choosing a branch.bottleneck— two bundles through a shared channel. All four endpoint pairings are equally consistent with the orientation field, so half of any connections found must be false. This is an information limit, not a resolution limit.
A high valid-connection ratio on those two indicates an informative prior, not a better algorithm. The phantoms are diagnostic instruments, not a leaderboard.
python benchmarks/run_phantom_benchmark.py — 60 seeds per bundle, SNR 30, lmax 8:
| phantom | SD_STREAM VC | iFOD1 VC | SD_STREAM IC | iFOD1 IC |
|---|---|---|---|---|
| bending (r=10 mm) | 0.92 | 0.35 | 0.00 | 0.00 |
| crossing 90° | 0.76 | 0.53 | 0.00 | 0.00 |
| crossing 45° | 0.26 | 0.22 | 0.41 | 0.21 |
| kissing | 0.67 | 0.33 | 0.09 | 0.06 |
| fanning | 0.40 | 0.26 | 0.00 | 0.00 |
| bottleneck | 0.38 | 0.12 | 0.05 | 0.09 |
VC = valid connections, IC = invalid connections, as fractions of all streamlines. The deterministic tracker wins on single-fibre geometry, where probabilistic sampling only adds dispersion; at an under-resolved 45° crossing it makes roughly twice as many false connections, while iFOD1 disperses into unscoreable streamlines instead. Neither result is a defect — they are the two failure modes the methods trade against.
Port fidelity. iFOD1 accepts a step after 8.16–8.24 FOD evaluations across all six
phantoms, against the calibrator's own independent prediction of 8.115. The port is
validated against MRtrix3 source pinned at commit 6735cfa, with file:line citations
throughout docs/mrtrix3_probabilistic_tractography.md.
python benchmarks/run_ismrm2015_benchmark.py --n-seeds 20000. The acquisition is 32
directions at b=1000, which caps the spherical-harmonic order at lmax=6 (lmax=8 would
need 45 directions), so these numbers are bounded by the data as much as by the trackers.
| bundles detected | mean OL | mean OR | mean F1 | unassigned | wall | |
|---|---|---|---|---|---|---|
| iFOD1 | 25 / 25 | 0.299 | 0.170 | 0.391 | 30% | 488 s |
| SD_STREAM | 25 / 25 | 0.293 | 0.183 | 0.380 | 38% | 233 s |
The Tier-A ordering inverts here, narrowly. On phantom single-fibre geometry the deterministic tracker wins decisively (bending F1 0.95 vs 0.85). On real data it does not: iFOD1 takes a slim lead on mean F1 (+0.011, better on 15 of 25 bundles), with noticeably lower overreach and a quarter fewer unassignable streamlines. The margin is small enough that the right reading is "probabilistic sampling stops being a liability once crossings are pervasive", not "iFOD1 is the better tracker".
Per-bundle results are the more useful output. Best recovered are the bundles with long unambiguous cores — UF_left (F1 0.565), ILF_left (0.533), CC (0.527), Cingulum_left (0.525). Worst are exactly the ones that must cross a dense region to reach their target: CA (0.191), POPT_right (0.194), CST_right (0.195), CST_left (0.207). The CST failure has the classic signature — overlap 0.11 with overreach only 0.04, meaning the tracker is not producing spurious streamlines, it is failing to reach the lateral projections at all.
Two honest limits. Bundle assignment here (metrics.assign_to_bundles) is a 20-line
auditable segmentation, not a reimplementation of the challenge's own scoring package, so
these figures are not comparable to published challenge leaderboards. And this
pure-Python iFOD1 at lmax=6 produces shorter streamlines than MRtrix3's default iFOD2
pipeline would (median 19 mm; terminations split roughly evenly between leaving the mask
and exhausting the 1000-trial budget where the peak FOD sits near the 0.10 cutoff), so
treat bundle recovery as a lower bound.
data/registry/datasets.yaml — 20 entries (3 Tier A, 11 Tier B, 6 Tier C) with DOIs,
licences, access paths and verification status. Downloaded and checksum-verified here:
the ISMRM 2015 Tractography Challenge (Zenodo 572345, CC-BY-4.0) — DWI at
90×108×90×33, 2 mm isotropic, 1 b0 + 32 directions at b=1000, with 25 ground-truth
bundles totalling 200,433 streamlines in MRtrix .tck format. All three archive md5s
match Zenodo's published values.
Tier C access was surveyed rather than downloaded: the Allen Mouse Connectivity API was
queried live (2,331 anterograde tracer experiments across 216 injection structures; Allen
publishes no mouse dMRI volume, a negative result worth knowing before planning around
it), and IronTract, BigMac, the marmoset atlases and PRIME-DE are recorded with DOI-backed
access paths. Ten registry entries are documented_only — reachable in principle but not
fetched here; docs/data_provenance.md records exactly which and why.
The lab asked whether quantum computing — QML, QAOA, annealing — could give either a
"superfast" tractography algorithm or one that "considers all possibilities in parallel."
tractobench.quantum exists to answer that with measurement instead of opinion. The full
argument is in docs/quantum_tractography_evaluation.md.
Headline findings, all reproducible from this repository:
- No prior work exists at this intersection. arXiv, OpenAlex and PubMed all return zero on-topic hits, against control queries returning 179 and 2,907. Quantum methods have reached brain connectomes, but only for community detection on an already built connectome.
- Streamline propagation is not a search problem. iFOD1 is a sequential Markov process, embarrassingly parallel across seeds and linear in streamline length. There is no exponential configuration space for a quantum optimiser to attack.
- On the parts that are combinatorial, QAOA loses. On 21 QUBO instances built from real tracked streamlines (n = 8–20, matched 100-shot/100-read budget), simulated annealing found the exact optimum in 21/21 cases in 16.8 ms; QAOA at p=3 found it in 4/21 and took 14.2 s. Greedy single-flip descent (0.1 ms) beat QAOA at p=1 and p=2. QAOA's gap to the optimum grows with problem size. This is exact noiseless statevector simulation with unlimited connectivity — strictly better than any real device.
- Binarising the objective costs more than any solver can recover. The continuous non-negative solution — what SIFT2 and COMMIT actually use — beats the exact binary optimum in 21/21 instances by a median factor of 78×. Making the problem a QUBO discards 98.7% of the achievable fit before a solver runs.
- The hardware gap is about seven orders of magnitude. Whole-brain global tractography involves ~10⁶ variables with dense couplings; quadratic clique embedding puts that at ~8×10¹⁰ physical qubits, against ~5,000 in the largest demonstrated programmable spin glass.
The size trend was then tested past the CPU wall on an NVIDIA GB10 (n up to 24,
benchmarks/gpu_qaoa_scaling.py), because a trend measured over 12 variables is worth
what its range is worth. QAOA's excess over the optimum keeps climbing — Spearman ρ = +0.93
to +0.98 against n, at every depth — while a greedy classical control stays flat
(ρ = +0.08), which is what rules out "the instances just got harder". At n ≥ 16,
simulated annealing solved 14/14 instances exactly and QAOA 0/14 at every depth.
A 2024–2026 follow-up survey screened 22 quantum-information-science methods from 437
retrieved records, selecting by which of the two walls a method could evade rather than by
algorithm novelty (docs/quantum_frontier_2026.md). Two candidates survived, five are
viable only as quantum-inspired classical methods, fifteen were screened out, and no
evasion route outside the three anticipated ones appeared in any record. The survey's own
conclusion is that the field is actively dequantizing its 2023 headlines — quantum-enhanced
MCMC now carries a rigorous no-go bound and a classical tensor-network surrogate that
reproduces its claimed advantage. Three research directions with pre-stated falsification
conditions are in docs/research_directions.md.
One methodological result came out of that work and is worth flagging on its own. On penalty-augmented QUBOs, which relaxation seeds a warm-started QAOA matters more than circuit depth: seeding from the full penalised QUBO halves the gap to the optimum at n=20 versus seeding from the objective term alone (28.6% vs 55.7%), never loses to standard QAOA at p=1 across 12 instances (9W/3T/0L), and cuts the size-trend slope from 5.03 to 2.39 %/variable. Seeding from the objective-only relaxation — the natural choice, and not a relaxation of the problem actually being solved — buys nothing: 6W/2T/4L at p=1 and net negative at p=2 and p=3.
It still does not make QAOA competitive: full-QUBO warm-started QAOA loses to simulated annealing in every one of the 9 instances where they differ.
What survives is worth doing: publishing the negative result with this testbed, contributing the QUBO instances to the quantum-optimisation community as a new structured family with known optima, and pursuing quantum-inspired classical methods (tensor networks) under that honest label. §8 of the evaluation lists, in advance, the four conditions that would change the assessment.
src/tractobench/
sh.py SH basis in the MRtrix3 convention + precomputed Legendre table
signal.py multi-tensor forward model, Rician noise
fod.py CSD via DIPY, with an explicit basis conversion and scale normalisation
calibrator.py port of MRtrix3's rejection-sampling calibrator
tracker.py DirectionGetter interface, FODField, propagation driver
algorithms/ iFOD1 (probabilistic), SD_STREAM (deterministic)
metrics.py Tractometer-style VC/IC/NC, overlap/overreach, tracer AUC
io.py .tck and NIfTI via nibabel
quantum/ QUBO formulations, solver suite, resource model
benchmarks/ reproducible scripts for every number in this README
docs/ MRtrix3 technical reference, literature review, evaluation, provenance
data/registry/ the 3-tier dataset registry
ci/ GitHub Actions workflow (see ci/README.md to activate it)
CI lives in ci/ rather than .github/workflows/ because the token used to publish this
repository lacked the workflow scope; ci/README.md gives the one-line git mv that
turns it on.
Two are worth stating because they cost days when missed:
SH basis. MRtrix3 uses even-order real SH with n = (l+1)(l+2)/2 and
index(l,m) = l(l+1)/2 + m, matching DIPY's real_sh_tournier(..., legacy=False). The
legacy=True variant is the MRtrix 0.2 convention and is not interchangeable. This
repository fixes legacy=False everywhere and verifies the basis against the analytic
zonal reproducing kernel to ~1e-15.
FOD amplitude scale. MRtrix3's default cutoff of 0.10 is absolute, and meaningful
only against the amplitude scale its own CSD produces on real data. A phantom's scale
depends on simulated volume fractions, so fod.normalise_fod_scale rescales to a fixed
peak and reports the factor. Copying 0.10 across without this is not conservative — it is
arbitrary.
Apache-2.0. MRtrix3 is MPL-2.0 and is not vendored here; the port cites it by
file:line at commit 6735cfa. Per MRtrix3's own reference block
(cmd/tckgen.cpp:148-155): iFOD1 and SD_STREAM are attributed to Tournier, Calamante &
Connelly (2012, IJIST), and iFOD2 to the second-order-integration paper (2010 ISMRM
abstract). The ISMRM 2015 challenge data is CC-BY-4.0; see docs/data_provenance.md for
per-dataset licences.