arXiv:1704.08562v1 [math.ST] 27 Apr 2017

ESTIMATING THRESHOLDING LEVELS FOR RANDOM
FIELDS VIA EULER CHARACTERISTICS
By Robert J. Adler∗,‡ Kevin Bartz∗
and Anthea Monod†

S. C. Kou∗,§

Technion, Renaissance Technologies, Harvard and Columbia.
We introduce Lipschitz-Killing curvature (LKC) regression, a new
method to produce (1 − α) thresholds for signal detection in random fields that does not require knowledge of the spatial correlation
structure. The idea is to fit observed empirical Euler characteristics to the Gaussian kinematic formula via generalized least squares,
which quickly and easily provides statistical estimates of the LKCs —
complex topological quantities that can be extremely challenging to
compute, both theoretically and numerically. With these estimates,
we can then make use of a powerful parametric approximation via
Euler characteristics for Gaussian random fields to generate accurate
(1 − α) thresholds and p-values. The main features of our proposed
LKC regression method are easy implementation, conceptual simplicity, and facilitated diagnostics, which we demonstrate in a variety of
simulations and applications.

1. Introduction.
1.1. Pre-History. The paper that you are (maybe) about to read was
started in 2008, when RJA spend a semester at Harvard and talked a lot
with SCK. During a lecture, RJA made the claim that estimating LipschitzKilling curvatures (defined below) was a hard problem, at which stage SCK
commented “It doesn’t seem all that hard. Why can’t you just treat it as
a regression problem?” At the time KB was a graduate student, listening
to the same lecture, and all three got together to write the first version
of this paper. Since RJA returned to the Technion, and KB graduated and
moved to Renaissance Technologies, it took until 2011 for a final first version
to be submitted for publication. A few months later it was rejected (from
∗

Research supported in part by US-Israel Binational Science Foundation, 2008262.
Research supported in part by FP7-ICT-318493-STREP.
‡
Research supported in part by ERC 2012 Advanced Grant 20120216.
§
Research supported in part by NIH/NIGMS,R01GM090202 and NSF, DMS-1510446.
AMS 2000 subject classifications: Primary 60G60, 62G32, 62E20; Secondary 60G15,
62M30, 62M40.
Keywords and phrases: Brain imaging, Gaussian kinematic formula, Lipschitz-Killing
curvature, design points, excursion set, regression, significance level, covariance.
†

1

2

ADLER, BARTZ, KOU AND MONOD

an excellent journal) but, since the authors still liked their paper they set
about preparing a revision, this time together with AM, who was then a
postdoctoral fellow at the Technion. In 2014 the new version was submitted
to a different (and even better) journal, only to be rejected since it was
too methodological and not applied enough, but with some nice ideas for
rewriting as long as the authors were prepared to rewrite everything in such
a way that it made the referees happy and them miserable.
Now that AM has left the Technion and is doing other things, RJA is
(happily) approaching retirement, KB is well entrenched in industry, and
AM is also busy with new problems, it seems unlikely that this paper will
ever get revised again. But, the four authors still like it, and it keeps getting
cited (seemingly implying that others share their opinion) and so the time
has come to put what there is up on arXiv. Which is why what you are
reading is a 2017 arXiv deposit of what is basically a 2011 paper.
1.2. Motivation. Random field models are widely used in many scientific
applications, including the description of spatial structures in environmental and epidemiological studies, the statistical analysis of brain images, and
the modeling of the cosmic microwave background radiation, as well as other
cosmological phenomena. An important problem of interest common to most
of these applications is the determination of threshold levels for the random
field, which indicate that regions with values above the level are, in some
sense, significant, while regions with values below are not. Accurate determination of threshold levels for random fields faces a major challenge that the
simpler setting of independent observations does not: the values of random
fields are correlated in space.
A concrete example that we will study in this paper is an experiment
involving functional magnetic resonance imaging (fMRI), which was a language priming experiment carried out by Dehaene-Lambertz et al. (2006)
that appeared in the functional image analysis contest (FIAC). fMRI responses were measured twice for 16 subjects after each heard a sentence
spoken under two different conditions: once with the same speaker both
times, and once with different speakers. After each repetition, hemodynamic
activity was measured at every point (voxel) in a 64 × 64 × 30 grid that
encompasses the brain for each subject. The fMRI scan for the first subject
is shown in Figure 1.1: the light and dark red domains represent mid- and
high-activation regions of the brain, respectively, under the same-speaker
(upper left) and different-speaker conditions (lower left).
This study inspires two important research questions: Firstly, are there
any significant differences at all? Secondly, if there are, where do these sig-

THRESHOLDS VIA EULER CHARACTERISTICS
Subject 1
Same-speaker

Subject 1
Difference

Subject 1
Studentized Residual

Subject 1
Different-speaker

All Subjects
Average Difference

All Subjects
t Statistic

3

Fig 1.1. Example of fMRI brain scans from the functional image analysis contest (FIAC).
Responses were measured for 16 subjects after hearing the same sentence spoken twice
under two conditions: both spoken by the same speaker (top left) and each by a different
speaker (lower left). The gray shell gives the outline of the brain, while the light and dark red
domains represent the contours at the hemodynamic activity levels (12,000, 15,000). The
difference brain (upper middle) shows the pointwise difference of these two scans. In red
regions, activity under the different-speaker condition exceeds that under the same-speaker
condition. In blue regions, the converse is true. The pointwise average of these differences
across the 16 experimental subjects is less noisy (lower middle). The residual brain (upper
right) gives the difference between the first subject and the experimental average. Finally,
the t brain (lower right) shows paired t statistics, where the light and dark colors represent
the nominal single t-test 95% and 99% thresholds of 2.13 and 2.95.

nificant differences lie? A natural way to look for differences is to compare
the two conditions at every voxel in the brain. Figure 1.1 (lower right) shows
the result of voxelwise paired t tests: in red regions, different-speaker activity exceeds same-speaker activity; in blue regions, the opposite is true. The
light and dark colored regions show areas where the t statistics exceed the
nominal single-test 95% and 99% levels of 2.13 and 2.95, respectively. The
multiplicity of comparisons of 24,759 voxels in the brain, each with its own t
test, poses a challenge in establishing significance. Using the 95% threshold,

4

ADLER, BARTZ, KOU AND MONOD

1,273 voxels indicate significant differences, which, as a proportion of 24,759,
amounts to approximately 5%, and therefore mainly false positives. In contrast to these excessive false positives using the t test 95% threshold, using
the 95% Bonferroni bound of 7.23 indicates there are no significant differences at all. This bound is known to be very conservative, and the resulting
conclusion of no significant results whatsoever is unconvincing, and moreover, unlikely to be inconsistent with neurological hypotheses. Furthermore,
in considering the nature of fMRI studies, the strength of spatial dependence depends on the resolution, among other factors, and should thus also
be taken into consideration. These preliminary exploratory results underline the importance of our motivating research questions and provide the
impetus for a sound method to determine and characterize significance.
1.3. Tail probabilities and Euler characteristics. As a starting point, let
the null hypothesis H0 assert that both brain scans, under both conditions,
are, on average, equivalent. For the moment, assume that rather than having
voxel-based data, we observe the fMRI images over a continuum. Issues of
resolution and discretization will be addressed and formalized later. Under
such an H0 , the grid of t-statistics is a smooth random field T over a region
S of the brain. Since high values of T usually indicate deviation from H0 , a
natural test statistic to consider is then the maximum
∆

MS = sup T (s).
s∈S

The problem with using MS as a test statistic, however, is that its null
distribution is required, but is virtually never known. In the FIAC data, MS
= 6.08; a test of H0 requires the p-value, P (MS ≥ 6.08|H0 ). In addition,
to identify the activated regions of T , we require a 5% threshold t such
that P (MS > t|H0 ) = 0.05. Obtaining these unknown quantities, the null
distribution and the threshold, is difficult because they both depend on the
correlation structure of T , which itself is also unknown.
To bypass this problem, we make use of a powerful parametric method
for determining the null tail probabilities of MS for a wide class of Gaussian and Gaussian-related random fields: the Euler characteristic heuristic
(ECH) (Adler, 2000), which provides an accurate approximation to the exceedence probability P (MS ≥ u) for large u (i.e. high levels or thresholds
u), where “large” refers to values of P (MS ≥ u) of the order of 5% or
smaller, and “accurate” refers to an error in approximation of the order
of 1% to 2% in the threshold level. (i.e. A true threshold of, say, 5, in
standardised units, will be estimated by a value in the range (4.99,5.01).)

THRESHOLDS VIA EULER CHARACTERISTICS

5

The ECH has been widely used in topological inference in fMRI studies (e.g., Cao and Worsley, 1999; Friston et al., 1994; Kilner and Friston,
2010; Shafie et al., 2003; Worsley, 1994, 1995a,b), cosmological data (e.g.,
Bardeen et al., 1986; Gott et al., 2008, 2007; Torres, 1994; Vogeley et al.,
1994), as well as a number of other areas. Moreover, for smooth Gaussian
random fields, the ECH is not merely a heuristic: it is a rigorous theory
(Taylor, Takemura and Adler, 2005; Adler and Taylor, 2007).
The ECH is based on the Euler characteristic (EC) ϕ, an important
topological quantity for many general classes of well-behaved sets. For a
3-dimensional Euclidean volume V , ϕ counts the number of each of the
three types of topological features of a manifold — (i) solid, simply connected regions of the manifold, or connected components, (ii) visible, open
holes or handles, and (iii) invisible, closed holes or voids — in an alternating sum. It is a topological invariant, i.e. a property of a topological space
that is invariant under homeomorphisms, or continuous deformations, of the
topological space. It is given by
(1.1)

ϕ(V ) = # connected components in V
−(# handles in V ) + (# voids in V ).

The ECH considers the ECs of a specific case, that of the excursion sets Au ,
(1.2)

∆

Au = {s ∈ S : T (s) ≥ u};

Figure 1.2 shows examples of Au and ϕ(Au ) for varying u; ϕ(Au ) is an
integer that can be negative (middle left), positive (middle right) or close to
zero (rightmost). Since topological features of the manifold lying below the
level u are effectively ignored, the manifold (and thus also the EC) becomes
simpler as u moves higher: At moderately high levels of u, Au typically takes
on the simplified form of a union of simply connected components, while at
the highest levels, it either disappears completely (i.e. it is empty) or all
that remains is a single simple component. Thus, at such high levels, the
EC is either zero or one. Based on this phenomenon, the ECH claims that
the expected EC for high u approximates the tail probability,
(1.3)

P (MS ≥ u) ≈ E[ϕ(Au )].

A further result based on this phenomenon of topological simplification of
excursion sets with increasing u is that the expected EC is guaranteed to lie
in [0, 1] for u high enough. From (1.3), we obtain the required p-value for H0 ,
which can then be inverted to find the uα , or the threshold, corresponding

6

ADLER, BARTZ, KOU AND MONOD

to a (1 − α) confidence level, via the approximation
(1.4)
ϕ(A−4 ) = 1

∆

u
cα = max{u : E[ϕ(Au )] ≥ α}.
ϕ(A0 ) = −493

ϕ(A2 ) = 98

ϕ(A3 ) = 3

Fig 1.2. Four excursion sets for a random field from the FIAC data. Each excursion set
∆
is defined as Au = {s ∈ S : T (s) ≥ u}. The light gray shell shows the outline of the brain,
while the dark gray regions show Au for u = -4, 0, 2, and 3.

When considering random fields (see Adler and Taylor, 2007) with constant mean and variance (which for the moment we take to be zero and one),
the ECH provides the additional practicality of a parametric closed form for
the expected EC E[ϕ(Au )], given by the Gaussian kinematic formula (GKF)
(Taylor, 2006),
dim(S)

(1.5)

E[ϕ(Au )] =

X
i=0

Li (S)ρi (u).

Here the ρi (u) are functions that usually take on a simple and explicit form.
For example, if the random field is Gaussian, then the ρi (u) are explicitly
expressible in terms of Hermite polynomials, as we shall see in (2.7) below.
Other examples can be found in Adler (1981); Adler and Taylor (2011);
Adler, Taylor and Worsley (2015?). The functions Li (S), which depend on
both the domain S and the covariance of T , are the Lipschitz-Killing curvatures (LKCs), or intrinsic volumes, of S — complex topological quantities
that are often extremely difficult to evaluate theoretically, even for dedicated
topologists, as well as numerically. When known, however, the LKCs Li (S)
provide the expected EC for all u via the GKF (1.5), which can in turn
be substituted into (1.3) provided by the ECH to obtain an approximate
p-value for H0 or into (1.4) to obtain a threshold level. The estimation of
the LKCs is the bottleneck of a successful application of the ECH and is the
central concern of this paper.

7

THRESHOLDS VIA EULER CHARACTERISTICS

1.4. Estimating Lipschitz-Killing curvatures. We introduce a new approach that estimates the LKCs by simply matching the expected and empirical ECs, which then directly yields the p-value and threshold level in a
straightforward manner via (1.3), (1.4), and (1.5). Figure 1.3 provides an
(i)
illustration of our method: the left panel shows the empirical ϕ(Au ) profiles for the 16 observed fields from the FIAC data (thin gray lines) and
their sample average (solid black line) for different values of u. We find the
best-fitting LKCs through a generalized least squares regression of (1.5).
The fitted LKCs then produce the best-fitting profile: the dashed black line.
The intersection of this fitted profile and the 0.05 line yields the 95% confidence threshold as u∗95% = 4.19 (black dot). This threshold can be applied
to the brain of t statistics to identify significantly activated regions. Doing
this indicates only 7 voxels activated beyond 4.19 — a striking contrast to
the 1,273 found using the naive t-test 95% threshold of 2.13, and also more
credible than the 0 found using the 95% Bonferroni bound of 7.23.
ϕ(Au), FIAC Residual Fields

ϕ(Au), FIAC Residual Fields, Zoomed

200
0.4

u*0.95 = 4.19

0.3

ϕ(Au)

ϕ(Au)

0

−200
Field 1
Field 2
Field 3
Field 4
Field 5
Field 6
Field 7
Field 8
Field 9
Field 10

−400

−600
−2

0

2

Level u

Field 11
Field 12
Field 13
Field 14
Field 15
Field 16
Average
Fitted
5%

4

0.2

0.1

u*0.95 = 4.19
2.5

3.0

3.5

4.0

4.5

Level u

Fig 1.3. Illustration of the proposed regression method to fitting Lipschitz-Killing curvatures (LKCs). The light grey lines show the observed Euler characteristic (EC) profiles for
each of the 16 observed fields from the FIAC data. The solid black line is their average.
The dashed black line is the expected EC estimated by our regression method. The black dot
is the point where the expected EC intersects with dotted grey line at 0.05, which denotes
the 5% threshold u∗95% . The plots are displayed both over a broad range of levels (left) and
zoomed in near the threshold (right).

One caveat of which to take note in matching expected and empirical
(i)
ECs is that the observed (empirical) ECs ϕ(Au ) at very high levels of u
are generally too noisy to directly estimate the threshold u∗95% ; calculating
a threshold from extreme empirical ECs is akin to estimating a tail quantile
with a small number of observations (16, for the FIAC data). The right panel
of Figure 1.3 illustrates this phenomenon by zooming in and showing that

8

ADLER, BARTZ, KOU AND MONOD

the empirical and expected ECs cross the 0.05 line at very different values of
u. In cases of very high u, for general p-value approximation and threshold
level determination, use of the expected EC yields significantly more reliable
and accurate results that do not depend on small sample sizes, but instead
derive from rigorous probability theory and topology.
Inference via LKC regression is parametric because it relies on the ECH,
which leverages on an assumed underlying Gaussianity. A commonly suggested nonparametric alternative is a permutation test on the group assignments (same-speaker, different-speaker) to approximate the null distribution
of MS . The main drawback of this nonparametric option is the tremendous
expense of computing the maximum t-statistic given a set of permuted labels,
which is severely impeded by the computation of many hundred thousand
voxelwise averages and standard deviations. For the FIAC data, a permutation test with just 100 permutation samples takes several hundred times
longer than LKC regression. Moreover, observed random fields are typically
Gaussian-related by construction, since they arise as the residual of some statistical procedure. Nonparametric methods generally underperform in cases
where such distributional properties are known.
To date, the only other existing estimation method for LKCs is the warping method of Taylor and Worsley (2007). Warping transforms the realizations T (1) , T (2) , . . . , T (F ) of the random field to a volume in a F -dimensional
space, where F is the number of random fields in the data; for the FIAC
data F =16. The intrinsic volumes of the transformed volume are shown to be
good estimates of the LKCs. Though accurate, understanding the method
requires a knowledge of Riemannian geometry, and can be computationally slow, particularly for high dimensions. Compared to Taylor and Worsley
(2007), our method has similar accuracy but is computationally faster. Furthermore, it requires only the ECs of the observed random fields. Since
observations are typically taken over a regular grid, efficient numerical optimization is available. The warping method, on the other hand, must perform
the slower computation of intrinsic volumes on a high-dimensional transformed space with irregularly spaced points.
1.5. Paper structure and acknowledgements. The remainder of this paper is organized as follows: In Section 2, we introduce our regression-based
estimation method, detailing on the various implementation factors and
combinations available, and their associated issues of optimality. Section
3 presents the results of a large scale simulation study, which compares
these various factors and their combinations, as well as our approach to the
warping method of Taylor and Worsley (2007). In Section 4 we extend this

THRESHOLDS VIA EULER CHARACTERISTICS

9

comparison from simulation studies to the concrete example of the FIAC
data. Section 5 concludes the paper with a summary of its main results and
proposals for other areas where our method would be applicable. The Appendix provides a brief treatment of the Gaussian kinematic formula (1.5).
Finally, we would like to thank Jonathan Taylor for helpful discussions at
various stages of this work, and Naor Alush for help with the programming.
2. Lipschitz-Killing Curvature Regression. The goal of our regression method is to obtain accurate estimates of the Lipschitz-Killing curvatures (LKCs) Li for use in p-value calculation and threshold level determination (see (1.2) and (1.4)). (While the LKCs are of interest in and of
themselves, in this paper we shall concentrate almost exclusively on their
importance for threshold determination.) The inputs are the observed realizations T (i) , i = 1, . . . , F , of a random field, assumed to be Gaussian or
Gaussian-related and normalized to mean zero and unit variance, over a
region S.
The Gaussian kinematic formula (GKF) (1.5), combined with the central limit theorem, suggests a linear model with heteroscedastic, possibly
correlated errors:
(2.6)

F
1 X
ϕ(A(i)
u ) =
F
i=1

dim(S)

X
i=0

Li ρi (u) + ε(u),

ε(u) ∼ N (0, σu2 ),

for some unknown σu2 = Var(ϕ Au ) /F . The response variable is the average
of the empirical Euler characteristics (ECs). The regressors ρi (u) take on a
simple and explicit form; when the random field is Gaussian,
(2.7)

ρi (u) = (2π)−(i+1)/2 Hi−1 (u)e−u

2 /2

,

where Hj is the j-th Hermite polynomial. We estimate the unknown LKCs
Li through generalized least squares (GLS).
Once the LKCs are estimated, we easily obtain an estimate of the expected
EC via the GKF (1.5) for any level u,
dim(S)

Ê[ϕ(Au )] =

X
i=0

L̂i (S)ρi (u),

which then substitutes into (1.3) to give an approximate p-value for hypothesis testing for large u, and into (1.4) to give a threshold level.

10

ADLER, BARTZ, KOU AND MONOD

The validity of the GKF (1.5) for all u forms the foundation of the effectiveness of our LKC regression method. For p-value and threshold calculation, one typically encounters large values of u, which correspond to
small-probability events, the direct estimation of which is unreliable. Our
regression approach, instead, is grounded on the observation that when u is
small or moderate, the expected EC E[ϕ(Au )] can be well estimated from
the data, since these cases do not correspond to small-probability events.
These reliable estimates then translate, in turn, through our regression, into
reliable estimates of the LKCs Li , which do not depend on u; the L̂i then
yield good approximations for p-values and threshold levels for large u. In
summary, LKC regression leverages the precision of estimation at low levels
of u to obtain accurate approximation at high levels of u.
Several issues need to be addressed in the construction of our LKC regres(i)
sion method: First, to evaluate the empirical EC, the region Au for fixed
u must be determined, which poses a potential difficulty because each field
T (i) is typically observed at a discrete grid of sample sites, yet, to apply
(i)
the GKF, the excursion sets Au must be described as smooth regions over
(i)
S. We detail the construction of Au in Section 2.1. Second, compared to
conventional linear regression, an interesting feature of our setting is that
arbitrary amounts of “data” can be generated at very little cost, simply by
exploiting (2.6) at a collection of specified levels u. Each u gives distinct
(i)
excursion sets Au and their corresponding average EC,
F
1 X
ϕ(A(i)
u ).
F
i=1

However, this supposed advantage in fact raises the problem of a tradeoff
between residual error in the regression and estimation error in the covariance matrix: additional data might help the regression, but the covariance
matrix becomes increasingly difficult to estimate with more data. We address this tradeoff in detail in Section 2.3. Third, after specifying the levels
u, the unknown error covariance matrix must be estimated before running
GLS. We present several estimation options in Section 2.2. In Section 3 we
find optimal settings for level selection of u and error covariance estimation.
As a tool for hypothesis testing, the LKC regression approach is powerful
because it overcomes the need for the complete specification of the covariance structure of the random field. Very few assumptions about covariance
are made; even isotropy is not required. The random fields only need be
Gaussian or Gaussian-related, which is often a natural consequence of the
data generating process.

THRESHOLDS VIA EULER CHARACTERISTICS

11

2.1. Observed Excursion Sets. Though assumed to be continuous, in
practice random fields are typically observed at a discrete set of sample
sites sk ∈ S, k = 1, . . . , K. An empirical realization of a random field T
is thus defined by the observations {T (sk ), k = 1, . . . , K}; Figure 2.4 provides an illustration. Sample sites commonly lie on a rectangular lattice for
fields on the square and cube (upper left and middle) or on a latitudinallongitudinal lattice for fields on the sphere (upper right). We measure the
resolution of the field by the grid size G, the number of sample sites in
each direction/dimension; for example, for a 3D cube G =20 corresponds to
K = G3 =8000.

Fig 2.4. Simulated examples of observed fields on the unit square (left), on the unit cube
(middle) and on the unit sphere (right). The two rows show the sample sites (upper row)
and approximated excursion sets (lower row). In each case darker points denote greater
values. The simulated random fields have Gaussian covariance and are simulated over a
discrete grid with G = 10 sample sites in each dimension for the square and cube and G =
50 for the sphere.
(i)

In theory, to compute the EC ϕ of an excursion set Au , the excursion
set needs to be represented as a smooth domain. However, as the fields are
observed only at discrete points sk , an approximation of the exact domain is
required. In the digital topology literature, there exist many constructions
(e.g., Imiya and Eckhardt, 1999) of domains created from a set of points,
known as quasi-objects. The methods are distinguished by their connectivity, which specifies how close two sample points must be for them to be
considered part of the same domain. For example, in considering points on
a planar lattice, 4-connectivity places vertically and horizontally adjacent

12

ADLER, BARTZ, KOU AND MONOD

points into the same domain, but diagonally adjacent points into separate
domains (unless they also share an immediate neighbor); 8-connectivity includes all diagonal adjacencies in the same domain.
To create excursion sets from observed random fields, we use planar (2D)
8-connectivity and spatial (3D) 26-connectivity, where each point touches
(i)
all its immediate and diagonal neighbors. This is equivalent to building Au
by introducing a small square or cube centered at sk whenever T (i) (sk ) > u,
and then joining them to form the excursion set. For spherical fields, we use
8-connectivity over the latitudinal/longitudinal grid of sample sites, with
the excursion sets made up of small plates. Figure 2.4 shows examples of
the resulting excursion sets for fields on the plane (lower left), cube (lower
middle) and sphere (lower right). We experimented with different connectivity structures, and concluded that they had little impact on the calculated
ECs. In the brain image analysis example (FIAC, 64 × 64 × 30), we observed
nearly identical EC profiles for spatial 6-connectivity, 18-connectivity, and
26-connectivity.
The final step in computing the response variable in the regression of (2.6)
is to evaluate the EC for excursion sets. (1.1) provides a simple formula,
provided we can count the number of connected components, handles, and
voids. This calculation is easy in the discrete setting, in which the excursion
set is made up of a large number of simple cells, such as squares, cubes or
plates; the formula thus reduces to an alternating sum over distinct vertices,
edges, faces and cubes of these cells:
(2.8)

ϕ = # vertices − # edges + # faces − # cubes

The same formula also holds in two dimensions without the final term. A key
advantage of our LKC regression method is the rapidity of this computation
for cells on a grid; the warping method of Taylor and Worsley (2007), in
contrast, must perform the slower computation of all intrinsic volumes (not
just the EC) over an irregular (nonrectangular), warped grid.
For an application with F realizations of the random field, LKC regression
calls for the computation of ECs at U different excursion levels. Taking
advantage of binary search and careful indexing, the computation cost of
our method is only O(F K log U ). This makes it extremely fast to evaluate
the EC for many levels U .
2.2. Covariance of the Error Terms. The covariance structure of the
error terms in LKC regression is defined by
∆

Cϕ (u, v) = Cov(ϕ(Au ), ϕ(Av )).

THRESHOLDS VIA EULER CHARACTERISTICS

13

We emphasize that Cϕ should not be mistaken for the covariance of the
random fields T (i) ; Cϕ describes the correlation of the Euler characteristics
across differing levels, and not of field values themselves. Cϕ can be viewed
(i)
as the covariance function of the stochastic process ϕ(Au ) in u.
The error covariance structure of LKC regression is nontrivial: First, they
are heteroskedastic; the observed EC tends to be less variable at higher levels
of u, as illustrated in Figure 2.5 (upper right) for the FIAC data. Second,
they are correlated; the sample correlation plot (lower left) and correlogram
(lower right) show pronounced negative correlation for levels separated by
about 1. This is due to the periodicity of the three swings in the EC profiles,
which are about one (standardised) unit apart.
Given U levels uj , LKC regression requires the covariance matrix Σ =
[Cϕ (uj1 , uj2 )]U ×U . Although a misspecified covariance matrix still gives consistent estimates under GLS, the efficiency of the estimate depends on the
estimation of Σ. There are five ways of producing a positive-definite estimates of the covariance matrix Σ̂:
1. Identity (I). As a baseline, heteroskedasticity and correlation are ignored under this option so that Σ̂ ∝ IU , yielding an ordinary least
squares (OLS) model.
2. Smoothed diagonal (SD). Heteroskedasticity is incorporated while cor(i)
d
relation is ignored. The sample variances V̂u = Var[{ϕ(A
u )}i ] are computed and then smoothed as a function of u; we employ local quadratic
smoothing with 10% nearest neighbors (Loader, 1999, 2010). To prevent the smoothing procedure from generating negative variances, the
smoothing is applied to the logarithm of the sample variances; values
of u with a sample variance of zero are dropped from the smoothing
procedure. Examples of original and smoothed variances are depicted
in Figure 2.6 (left panel). The covariance matrix is Σ̂ = Diag(Ṽ ),
where Ṽ denotes the vector of smoothed variances.
3. Smoothed correlogram (SC). This option uses the same smoothed variances of SD but assumes that the correlation between ϕ(Au ) and ϕ(Av )
is stationary, i.e. a function of the separation |u − v|. Following a procedure in Hall, Fisher and Hoffman (1994), the sample correlogram is
smoothed and then its negative Fourier frequencies are truncated to ensure positive-definiteness. Figure 2.6 shows the sample and smoothed
correlograms (right) for the FIAC data.
4. Sampson-Guttorp Warping (SGW). Warping (Sampson and Guttorp,
1992) does not assume stationarity in its smoothing of covariance.
Multi-dimensional scaling (MDS) is applied to the sample semivari-

14

ADLER, BARTZ, KOU AND MONOD
ϕ(Au), FIAC Residual Fields

Var(ϕ(Au))
10000

200

8000

Var(ϕ(Au))

ϕ(Au)

0

−200

6000

4000

−400
2000

−600

0
−3

−2

−1

0

1

2

3

−3

−2

−1

Level u

0

1

2

3

Level u

Sample Cor(ϕ(Au), ϕ(Av))

ϕ(Au) Correlogram

1.0
1.0

2

Level v

0

0.0

−1
−0.5

Cor(ϕ(Au), ϕ(Av))

0.5
1

0.5

0.0

−0.5

−2

−1.0
−2

−1

0

Level u

1

2

0

1

2

3

4

Level separation: u − v

5

Fig 2.5. Illustrations of the profiles (upper left), variance (upper right), correlation (lower
left) and correlogram (lower right) of the observed EC ϕ(Au ) for the 16 FIAC fields. In
the correlation plots, the grayness reflects the strength of the correlation between two levels
u and v. The correlogram reports average correlation as a function of |u − v|.

ogram, and then all but the first k MDS components are dropped.
The top k are re-expanded to produce a smoothed semivariogram,
which is then translated to a covariance matrix. In our experience, a
minimum of k = 5 components are necessary; we use k = 10.
5. Pseudo-inverse (PI). The sample covariance matrix cannot be directly
used as Σ̂ because it is singular when F < U ; Σ−1 can be replaced with
the pseudoinverse of Σ (Rao, 1962), obtained by dropping the eigenvectors corresponding to eigenvalue λ = 0 in the sample covariance
matrix.

THRESHOLDS VIA EULER CHARACTERISTICS

15

The first four options — I, SD, SC and SGW — represent a complexity
continuum in modeling C. Figure 2.7 illustrates SC and SGW, comparing the
results of SC (middle) and SGW (right) to the original sample covariance.
SC depends on the assumption that ϕ(Au ) is stationary, which accounts for
the diagonal stripes in the SC covariance plot. SGW is more flexible, yet we
shall see that SGW tends to overfit Σ, which tends to increase the variance
of the final regression estimates.
The fifth option, PI, unlike the others, does not estimate the function
Cϕ (u, v) for arbitrary u and v; it produces only Σ̂, and only for a prescribed
set of uj . While PI does incorporate both heteroskedasticity and correlation,
we shall see that it too suffers from overfitting and yields highly variable
regression estimates.
Our final recommendation is the SD covariance in the regression, as it provides the best balance between computation speed, flexibility and stability.
SC, SGW, and PI are flexible, however tend to overfit the sample covariance matrix, producing estimates of high variance. We conduct a detailed
comparison of the five approaches in Section 3.
Sample, Smoothed Var(ϕ(Au))

Sample, Smoothed ϕ(Au) Correlogram
Sample
Smoothed

1.0
10000

Cor(ϕ(Au), ϕ(Av))

V(ϕ(Au))

8000

6000

4000

0.5

0.0

−0.5

2000

0
−3

−2

−1

0

Level u

1

2

3

0

1

2

3

4

Level separation: u − v

5

Fig 2.6. Smoothing applied to the variogram (left) and correlogram (right) by the smoothed
diagonal (SD) and smoothed covariance (SC) methods. Each plot shows both sample
(dashed) and smoothed (solid) estimates of the variance and correlation of ϕ(Au ).

2.3. Design Selection. The levels u = {uj } are free to be chosen in LKC
regression; both the number U of levels and their locations (placement) need
to be specified. The primary tradeoff is between accuracy and speed. If the
levels are too few or too coarsely spaced, then estimation is fast but accuracy
suffers. On the other hand, if the levels are too many or too finely spaced,
computational time is wasted with little marginal gain because the resulting
ECs are highly correlated with one another.
A guiding principle is to minimize the theoretical variance of L̂i , the
regression-estimated LKCs. Referring to (2.6), let X(u) be the U × dim(S)

16

ADLER, BARTZ, KOU AND MONOD
Sample Cov(ϕ(Au), ϕ(Av))

SC Cov(ϕ(Au), ϕ(Av))

10000
2

SW Cov(ϕ(Au), ϕ(Av))

2

2

1

1

0

−1
−5000
−2

Level v

0

Level v

Level v

5000
1

0

0

−1

−1

−2

−2

−10000
−2 −1

0

1

2

−2

Level u

−1

0

1

2

Level u

−2

−1

0

1

2

Level u

Fig 2.7. Covariance matrix of ϕ(Au ) estimated in three ways: the sample covariance (left),
the SC method (middle), and the SGW method (right).

matrix of regressors, with Xij = ρi (uj ), and denote the error covariance as
Σ(u) to emphasize its dependence on u. This gives
V [L̂] = X ′ (u)Σ(u)−1 X(u).
Since V [L̂] is a matrix, we require a real-valued function g of V [L̂] to perform optimization. We consider A-optimality, where g(M ) = tr(M ), and
D-optimality, where g(M ) = det(M ) (Pukelsheim, 2006), and substitute
an estimated covariance matrix from Section 2.2 for Σ(u). Note that for
a diagonal Σ(u), the optimizing uj lie at the extremes of the permissible
range of u, which requires Σ(u) to have nonzero off-diagonal terms for the
optimization of u to be sensible.
We first examine the question of how many levels to use. Using the FIAC
data, we tried possible values of U between 5 and 200. For illustration,
we fix u to be U equally-spaced levels between the field minimum and the
maximum, though our results are not sensitive to this mode of spacing.
Figure 2.8 (left) shows the result for the SC covariance smoothing method
as U varies. Although a larger U may be better, the improvement levels off
by about U = 50. Thus, we fix U = 50 throughout our subsequent analysis.
The SGW and PI covariance options gave nearly identical results.
We next attempt to find an optimal distribution of levels. Ideally, we could
solve the optimization problem
uoptimal = arg min g(X ′ (u)Σ(u)−1 X(u)).
u

However, this search is challenging because the domain is U -dimensional.
In addition, the objective is discontinuous in u due to a variety of instabilities in the covariance smoothing function. Rather than attempting a direct
optimization, we instead compare the objective values for a few placement
heuristics:

17

THRESHOLDS VIA EULER CHARACTERISTICS

1. Equal spacing of u from mini,s T (i) (s) to maxi,s T (i) (s).
2. Quantile spacing of u at the 1%, 3%, . . ., 99% quantiles of the values
in T (i) , which places more levels where there are more field values.
3. Variance spacing of u, where placement density is proportional to V̂u ,
which places more levels where the EC is more variable.
Figure 2.8 (right) shows the results for each of the options SC, SGW
and PI. Using SC, all design selection schemes perform similarly. However,
under the SGW and PI covariance matrices, equal spacing outperforms by
a substantial margin. We also experimented with the design styles in the
simulation experiment of Section 3 and found a similar advantage for equal
spacing. We henceforth use equally spaced levels exclusively.
^
g(V(β)) by U
Determinant

^
det(V(β)) by u Spacing

Trace
10^7.0

10^12 10^13 10^14 10^15 10^16 10^17

Quantile

10^6.0

10^6.5

u spacing style

Equal

SW Covariance
Variance
Quantile
Equal

PI Covariance
Variance

10^5.5

^
g(V(β))

SC Covariance
Variance

20

50

100 200

Quantile
Equal

20

U

50

100 200

10^13.5

10^14.5

10^15.5

Level u

Fig 2.8. Factors of optimal design in the LKC regression: the number of levels U (left)
and the spacing mode of the levels uj (right). The plots show how these factors impact the
optimal design criterion g(V [β̂]), for g(M ) = det(M ) and g(M ) = tr(M ). Level counts
ranging from U = 5 to 200 are considered along with equal, quantile and variance-based
spacing. The U plot uses equal spacing and the spacing plot uses U = 50 levels.

2.4. Comparison to Warping. The only published alternative method for
LKC estimation is the warping approach of Taylor and Worsley (2007). The
main idea of warping is to transform the input random fields such that they
become roughly isotropic in a higher-dimensional space while their LKCs
are held constant.
For comparison purposes, we provide a brief procedural review of the
warping method. Suppose the input random fields T (i) lie on a domain S
and have sample locations at sk ∈ S, k = 1, . . . , K, and we are given a triangulation ∆ of S with a vertex at each sk . Warping transforms ∆ to ∆′ ∈ RF ,
a mesh in F -dimensional space also with K vertices. The coordinates of each
vertex of ∆′ are the observed field values T (1) (sk ), T (2) (sk ), . . . T (F ) (sk )

18

ADLER, BARTZ, KOU AND MONOD

from a particular location sk . Formally,
T (j) (sk )
j th coordinate of kth vertex of ∆′ = qP
F
(i)
2
i=1 T (sk )

The (Euclidean) intrinsic volumes of ∆′ are then the estimates of the LKCs
L1 , . . . , Ldim(S) . The precise formulation is provided in Taylor and Worsley
(2007), and involves the edge lengths, interior angles, surface areas and volumes of the triangulation in the transformed space.
An apparent downside of warping is its conceptual as well as computational complexity. The geometric calculations in the transformed space can
be time-consuming for a large triangulation. Unlike ∆, ∆′ does not lie along
a regular grid, so there are no cost-saving computational measures available and runtimes can be up to an order of magnitude longer, which makes
the method time-prohibitive when multiple sets of LKCs need to be found
(e.g. for a time series of repeated experiments). In contrast, our LKC regression requires the computation of only the EC instead of all the intrinsic
volumes. It also works directly in the space of the input fields, where the
regular layout of sample locations permits speedy computation.
LKC regression is fundamentally more transparent and conceptually accessible, while the warping procedure may be more challenging for nonspecialists to understand. Diagnostics are difficult under warping because it
can be unclear to detect the cause of an erroneous estimate, while LKC
regression is interpreted simply as the curve that best matches expected
with observed ECs. Any necessary diagnostics — such as outliers or data
integrity issues — are revealed in a straightforward manner by examining
the EC profile.
3. Simulation Study. To test the accuracy of the Lipschitz-Killing
curvature (LKC) estimates resulting from our regression procedure proposed
in Section 2, a numerical experiment using simulated random fields was
performed. Our aims are twofold: first, to determine optimal settings for
our LKC regression; and second, to compare it to the main alternative, the
warping method (Taylor and Worsley, 2007).
Rather than on the LKC estimates themselves, our simulation study focuses on the estimated 95% threshold for the field maximum, which are
determined by the LKC estimates via (1.4). This threshold is the driver for
95% hypothesis tests, and with practitioners’ needs in mind, the primary
purpose and application for which our LKC method was developed. Our
simulation study computes the standard deviation and bias in this thresh-

THRESHOLDS VIA EULER CHARACTERISTICS

19

old. To assess the generality of the method, we test a range of random fields
commonly found in practice.
3.1. Experimental Metrics. We compare the LKC estimation methods
by their computational runtime, error and bias. Care has been taken in all
cases to employ the optimal implementation of each method and to exclude
fixed costs (e.g. loading the data). Standard deviation describes how much
an estimated threshold u[
95% varies over multiple simulations.
Bias reveals how far, on average, u[
95% is from the true threshold. Bias is
harder to measure as it requires a specification for the input random fields
— a domain and covariance function — for which the true threshold is
known. For stationary and isotropic fields, it is known that the LKCs can be
comparatively simple to compute, particularly if S has a simple topology. We
2
thus base our simulations on the Gaussian covariance (C(x, y) = e−α|x−y| )
on the square, cube and sphere. The true LKCs and thresholds, calculated
via the ECH, are shown in Table 1 (first row). Our experiment is made up of
simulations from these three field types with suitable values of α. The first
LKC, L0 , is not included, Since L0 is the Euler characteristic, it is known
deterministically, and, for the three examples in the table it is 1, 1, and 2,
respectively.
L
√1
2 2α
28.3

Square
L2
u∗95%
2α
−
200 3.72

L
√1
3 2α
19.0

Cube
L2
L3
3(2α) (2α)3/2
120
253

u∗95%
−
3.96

L1
0
0

Sphere
L2
u∗95%
4π(2α)
−
503
3.96

Table 1
True continuous LKCs and 95% thresholds for continuous random fields in the square
(α = 100), in the cube (α = 20), and on the surface of the sphere (α = 20). Both
analytical expressions (upper row) and simulation-specific numerical values (lower row)
are provided. In each case, L0 has an easy known form as the EC of the domain S. The
other LKC formulae are based on analytical results from Adler and Taylor (2007). The
95% thresholds do not have easy analytical forms, but are calculated using the ECH.
These expressions act as the “truth” in our experimental assessment of bias.

A complication for bias computation is that the known thresholds are for
continuous random fields, whereas our experimental fields are generated at
discrete sample sites (see Section 2.1). Discrete and continuous random fields
have different thresholds even when they share the same covariance function
and domain. There are no known analytical forms for the LKCs of discrete
fields (Adler and Taylor, 2007), but it is not hard to show that, for the random fields under consideration, the threshold for discrete fields converges to
that for the continuous limit as the simulation resolution increases. The convergence of discrete to continuous thresholds provides a basis for computing

20

ADLER, BARTZ, KOU AND MONOD

bias; for each estimation method, the average u[
95% at varying grid sizes thus
forms a convergent pattern. The approximate limit yields an estimate of the
continuous threshold, which can be compared to the known true value to
measure bias. (Figure 3.11 provides an illustration.)
3.2. Experimental Factors. Observed random fields are commonly classified by the first four factors listed in Table 2: number of fields, domain,
covariance function, and resolution. Our experiment is factorial and tests
10,000 replicates of each combination of these factors. The domains are the
unit square, the unit cube, and the surface of a unit sphere. The grid sizes
range from G=5 to 200 sample locations in each direction, which spans most
applications and reveals trends in the estimates as the resolution increases.
We simulate F =15 Gaussian random fields in each run, representing a typical
thresholding experiment. Each field is equipped with Gaussian covariance
for a suitable value of α. (All fields are simulated via the turning bands
method (Mantoglou and Wilson, 1982): we use 1,000 lines and 4G points
along each line, twice as many as the authors recommend to help ensure
that simulation error is small compared to estimation error. We adapt the
implementation of Schlather (2009).)
To determine the best settings for LKC regression, we also test the five
positive-definite covariance construction options discussed in Section 2.2.
All estimations use U = 50 equally-spaced levels uj . For the spherical example L1 is always zero for topological reasons, and so neither the LKC
regression nor warping attempt to estimate it. We employ the technique of
Taylor and Worsley (2007) for warping (for the triangulation ∆, we input
a mesh of standard right triangles for the cube and plane and latitudinallongitudinal triangles for the sphere).
Factor
# fields
Domain
Covariance
Grid size
Method
Σ (LKCR only)
# levels (LKCR only)

Values
F = 15
Square, cube, sphere
2
2
c(x, y) = e−100|x−y| (square), e−20|x−y| (cube and sphere)
G = 5, 10, 20, 50, 100, 200 (square only)
Lipschitz-Killing curvature regression (LKCR), warping
From Section 2.2: I, SD, SC, SGW, PI
U = 50

Table 2
Factors for the simulation experiment of Section 3. Factors describe both the nature of the
input random fields (first four rows) and the LKC estimation methods (last three rows).

21

THRESHOLDS VIA EULER CHARACTERISTICS

3.3. Results. Overall, we find that LKC regression is much faster than
warping, and maintains similar accuracy. The increase in speed is roughly
of a factor of up to eight for high-resolution random fields on the square and
cube, which are the most common domains in applications. This gain comes
with negligible loss in standard deviation or bias. Among LKC regression
covariance styles, the smoothed diagonal (SD) method is a consistent winner,
providing the best combination of accuracy and speed.
3.3.1. Runtime. Computational runtime provides a compelling case for
LKC regression. The averages shown in Figure 3.9 show a dramatic advantage for most LKC regression methods (dashed and dotted lines) against
warping (solid line). This advantage is most pronounced for large grids that
are commonly found in practice. For cubic fields with G = 100, warping
takes an average 40 sec compared to 6 sec for LKC regression methods I,
SD and pseudoinverse (PI).
Among LKC regression styles, the slowest are smoothed covariance (SC)
and Sampson-Guttorp warping (SGW). These two options involve complex
estimation of the off-diagonal terms in the regression covariance matrix,
which thus incurs a heavy fixed cost in covariance estimation that dominates
the runtime.

Median computational runtime in seconds

Average Runtime by Method
LKCR: I
LKCR: SD
Square

LKCR: SC
LKCR: SW
Cube

LKCR: PI
Warping
Sphere

10^1

10^0

10^−1

5

10

20

50

100

200

5

10

20

50

100

5

10

20

50

100

Grid size

Fig 3.9. Median runtimes for estimating LKCs in the simulation experiment. Results are
given by method (line styles), random field domain (square, cube or sphere) and simulation
grid size (horizontal axis). Medians are over 10,000 independent replications of the input
data.

3.3.2. Standard Deviation. Comparing the standard deviation of LKC
regression and warping method depends on simulation resolution. Figure

22

ADLER, BARTZ, KOU AND MONOD

3.10 shows generally similar standard deviations between warping and the
top-performing LKC regression styles, identity and SD covariances, for large
grid sizes. For fields with G = 50, the differences are slight: 0.0096 for warping versus 0.0103 for LKC regression using SD covariance for square fields,
0.0056 versus 0.0059 for spherical fields, and 0.0127 versus 0.0158 for cubic
fields. Moreover, both standard deviations are negligible when viewed with
regard to the underlying application. Applied random fields usually have
thresholds of around 5 (see Section 4 for an example). Standard deviations
on the scale of those we observed have little impact on the threshold. They
essentially do not affect the task of identifying significant regions in highresolution fields. Warping has smaller standard deviations for fields with low
resolutions, particularly G ≤ 10.
The best LKC regression covariance structure is shown to be SD, although
I and PI also perform well. The poor performance of SC and SGW is likely
due to overfitting; both fit nonparametric models of the regression covariance
matrix. Errors in these procedures accumulate in the ultimate regression,
causing higher standard deviations. The I and PI covariance models are less
flexible but more parsimonious, stabilizing the estimates. SD is similar to
I, but comprises a relatively conservative smoothing procedure (see Figure
2.6) to model heteroskedasticity.
Standard Error of 95% Threshold by Method

Standard error of 95% threshold

LKCR: I
LKCR: SD
Square

LKCR: SC
LKCR: SW
Cube

LKCR: PI
Warping
Sphere

10^−0.5

10^−1.0

10^−1.5

10^−2.0

5

10

20

50

100

200

5

10

20

50

100

5

10

20

50

100

Grid size

Fig 3.10. Standard deviations of 95% threshold estimates.

3.3.3. Bias. As mentioned in Section 3.1, a direct bias calculation is difficult because true values of LKCs for discrete random fields are unknown;
true values are known only for continuous random fields with Gaussian covariance. However, the continuous LKCs can be approximated by extrapo-

23

THRESHOLDS VIA EULER CHARACTERISTICS

Average 95% Threshold by Method
LKCR: I
LKCR: SD

LKCR: SC
LKCR: SW

LKCR: PI
Warping
Cube

3.9

3.9

3.8

3.8

3.6

3.7

3.7

3.5
5

10

20

50 100 200

3.6

3.6

3.4

Median 95% threshold

3.7

Square

Empirical
Continuous
Sphere

5

10

20

50 100 200

5

10

20

50 100 200

Grid size

Fig 3.11. Medians of 95% threshold estimates. The extrapolated curves are fitted using
nonlinear regression, Equation (3.9). The unlinked empirical dots are the “true” thresholds
approximated using 150,000 simulated discrete random fields.

lating from a pattern of discrete LKCs at varying resolutions. Similarly, the
median estimated 95% thresholds in Figure 3.11 show clear convergent patterns. Assuming convergence of order ς, the limit is found using nonlinear
regression:
(3.9)

u95% (g) = u∗95% + βg ς + ε(g).

Here, u∗95% denotes the implied continuous threshold, while u95% (g) is the
estimated 95% threshold at grid size g. The unknown parameters β and
ς govern the rate of approach to the continuous value, for which typical
estimates are β̂ = −7.0 and ςˆ = −1.5; ε is an error term that we take as
normal for simplicity.
∗ , an approximate continuous
The end result of regression (3.9) is u[
95%
threshold. One regression is needed for each method and random field type,
representing a single pattern of dots in Figure 3.11. The fitted regression lines
are also shown; bias is measured as the difference between approximate and
true continuous thresholds.
Neither LKC regression nor warping emerges as the clear winner in bias.
Figure 3.12 (left column) indicate small, similar levels of bias across all
methods. The only outlier is the SGW covariance, which again produces
poor estimates, likely to be the result of an overfitting of covariance.
An alternative approach to measuring bias is to approximate the true
95% thresholds for discrete random fields with intensive simulation. Given
B simulated fields, the 95% sample quantile of the field maxima provides an

24

ADLER, BARTZ, KOU AND MONOD

approximation (we use B =150,000). These estimates appear in Figure 3.11
(“Empirical”), and are used to measure bias in Figure 3.12 (right column).
The result is the same: all methods demonstrate a similar level of bias,
however Figure 3.11 also reveals an interesting self-regulatory mechanism
in LKC estimation. The estimated discrete thresholds always lie between
the continuous thresholds computed via the ECH (1.4) and the discrete
thresholds computed via simulation. Despite their basis in continuous theory,
the estimated thresholds of the LKC regression approach the discrete truth.
Bias in 95% Threshold by Method
Square
Estimated vs. Continuous

Square
Estimated vs. Empirical (G = 50)

Cube
Estimated vs. Continuous

Cube
Estimated vs. Empirical (G = 50)

Sphere
Estimated vs. Continuous

Sphere
Estimated vs. Empirical (G = 50)

Estimation method

LKCR: I
LKCR: SD
LKCR: SC
LKCR: SW
LKCR: PI
Warping

LKCR: I
LKCR: SD
LKCR: SC
LKCR: SW
LKCR: PI
Warping

LKCR: I
LKCR: SD
LKCR: SC
LKCR: SW
LKCR: PI
Warping
−0.02

0.00

0.02

0.04

0.06

0.08

−0.02

0.00

0.02

0.04

0.06

0.08

Bias in 95% threshold

Fig 3.12. Biases of two types of 95% threshold estimates: the continuous threshold estimated via extrapolation (left); and the discrete threshold for resolution G = 50.

Fig 4.13. Examples of three residual fields from a voxel-wise AR(1) model fit to the FIAC
brain data. The scans shown come from the first experimental subject at times 1, 50 and
100. The outer gray layer shows the surface of the brain, while the red and blue show
domains with positive and negative residuals, respectively. Light and dark colors represent
the data thresholded at the single t-test 95% and 99% thresholds of 1.96 and 2.58.

THRESHOLDS VIA EULER CHARACTERISTICS

25

4. Applications. In this section we discuss two applications. One, on
the analysis of CMB data, is discussed briefly at the end of the section.
For a more detailed example, we return to the fMRI data discussed in the
Introduction.
4.1. A fMRI example. We return to the language priming experiment of
Dehaene-Lambertz et al. (2006) and apply our LKC regression method.
Contrary to the group-level comparison of Section 1, we conduct a withinsubject study to identify stimulated regions for each specific individual. Here
the stimulus refers to the difference in activation between the “different
speaker” and “same speaker” conditions, under which each subject hears
a sentence repeated twice either by two different speakers or by the same
speaker. The subject is then scanned 191 times, alternatively undergoing
the two stimuli. There are 14 such data sets, one for each individual, which
are analyzed by Taylor and Worsley (2006) using the warping thresholding
method. We reanalyze the data using our LKC regression method.
Prior to the analysis we conducted a standard preprocessing (as detailed
in Section 2 of Taylor and Worsley, 2006) of the data: First, a subject’s scans
are aligned to correct for motion drift (Smith et al., 2004). An AR(1) time
series model is then fitted to every voxel of the brain. The response variable
is the fMRI activation level, which is a time series indexed by each of the 191
time points; in the AR(1) model, the stimulus effect is treated as a regression
coefficient in front of a time-varying contrast that represents the level of
stimulus administered at each time. There is one such model for every voxel
in the brain; fitting all such models produces a random field of estimated
stimulus effects across the brain: this is the field of interest, analogous to the
field T in Section 1. The random fields T (i) , i = 1, . . . , 191 correspond to
the AR(1) residuals at every voxel and time. As in Section 1, the threshold
is obtained from the T (i) and then applied to T to identify regions where
the stimulus effect is significant. There is one analysis — one T and one set
of T (i) — for every subject. As an example, three residual fields for the first
subject in the study, corresponding to three separate time points, are shown
in Figure 4.13. In Taylor and Worsley (2006), warping is used to estimate
the threshold, which found values of around 5 for each subject.
We are interested in analyzing how these thresholds, found using warping (Taylor and Worsley, 2007), compare to those found using our proposed
LKC regression method. For the purpose of such a comparison, the thresholds are estimated using warping and our LKC regression with smoothed
diagonal (SD) covariance (which is found to be the best-performing option
found in Section 3).

26

1
2
3
4
5
6
7
8
9
10
11
12
13
14

ADLER, BARTZ, KOU AND MONOD
95% thres.
LKC Warp
5.01
4.94
4.93
4.75
5.10
4.98
5.03
5.14
5.06
5.11
4.99
4.98
5.03
5.09
4.94
5.05
4.97
5.02
5.15
5.28
4.98
4.91
5.01
4.89
5.01
5.01
5.01
5.13

L1
LKC
1,686
1,237
943
1,252
1,961
1,158
1,263
1,737
2,017
1,437
2,394
1,333
1,933
1,757

L2
Warp
274
-32
201
249
49
22
419
382
191
386
320
307
255
105

LKC
7,774
7,801
4,694
5,584
9,484
7,712
4,702
8,860
8,782
4,857
8,967
6,528
7,658
7,724

L3
Warp
2,343
2,351
2,242
2,169
2,277
2,427
2,091
2,162
2,235
2,245
2,274
2,113
2,063
2,262

LKC
26,706
23,878
15,234
19,886
30,229
23,317
17,543
28,163
31,009
17,988
30,473
21,540
27,509
26,367

Warp
29,382
32,334
30,659
29,537
30,683
34,124
28,103
26,701
31,327
30,403
31,133
28,898
29,670
32,376

Runtime
LKC Warp
18.3
38.7
12.5
34.9
15.0
38.7
15.8
33.2
12.7
33.8
14.8
32.7
16.7
37.5
13.9
34.1
20.9
35.4
18.5
42.6
15.3
37.3
15.3
36.1
13.6
34.5
12.5
34.1

Table 3
Per-subject LKC estimation results for LKC regression (LKCR) and warping. The data
are residual random fields from a voxel-wise AR(1) model fit to the FIAC brain scan
data, with F = 191 fields per subject. The estimated 95% thresholds (left), LKCs
(middle) and runtime (in seconds, right) are displayed for each subject and method.

Table 3 compares the LKC estimates and the computational time of each
method. The major difference is in the latter (rightmost column): the LKC
regression is over twice as fast as warping. At first, while any improvement in
computational speed is desirable, it would not seem that this improvement is
all that important. Statistical analysis of fMRI data is typically done offline,
as opposed to the actual production of fMRI data which is online and must
be minimised for subject comfort and cost considerations.
Nevertheless, computation time is still an issue, as evidence by the fact
that practitioners have developed data structures optimized for 3D fields
(e.g., Theis and Tanaka, 2005) and schemes for parallel processing (e.g.,
Zhao et al., 2007; Wilde et al., 2009). Computational firepower is necessary
because over the course of an fMRI study, the number of hypothesis tests
conducted and thresholds found easily reaches into the hundreds. These computational demands occur for two reasons: first, there are usually multiple
experimental stimuli, each with a categorical set of possible conditions (for
instance, Taylor and Worsley (2006) consider three factors with two or three
conditions each, yielding 12 significance analyses); and second, it is also common to perform subject-specific tests (one test per subject), or time-specific
tests (one test per unit time) (see Beckmann, Jenkinson and Smith, 2003).
In light of these repeated computational demands that are inherent to various applications, our LKC regression can save substantial development time.

THRESHOLDS VIA EULER CHARACTERISTICS

27

In addition to speed, accuracy also carries great importance for practictioners. As seen from the spread of the estimated 95% thresholds (leftmost
column), LKC regression and warping produce very similar thresholds, with
the mean values of the 14 subjects being 5.01 and 5.02, respectively. Thus,
as to be expected, the activated regions change only slightly (not shown)
and we are still able to reproduce the activation in the left and right midtemporal gyri seen in Taylor and Worsley (2007).
What we have not yet checked in detail, but, from preliminary studies
are quite certain is true, is that in higher dimensional problems the gain
in speed of LKC regression over warping will be quite significant. What is
definitely true is that the LKC approach, even in this three dimensional
example, involves much less coding, since, following Section 3, it is basically
just a regression analysis.
As a side note, it is interesting that the similarity in thresholds between
LKC regression and warping belies often large differences in the estimated
LKCs (middle three columns). These phenomena underscore the fact that
different LKC combinations can lead to the same 95% thresholds, which
are what matter in practice. It is impossible to tell which is correct because only the thresholds are observable, not the LKCs themselves. Interestingly, the discrepancies in L1 and L2 are not apparent in the cubic random
field simulations of Section 3, where LKC regression and warping produce
approximately the same LKC estimates. We suspect that isotropy plays a
role: the simulated cubic fields are isotropic, while the FIAC data are likely
anisotropic.
4.2. A cosmic microwave background radiation example. Working from
an earlier version of the current paper, Fantaye et al. (2014) have used LKC
regression to estimate the LKCs (‘Minkowski functionals’ in their language)
for CMB fluctuation models.
To quote from their paper: “A general trend in modern cosmological research is the implementation of more and more sophisticated statistical tools
to perform data analysis. Indeed, as well-known cosmological data have
reached over the last decade an un- precedented accuracy, so that it has
become customary to speak about a golden era for Cosmology, featuring
a data deluge from a bunch of satellite - and ground based-experiments.
As the data grow in size and precision, more and more detailed questions
can be addressed, and exploiting techniques at the frontier of statistical and
mathematical research becomes mandatory to warrant a full exploration of
the available evidence.
Among these techniques, stochastic geometry tools have now become very

28

ADLER, BARTZ, KOU AND MONOD

well established, especially in the field of Cosmic Microwave Background radiation experiments. In this area, one of the most popular geometric tools for
data analysis are certainly the so-called Minkowski functionals (MFs), which
have been extensively exploited as tools to search for non-Gaussianities,
anisotropies, asymmetries and other features of CMB data. The use of MFs
in Cosmology goes back at least to [1, 2]; a complete bibliography would
certainly include hundreds of entries, so we refer only to the earlier works
by [3–10] and to the more recent ones by [11–16].” (References not included
here.)
We recommend their paper, which contains detailed calculations and scientific discussions, as an excellent application of the techniques developed
here (and, of course, for the many other things there).
5. Conclusion. We have presented a new method, the Lipschitz-Killing
curvature regression, for estimating tail probabilities of a Gaussian random
field. The LKC regression procedure is to estimate the Lipschitz-Killing curvatures and then substitute them into the Euler characteristic heuristic (1.4)
to generate a (1 − α) threshold. The estimation does not require knowledge
of the covariance structure of the field or even whether it is isotropic. It need
not be Gaussian, although it must be Gaussian related. The high efficiency
of our method hinges on the fact that the Euler characteristic can be well
estimated from the data when the exceedence level u is small or moderate.
The LKC regression allows us to leverage the estimation strength at a low
level u for accurate approximation at high exceedence level.
The primary advantage of the LKC regression is its simplicity and speed.
It is transparent, which offers straightforward interpretation, implementation, and diagnosis. It runs faster than its chief competitor, warping (Taylor and Worsley,
2007), with comparable accuracy. The key procedural difference is that LKC
regression computes Euler characteristics of the input random fields, which
is a fast computation when the sample sites lie on a grid. In practical terms,
for the two- and three-dimensional random fields we considered, the gain in
speed ranges from a factor of two to a factor of eight. It is most pronounced
for high resolution fields, which are commonly found in practice. We believe that the gain rate will be considerably higher in higher dimensional
cases, such as the scale-space (5D) of Siegmund and Worsley (1995) and the
rotation-space fields (8D) of Shafie et al. (2003), for both of which the analytic evaluation of the LKCs is orders of magnitude harder than the isotropic
cases we have considered in this paper (cf. Adler, Subag and Taylor (2012)).
Appendix. In this section we shall briefly discuss the Gaussian kinematic formula, which gives an exact formula for the expected Euler charac-

THRESHOLDS VIA EULER CHARACTERISTICS

29

teristic E[ϕ(Au )]. We shall explain its general structure without going into
technical details.
A.1. Gaussian and Gaussian related random fields. The basic building
blocks of all the random fields we consider are smooth Gaussian random
fields, or processes,


g = g1 , . . . , gk : S ∈ RN → Rk .
The first assumption that we place on g, and the most important one,
is that each of its k components are twice differentiable, and that these
derivatives are themselves continuous. Some additional minor assumptions
of non-degeneracy also need to be made, but since these almost always hold
in practice we direct the interested reader to Chapter 11 of Adler and Taylor
(2007) or Chapter 4 of Adler, Taylor and Worsley (2015?) for details. For
the rest of this paper we shall assume these conditions are met. The second
assumption is a minor one, that all means be fixed at zero. A little more
restrictively, we also assume that the gj have constant variance throughout
S. Note that this is a much weaker assumption than either isotropy or stationarity, which we do not require, and is achievable in general by replacing
a random field that does not have constant variance by a normalized version of itself. Unless stated otherwise, for the remainder of this section we
shall assume that this constant variance be 1. Finally, we assume that the
components of g are all independent and identically distributed.
Real valued Gaussian related random fields f are defined by taking a
smooth function F : Rk → R and setting f (s) = F (g(s)) for all s ∈ S,
where g is as above. Gaussian related fields are typically quite different to
Gaussian ones (for which k = 1 and F is the identity function). Three useful
examples are given by the following choices for F , where in the third we set
k = n + m.
(A.1)

k
X
1

x2i ,

√
x1 k − 1
,
P
( k2 x2i )1/2

P
m n1 x2i
Pn+m 2 .
n n+1
xi

The corresponding random fields are known as χ2 fields with k degrees of
freedom, the t field with k − 1 degrees of freedom, and the F field with n
and m degrees of freedom. These three random fields all have very different
spatial behavior, and each is as fundamental to the statistical applications
of random field theory as is its corresponding univariate distribution to
standard statistical theory.

30

ADLER, BARTZ, KOU AND MONOD

Throughout we shall assume that F is twice continuously differentiable,
and that for all real u the sets F −1 [u, ∞) ⊂ Rk are well behaved, as described
in the following subsection. Under these conditions, Gaussian related fields
are also twice continuously differentiable, with well behaved excursion sets.
A.2. Parameter spaces. The modern general theory of Gaussian fields, as
developed in Adler and Taylor (2007), allows the parameter S to be what is
known as a Whitney stratified manifold, satisfying some mild side conditions.
Roughly speaking, these are compact subsets of RN which can be written
as a disjoint finite union
S=

N
G

∂i S,

i=1

where ∂i S is an i-dimensional manifold (the open i-dimensional ‘boundary’
of S).
An easy example is given by a N -dimensional rectangle, in which ∂N S
is its interior, ∂N −1 S the collection of its (open) (N − 1)-dimensional faces,
etc., down to ∂0 S which is the collection of its corners. If S were a ball, we
would decompose it into its interior, again ∂N S, and its surface, the sphere
∂N −1 S.
The ‘Whitney’ part of the definition contains rules about how all the various pieces must be put together, and details can be found in Adler and Taylor
(2007). More or less every (non-fractal) parameter space arising in statistical
practice satisfies these rules.
A.3. The expected Euler characteristic of excursion sets. The Euler characteristic ϕ(S) of a well behaved set S ⊂ RN is a topological invariant, that,
typically, is easy to compute. For example, if S is two-dimensional, then
ϕ(S) is simply the number of connected components in S minus the number of holes. If S is three-dimensional, then it is the number of components
minus the number of handles plus the number of internal holes.
If S is made up of a union of cubes, all of whose corners sit at the points
of a rectangular lattice, then
X
ϕ(S) =
(−1)n µn ,
n=1

where µn is the number of distinct facets of dimension n in S.
The central result behind this paper is then the following formula, which
is a special case of a far more general result known as the Gaussian kinematic
formula. (cf. Adler and Taylor (2007); Taylor (2006); Taylor and Adler (2003)

THRESHOLDS VIA EULER CHARACTERISTICS

31

for details.) It says that for Gaussian related f = F (g) field defined over
well-behaving sets S,
(A.2)

E {ϕ ({t ∈ S : f (t) ≥ u})} =

dim
XS
j=0

Lj (S) ρFj (u).

The Lj (S) are the Lipschitz-Killing curvatures, or intrinsic volumes, of S,
and estimating these in statistical practice was one of the aims of this paper. If the gi are stationary and isotropic, and the variance of their first
order partial derivatives is given by the second spectral moment λ2 , then
j/2
E
Lj (S) = λ2 LE
j (S) where the Lj (S) are the Euclidean Lipschitz-Killing
curvatures . The Euclidean Lipschitz-Killing curvatures are known under a
variety of names, including Quermassintegrales, Minkowski or Steiner functionals, integral curvatures, and intrinsic volumes, the differences between
them generally being of ordering and scaling. In general, LE
j (S) can be
thought of as a measure of the ‘j-dimensional size’ of S. For example, when
E
N = 2, LE
2 (S) is the two dimensional area of S, L1 (S) is half its boundary
E
length, and L0 (S) its Euler characteristic. When N = 3, LE
3 (S) is the threeE (S) is twice the
(S)
is
half
the
surface
area,
L
dimensional volume of S, LE
1
2
caliper diameter of S, (where the caliper diameter of a convex S is defined
by placing the solid between two parallel planes (or calipers), measuring
the distance between the planes, and averaging over all rotations of S) and
E
LE
0 (S) is again the Euler characteristic . In fact L0 (S) = L0 (S) = ϕ(S) in
all cases.
While these are simple examples of Lipschitz-Killing curvatures, when
the underlying gi are not isotropic their covariance induces a Riemannian
metric on S, and the Lj (S) become quite complicated. To be a little more
specific, they involve integrals of traces of powers of the curvature tensor
and second fundamental forms, with respect to the volume form determined
by the induced metric. In brief, these are not typically things that a practitioner wants to compute analytically, particularly since computing them
also requires knowledge of the covariance function of the gi , something which
typically needs to be estimated from data.
On the other hand, the functions ρFj in (A.2) are much easier to compute,
and in many cases are given by
ρFj (u) = (−1)j (2π)−j/2

dj
P {F (Z) ≥ x}
dxj

x=g(u)

,

where Z ∼ N (0, Ik×k ) and g is a function determined by F . In such cases,
computing the ρk is thus simply a matter of calculus, and for most interesting

32

ADLER, BARTZ, KOU AND MONOD

F this has already been done. (Adler and Taylor (2007); Adler, Taylor and Worsley
(2015?) have a number of useful examples with references to others.) A particularly simple example is given in the Gaussian case, in which, as mentioned above, k = 1 and F is the identity function. Then the ρj are then
given by
ρj (u) = (2π)−(j+1)/2 Hj−1 (u)e−u

2 /2

j = 0, 1, . . . N,

where Hj , j ≥ 1 is the j-th Hermite polynomial and
∆

H−1 (u) =

√
2
2πeu /2 P (N (0, 1) ≥ u).

To give an example in a more complicated case, we consider the χ2 case,
with k degrees of freedom, so that F is given by the leftmost function in
(A.1). In this case, the ρj are, for j ≥ 1 and u > 0,
⌊ j−1 ⌋

j−1−2l
2
X
X
u(k−j)/2 e−u/2
ρj (u) =
j/2
(k−2)/2
(2π) Γ(k/2)2
l=0 m=0


k−1
(−1)j−1+m+l (j − 1)! m+l
u
.
×1{k≥j−m−2l}
m!l!2l
j − 1 − m − 2l

When j = 0,

ρ0 (u) = P χ2k ≥ u .

While this may look a little complicated, it is trivial to code, and the LKC
regression procedure then continues as in the Gaussian case.
REFERENCES
Adler, R. J. (1981). The Geometry of Random Fields. John Wiley & Sons Ltd., Chichester. Reprinted in 2010 by SIAM. MR82h:60103
Adler, R. J. (2000). On Excursion Sets, Tube Formulae, and Maxima of Random Fields.
The Annals of Applied Probability 10 1–74.
Adler, R. J., Subag, E. and Taylor, J. E. (2012). Rotation and Scale Space Random
Fields and the Gaussian Kinematic Formula. Annals of Statistics 40 2910–2942.
Adler, R. J. and Taylor, J. E. (2007). Random Fields and Geometry. Springer.
Adler, R. J. and Taylor, J. E. (2011). Topological Complexity of Smooth Random
Functions. Lecture Notes in Mathematics 2019. Springer, Heidelberg. Lectures from the
39th Probability Summer School held in Saint-Flour, 2009, École d’Été de Probabilités
de Saint-Flour. [Saint-Flour Probability Summer School]. MR2768175
Adler, R. J., Taylor, J. E. and Worsley, K. J. (2015?). Applications of Random
Fields and Geometry: Foundations and Case Studies. Springer-Verlag In preparation,
early chapters available at http://webee.technion.ac.il/people/adler/publications.html.

THRESHOLDS VIA EULER CHARACTERISTICS

33

Bardeen, J. M., Bond, J. R., Kaiser, N. and Szalay, A. S. (1986). The Statistics of
Peaks of Gaussian Random Fields. Astrophysical Journal 304 15–61.
Beckmann, C. F., Jenkinson, M. and Smith, S. M. (2003). General Multilevel Linear
Modeling for Group Analysis in FMRI. Neuroimage 20 1052-1063.
Cao, J. and Worsley, K. J. (1999). The Geometry of Correlation Fields With an Application to Functional Connectivity of the Brain. Annals of Applied Probability 9 1021–
1057. MR2000k:60099
Dehaene-Lambertz, G., Dehaene, S., Anton, J. L., Campagne, A., Ciuciu, P., Dehaene, G. P., Denghien, I., Jobert, A., LeBihan, D., Sigman, M., Pallier, C. and
Poline, J. B. (2006). Functional Segregation of Cortical Language Areas by Sentence
Repetition. Human Brain Mapping 27 360-371.
Fantaye, Y., Hansen, F., Maino, D. and Marinucci, D. (2014). Cosmological Applications of the Gaussian Kinematic Formula. arxiv.org/pdf/1406.5420v1.pdf.
Friston, K. J., Worsley, K. J., Frackowiak, R. S. J., Mazziotta, J. C. and
Evans, A. C. (1994). Assessing the Significance of Focal Activations Using Their Spatial Extent. Human Brain Mapping 1 214–220.
Gott, J. R., Colley, W. N., Park, C. G., Park, C. and Mugnolo, C. (2007). Genus
Topology of the Cosmic Microwave Background From the WMAP 3-year Data. Monthly
Notices of the Royal Astronomical Society 377 1668–1678.
Gott, J. R., Hambrick, D. C., Vogeley, M. S., Kim, J., Park, C., Choi, Y. Y.,
Cen, R. and Ostriker, K. J. P. Nagamine (2008). Genus Topology of Structure in
the Sloan Digital Sky Survey: Model Testing. Astrophysical Journal 675 16–28.
Hall, P., Fisher, N. I. and Hoffman, B. (1994). On the Nonparametric Estimation of
Covariance Functions. Annals of Statistics 22 2115-2134.
Imiya, A. and Eckhardt, U. (1999). The Euler Characteristics of Discrete Objects and
Discrete Quasi-objects. Computer Vision and Image Understanding 75 307-318.
Kilner, J. M. and Friston, K. J. (2010). Topological Inference for EEG and MEG
Data. Annals of Applied Statistics 4 1272-1290.
Loader, C. (1999). Local Regression and Likelihood. Springer, New York.
Loader, C. (2010). locfit: Local Regression, Likelihood and Density Estimation R package
version 1.5-6.
Mantoglou, A. and Wilson, J. L. (1982). The Turning Bands Method for Simulation of
Random Fields Using Line Generation by a Spectral Method. Water Resources Research
18 1379-1394.
Pukelsheim, F. (2006). Optimal Design of Experiments. Society for Industrial and Applied Mathematics.
Rao, C. R. (1962). A Note on a Generalized Inverse of a Matrix With Applications
to Problems in Mathematical-statistics. Journal of the Royal Statistical Society Series
B-Statistical Methodology 24 152-158.
Sampson, P. and Guttorp, P. (1992). Nonparametric Estimation of Nonstationary Spatial Covariance Structure. Journal of the American Statistical Association 87 108-119.
Schlather, M. (2009). RandomFields: Simulation and Analysis of Random Fields R
package version 1.3.41.
Shafie, K., Sigal, B., Siegmund, D. and Worsley, K. J. (2003). Rotation Space
Random Fields With an Application to fMRI Data. Annals of Statistics 31 1732–1771.
Siegmund, D. O. and Worsley, K. J. (1995). Testing for a signal with unknown location
and scale in a stationary Gaussian random field. Ann. Statist 23 608–639.
Smith, S. M., Jenkinson, M., Woolrich, M. W., Beckmann, C. F.,
Behrens, T. E. J., Johansen-Berg, H., Bannister, P. R., De Luca, M., Drobnjak, I., Flitney, D. E., Niazy, R. K., Saunders, J., Vickers, J., Zhang, Y. Y.,

34

ADLER, BARTZ, KOU AND MONOD

De Stefano, N., Brady, J. M. and Matthews, P. M. (2004). Advances in Functional and Structural MR Image Analysis and Implementation as FSL. Neuroimage 23
S208-S219.
Taylor, J. E. (2006). A Gaussian kinematic formula. Annals of Probability 34 122-158.
Taylor, J. E. and Adler, R. J. (2003). Euler characteristics for Gaussian fields on
manifolds. Annals of Probability 31 533–563.
Taylor, J. E., Takemura, A. and Adler, R. J. (2005). Validity of the
Expected Euler Characteristic Heuristic. Annals of Probability 33 1362–1396.
MR2150192 (2006b:60073)
Taylor, J. E. and Worsley, K. J. (2006). Inference for Magnitudes and Delays of
Responses in the FIAC Data Using BRAINSTAT/FMRISTAT. Human Brain Mapping
27 434-441.
Taylor, J. E. and Worsley, K. J. (2007). Detecting Sparse Signals in Random Fields,
With an Application to Brain Mapping. Journal of the American Statistical Association
102 913-928.
Theis, F. J. and Tanaka, T. (2005). A fast and efficient method for compressing
fMRI data sets. In Artificial Neural Networks: Formal Models and Their Applications ICANN 2005, Pt. 2, Proceedings, ( Duch, W and Kacprzyk, J and Oja, E and Zadrozny,
S, ed.). Lecture Notes in Computer Science 3697 769-777.
Torres, S. (1994). Topological Analysis of COBE-DMR Cosmic Microwave Background
Maps. Astrophysical Journal 423 9–12.
Vogeley, M. S., Park, C., Geller, M. J., Huchin, J. P. and Gott, J. R. (1994).
Topological Analysis of the CfA Redshift Survey. Astrophysical Journal 420 525–544.
Wilde, M., Foster, I., Iskra, K., Beckman, P., Zhang, Z., Espinosa, A., Hategan, M., Clifford, B. and Raicu, I. (2009). Parallel Scripting for Applications at
the Petascale and Beyond. Computer 42 50-60.
Worsley, K. J. (1994). Local Maxima and the Expected Euler Characteristic of Excursion Sets of χ2 , F and t Fields. Advances in Applied Probability 26 13–42. MR94i:60064
Worsley, K. J. (1995a). Boundary Corrections for the Expected Euler Characteristic of
Excursion Sets of Random Fields, With an Application to Astrophysics. Advances in
Applied Probability 27 943–959. MR97b:60068
Worsley, K. J. (1995b). Estimating the Number of Peaks in a Random Field Using
the Hadwiger Characteristic of Excursion Sets, With Applications to Medical Images.
Annals of Statistics 23 640–669.
Zhao, Y., Hategan, M., Clifford, B., Foster, I., von Laszewski, G., Nefedova, V.,
Raicu, I., Stef-Praun, T. and Wilde, M. (2007). Swift: Fast, Reliable, Loosely
Coupled Parallel Computation. IEEE Congress on Services 0 199-206.

Robert J. Adler
Electrical Engineering
Technion, Haifa, Israel 32000
E-mail: robert@ee.technion.ac.il
URL: robert.net.technion.ac.il

Kevin Bartz
Renaissance Technologies LLC
600 Route 25A, East Setauket, NY 11733
E-mail: kevin@kevinbartz.com

THRESHOLDS VIA EULER CHARACTERISTICS
S. C. Kou
Department of Statistics
Harvard University
1 Oxford Street Cambridge, MA 02138
E-mail: kou@stat.harvard.edu
URL: www.people.fas.harvard.edu/∼skou

35

Anthea Monod
Department of Systems Biology
Columbia University
1130 St. Nicholas Avenue
New York, NY 10032
E-mail: am4691@cumc.columbia.edu
URL: systemsbiology.columbia.edu/people/anthea-monod

