Smoothing in SPM: The Do's, Dont's, and Maybes

Different sized smoothing kernels applied to a functional dataset. Note that larger smoothing kernels cause a loss of spatial resolution by turning the relatively high resolution, jagged-edged dataset in the upper left, into the soft, puffy, amorphous cotton ball in the lower right.

Smoothing is one of the most straightforward processing steps, simply involving the application of spatial filtering to your data. Signal is averaged over a range of nearby voxels in order to produce a new estimate of the signal at each voxel, and the range can be narrowed or extended to whatever range suits the researcher's delectation. It is rare for this step to fail, as it is not contingent on overlapping modalities; nor is it susceptible to typical neuroimaging landmines such as entrapment in local minima. Furthermore, the benefits are several: True signal tends to be amplified while noise is canceled out, and power is therefore increased. As a result, often this step is thrown in almost as an afterthought, the defaults left flicked into the "On" position, and quickly forgotten about, as the researcher scampers out of the lab and into his Prius for a quick connection before dinner.

However, smoothing can also be deceptively treacherous. For those researchers intending to tease apart discrete cortical or subcortical regions - for example, the amygdala, if you're into that kind of thing - will find that smoothing tends to smear signal across a wide area, leading to a reduction in spatial specificity. Furthermore, ridiculously large smoothing kernels can actually lead to lower t-values in peak voxels. This may appear to be counterintuitive at first; however, note that increasing the range of voxels can begin to recruit voxels which have nothing to do with the signal you are looking at, and can even begin to average signal from voxels which have an opposite deflection to the signal you are interested in.

Effect of smoothing kernels on statistical results. Here, a contrast of left-right was performed on datasets smoothed with a 4mm kernel. Note that as the smoothing kernel increases, the peak t-value decreases, as depicted by the thermometer bar.

8mm kernel

15mm kernel

For example, let's say we are interested in the contrast of left button presses minus right button presses, as pictured above; as we increased the smoothing kernel, more and more voxels become part of the big blog - I mean, blob! - and it appears that our power increases as well. However, as we extend our averaging over a wider expanse over the fields and prairies of voxels, we risk beginning to smooth in signal from white matter and increasingly unrelated areas. At the most extreme, one can imagine smoothing in signal from the opposite motor cortex, which, for this contrast, will have strongly negative beta estimates.

Your sensitive FMRI antennae should also be attuned to the fact that smoothing can be applied at different magnitudes in the x-, y-, and z-directions. For example, if you are particularly twisted, you could smooth eight millimeters in the y- and z-directions, but only six millimeters in the x-direction. This also comes into play when estimating the smoothness of a first- or second-level analysis, as the smoothing extent may differ along all three coordinates.

For more details, along with a sample of my writing style as a younger man, see the following posts:

Group Level Smoothness Estimation in SPM

Smoothing in AFNI

Group-Level Smoothness Estimation in SPM

Fathers and teachers, I ponder: What is hell? I maintain that it is the suffering of being unable to love. Oh, and it also might be the fact that you can't get an accurate smoothness estimate on second-level results. I mean really, what's up with that?

-Father Zossima

Last week I received an email from a fellow neuroimager asking about smoothness estimates for second-level results. As I discussed in a previous post, this is an important question when calculating cluster correction thresholds, as the smoothing kernel applied to FMRI data is not the same smoothness that should go into cluster correction estimates; and failing to account for this will lead to shame and ridicule. (The exception to this would be using AFNI's 3dBlurInMask, but I assume that most people use either AFNI's 3dmerge or SPM's spm_smooth function.)

To highlight this distinction, imagine a cake - the cake itself is the smoothness on the images that come out of the scanner. Now, pipe some frosting onto the cake. That frosting is the extra smoothing applied on top of those images. The baker piping the frosting is the analysis package you are using, and you can either tell him to add more frosting or less. Sometimes he will malfunction and give you an error message, and you don't understand why he is not working, but you cannot fire him, as you have no other alternative. In any case, when someone asks you how many calories are in the cake, you do not lie and tell them only the calories that are in the frosting; you must account for the underlying cake as well. And then there are times where you will dip your finger into the frosting, just to get a little taste, but your pudgy fingers leave a noticeable divot; and you attempt to fill in the gap with frosting from the edges, but it only ends up looking worse.

In sum, if you are using SPM to estimate smoothness for second-level contrasts or ANOVAs, do not simply use spm_est_smoothness on the ResMS.hdr images output from that analysis. I have tried and I have failed to get any meaningful results, and I have also been unable to find an answer on the message boards or listservs. Therefore, instead of focusing on the second-level analyses, I recommend averaging the smoothness estimates across each subject, and using those averages for your cluster correction estimates. This can be done by looping over each subject, applying spm_est_smoothness to each individual's ResMS.hdr file, and storing the results in a matrix, after which the results are averaged in the x-, y-, and z-directions. The following .m file allows you to do this.

I confess that I know little else about smoothing; this was the only reasonable approach that I could find, and I have little intuition about why spm_est_smoothness gives such wonky results on second-level residuals datasets. It may be that these are to remain mysteries, too profound for our common clay.

Thanks to Jillian Hardee, who never smoothed a dataset she didn't like.

FSL Tutorial 2: FEAT (Part 1)

A new tutorial about FEAT is now up; depending on how long it takes to get through all of the different tabs in the interface, this may be a three-part series. In any case, this will serve as a basic overview of the preprocessing steps of FEAT, most of which can be left as a default.

The next couple of tutorials will cover the set up of models and timing files within FSL, which can be a little tricky. For those of you who have stuck with it from the beginning (and I have heard that there are a few of you out there: Hello), there will be some more useful features coming up, aside from reviewing the basics.

Eventually we will get around to batch scripting FEAT analyses, which can save you several hours of mindless pointing and clicking, and leave you with plenty of time to watch Starcraft 2 replays, or Breaking Bad, or whatever it is kids watch these days.

Smoothing in AFNI

What Smoothing Is

Smoothing, the process of introducing spatial correlation into your data by averaging across nearby data points (in fMRI datasets, nearby voxels) is an essential step in any preprocessing stream. Traditionally, smoothing applies a Gaussian filter to average signal over nearby voxels in order to summate any signal that is present, while attenuating and cancelling out noise. This step also has the advantage of spreading signal over brain areas that have some variation across subjects, leading to increased signal-to-noise ratio over larger cortical areas at the loss of some spatial specificity. Lastly, smoothing mitigates the fiendish multiple comparisons problem, since voxels are not truly independent in the first place (neither are timepoints, but more on this at a later time).

Before (left) and after (right) smoothing. Note that there is much more spatial detail in the image before smoothing is applied.

Usually this smoothing step is taken care of by your fMRI software analysis package's built-in smoothing program, but something I learned back at the AFNI bootcamp is that images coming off the scanner are already smoothed, i.e., there is already some inherent spatial correlation in the data. This makes sense, given that functional imaging datasets are acquired at relatively low resolution, and that voxels with a given intensity will most likely be immediately surrounded by other voxels with similar intensities. In the above image on the left which has not been smoothed, voxels that are brighter tend to cluster together, while voxels that are darker tend to cluster together. This spatial correlation will be increased with any additional step requiring spatial interpolation, including coregistration (aligning all volumes in a run to a reference volume), normalization (warping an image to a template space), and smoothing.

Smoothing Programs in AFNI: 3dmerge and 3dBlurToFWHM

In the good old days, the staple program for blurring in AFNI was 3dmerge with the -blur option. Since 2006, there has been a new program, 3dBlurToFWHM, which estimates the smoothness of the dataset, and attempts to smooth only to specified blur amounts in the x, y, and z directions. (Another thing I didn't recognize, in my callow youth, was that there can indeed be different smoothness amounts along the axes of your dataset, which in turn should be accounted for by your smoothing algorithm.) 3dmerge will add on the specified smoothness to whatever smoothness is already present in the data, so the actual smoothness of the data is not what is specified on the command line with 3dmerge. For example, if 3dmerge with a 6mm Full-Width at Half-Maximum (FWHM) is applied to a dataset that already has smoothness of about 3mm in the x, y, and z directions, the resulting smoothness will be about 8-9mm. (This example, similar to many simulations carried out in the lab, is not a purely linear addition, but is relatively close.)


3dBlurToFWHM, on the other hand, will begin with estimating the smoothness of your dataset, and then iteratively blur (i.e., smooth) the data until it reaches the desired level of smoothness in the x, y, and z directions. These can be specified separately, but often it is more useful to use the same value for all three. For example, the following command will blur a dataset until it reaches approximately 6mm of smoothness in each direction:

3dBlurToFWHM -FWHM 6 -automask -prefix outputDataset -input inputDataset

Applying this to a dataset that has been coregistered and normalized will produce output that looks like this:

Note that the first set of numbers highlighted in red are the estimates produced by 3dBlurToFWHM, and the last set of numbers of red are estimates from a separate program, 3dFWHMx. The reason for the discrepancy is that 3dBlurToFWHM uses a slightly faster and less accurate estimator when blurring; however, the two sets of numbers are functionally equivalent. In any case, the resulting blur estimates are very close to the specified FWHM kernel, typically within 10% or less.

What Mask to Use?

Often using a mask generated after the coregistration step (i.e., the 3dvolreg command) is sufficient, such as full_mask+tlrc or mask_epi_extents+tlrc. The -automask option can also be used, which usually works just as well. The only considerations to take into account are that masking out non-brain voxels will limit your field of view; some experienced users advise not masking until after performing statistical tests, as artifacts or consistent signal outside of the brain can point to hardware problems or issues with your experimental design.

What to Use for the -blurmaster Option?

According to the AFNI help for 3dBlurToFWHM, the -blurmaster option can be used to control the smoothing process, which can lead to slightly better smoothing results. As the FWHM is supposed to represent the underlying spatial structure of the noise, and not the spatial structure of the brain's anatomy, using a dataset of residuals is often a good choice here; for example, the output of -errts in 3dDeconvolve. However, this requires running 3dDeconvolve first, and then using the errts dataset in the blurmaster option all over again. This is usually more trouble than it's worth, as I haven't found significant differences between having a blurmaster dataset and not having one.

Estimating Smoothness with 3dFWHMx

After blurring your data, 3dFWHMx can be used to estimate the smoothness of a given dataset. Useful options are -geom, to estimate the geometric mean of the smoothness, and -detrend, which accounts for outlier voxels. For example,

3dFWHMx -geom -detrend -input inputDataset

What about this 3dBlurInMask Thing?

3dBlurInMask does much the same thing as 3dBlurToFWHM, except that you can use multiple masks if you desire. This procedure may also be useful if you plan to only blur within, say, the gray matter of a dataset, although this can also lead to problems such as biasing signal in certain locations after the grey matter masks have all been warped to a template space, as there is far greater variability in smoothing across masks than there is applying a smoothing kernel to the same averaged mask that incorporates all of the brain voxels across all subjects.

How to Estimate Smoothness in SPM

SPM comes with it's own smoothness estimation program, spm_est_smoothness. It requires a residuals dataset and a mask for which voxels to analyze. For example,

[fwhm, res] = spm_est_smoothness('ResMS.hdr', 'mask.hdr')

Will return the estimated FWHM in the x, y, and z directions, as well as a structure containing information about the resels of the dataset.

Why it Matters

I tuned into this problem a couple months ago when one of my colleagues got nailed by a reviewer for potentially not estimating the smoothness in his dataset appropriately. Evidently, he had used the smoothing kernel specified in his preprocessing stream when calculating cluster correction thresholds, when in fact the actual smoothness of the data was higher than that. It turned out to be a non-significant difference, and his results held in either case. However, to be on the safe side, it always helps to report estimating your smoothness with an appropriate tool. If your results do significantly change as a result of your smoothness estimation approach, then your result may not have been as robust to begin with, and may require some more tweaking.

(By tweaking, of course, I mean screwing around with your data until you get the result that you want.)