SurfClustSim - native surface cluster simulation and alpha tables
========================================================================
PURPOSE
-------
SurfClustSim estimates surface-cluster area thresholds under a Gaussian
noise null hypothesis. It is a native replacement for the simulation,
smoothing, SurfClust, and quick.alpha.vals.py stages generated by
slow_surf_clustsim.py. Surfaces, smoothing weights, and cluster topology
are loaded once; simulations run in memory and may use OpenMP.
For each simulated field the program generates noise, applies HEAT_07
surface smoothing, rescales the field to unit standard deviation, and saves
the largest cluster area at each nodewise p threshold. Quantiles of those
maxima form a C(p,alpha) table: a row is an uncorrected nodewise p threshold
and a column is a whole-surface corrected alpha level.
SYNOPSIS
--------
Surface-domain noise:
SurfClustSim -spec SPEC -surf_A SURF -on_surface \
SMOOTHING_OPTIONS -niter N -prefix PREFIX [options]
Volume-domain noise mapped between two surfaces:
SurfClustSim -spec SPEC -surf_A SURF_A -surf_B SURF_B -sv SURFVOL \
-vol_mask MASK SMOOTHING_OPTIONS -niter N -prefix PREFIX [options]
QUICK RECIPE: REPEAT A slow_surf_clustsim.py ANALYSIS
----------------------------------------------------
The closest native equivalent to an old surface-domain run with blur=10,
itersize=10, and pthr_list=(0.05 0.01 0.001) is:
mkdir -p clust.lh.results.10000
SurfClustSim \
-spec fsaverage_SUMA/fsaverage_lh.spec \
-surf_A smoothwm -on_surface \
-compat -target_fwhm 10 -itersize 10 \
-niter 10000 -pthr 0.05 0.01 0.001 \
-legacy_sided -athr 0.10 0.05 0.02 0.01 \
-prefix clust.lh.results.10000/surf.lh.clustsim -niml
The historical surf_vol input is not needed when noise is generated directly
on the surface. The output of interest is surf.lh.clustsim.legacy.1D
(and .legacy.niml with -niml).
This is statistical compatibility, not iteration-for-iteration identity.
The native program uses a different Gaussian random stream. Also, the old
quick.alpha.vals.py rounded simulated areas to integers and printed an entire
area/alpha curve; this program retains continuous areas and directly reports
the requested alpha quantiles. Results should agree within Monte Carlo and
small tabulation differences, not necessarily character for character.
SMOOTHING MODES
---------------
Both modes use the same SUMA HEAT_07 weights and iterative smoothing kernel.
They differ in how the kernel bandwidth and number of passes are selected.
-compat
Adaptive, block-level target-FWHM smoothing. This reproduces the path
used by slow_surf_clustsim.py and SurfSmooth with -target_fwhm,
-blurmaster, and -detrend_master.
For every block of B=-itersize fields, the program:
1. copies the B noise fields to make a smoothing master;
2. detrends the master across its B columns;
3. repeatedly estimates mean master FWHM and applies one HEAT_07 pass
until that mean exceeds -target_fwhm; and
4. applies the selected number of passes to the original fields.
One Niter is therefore selected for the entire block, not independently
for every field. Different blocks can select different values.
-target_fwhm is the requested FINAL smoothness, not additional blur.
If -sigma is omitted, a mesh-aware bandwidth is estimated from target
FWHM and mean edge length, as SurfSmooth did. Supplying -sigma fixes the
bandwidth of each pass, but the number of passes remains adaptive.
Do not give -Niter with -compat.
Use -compat when comparing with old slow_surf_clustsim.py results, or
whenever reaching a measured target FWHM is more important than speed.
It is slower because FWHM is repeatedly re-estimated and the master is
smoothed before the original fields are smoothed.
-fixed
Fixed-parameter smoothing. Every field receives exactly K HEAT_07
passes with bandwidth S:
-fixed -sigma S -Niter K
-target_fwhm is not used in this mode. Fixed mode does not promise that
the final measured smoothness equals a particular value unless S and K
have first been calibrated for the same surface, mask, and noise domain.
It is faster and its results are independent of -itersize.
A practical calibration workflow is to run a smaller -compat pilot using
the intended surface and mask, inspect the reported sigma and block Niter
values, choose a representative Niter (often the median), and then use
that sigma/Niter pair in a large -fixed run. Validate the achieved FWHM
before treating the fixed run as interchangeable with adaptive smoothing.
-sigma S
Bandwidth of the HEAT_07 kernel for ONE iteration, in surface coordinate
units (normally mm). It is not the final FWHM. Repeated graph smoothing
and surface geometry make the relation between sigma, Niter, and final
FWHM nontrivial. Very small sigma can require thousands of passes and
can cause numerical precision warnings on a coarse mesh.
-target_fwhm F, -blur F
Requested final FWHM for -compat (default 4). FWHM is an equivalent
Gaussian width in surface coordinate units. On a curved mesh it is an
estimate based on neighboring-node differences.
-Niter K
Nonnegative number of HEAT_07 passes for -fixed. If neither smoothing
mode is named, supplying -Niter implies -fixed; otherwise the default is
-compat. -compat and -fixed are mutually exclusive.
-max_Niter K, -max_smooth_niter K
Safety ceiling for an adaptive -compat search (default 3000). A warning
is issued if a block reaches this ceiling below the target FWHM. Raising
sigma is often more useful than merely raising this ceiling.
-acf a b c (or -ACF)
-acf_file FILE
A third smoothing mode, alternative to -compat and -fixed. Instead of
smoothing white noise to a target FWHM, generate noise whose spatial
AUTOCORRELATION matches a requested curve:
ACF(r) = a*exp(-r*r/(2*b*b)) + (1-a)*exp(-r/c)
with r a geodesic distance. -acf_file reads a, b and c from the output
of 'SurfFWHM -acf' run on your own data, which is the intended workflow.
WHY: a FWHM pins the WIDTH of the autocorrelation, not its SHAPE, and
cluster extent at high thresholds is governed by the tail. Smoothing
white noise gives a nearly Gaussian tail (a near 1); real data has a
much heavier one (a nearer 0.5). A too-light tail makes simulated
clusters too small, so the area threshold comes out too lenient. This
is the same failure that inflated cluster false positives in the volume
and the reason 3dClustSim gained -acf. Measure your data first: if it
reports a near 0.9, the ordinary -fixed/-compat noise is already close
and this mode buys little.
HOW: a single heat kernel cannot produce an exponential tail at all, so
the noise is built as a weighted sum of INDEPENDENT fields smoothed by
differing numbers of passes. Independent unit-variance components make
the mixture's ACF the weighted sum of theirs, so the weights follow from
an ordinary least-squares fit to the requested curve, solved once at
startup. The program then GENERATES a few fields, measures their ACF,
and reports what it actually achieved -- read that line rather than
trusting the fit. It warns if the generated 'a' misses the request by
more than 0.10.
ACCURACY: expect the generated a within roughly 0.15 of the request, and
b and c looser -- both are weakly identified whenever their component
carries little weight, so do not read much into them on their own. The
achieved value is printed on the 'ACF verify' line every run: READ IT
rather than assuming the request was met, and treat the request as a
target the program aims at, not a guarantee. The generated shape tends
to land on the heavy-tailed side of what was asked for (a lower than
requested), which is the conservative direction for cluster inference:
a heavier tail yields larger null clusters and so a stricter threshold.
The
generated curve tends to sit slightly below the requested one at large
r, because an empirical autocorrelation is biased low at long lags while
the target is an exact analytic curve; the true ACF of the noise is
likely closer to the request than the report suggests.
COST: calibration is a fixed startup cost, dominated by measuring the
basis autocorrelations, and takes roughly half a minute on a 10k-node
surface. After that the per-simulation cost is small, so this mode
suits large -niter: 1000 simulations added under a second to a 27-second
calibration in testing.
-acf_nbasis K
Number of basis fields in the mixture (default 5, range 2-16). Raise it
if the achieved 'a' misses the request; more components can follow a
more awkward curve, at proportionally more smoothing per simulation.
-acf_radius R
Largest geodesic distance used when fitting and verifying the ACF. The
default reaches a couple of decay lengths of c, which is what lets the
fit see the tail at all, capped to keep startup affordable.
-itersize B
Number of fields per processing block (default 10). In -compat this is
also the blurmaster length, matching the old script default; changing it
can change detrending, selected Niter, and results. In -fixed it is only
a memory/parallel-work batching parameter and does not change results.
Blocks are distributed across threads, so -niter/B (the block COUNT) is
what bounds parallelism -- a larger B means fewer, bigger blocks and less
parallelism, not more. The default suits both roles; prefer to leave it.
NOISE DOMAIN AND MASKING
------------------------
-on_surface [yes|no]
Generate independent N(0,1) noise at surface nodes before smoothing.
A bare -on_surface means yes. This matches '-on_surface yes' in the old
script and does not require -sv or a surface volume.
-surf_mask DSET
Restrict generated noise, FWHM estimation, smoothing, rescaling, and
clustering to nonzero surface nodes. Use the same analysis mask that
will be used on real data. This option can also restrict nodes after
volume-to-surface mapping.
-vol_mask MASK
Select volume-domain simulation and use MASK both as the noise grid and
nonzero voxel mask. Volume noise is mapped between -surf_A and -surf_B
using the coordinate system supplied by -sv SURFVOL.
-map_func FUNC 3dVol2Surf mapping function (default ave).
-f_steps N Samples between surfaces (default 10).
-nsteps is an alias.
SIMULATION AND CLUSTER OPTIONS
------------------------------
-niter N
Total number of independent simulated fields (default 1000). Unlike the
generated legacy script, this is not converted to a number of blocks and
is not rounded upward. For a final alpha=0.01 table, 10000 or more runs
are advisable: 10000 simulations provide only about 100 observations in
a 1% tail. More simulations reduce Monte Carlo variation.
-pthr p1 ... pn, -pval p1 ... pn
Uncorrected nodewise probability thresholds. Values must be in (0,.5].
Defaults are .05 .02 .01 .005 .002 .001 .0005 .0002 .0001. Each value
is converted internally to a standard-normal threshold according to the
selected sidedness.
-athr a1 ... an
Corrected whole-surface alpha levels to report. Defaults are
.10 .05 .02 .01. A table entry C(p,alpha) is the cluster area threshold
obtained from the empirical distribution of maximum cluster area.
-LOTS
Use dense 3dClustSim-style pthr and athr grids (29 p thresholds and
10 alpha levels). Explicit -pthr or -athr options appearing later can
replace the corresponding list.
-rmm R
Cluster connection rule (default -1). A negative integer -N connects
nodes within N mesh-edge steps; -1 means immediate neighbors. A positive
R uses SUMA's surface-offset radius in surface coordinate units. Use the
same rule when clustering the real statistical map.
-seed S Nonnegative deterministic seed (default 123456789).
-nthreads N OpenMP thread limit; 0 means the runtime maximum.
Results are invariant to thread count.
-quiet Suppress routine information and progress messages.
Warnings and errors are still printed.
-verbose, -verb Show per-block smoothing messages and progress
pacifiers. By default, routine setup/output messages
and one global progress update every 10 blocks are
printed.
SIDEDNESS: WHAT pthr MEANS
-------------------------
By default the program writes 1-sided, 2-sided, and bi-sided tables. In
-compat it additionally writes a legacy table. If any selector below is
given, only explicitly selected modes are written; selectors may be repeated.
-1sided, -onesided
Positive clusters only. Threshold z=Phi^-1(1-pthr), so pthr is the
positive upper-tail probability.
-2sided, -twosided
Threshold |z| at Phi^-1(1-pthr/2), so total nodewise probability is pthr.
Positive and negative active nodes belong to one binary excursion set and
may join the same connected component. This follows 3dClustSim's
2-sided convention.
-bisided
Use the same pthr/2 threshold in each tail, but cluster positive and
negative nodes separately. The simulation records the larger maximum
from the two signs.
-legacy_sided
Match slow_surf_clustsim.py/SurfClust thresholding: threshold |z| at
Phi^-1(1-pthr). Thus EACH tail has probability pthr and the total
nodewise false-positive probability is approximately 2*pthr. Opposite
signs can join. This is intentionally not the modern meaning of a
2-sided pthr; use it only for historical compatibility.
OUTPUT
------
-prefix P
Write P.1sided.1D, P.2sided.1D, P.bisided.1D, and/or P.legacy.1D.
The parent directory must already exist. Each file records the command,
surface/mask sizes, smoothing parameters, seed, threads, pthr rows, and
athr columns.
-niml Also write matching .niml table files.
-niml_only Write only NIML tables.
-both Write both 1D and NIML tables (same as -niml).
-maxarea_1D
Also write the raw per-simulation maximum cluster areas, one file per
p-value, named PREFIX.MODE.max.area.PTHR. Each file holds -niter lines,
one area per line, unsorted -- the same layout slow_surf_clustsim.py
wrote as z.max.area.PTHR. Use it to compare against the old pipeline:
quick.alpha.vals.py -niter N PREFIX.legacy.max.area.0.001
The built-in alpha table is NOT directly comparable to that tool, which
rounds areas to whole units and builds an empirical survival curve; the
table here uses an interpolated quantile instead.
One format difference from slow_surf_clustsim.py: when a simulation has
no suprathreshold cluster, that file wrote nothing, so it held FEWER
than niter lines. This writes an explicit 0, always niter lines. Alpha
values still agree because quick.alpha.vals.py normalizes by its -niter
argument -- but pass -niter, or it will normalize by the line count and
the two pipelines will disagree by the fraction of empty simulations.
-multithresh (or -mthresh)
Also write a jointly calibrated multi-threshold table, named
PREFIX.MODE.mthresh.1D.
THE IDEA, AND WHERE IT COMES FROM
This option is inspired by ETAC (Equitable Thresholding And Clustering),
available in the volume via 3dttest++ -ETAC and 3dXClustSim. ETAC's
insight is that picking a single per-node p-threshold before you know
what the signal looks like is arbitrary, and the choice changes the
answer: a low threshold favours large diffuse clusters, a high one
favours small focal ones. ETAC's answer is to stop choosing. It tests
several thresholds at once and calibrates them together, 'equitably' --
meaning each member of the family is tuned to contribute about equally
to the overall false positive rate, so no single one dominates.
That equitable multi-threshold principle is what this option borrows.
A result counts as significant if a cluster survives at ANY of the
p-thresholds, and all the cutoffs are chosen together so the probability
of that happening under the null is alpha. Concretely, one shared tail
probability q sets every threshold's cutoff, which is what makes the
contributions equitable, and q is solved for by bisection. The
arbitrary choice of threshold goes away; the cost is larger cutoffs,
typically 15-25% above the per-threshold table, which is the honest
price of testing a whole family.
HOW TO READ THE TABLE
Each COLUMN is one calibrated family, not a menu of alternatives. To
use the alpha column, test EVERY p-threshold listed in it and accept a
cluster surviving at any of them; the family then carries the stated
false positive rate. Using a single row on its own is still valid, just
conservative -- that cutoff is larger than the per-threshold table's.
The achieved null FPR is written into the header so you can see how well
the calibration converged. It needs a few thousand simulations to mean
much, and a warning is issued if it lands far from what was requested.
WHAT THIS DOES *NOT* INHERIT FROM ETAC
Real ETAC is substantially more than multi-threshold calibration, and
this option should not be mistaken for it. ETAC additionally offers:
* a null built by randomizing/permuting the residuals of your ACTUAL
data, rather than from simulated noise. This is the big one: it is
where most of ETAC's advantage over a simulated null comes from,
because it assumes nothing about the shape of the spatial
autocorrelation. This program calibrates over the noise it
SIMULATES, and so inherits every assumption in that noise model --
in particular that smoothing white noise to a target FWHM yields a
realistic autocorrelation, which is exactly the assumption that
could inflate cluster false positives in the volume literature.
* a more general cluster figure of merit than area alone; ETAC can
also use sums of powers of |z| within a cluster. Here it is area.
* simultaneous calibration across multiple blur levels. Here there is
one smoothing level per run.
* spatially varying thresholds, so that ETAC produces a MAP of cutoffs
adapted to local smoothness. This produces one number per
(p, alpha) pair, applied everywhere, and so still assumes the
surface is statistically stationary.
In short: this borrows ETAC's equitable multi-threshold idea, and with
it the benefit of not having to pick a threshold. It does not deliver
ETAC's data-driven null, its generalized merit, its multi-blur
calibration, or its spatial adaptivity. If you can run a permutation
test on your data, that remains the better tool.
-selfcheck
Recompute every maximum cluster area with an independent breadth-first
search and abort if it disagrees with the threshold-sweep result. The
sweep activates nodes in descending order and merges components with
union-find to get every threshold from one pass; the BFS simply floods
each thresholded field separately. They must agree exactly. This is a
correctness check, not a sampling option -- it is roughly an order of
magnitude slower, so use it on a small -niter, not a production run.
Table entries are areas in the SQUARED coordinate units of the surface.
They are mm^2 only when the surface coordinates are expressed in mm. Apply
a table only to the same surface geometry, mask, smoothing procedure, noise
domain, rmm, and sidedness used to generate it. A cluster is significant at
the requested (pthr,alpha) when its area meets the corresponding threshold.
The empirical quantiles retain continuous node areas and are monotonized so
that stricter pthr or alpha settings cannot produce a smaller threshold merely
from Monte Carlo jitter.
REPRODUCIBILITY AND PERFORMANCE
-------------------------------
The random stream is indexed by global simulation number, so results do not
depend on -nthreads. Fixed-mode results also do not depend on -itersize.
Compatibility-mode results can depend on -itersize because the block is the
adaptive smoothing master.
OpenMP distributes whole blocks across threads, so parallelism is limited by
the BLOCK COUNT (-niter / -itersize), not by -itersize itself. Raising
-itersize therefore REDUCES available parallelism, and in -compat it also
changes the results; leave it alone unless you mean to. Memory scales as
-nthreads * -itersize * nodes, since each thread holds a whole block.
Volume-mode noise generation is serialized (the vol2surf routines beneath it
are not reentrant); everything after it runs in parallel. Final table output
is serial. With -verbose, pacifiers are capped at 50 marks per phase.
EXAMPLES
--------
1. Historical surface-domain compatibility, one old-style table:
SurfClustSim -spec std.lh.spec -surf_A smoothwm -on_surface \
-compat -target_fwhm 10 -itersize 10 -niter 10000 \
-pthr .05 .01 .001 -legacy_sided -prefix lh.old -niml
2. New analysis with adaptive FWHM and modern sidedness tables:
SurfClustSim -spec std.lh.spec -surf_A smoothwm -on_surface \
-compat -target_fwhm 10 -niter 10000 \
-1sided -2sided -bisided -prefix lh.adaptive -niml
3. Calibrated fixed smoothing:
SurfClustSim -spec std.lh.spec -surf_A smoothwm -on_surface \
-fixed -sigma S -Niter K -niter 10000 -itersize 32 \
-nthreads 16 -prefix lh.fixed -niml
4. Volume noise mapped between white and pial surfaces:
SurfClustSim -spec std.lh.spec -surf_A smoothwm -surf_B pial \
-sv SurfVol.nii -vol_mask epi.mask.nii \
-compat -target_fwhm 10 -niter 10000 -prefix lh.volume
See also: SurfSmooth, SurfClust, 3dClustSim, slow_surf_clustsim.py
Author: P. Molfese
Date: 20 Aug 2026
Version: 0.3.0
=========================================================================
* This binary version of SurfClustSim is compiled using OpenMP, a semi-
automatic parallelizer software toolkit, which splits the work across
multiple CPUs/cores on the same shared memory computer.
* OpenMP is NOT like MPI -- it does not work with CPUs connected only
by a network (e.g., OpenMP doesn't work across cluster nodes).
* For some implementation and compilation details, please see
https://afni.nimh.nih.gov/pub/dist/doc/misc/OpenMP.html
* The number of CPU threads used will default to the maximum number on
your system. You can control this value by setting environment variable
OMP_NUM_THREADS to some smaller value (including 1).
* Un-setting OMP_NUM_THREADS resets OpenMP back to its default state of
using all CPUs available.
++ However, on some systems, it seems to be necessary to set variable
OMP_NUM_THREADS explicitly, or you only get one CPU.
++ On other systems with many CPUS, you probably want to limit the CPU
count, since using more than (say) 16 threads is probably useless.
* You must set OMP_NUM_THREADS in the shell BEFORE running the program,
since OpenMP queries this variable BEFORE the program actually starts.
++ You can't usefully set this variable in your ~/.afnirc file or on the
command line with the '-D' option.
* How many threads are useful? That varies with the program, and how well
it was coded. You'll have to experiment on your own systems!
* The number of CPUs on this particular computer system is ...... 1.
* The maximum number of CPUs that will be used is now set to .... 1.
* SurfClustSim also accepts -nthreads N to set a smaller thread limit
for one invocation. Parallel simulation workers are also limited by
-itersize (default 10).
