1- # tfce (Python)
1+ # tfce
22
3- Exact threshold-free cluster enhancement, and the permutation machinery around it,
4- as a plain Python package.
3+ ** Exact threshold-free cluster enhancement, and the permutation inference around it.**
54
6- This is ** not a reimplementation** . The transform is the same C max-tree the
7- MATLAB toolbox runs — compiled straight from the repository root into a Cython
8- extension — so the two give bit-identical answers, and there is only ever one
9- implementation of TFCE in this repository to get right.
5+ [ ![ PyPI] ( https://img.shields.io/pypi/v/tfce.svg )] ( https://pypi.org/project/tfce/ )
6+ [ ![ Python] ( https://img.shields.io/pypi/pyversions/tfce.svg )] ( https://pypi.org/project/tfce/ )
7+ [ ![ License] ( https://img.shields.io/pypi/l/tfce.svg )] ( https://github.com/ChristianGaser/tfce/blob/master/python/LICENSE )
108
11- ``` python
12- import numpy as np
13- import tfce
9+ TFCE combines focal effects of large height with broad effects of large extent, and needs ** no
10+ cluster-forming threshold ** — the arbitrary choice that cluster-based inference forces on you, and
11+ that the result can depend on heavily.
1412
15- # volume: (nx, ny, nz) or (nx, ny, nz, n_permutations)
16- t = tfce.tfce(stat_map, E = 0.5 , H = 2.0 )
13+ ``` bash
14+ pip install tfce
15+ ```
1716
18- # surface: faces are 1-based, as GIFTI stores them
19- adj = tfce.adjacency_from_faces(faces, n_vertices)
20- t = tfce.tfce(surf_map, adjacency = adj, E = 1.0 , H = 2.0 )
17+ Wheels for Linux, macOS and Windows. No compiler needed.
2118
22- # a block of permutations, one per thread
23- t = tfce.tfce(perms, E = 0.5 , H = 2.0 , n_jobs = - 1 )
24- ```
19+ ---
2520
26- ## Why exact matters
21+ ## Exact, not stepped
2722
28- The TFCE of an element is an integral,
23+ The TFCE of an element is an integral:
2924
30- ```
25+ ``` text
3126TFCE(v) = ∫ e_v(h)^E · h^H dh
3227```
3328
34- and implementations normally approximate it by stepping ` h ` over a grid and
35- summing. That costs a step size, a discretisation error that depends on it, and an
36- accuracy parameter the caller has to guess. This one builds the ** max-tree** (the
37- component tree) with union-find, and because the extent function ` e_v(h) ` is
38- piecewise constant it integrates each piece in closed form. There is nothing to
39- tune.
29+ over the extent ` e_v(h) ` of the cluster containing ` v ` at height ` h ` . Implementations normally
30+ approximate it by stepping ` h ` over a grid and summing. That costs a step size ` dh ` , a discretisation
31+ error that depends on it, and an accuracy parameter you have to guess.
32+
33+ This one doesn't. It builds the ** max-tree** (the component tree) with union-find, and because
34+ ` e_v(h) ` is piecewise constant, integrates each piece in closed form. The answer is the integral, not
35+ a sample of it. ** There is nothing to tune.**
4036
41- The difference is not academic. Against nilearn's stepped transform on a 60×72×60
42- volume:
37+ The difference is not academic. Against a stepped implementation on a 60×72×60 volume:
4338
44- | n_steps | max error vs exact |
39+ | steps | error vs exact |
4540| --- | --- |
4641| 50 | 3.4% |
47- | 100 * (nilearn's default) * | ** 1.7%** |
42+ | 100 | ** 1.7%** |
4843| 200 | 0.8% |
4944| 400 | 0.4% |
5045
51- The error halves every time the steps double — first order, exactly as a step-size
52- approximation must. The exact transform has no such term, and is also faster:
46+ The error halves every time the steps double — first order, exactly as a step-size approximation must
47+ behave. The exact transform has no such term.
48+
49+ ## Quick start
5350
54- | | one volume | 16 permutations |
51+ ``` python
52+ import numpy as np
53+ import tfce
54+
55+ # a volume: (nx, ny, nz)
56+ t = tfce.tfce(stat_map, E = 0.5 , H = 2.0 )
57+
58+ # a surface: faces are 1-based, as GIFTI stores them
59+ adj = tfce.adjacency_from_faces(faces, n_vertices)
60+ t = tfce.tfce(surf_map, adjacency = adj, E = 1.0 , H = 2.0 )
61+
62+ # a block of permutations, one per thread
63+ # (nx, ny, nz, n_perm) -> (nx, ny, nz, n_perm)
64+ t = tfce.tfce(perms, E = 0.5 , H = 2.0 , n_jobs = - 1 )
65+ ```
66+
67+ Volumes take ` connectivity=6 | 18 | 26 ` (26 is the default, and what fslmaths and the MATLAB toolbox
68+ use). Surfaces take the mesh, so the neighbourhood is whatever the mesh says it is.
69+
70+ Everything is arrays in, arrays out. No image objects, no file I/O, no design parsing — those belong
71+ in a layer above, so the core can be dropped anywhere.
72+
73+ ## Speed
74+
75+ Permutations are independent, so they are TFCE'd a block at a time, one per thread, with the GIL
76+ released. On the same 60×72×60 volume:
77+
78+ | | one map | 16 permutations |
5579| --- | --- | --- |
56- | nilearn (100 steps) | 0.38 s | 9.1 s ( 569 ms/perm) |
57- | exact max-tree | ** 0.03 s** (14×) | ** 0.13 s** ( 8 ms/perm, ** 72×** ) |
80+ | stepped (100 steps) | 0.38 s | 9.1 s — 569 ms/perm |
81+ | ** tfce ** | ** 0.03 s** (14×) | ** 0.13 s** — 8 ms/perm ( ** 72×** ) |
5882
59- ## Using it inside nilearn
83+ ## Using it with nilearn
6084
61- nilearn already has ` permuted_ols(..., tfce=True) ` . Its TFCE is
62- ` nilearn.mass_univariate._utils.calculate_tfce ` , and ` tfce.nilearn_compat `
63- provides a drop-in with the same signature:
85+ nilearn has ` permuted_ols(..., tfce=True) ` , and its TFCE is a stepped approximation. ` tfce ` ships a
86+ drop-in with the same signature:
6487
6588``` python
6689from nilearn.mass_univariate import _utils
@@ -69,49 +92,73 @@ import tfce.nilearn_compat as tc
6992_utils.calculate_tfce = tc.calculate_tfce # now exact, and much faster
7093```
7194
72- One thing to know: nilearn builds its neighbourhood with
73- ` generate_binary_structure(3, 1) ` , which is ** 6-connectivity** . The MATLAB toolbox
74- and fslmaths use ** 26** . The drop-in reads the neighbourhood out of the
75- ` bin_struct ` it is handed, so it reproduces whichever one the caller meant, and
76- ` tfce.tfce(..., connectivity=6|18|26) ` lets you say so directly.
95+ It reads the neighbourhood out of the ` bin_struct ` it is handed, so it reproduces whichever
96+ connectivity the caller meant.
97+
98+ ## Fewer permutations
99+
100+ Counting exceedances cannot report a p-value below ` 1/n_perm ` . That floor — not the statistic — is
101+ what forces a permutation test to run many thousands of permutations, and it caps FDR too, since FDR
102+ is computed from the uncorrected p-values. ` tfce.tails ` removes it:
77103
78- ## What's here
104+ - ** Gamma** fit to the null of the * maximum* , for FWE-corrected p-values.
105+ - ** Generalised Pareto** fit to each element's own tail, for the uncorrected ones — with the * shape*
106+ pooled across elements, since they all carry the same statistic under the same design and so differ
107+ in scale, not in shape.
79108
80- | module | what it is |
109+ From 1000 permutations this recovers ` p ≈ 1e-4 ` with the median right, where plain counting returns
110+ zero for 91% of elements and tells you nothing at all.
111+
112+ ## What's in the box
113+
114+ | | |
81115| --- | --- |
82- | ` tfce.core ` | the exact transform: arrays in, arrays out |
116+ | ` tfce.core ` | the exact transform — arrays in, arrays out |
83117| ` tfce.tails ` | Gamma and Generalised Pareto tail approximations |
84118| ` tfce.glm ` | the permuted GLM, which never forms the permuted data |
85- | ` tfce.nilearn_compat ` | drop-in for ` calculate_tfce ` |
119+ | ` tfce.nilearn_compat ` | drop-in for nilearn's ` calculate_tfce ` |
86120
87- The core is deliberately framework-free — numpy and scipy, nothing else. No image
88- objects, no file I/O, no design parsing. Those belong in a layer above, so that
89- the core can be vendored anywhere (nilearn included) without dragging a dependency
90- tree behind it. ` nibabel ` is an optional extra, wanted only by I/O.
121+ Dependencies are deliberately thin: ** numpy and scipy** , nothing else. ` nibabel ` is an optional extra,
122+ wanted only by I/O.
91123
92- ### Fewer permutations
124+ ## Trust
93125
94- Counting exceedances cannot report a p-value below ` 1/n_perm ` . That floor — not
95- the statistic — is what forces a permutation test to run many thousands of
96- permutations, and it caps FDR too, since FDR is computed from the uncorrected
97- p-values. ` tfce.tails ` removes it: a ** Gamma** fit to the null of the maximum for
98- the FWE-corrected p-values, and a ** Generalised Pareto** fit to each element's own
99- tail for the uncorrected ones, with the * shape* pooled across elements (they all
100- carry the same statistic under the same design, so they differ in scale, not in
101- shape).
126+ This is not a reimplementation. It is the same C core that the
127+ [ MATLAB TFCE toolbox] ( https://github.com/ChristianGaser/tfce ) has been running for years, with a
128+ Cython binding instead of a MEX one — so the two give ** bit-identical** results, and the test suite
129+ holds them to it. It ships a validation suite of its own: the exactness of the max-tree is established
130+ against an * independent* stepped implementation, which must converge onto it at first order.
102131
103- ## Install
132+ ## Citing
104133
105- ``` bash
106- pip install -e . # needs a C compiler; the C core sits at ../
107- pytest # 49 checks
108- ```
134+ If you use this, please cite the method:
135+
136+ > Smith SM, Nichols TE (2009). * Threshold-free cluster enhancement: addressing problems of smoothing,
137+ > threshold dependence and localisation in cluster inference.* NeuroImage 44:83–98.
138+ > [ doi:10.1016/j.neuroimage.2008.03.061] ( https://doi.org/10.1016/j.neuroimage.2008.03.061 )
139+
140+ The tail approximations are from:
141+
142+ > Winkler AM, Ridgway GR, Douaud G, Nichols TE, Smith SM (2016). * Faster permutation inference in
143+ > brain imaging.* NeuroImage 141:502–516.
144+ > [ doi:10.1016/j.neuroimage.2016.05.068] ( https://doi.org/10.1016/j.neuroimage.2016.05.068 )
109145
110146## Status
111147
112- Early. The transform, the tails and the permuted GLM are here and tested. The
113- design layer — turning a ** BIDS Stats Model** into a design matrix, contrasts and,
114- crucially, an * exchangeability* structure — is not, and is the real work: BIDS-SM
115- specifies ` X ` , ` Formula ` , ` Contrasts ` and ` GroupBy ` , but has ** no** concept of
116- exchangeability blocks or variance groups, so what may be permuted with what has
117- to be defined on top of it rather than read out of it.
148+ Early, but the parts that are here are tested and are the parts that matter. The transform, the tail
149+ approximations and the permuted GLM are done. Not yet done: the ** design layer** — turning a
150+ [ BIDS Stats Model] ( https://bids-standard.github.io/stats-models/ ) into a design matrix, contrasts and,
151+ crucially, an * exchangeability* structure. That last one is the real work, because BIDS-SM specifies
152+ ` X ` , ` Formula ` , ` Contrasts ` and ` GroupBy ` but has ** no** concept of exchangeability blocks, so what may
153+ be permuted with what has to be defined on top of it rather than read out of it.
154+
155+ ## Links
156+
157+ - ** Source, issues, MATLAB toolbox:** < https://github.com/ChristianGaser/tfce >
158+ - ** Licence:** BSD-3-Clause. Permissive on purpose: ` nilearn ` is BSD-3 and
159+ ` nipreps ` is Apache-2.0, and neither can take on a GPL dependency. Nothing here
160+ is derived from SPM, so nothing here has to be GPL. (The MATLAB/SPM toolbox in
161+ the same repository * is* GPL, because it genuinely is derived from SPM — see
162+ [ LICENSE.md] ( https://github.com/ChristianGaser/tfce/blob/master/LICENSE.md ) .)
163+
164+ Developed by Christian Gaser, Structural Brain Mapping Group, Jena University Hospital.
0 commit comments