1

FAST Adaptive Smoothing and Thresholding for
Improved Activation Detection in Low-Signal fMRI
Israel Almodóvar-Rivera and Ranjan Maitra

arXiv:1702.00111v2 [stat.ME] 2 Feb 2017

Abstract
Functional Magnetic Resonance Imaging is a noninvasive tool used to study brain function. Detecting activation
is challenged by many factors, and even more so in low-signal scenarios that arise in the performance of high-level
cognitive tasks. We provide a fully automated and fast adaptive smoothing and thresholding (FAST) algorithm that
uses smoothing and extreme value theory on correlated statistical parametric maps for thresholding. Performance on
simulation experiments spanning a range of low-signal settings is very encouraging. The methodology also performs
well in a study to identify the cerebral regions that perceive only-auditory-reliable and only-visual-reliable speech
stimuli as well as those that perceive one but not the other.
Index Terms
AM-FAST, AR-FAST, Adaptive Segmentation, AFNI, BIC, circulant matrix, CNR, Cluster Thresholding,
Gumbel distribution, SNR, SPM, SUMA, reverse Weibull distribution

I. I NTRODUCTION
UNCTIONAL Magnetic Resonance Imaging (fMRI) [1]–[9] studies the spatial characteristics and
extent of brain function while at rest or, more commonly, while performing tasks or responding
to external stimuli. In the latter scenario, which is the setting for this paper, the imaging modality
acquires voxel-wise Blood-Oxygen-Level-Dependent (BOLD) measurements [10], [11] at rest and during
stimulation or performance of a task. After pre-processing, a statistical model such as a general linear
model [4], [12] is fit to the time course sequence against the expected BOLD response [13]–[15]. Statistical
Parametric Mapping (SPM) [16] techniques provide voxel-wise test statistics summarizing the association
between the time series response at that voxel and the expected BOLD response [3]. The map of test
statistics is then thresholded to identify significantly activated voxels [17]–[19].
The analysis of fMRI datasets is challenged [20]–[23] by factors such as scanner, inter- and intrasubject variability, voluntary/involuntary or stimulus-correlated motion and also the several-seconds delay
in the BOLD response while the neural stimulus passes through the hemodynamic filter [23]–[25]. Preprocessing steps [22], [26] mitigate some of these effects, but additional challenges are presented by
the fact that an fMRI study is expected to have no more than 1-3% activated voxels [8], [27]. Also,
many activation studies involving high-level cognitive processes have low contrast-to-noise ratios (CNR),
throwing up more issues as illustrated next.

F

A. Activation Detection during Noisy Audiovisual Speech
The most important way in which humans communicate is through speech [28]–[30], which the brain
is particularly adept at understanding, even when this form of communication is degraded by noise. It is
I. Almodóvar-Rivera is with the Department of Biostatistics and Epidemiology at the University of Puerto Rico, Medical Science Campus,
San Juan, Puerto Rico, USA.
R.Maitra is with the Department of Statistics, Iowa State University, Ames, Iowa, USA.
This research was supported in part by the National Institute of Biomedical Imaging and Bioengineering (NIBIB) of the National Institutes
of Health (NIH) under its Award No. R21EB016212, I. Almodóvar-Rivera also acknowledges receipt of a fellowship from Iowa State
University’s Alliance for Graduate Education and the Professoriate (AGEP) program for underrepresented graduate students in STEM fields.
The content of this paper however is solely the responsibility of the authors and does not represent the official views of either the NIBIB
or the NIH.

believed [31] that the ability to understand noisy speech may be explained by the capacity of the brain
for multisensory integration of input information, with the information coming in independently from the
visual and the auditory modalities integrated to reduce noise and allow for more accurate perception [32],
[33]. Recently, [31] studied the role of the superior temporal sulcus (STS) in the perception of noisy
speech, through fMRI and behavioral experiments and established increased connectivity between the
STS and the auditory or the visual cortex depending on whichever modality was more reliable, that is,
less noisy.
[31] reported their fMRI analyses at the level of regions of interest (ROIs) drawn on the STS and the
auditory and visual cortices. However, the full benefit of fMRI can be realized only if we are able to
go beyond statements on functionality at the ROI level and delve into understanding cerebral function at
the voxel level. Further, for using fMRI in a clinical setting, it may also be worthwhile to reliably detect
activation at the voxel level and in single subjects. All these are potential scenarios with low CNRs where
reliable activation detection methods are sorely needed. We return to this application in Section IV.
B. Background and Current Practices
Many thresholding methods [34]–[39] address the issue of multiple testing in determining significance
of test statistics but generally ignore the spatial resolution of fMRI. This shortcoming is often accounted
for by spatially smoothing the fMR images in the pre-processing stage, but such non-adaptive smoothing
reduces both the adaptive spatial resolution and the number of available independent tests for detecting
activation [40]. Iterative adaptive smoothing and segmentation methods in the form of propagationseparation (PS) [40] and adaptive-segmentation (AS) [41] have been developed. While operationally
different, PS and AS are broadly similar, essentially yielding a segmentation of the SPM into activated
and inactivated voxels. PS results (approximately) in a random t-field and uses Random Field Theory to
segment the volume while AS uses multi-scale testing to determine activated segments. [41] argue for
AS, because of its more general development and fewer model assumptions. AS also requires no heuristic
corrections for spatial correlation, and attractively, provides decisions at prescribed significance levels.
[41] demonstrated better performance of AS over PS in an auditory experiment. However, AS requires
pre-specified bandwidth sequences and ignores correlation within the SPM. In this paper, we therefore
develop a fully automated and computationally speedy adaptive smoothing and thresholding approach that
accounts for correlation and obviates the need for setting any parameter beyond the significance level.
The remainder of this paper is organized as follows. In Section II, we provide theoretical and implementation details of our algorithm. We evaluate performance in a large-scale simulation experimental
study in Section III. Both visual and numerical summaries of performance are presented. In Section IV,
we illustrate applicability of our FAST algorithms on the dataset of Section I-A. We end with some
discussion.
II. T HEORY AND M ETHODS
A. Preliminaries
Let Y i be the time series vector of the observed BOLD response at the ith voxel obtained after
preprocessing for registration and motion correction. A common approach relates Y i to the expected
BOLD response via the general linear model
Y i = Xβ i + i ,

(1)

where i is a pth-order auto-regressive (AR) Gaussian error vector with AR coefficients φi=(φi1 , φi2 , . . . , φip )
and marginal variance σi2 . Without loss of generality (w.l.o.g.), assume that the design matrix X has the
intercept in the first column, the expected BOLD response for the k stimulus levels in the next k columns,
and polynomial terms for the drift parameter in the remaining m columns. Thus, β is a coefficient vector
of length d = k + m + 1. We assume that the image volume has n voxels: thus i = 1, 2, . . . , n in (1).

The parameters (β̂ i , σ̂i2 , φ̂i )s are usually estimated via generalized least squares or restricted maximum
likelihood. A typical analysis approach then applies (voxel-wise) hypothesis tests with the null hypothesis
specifying no activation owing to the stimulus or task. SPMs of the form Γ = {c0 β̂ i }i∈V with appropriate
contrasts c0 β i are then formulated at each voxel.
Many researchers use models that assume independence or AR(1) errors: others pre-whiten the time series before fitting (1) under independence. Incorrect model specifications can yield less precise SPMs [42]–
[45]. In this paper, we used Bayes Information Criterion (BIC) [46], [47] which measures a fitted model’s
quality by trading its complexity against its fidelity to the data. Tests on the SPM Γ identify voxels that
are activated related to the experiment. Specifically, voxel i ∈ V is declated to be activated if a suitable
test rejects the hypothesis H0 : c0 β i = 0 i.e., if c0 β̂ i significantly deviates from zero. Our objective in this
paper is to develop an approach that adaptively and automatically smooths and thresholds the SPM while
accounting for the spatial correlation in the SPM and the fact that these sequence of thresholds yield SPMs
from truncated distributions. Before detailing our algorithm, we provide some theoretical development.
B. Theoretical Development
We assume that the SPM is t-distributed with degrees of freedom large enough to be well-approximated
via the standard normal distribution. We also assume a homogeneous correlation structure for the SPM,
which is reasonable given our choice of a Gaussian kernel for smoothing. We have the
Theorem 1. Let X ∼ Nn (0, R) where X = (X1 , . . . , Xn )0 and R is a circulant correlation matrix.
1
1
Writing 1 = (1, 1, . . . , 1)0 , we let ρ = R− 2 1 be the sum of the elements in the first row of R− 2 . Further, let
X(n) be the maximum value of X, that is, X(n) = max {X1 , X2 , . . . , Xn }. Then the cumulative distribution
function (c.d.f.) Fn (x) of X(n) is given by Fn (x) = P (X(n) ≤ x) = [Φ(ρx)]n , where Φ(·) is the c.d.f. of
the standard normal random variable.
Proof: From the definition of X(n) , we have
Fn (x) = P (X ≤ x1)
= P (R−1/2 X ≤ xR−1/2 1)
= P (Y ≤ ρx1), where Y ∼ Nn (0, I n )
= P (Yi ≤ ρx ∀ i = 1, 2, . . . , n) = [Φ(ρx)]n .

(2)

In the limit, we are led to the following
Corollary 2. Let X and X(n) be as in Theorem 1. Then the limiting distribution of X(n) belongs to the
domain of attraction of the Gumbel distribution, and satisfies:
lim [Fn (an x + bn )] = exp{e−x },

n→∞

(3)

where an = ρ/[nφ(bn )] and bn = Φ−1 (1 − 1/n)ρ, with φ(·) being the standard normal probability density
function (p.d.f.).
Proof: For Y ∼ Nn (0, I n ), the limiting distribution of Y(n) satisfies limn→∞ [Φ(an x+bn )]n = exp{e−x },
with an = 1/[nφ(bn )] and bn = Φ−1 (1 − 1/n) [48]. The result follows from Theorem 1.
Theorem 1 and Corollary 2 provide the wherewithal for choosing a threshold using the limiting
distribution of the maximum of correlated normal random variables with circulant correlation structure.
Note however, that the first thresholding step results in truncated (and correlated) random variables that
are under consideration for thresholding in subsequent steps. We account for this added complication by
deriving the limiting distribution of the maximum of a correlated sample from a right-truncated normal

distribution. We do so by first noting that if Y1 , Y2 , . . . Yn are independent identically distributed random
variables from the standard normal p.d.f truncated at η, then each Yi has p.d.f. φ•η (·) and c.d.f. Φ•η (·):
φ(x)
Φ(x)
I(x < η); Φ•η (x) =
I(x < η).
(4)
Φ(η)
Φ(η)
= max {Yi , i = 1, 2, . . . , n} has c.d.f. G•η (x) = [Φ(x)/Φ(η)]n I(x < η) with limiting distribution
φ•η (x) =

Then Y(n)
given by

Theorem 3. Let Y1 , Y2 , . . . , Yn be a sample from (4). Then the limiting distribution of Y(n) belongs to the
domain of attraction of the reverse Weibull distribution and satisfies
lim [G•η (a•n x + b•n )] = exp {−(−x−τ )}I(x ≤ 0).

n→∞

(5)

for some τ > 0. Here a•n = η − Φ•η −1 (1 − 1/n) and b•n = 0.
Proof: Note that η = sup{x | Φ•η (x) < 1}. From Theorem 10.5.2 in [49], a sufficient condition for
Y(n) to be the domain of attraction is to demonstrate that
(η − x)φ• (x)
=τ
lim
x→η 1 − Φ• (x)
for some τ > 0 [50]. In our case, the limit holds because η < ∞. Then, using L’Hôpital’s rule, we have
d •
(η − x) dx
φη (x) − φ•η (x)
(η − x)φ•η (x)
=
lim
lim
x→η
x→η 1 − Φ•
−φ•η (x)
η (x)
(η − x)φ0 (x)/Φ(η) − φ(x)/Φ(η)
= lim
x→η
−φ(x)/Φ(η)
=1
Thus the right-truncated normal distribution satisfies the reverse Weibull condition and converges to the
reverse Weibull distribution with τ ≡ 1 in (5). The constants in the theorem are as per extreme value
theory [48], [49].
Theorem 4. Let X be a random vector from the Nn (0, R) density but that is right-truncated in each
1
coordinate at η, with R being a circulant correlation matrix. Let ρ = R− 2 1 be the sum of the elements
1
in the first row of R− 2 . Also, let X(n) be defined as before. Then the c.d.f. Fη• (x) of X(n) is given by
Fη• (x) = P (X(n) ≤ x) = [Φ(ρx)/Φ(η)]n I(x < η).
Proof: Proceeding as in the proof of Theorem 1 yields
Fη• (x) = P (X ≤ x1)
= P (R−1/2 X ≤ xR−1/2 1)
= P (Y ≤ ρx1), where Yi s are i.i.d. Φ•ρη (y)
= G•ρη (ρx) = [Φ(ρx)/Φ(η)]n I(x < η).

(6)

Corollary 5. Let X and X(n) be as in Theorem 4. Then the limiting distribution of X(n) belongs to the
domain of attraction of the reverse Weibull distribution, and satisfies:
lim [Fη• (a•n x + b•n )] = exp {−(−x−α )}I(x ≤ 0).

n→∞

(7)

where a•n = (η − Φ•ρη −1 (1 − 1/n))ρ and b•n = 0.
Proof: From Theorem 3, the limiting distribution of Ynm) (5). Then Fη• (x) = G•ρη (ρx) from Themrem
4 and the result follows.
Having developed our results, we are now ready to propose our fast adaptive smoothing and thresholding
algorithm.

C. Fast Adaptive Smoothing and Thresholding
Our proposed algorithm adaptively and, in sequence, smooths and identifies activated regions by
thresholding. We estimate the amount of smoothing robustly or from the correlation structure, which we
assume is well-approximated by a 3D Gaussian kernel. Thus, under the null hypothesis (of no activation
anywhere), we assume that the SPM Γ ∼ N (0, R) where R = S h with h the smoothing parameter of S h
−1
given in terms of its full-width-at-half-maximum. Let Γ(−h) ∼ S h 2 Γ. Then estimate h by maximizing
the loglikelihood function
1
1
n
(8)
`(h | Γ(−h) ) = − log(2π) − log |S h | − Γ0(−h) Γ(−h) .
2
2
2
Note that Γ(−h) and |Sh | are speedily computers using Fast Fourier transforms. Starting with the SPM Γ,
obtained as discussed in Section II-A, we propose the algorithm:
1) Initial Setup. At this point, assume that ζi ≡ 0 ∀ i, where ζi is the activation status of the ith voxel.
(0)
That is, assume that all voxels are inactive. Set ζi ≡ ζi . Also denote Γ(0) = Γ, and n0 = n, where
(k)
nk denotes the number of voxels for which ζi = 0.
2) Iterative Steps. For k = 1, 2, . . . , iterate as follows:
a) Smoothing. We smooth Γ(k−1) in one of two ways:
i) Adaptive Maximum Likelihood (AM-FAST, pronounced ăm-fast): Maximize (8) given Γ(k−1)
1
to obtain hk . Obtain Γ(k) by smoothing Γ(k−1) with S hk . Also obtain ρk = S hk − 2 1.
ii) Adaptive Robust Smoothing (AR-FAST, pronounced ahr-fast): Use robust methods [51] on
1
Γ(k−1) to get Γ(k) . Now, maximize (8) given Γ(k) to obtain hk and ρk = S hk − 2 1.
b) Adaptive Thresholding. This consists of two steps:
i) For k = 1, use Corollary 2 to obtain an and bn , otherwise (i.e. for k > 1) use Corollary 5
to obtain a•nk−1 . In both cases, we use the spatial correlation estimated by S hk−1 .
ii) From the Gumbel (for k = 1) or reverse Weibull distributions (for k > 1), get
(
an0 ιGα + bn0 for k = 1
ηk =
(9)
a•nk−1 ιW
otherwise.
α
where ιGα and ιW
α are the upper-tail α-values for the Gumbel and the reverse Weibull (with
τ = 1) distributions, respectively.
(k)
(k−1)
iii) Set all ζi = 1 if ζi
= 0 and if the ith coordinate of Γ(k) exceeds ηk . Also, update
Pn (k)
nk = i=1 ζi .
3) Stopping criterion. Let J(ζ (k) , ζ (k−1) ), be the Jaccard Index [52], [53] of the activation maps in successive iterations (i.e. between the (k − 1)th and kth iterations). If J(ζ (k) , ζ (k−1) ) ≤ J(ζ (k+1) , ζ (k) ),
then the algorithm terminates and ζ (k) is the final activation map.
Comments: A few comments are in order:
a) Two-sided tests: Our development here assumes one-sided tests where (w.l.o.g.) large values are
the extremal values of the SPM. Two-sided tests can be easily accommodated by simply replacing the α
in (9) by α/2 and by using the absolute value of the components of Γ(k) to decide on activation.
b) Comparison with AS: AS [41] also provides an adaptive weighted smoothing approach based on
multi-scale testing. Both AM-FAST and AR-FAST have similarities with AS, in that they also smooth
and threshold iteratively. However, there are a few fundamental differences. For one, the AS approach has
a set user-specified sequence of bandwidths that smooths Γ(k) at each step. In contrast, we use likelihood
maximization (for AM-FAST) or robust (for AR-FAST) methods to optimally determine h at each step.
AS also performs thresholding in a similar manner as us, but they use a general Fréchet extreme value
distribution. Their development ignores both the spatial context and also the fact that thresholding results
in subsequent truncated decisions. Our development represents the procedure more accurately because

we account for both the correlation structure (resulting in the initial cut-off decided as per the Gumbel
distribution) and the truncation (with the subsequent cut-offs being determined by the reverse Weibull
distribution). Finally, our method is completely data-driven, with termination declared only when the
Jaccard index indicates no substantial changes in the detected activation.
III. P ERFORMANCE E VALUATIONS
This section reports performance of the FAST algorithm relative to the AS and also the more popular
cluster-thresholding algorithms. We use a 2D setup to make it computationally possible to evaluate
performance across many settings.
A. Experimental Setup
Our simulation setup was a modified version of the digitized 128×128 2D Hoffman phantom [54] often
used in Positron Emission Tomography. The phantom (Figure 1) has 3465 in-brain pixels, representing
two types (say, A and B) of anatomic structures. 138 of the B-type pixels from two contiguous regions
were deemed to be truly activated. The phantom provided the values (see Figure 1) for β i = (βi0 , βi1 , βi2 )
Region
Background
Brain A
Brain B
Activated

βi0
0
4500
6000
6000

βi1
0
0
0
600

βi2
0
-155.32
-155.32
-155.32

Fig. 1: The phantom used in our simulation experiments. Putative anatomic regions (A and B) are in
shades of grey, while truly activated pixels are in red. The table on the right lists the β i s used in our
simulations.
in (1) as per the location of the ith pixel in the phantom.
As in (1), the design matrix X had the intercept in the first column. For the second column, the input
stimulus time series was set to alternate as 16 on-off blocks of 6 time-points each. The first block signaled
rest, with input stimulus 0, while the second block was 1 to signify activation. These blocks alternated,
yielding values at 96 time-points. These temporal values were convolved with the hemodynamic response
function [7] to provide the second column of X. The third column of X represented linear drift and was
set to t (t = 1, 2, . . . , 96). Also, as per (1), AR(p) errors were simulated for different p and for different
structures. Specifically, for each p, we considered AR coefficients for a range of φ ≡ φi s with coefficients
(φ1 , φ2 , . . . , φp ) that were, with lag, (a) all equal, (b) decreasing, (c) increasing, (d) first decreasing and
then
Pp increasing and (e) first increasing and then decreasing. To ensure stationary solutions, we restricted
j=1 φj = 0.9. Thus, for AR(1), we have φ1 ≡ 0.9 for all cases. For p = 2, 3, 4, we have φi ≡ 0.9/p
for the equal AR coefficients scenario. Table I contains the φ-values for the other cases. Finally, we
chose σ to be such that they corresponded to very low to moderate CNR settings. Specifically, we set
σ0 = 240, 300, 400 yielding CNR = 1.0, 1.5, and 2.0. By design, our SNRs were 10 times our CNRs.
We simulated realizations of time series images using (1) and the setup of Figure 1 and the AR
structures in Table I. For each pixel, we fitted (1) with different AR(p), p = 0, 1, 2, 3, 4, 5, and chose
the model that provided the minimum BIC. SPMs were generated as discussed in Section II-A. Figure 2
provides sample realizations of SPMs obtained using our simulation setup and using AR(4) time series.
We also performed a larger-scale simulation study, where we replicated our experiment 25 times for each
simulation setting. Activation detection was done with AM-FAST and AR-FAST and also the commonlyused (in fMRI) cluster-thresholding (CT) at [39]’s suggested significance level (α = 0.001) with a secondorder neighborhood and number of voxels determined by [55]’s 3dClustSim function. The R package

TABLE I: φs for the AR(p) scenarios used in our simulations.
p
2
3
4
5
p
2
3
4
5

Decreasing
(0.6, 0.3)
(0.4, 0.3, 0.2)
(0.3, 0.25, 0.20, 0.15)
(0.3, 0.25, 0.20, 0.10, 0.05)
Decreasing-Increasing
(0.6,0.3)
(0.4,0.1,0.4)
(0.4,0.05,0.05,0.4)
(0.25,0.15,0.1,0.15,0.25)

Increasing
(0.3, 0.6)
(0.2, 0.3, 0.4)
(0.15, 0.20, 0.25, 0.3)
(0.05, 0.10, 0.20, 0.25, 0.30)
Increasing-Decreasing
(0.3,0.6)
(0.1,0.7,0.1)
(0.05,0.4,0.4,0.05)
(0.1,0.15,0.4,0.15,0.1)

25.17

22.17

19.17

16.17

13.17

10.17

7.17

4.17

1.17

−1.83

−4.83

Fig. 2: SPM for a sample simulation setting, for the five settings (from left to right) of equal, increasing,
decreasing, increasing-decreasing and decreasing-increasing AR coefficients. The three rows (from top to
bottom) correspond to simulation settings with CNR of 1.0, 1.5 and 2.0, respectively.

AnalyzeFMRI [56] was used to perform cluster-wise thresholding. We also performed AS [41] (henceforth,
TP11) using the R package fmri [57]. In each case, performance was evaluated in terms of the Jaccard
Index [52], [53] between the estimated activation map and the truth as per Figure 1.
B. Results
Figure 3 provides activation maps detected using AM-FAST (with α = 0.025), CT and TP11 for each
case of Figure 2. AR-FAST results were essentially indistinguishable from AM-FAST and so are not
displayed. It is encouraging to note that AM-FAST visually tracks very closely to the activation regions
of Figure 1. Both CT and TP11 also identify the true regions as activated but they also identify several
other pixels as activated, with more false positives detected at lower CNR levels. At moderate CNR=2.0,
CT does very well under the equal AR coefficients case but not as well as for some of the other cases.
TP11 however has more false positives even for CNR=2.0. In summary, however, AM-FAST performs
quite creditably at all CNR levels.
Figure 4 displays performance of our algorithms and its competitors in the large-scale simulation study.
For low CNRs, both CT and TP11 perform poorly, but AM-FAST and AR-FAST perform quite well. As
with Figure 3, both CT and TP11 improve with increased CNR: indeed CT is marginally the best performer
for when CNR=2.0 and with no autocorrelation (p = 0). Interestingly, for CNR=2.0, CT does worse than
AM-FAST or AR-FAST for the decreasing AR coefficients case but not necessarily for the other cases.
AM-FAST and AR-FAST perform very similarly at all settings, however AM-FAST is computationally
faster than AR-FAST. The poorer performance of CT relative to the FAST algorithms is not surprising
because 3dClustSim and other such functions are not very adaptive in their execution. TP11 is however,

(a) CNR=1.0

(b) CNR=1.5

(c) CNR=2.0

Fig. 3: Activation maps obtained using AM-FAST (top row), CT (middle) and TP11 (bottom row) from
the SPMs of Figure 2 at each CNR. Columns are as in Figure 2.

somewhat more adaptive, through the choice of the sequence of smoothing parameters that is left to the
user. Specifying this sequence appropriately may, however, require considerable skill, and the default
values provided in the fmri package do not appear to be adequate enough for lower CNR situations. On
the other hand, our approach determines the optimal smoothing from the SPM at each iteration through
likelihood maximization (AM-FAST) or robust methods (AR-FAST). Interestingly though, we have found
that α plays a role in the FAST algorithms with smaller values performing marginally better with higher
CNRs and higher values performing substantially better for lower CNRs. As a via media, we use α = 0.025
in our application. In summary, however, both FAST algorithms perform very well and over a wide range
of low-to-moderate CNR settings.
IV. ACTIVATION D URING P ERCEPTION OF N OISY S PEECH
The dataset is from [31] and also provided by [55] as data6 to illustrate the many features of the
AFNI software. In this experiment the subject was presented with audiovisual speech in both auditory and
visual formats. There were two types of stimuli: (1) audio-reliable where the subject could clearly hear the

CNR = 1.00

CNR = 1.50

CNR = 2.00

1.00

0.75

Equal

0.50

0.25

0.00
1.00

Decreasing

0.75

0.50

0.25

0.00

0.75

Increasing

Jaccard Index

1.00

0.50

0.25

0.00

Increasing-Decreasing

1.00

0.75

0.50

0.25

0.00

Decreasing-Increasing

1.00

0.75

0.50

0.25

0.00
AR(0)

AR(1)

AR(2)

AR(3)

AM-FAST0.001

AR(4)

AR(5)

AR-FAST0.001

AR(0)

AR(1)

AM-FAST0.01

AR(2)

AR(3)

AR-FAST0.01

AR(4)

AR(5)

AM-FAST0.05

AR(0)

AR(1)

AR(2)

AR-FAST0.05

AR(3)

CT

AR(4)

AR(5)

TP11

Fig. 4: Summary of performance of the AM-FAST, AR-FAST, CT and TP11 algorithms for the different
simulation settings. The subscript on the FAST legends represents the value of α = 0.001, 0.01, 0.05 used
in the experiments.

spoken word but the visual of the female speaker was degraded and (2) visual-reliable where the subject
could clearly see the female speaker vocalizing the word but the audio was of reduced quality. There
were three experimental runs, each consisting of a randomized design of 10 blocks, equally divided into
blocks of audio-reliable and visual-reliable stimuli. An echo-planar imaging sequence (TR=2s) was used
to obtain T∗2 -weighted images with volumes of 80 × 80 × 33 (voxel dimensions = 2.75 × 2.75 × 3.0 mm3 )
over 152 time-points. Our goal was to determine activation in areas that responded to the audio task
(H0 : βa = 0), areas that responded to the visual task (H0 : βv = 0) and finally the contrast of visual

(a)

(b)

(c)

Fig. 5: Activation regions in AFNI data6 using SPM based on AR(p̂) of (a) Visual-reliable stimulus (b)
Audio-reliable stimulus and (c) the difference contrast between Visual-reliable and Audio-reliable stimuli.
against audio (H0 : βv − βa = 0). We fitted AR models at each voxel using p = 0, 1, 2, 3, 4, 5 and chose
p minimizing the BIC. After obtaining the SPM, we applied AM-FAST with results displayed via AFNI
and Surface Mapping (SUMA) in Figure 5. Note that most of the activation occurs in Brodmann areas 18
and 19 (BA18 and BA19), which comprise the occipital cortex and the extrastriate (or peristriate) cortex.
In humans with normal sight, extrastriate cortex is a visual association area, where feature-extraction,
shape recognition, attentional, and multimodal integrating functions take place. We also observed activated
regions in superior temporal sulcus (STS). Recent studies have shown increased activation in this area
which are related to: voices versus environmental sounds, stories versus nonsensical speech, moving faces
versus moving objects, biological motion and so on [58]. Although a detailed analysis of the results of
this study is beyond the purview of this paper, we note that AM-FAST finds interpretable results even
when applied to a single subject high-level cognition experiment.
V. D ISCUSSION
In this paper, we have proposed a new fully automated and adaptive smoothing and thresholding
algorithm called FAST that has the ability to perform activation detection in a range of low-signal
regions. Two variants AM-FAST and AR-FAST, have been proposed that show similar performance. We
recommend AM-FAST because it is computationally faster than AR-FAST. Our methodology accounts
for both spatial correlation structure and uses more accurate extreme value theory in its development.
Simulation experiments indicate good performance over a range of low-SNR and low-CNR settings.
Developing FAST for more sophisticated time series and spatial models, as well as increased use of
diagnostics in understanding activation and cognition are important research areas and directions that
would benefit from further attention.
S UPPLEMENTARY M ATERIALS
The datasets and R code used in this paper and minimal examples for using the code are publicly
available at https://github.com/maitra/FAST-fMRI.
R EFERENCES
[1] J. W. Belliveau, D. N. Kennedy, R. C. McKinstry, B. R. Buchbinder, R. M. Weisskoff, M. S. Cohen, J. M. Vevea, T. J. Brady, and
B. R. Rosen, “Functional mapping of the human visual cortex by magnetic resonance imaging,” Science, vol. 254, pp. 716–719, 1991.
[2] K. K. Kwong, J. W. Belliveau, D. A. Chesler, I. E. Goldberg, R. M. Weisskoff, B. P. Poncelet, D. N. Kennedy, B. E. Hoppel, M. S.
Cohen, R. Turner, H.-M. Cheng, T. J. Brady, and B. R. Rosen, “Dynamic magnetic resonance imaging of human brain activity during
primary sensory stimulation,” Proceedings of the National Academy of Sciences of the United States of America, vol. 89, pp. 5675–5679,
1992.
[3] P. A. Bandettini, A. Jesmanowicz, E. C. Wong, and J. S. Hyde, “Processing strategies for time-course data sets in functional MRI of
the human brain,” Magnetic Resonance in Medicine, vol. 30, pp. 161–173, 1993.

[4] K. J. Friston, A. P. Holmes, K. J. Worsley, J.-B. Poline, C. D. Frith, and R. S. J. Frackowiak, “Statistical parametric maps in functional
imaging: A general linear approach,” Human Brain Mapping, vol. 2, pp. 189–210, 1995.
[5] A. M. Howseman and R. W. Bowtell, “Functional magnetic resonance imaging: imaging techniques and contrast mechanisms,”
Philosophical Transactional of the Royal Society, London, vol. 354, pp. 1179–94, 1999.
[6] W. D. Penny, K. J. Friston, J. T. Ashburner, S. J. Kiebel, and T. E. Nichols, Eds., Statistical Parametric Mapping: The Analysis of
Functional Brain Images, 1st ed. Academic Press, 2006.
[7] M. A. Lindquist, “The statistical analysis of fMRI data,” Statistical Science, vol. 23, no. 4, pp. 439–464, 2008.
[8] N. A. Lazar, The Statistical Analysis of Functional MRI Data. Springer, 2008.
[9] F. G. Ashby, Statistical Analysis of fMRI Data. MIT Press, 2011.
[10] S. Ogawa, T. M. Lee, A. S. Nayak, and P. Glynn, “Oxygenation-sensitive contrast in magnetic resonance image of rodent brain at high
magnetic fields,” Magnetic Resonance in Medicine, vol. 14, pp. 68–78, 1990.
[11] S. Ogawa, T. M. Lee, A. R. Kay, and D. W. Tank, “Brain magnetic resonance imaging with contrast dependent on blood oxygenation,”
Proceedings of the National Academy of Sciences, USA, vol. 87, no. 24, pp. 9868–9872, 1990.
[12] K. J. Worsley, C. H. Liao, J. Aston, V. Petre, G. H. Duncan, F. Morales, and A. C. Evans, “A general statistical analysis for fmri data,”
NeuroImage, vol. 15, pp. 1–15, 2002.
[13] K. J. Friston, P. Fletcher, O. Josephs, A. Holmes, M. Rugg, and R. Turner, “Event-related fMRI: characterizing differential responses,”
Neuroimage, vol. 7, no. 1, pp. 30–40, 1998.
[14] G. H. Glover, “Deconvolution of impulse response in event-related bold fmri,” NeuroImage, vol. 9, pp. 416–429, 1999.
[15] R. B. Buxton, K. Uludağ, D. J. Dubowitz, and T. T. Liu, “Modeling the hemodynamic response to brain activation,” Neuroimage,
vol. 23, pp. S220–S233, 2004.
[16] K. J. Friston, C. D. Frith, P. F. Liddle, R. J. Dolan, A. A. Lammertsma, and R. S. J. Frackowiak, “The relationship between global
and local changes in PET scans,” Journal of Cerebral Blood Flow and Metabolism, vol. 10, pp. 458–466, 1990.
[17] K. J. Friston, A. P. Holmes, K. J. Worsley, J.-P. Poline, C. D. Frith, and R. S. Frackowiak, “Statistical parametric maps in functional
imaging: a general linear approach,” Human brain mapping, vol. 2, no. 4, pp. 189–210, 1994.
[18] K. J. Worsley and K. J. Friston, “Analysis of fMRI time-series revisitedagain,” Neuroimage, vol. 2, no. 3, pp. 173–181, 1995.
[19] C. R. Genovese, N. A. Lazar, and T. Nichols, “Thresholding of statistical maps in functional neuroimaging using the false discovery
rate,” Neuroimage, vol. 15, no. 4, pp. 870–878, 2002.
[20] J. V. Hajnal, R. Myers, A. Oatridge, J. E. Schweiso, J. R. Young, and G. M. Bydder, “Artifacts due to stimulus-correlated motion in
functional imaging of the brain,” Magnetic Resonance in Medicine, vol. 31, pp. 283–291, 1994.
[21] B. Biswal, A. E. DeYoe, and J. S. Hyde, “Reduction of physiological fluctuations in fMRI using digital filters.” Magnetic Resonance
in Medicine, vol. 35, no. 1, pp. 107–113, January 1996. [Online]. Available: http://view.ncbi.nlm.nih.gov/pubmed/8771028
[22] R. P. Wood, S. T. Grafton, J. D. G. Watson, N. L. Sicotte, and J. C. Mazziotta, “Automated image registration. ii. intersubject validation
of linear and non-linear models,” Journal of Computed Assisted Tomography, vol. 22, pp. 253–265, 1998.
[23] R. P. Gullapalli, R. Maitra, S. Roys, G. Smith, G. Alon, and J. Greenspan, “Reliability estimation of grouped functional imaging data
using penalized maximum likelihood,” Magnetic Resonance in Medicine, vol. 53, pp. 1126–1134, 2005.
[24] R. Maitra, S. R. Roys, and R. P. Gullapalli, “Test-retest reliability estimation of functional mri data,” Magnetic Resonance in Medicine,
vol. 48, pp. 62–70, 2002.
[25] R. Maitra, “Initializing partition-optimization algorithms,” IEEE/ACM Transactions on Computational Biology and Bioinformatics,
vol. 6, pp. 144–157, 2009. [Online]. Available: http://doi.ieeecomputersociety.org/10.1109/TCBB.2007.70244
[26] Z. S. Saad, D. R. Glen, G. Chen, M. S. Beauchamp, R. Desai, and R. W. Cox, “A new method for improving functional-to-structural
mri alignment using local pearson correlation,” NeuroImage, vol. 44, pp. 839–848, 2009.
[27] E. E. Chen and S. L. Small, “Test-retest reliability in fMRI of language: Group and task effects,” Brain and Language, vol. 102,
no. 2, pp. 176–85, 2007. [Online]. Available: http://dx.doi.org/10.1016/j.bandl.2006.04.015
[28] K. D. Kryter, The handbook of hearing and the effects of noise: Physiology, psychology, and public health. Academic Press, 1994.
[29] M. D. Hauser, The evolution of communication. MIT press, 1996.
[30] S. Dupont and J. Luettin, “Audio-visual speech modeling for continuous speech recognition,” IEEE transactions on multimedia, vol. 2,
no. 3, pp. 141–151, 2000.
[31] A. R. Nath and M. S. Beauchamp, “Dynamic changes in superior temporal sulcus connectivity during perception of noisy audiovisual
speech,” The Journal of Neuroscience, vol. 31, no. 5, p. 1704 1714, 2011.
[32] W. H. Sumby and I. Pollack, “Visual contribution to speech intelligibility in noise,” The journal of the acoustical society of america,
vol. 26, no. 2, pp. 212–215, 1954.
[33] B. E. Stein and M. A. Meredith, The merging of the senses. The MIT Press, 1993.
[34] Y. Benjamini, A. M. Krieger, and D. Yekutieli, “Adaptive linear step-up procedures that control the false discovery rate,” Biometrika,
vol. 93, no. 3, pp. 491–507, 2006.
[35] R. Heller, D. Stanley, D. Yekutieli, N. Rubin, and Y. Benjamini, “Cluster-based analysis of fMRI data,” NeuroImage, vol. 33, no. 2,
pp. 599–608, Nov. 2006. [Online]. Available: http://dx.doi.org/10.1016/j.neuroimage.2006.04.233
[36] Y. Benjamini and R. Heller, “False discovery rates for spatial signals,” Journal of the American Statistical Association, vol. 102, no.
480, pp. 1272–1281, 2007.
[37] M. Smith and L. Fahrmeir, “Spatial Bayesian variable selection with application to functional Magnetic Resonance
Imaging,” Journal of the American Statistical Association, vol. 102, no. 478, pp. 417–431, 2007. [Online]. Available:
http://pubs.amstat.org/doi/abs/10.1198/016214506000001031
[38] S. M. Smith and T. E. Nichols, “Threshold-free cluster enhancement: Addressing problems of smoothing, threshold dependence and
localisation in cluster inference,” Neuroimage, vol. 44, pp. 83–98, 2009.
[39] C.-W. Woo, A. Krishnan, and T. D. Wager, “Cluster-extent based thresholding in fMRI analyses: Pitfalls and recommendations,”
Neuroimage, vol. 91, p. 412419, 2014.

[40] K. Tabelow, J. Polzehl, H. U. Voss, and V. Spokoiny, “Analyzing fMRIexperiments with structural adaptive smoothing procedures,”
NeuroImage, vol. 33, no. 1, pp. 55–62, 2006.
[41] J. Polzehl, H. U. Voss, and K. Tabelow, “Structural adaptive segmentation for statistical parametric mapping,” NeuroImage, vol. 52,
no. 2, pp. 515–523, 2010.
[42] M. M. Monti, “Statistical analysis of fMRI time-series: A critical review of the glm approach,” Frontiers in Human Neuroscience,
vol. 5, no. 00028, pp. 1–13, 2011. [Online]. Available: http://www.frontiersin.org/Journal/Abstract.aspx?s=537&name=human
neuroscience&ART DOI=10.3389/fnhum.2011.00028
[43] W.-L. Luo and T. E. Nichols, “Diagnosis and exploration of massively univariate neuroimaging models,” Neuroimage, vol. 19, no. 3,
pp. 1014–1032, 2003.
[44] J. M. Loh, M. A. Lindquist, and T. D. Wager, “Residual analysis for detecting mis-modeling in fMRI,” Statistica Sinica, pp. 1421–1448,
2008.
[45] M. Lindquist, J. Loh, L. Atlas, and T. Wager, “Modeling the hemodynamic response function in fMRI: Efficiency, bias and mismodeling,” Neuroimage, vol. 45, no. 1, p. S187S196, 2009.
[46] G. Schwarz, “Estimating the dimensions of a model,” Annals of Statistics, vol. 6, pp. 461–464, 1978.
[47] R. H. Shumway and D. S. Stoffer, Time Series Analysis and Its Applications, second edition ed. Springer, 2006.
[48] S. I. Resnick, Extreme values, regular variation and point processes. Springer, 2013.
[49] H. A. David and H. N. Nagaraja, Order Statistics. Hoboken, New Jersey: John Wiley and Sons, Inc., 2003.
[50] R. Von Mises, “La distribution de la plus grande de n valeurs,” Rev. math. Union interbalcanique, vol. 1, no. 1, 1936.
[51] D. Garcia, “Robust smoothing of gridded data in one and higher dimensions with missing values,” Computational statistics & data
analysis, vol. 54, no. 4, pp. 1167–1178, 2010.
[52] P. Jaccard, “Ètude comparative de la distribution florale dans une portion des alpes et des jura,” Bulletin del la Sociètè Vaudoise des
Sciences Naturelles, vol. 37, p. 547579, 1901.
[53] R. Maitra, “A re-defined and generalized percent-overlap-of-activation measure for studies of fMRI reproducibility and its use in
identifying outlier activation maps,” Neuroimage, vol. 50, no. 1, pp. 124–135, 2010.
[54] E. J. Hoffman, P. D. Cutler, W. M. Digby, and J. C. Mazziotta, “3-d phantom to simulate cerebral blood flow and metabolic images
for pet,” IEEE Transactions on Nuclear Science, vol. 37, pp. 616–620, 1990.
[55] R. W. Cox, “Afni: software for analysis and visualization of functional magnetic resonance neuroimages,” Computers and Biomedical
research, vol. 29, no. 3, pp. 162–173, 1996.
[56] C. Bordier, M. Dojat, and P. L. de Micheaux, “Temporal and spatial independent component analysis for fmri data sets
embedded in the AnalyzeFMRI R package,” Journal of Statistical Software, vol. 44, no. 9, pp. 1–24, 2011. [Online]. Available:
http://www.jstatsoft.org/v44/i09/
[57] K. Tabelow and J. Polzehl, “Statistical parametric maps for functional mri experiments in R: The package fmri,” Journal of Statistical
Software, vol. 44, no. 11, pp. 1–21, 2011. [Online]. Available: http://www.jstatsoft.org/v44/i11/
[58] E. Grossman and R. Blake, “Brain activity evoked by inverted and imagined biological motion,” Vision research, vol. 41, no. 10, pp.
1475–1482, 2001.

