Showing posts with label FSL. Show all posts
Showing posts with label FSL. Show all posts

Tuesday, February 17, 2015

Notes on Resting-state fMRI Analyses (Part II)

Fieldmap Correction and Coregistration
What is fieldmap correction?
Apart from having a poor resolution, functional images acquired using common EPI sequences suffer distortion due to magnetic field (B0) inhomogeneity introduced by different tissue types in our heads. Such a thing occurs due to the existence of non-homogeneity in RF receive and transmit of the head coils. The more channels you have, the more inhomogeneity the image may have. The most severe inhomogeneity includes the air-bone or air-brain tissue interfaces in the sinuses in the inferior frontal gyrus and medial temporal lobes. This poses a serious effect on our data, a geometrical distortion and signal loss as depicted in Fig-1.

Magnetic field inhomogeneity can be measured with fieldmap images; which can give us a geometric distortion and signal loss. These values can then be used to compensate for the loss by geometrically unwarping the EPI images, and applying cost-function masking in registrations to ignore areas of signal loss. The correction is most useful during image co-registration as it dramatically improves the registration accuracy. Areas where signal loss has occurred unfortunately cannot be restored with any form of post-processing. In other words, it is impossible to recover time-series data in those locations.

There is no separate sequence for acquiring the fieldmap and different scanners give different images. The sequence can be EPI, Spin-echo, or Gradient-echo sequences, but it isn't recommended to use the EPI-based sequence since it will suffer the same problem. There exist 2 different methods of acquiring fieldmap images for the purpose of correction. 

When you do the fieldmap acquisition, you usually acquire two different images: a pair of magnitude images captured with different echo times, and a phase difference image (Fig-1). The acquisition can also be controlled either in the AP (j+) or PA (j-) direction. These images should be acquired in the same orientation as the target EPIs. The phase difference between the two images is proportional to the difference in echo time (ΔTE) and the B0 inhomogeneity observed. The fieldmap is calculated by taking the difference between the two-phase images, and dividing that by the echo time difference.

Method-2 is called the blip-up blip-down method, which calculates the fieldmap based on the difference in distortion between the two consecutive acquisitions. This method acquires two diffusion-weighted images (DWI) with opposite phase encoding directions, that is, the AP and AP directions. It is assumed that there is no change in the magnetic field and sudden motion during the two acquisitions. You can use TOPUP in FSL to help you do fieldmap processing using this method.

Fig-1: Images obtained from the scanner (left) are converted to get a fieldmap image (right). Red circles show distorted regions that require correction. This is a standard procedure of double gradient-echo performed in Siemens 3T scanner.


How to process this in FSL?

At the MNI, our brain imaging center uses Siemens 3T scanner, which is a good thing as FSL provides a ready-to-use tool, fsl_prepare_fieldmap, to obtain a fieldmap phase image in rad/sec. The magnitude image resembles a lower resolution version of the T1 structural (anatomical) image. FSL FUGUE, which is incorporated in FEAT, helps us to do distortion correction using this method. Both the complete and skull-stripped versions of the magnitude image and the processed phase image (rad/sec) should be defined in FEAT. FSL will then attempt to unwarp the distorted EPI image before mapping it to the structural image. The unwarp direction has to be specified and is typically given by the scanner operator depending on how the fieldmap acquisition is set.

In FSL, fieldmap correction is incorporated as part of the registration (preprocessing) pipeline. The highly accurate functional-to-structural coregistration is also called boundary-based registration or BBR  (Greve and Fischl, 2009). The method is based on changes in the intensity along the white matter boundaries instead of the less reliable grey matter boundaries. This means that an accurate segmentation of the structural image is required and bias-field correction reliably improves the accuracy. Performing BBR registration without a fieldmap correction doesn't give many benefits than the usual 6DOF method with FLIRT (Fig 2-3).

Also, there must be some grey-white intensity contrast in the EPI, though it doesn't have to be good enough for segmentation. The FSL website said since only intensities near the white-matter boundary are used by BBR, it is likely to be more robust to a range of pathologies and artefacts in the EPI or the structural.
Fig-2: Comparison of 3 situations using BBR coregistration with fieldmap correction: when there full magnitude image with the skull wasn't supplied to FSL(left); when the correct magnitude image was used but the unwarping direction was the opposite (middle); the correct BBR-registered image with a superior accuracy (right).
Fig-3: Performing coregistration of a functional image to a structural image using BBR is superior than the usual linear 6DOF registration in FSL. Note that the asterisk ( * ) sign indicates the region with severe signal loss. Without using the fieldmap correction, the corpus callosum mapping becomes inaccurate as denoted by a hex sign (#).

 




Slice Timing Correction?
Scientists more or less agree that the slice timing correction is important.  For the more recent multiband sequence, some experts said that slice timing misalignment may not have a huge impact on the analysis. In the earlier version of the Siemens WIP, I was told that the slice timing information contained in the DICOM files was not correct. This can be retrieved easily with a Matlab function. As a result, I didn't perform this correction in my fMRI paper (MB3, TR=1690 msec).

Until recently, one can deduce the slice timing information based on the CMRR Multiband protocol here. For comparison, I have included the effect of slice timing correction to my resting state data with an MB 3x acceleration measured on a single voxel.
Fig-4: Time series with and without the slice timing correction measured on a single voxel @ MNI coordinate (67,41,49).

Adding additional EVs to GLM
Additional regressors (EVs) can be added to the GLM in FSL FEAT. I find that the GUI is a bit tricky, better write a script for that. First, we have to recreate the design matrix by adding the extra regressors, using either:
       ⁍ Pointing to a file for each regressor by constructing a full model design    
       ⁍ Creating a space-delimited text-file comprising all confound EVs
Fig-5: If you click the "Full model setup", a new GUI will appear as shown on the left. Select an appropriate setting (number 1-3). You don't have to perform another temporal filtering. The temporal derivative is optional too. Another way is to construct a text file and select "Add additional confound EVs" (number 4). 

Refer to the Fig-5 above. If you click "Full model setup", a new GUI will appear. Choose the input file as 1-entry per volume (see number 1), no need to convolve it with the HRF anymore (number 2). Also, you do not have to perform another temporal filtering if the regressors are derived from the prefiltered data (number 3). This method is time-consuming. A better option is to use the additional confound EVs (number 4), just that the text file has to be space-delimited, not a comma-separated file! Using a wrong delimiter will cause FEAT to ignore these additional regressors! 

It is also safer to write a code or script rather than getting restricted with the GUI features. To create a design matrix, use the feat_model command. To manually perform GLM according to the design matrix with prewhitening, use the film_gls command. The command is equipped with sophisticated estimations of autocorrelation and Tukey tapering, which is very important for making statistical inferences in task-based fMRI.

One last note about the design matrix is about Orthogonalization. Most EVs generally are almost orthogonal, so enforcing orthogonality is not gonna help much in the results.

Denoising nuisance components with ICA
As mentioned often, rs-fMRI has one major drawback: the data is recorded at rest so it is prone to noise or artefacts. This is so because we don't have any reference pattern as we do when we perform task-based functional imaging. Hence. proper cleanup is paramount to getting a correct deduction or conclusion. This is even more serious for me who is doing learning-related rs-fMRI. One of the things I'm struggling with is choosing the best cleaning method! Technically for my project, I plan to try out the ICA denoising method since the work of our previous postdoc used a different technique.

As mentioned in the previous post, the subject-level data cleaning steps cover the following:
  1. First, you perform ICA, e.g. using FSL MELODIC and identify the nuisance components. From my experience, the tool performs a pretty good job. Let the algorithm choose the best number of ICs.
  2. Identify the nuisance components and run fsl_regfilt script in FSL, producing a so-called clean or denoised fMRI dataset. 
  3. Now, to assess its performance after cleanup, you can either conduct another round of ICA on the residual image or compute the temporal standard deviation of the GM regions (or a specific ROI in the motor cortex).


The script above effectively removes nuisance components identified from the original dataset. If I conducted another ICA on the denoised data, the resulting components are much cleaner!

How does the script differ from the usual GLM function in FEAT? They are not quite the same. The core function of the GLM in FEAT is film_gls. The fsl_regfilt, on the other hand, uses a simple time-varying GLM function called fsl_glm, which does not perform any sophisticated modelling of temporal autocorrelation and whitening. I finally found this difference after struggling for so long!    
Fig-6: Demeaned time-series of the same voxel obtained from FEAT and ICA denoising tool.

The output of regression is called the residual image, res4d.nii.gz, which is supposed to be as clean as the denoised image, but having the mean removed (just to add back). See Fig-6 for representative time-series of a voxel from the residual outputs of each FEAT and fsl_regfilt. There are some minor differences. I think pre-whitening is not necessary because the nuisance EVs are all spatially independent from the MELODIC. 
On the other hand, I finally found out that the residual output of the fsl_regfilt is in fact the same as the one produced by the fsl_glm. The data has to be demeaned, and the design matrix regressors des_norm has to be normalized into unit variance. The FSL gurus in their forum claimed that both methods use the same GLM methods basically, but I'm not sure which GLM it was. Knowing this similarity is essential because now I can compare different methods of denoising, e.g. using WM/CSF average time-series.
To regress out nuisance components following either ICA or other noise modelling tool (e.g. RETROICOR), you can just use the simple a GLM function. Sophisticated estimation of temporal autocorrelation becomes important when you want to do statistical inference of neural activity or connectivity.



Saturday, October 18, 2014

Basic Independent Component Analysis (ICA)

What is ICA?
As an exploratory method, Independent Component Analysis (ICA) provides an alternative to seed-based resting-state fMRI analysis. Traditionally, ICA attempts to solve the cocktail party problem. The story is this. Suppose I have five people with a microphone talking at the same time. How can I separate my speaker output into voices of the first person, second, third, and so on? Indeed, ICA is a powerful method to discover hidden independent patterns or features from a set of data, thus, an exploratory analysis. Unlike GLM, it is model-free because no assumption is made on the shape/pattern of the actual BOLD response.

There are two types of ICA, Temporal and Spatial ICA. The cocktail party problem is an example of Temporal ICA. The data to be separated contains temporal information, the mixing coefficients or weights of which vary uniquely across different locations. In contrast, the dataset in the Spatial ICA contains spatial maps or networks that get activated altogether at a given time. The weights vary uniquely across time for different networks. Spatial ICA is more popular for fMRI data analysis simply because there are more voxels than time points. An image of 64 x 64 x 64 acquired for 250 volumes or TRs contains 262,144 voxels but 250 time points.
Running ICA on the resting-state data will ideally yield to a set of ICs, some of which are clearly related to activation or network, physiological processes (heart rate, breathing), or even some artefacts (e.g. motion, ghosting, slice dropout, noise, etc). Since rs-fMRI is not task-based, we don't have prior knowledge of the temporal waveform of our resting-state networks. This is where data exploration is useful. Spatial ICA attempts to split the data into a set of spatial maps, each with an associated time course. No assumption is made, no experimental paradigm to be specified. The diagram from the FSL website below summarizes this. The content below is thus based on FSL-MELODIC tool.



ICA for Resting-state fMRI
Our observation is the fMRI data. Like standard fMRI analysis, the pipeline starts with the common fMRI data preprocessing steps such as motion and slice timing correction, registration, spatial smoothing, etc. We have to also prepare the preprocessed dataset before ICA. First, simplify the 4D dataset into time vs space. Consider Yt´m as the whole-brain fMRI image of m voxels and t number of volumes, where t < m. We then mean-center the data by removing the mean spatial map from each row of Y. We normalize the time series variance such that each column of Y will have unit variance, a step called variance normalization. Next, we remove the mean time series from each column of Y. Without variance normalisation, the PCA step below will be biased towards tissue exhibiting high temporal variability.

Fig-1: Temporal std deviation of a sample dataset together with its PCA components (FSL slides).


Two typical yet strong assumptions are made when doing ICA in FSL: (a) the observed source signals are statistically independent; and (b) they are non-Gaussian or not normally distributed. So, ICA is performed on observations that are assumed to be a linear combination of independent sources. Our fMRI image Y is a mixture of independent sources such that Y = A.s , where the A is a t ´ t mixing matrix, and s is a t ´ m spatially independent components.

The first part of ICA is carrying out dimensionality reduction to simplify the problem or computation. This can be achieved with the help of PCA decomposition using SVD. The criterion is that only components that represent a large amount of the dataset remain. The original dataset is now whitened with unit variance and reduced dimension. Rewrite the equation, Yw = P.Y = P. A.s ; so we have the new mixing matrix Aw = P.A. Note: removing autocorrelation in the fMRI dataset is called whitening, a necessary step in the GLM parametric analyses to make the results more accurate! (Woolrich, et al., 2001). Different statistical software tools (AFNI, FSL, and SPM) have their own steps to accurately model this temporal correlation.

Because Yw is whitened, its variance-covariance matrix Σw is equal to identity matrix I. In other words:  Σw = (Yw YwT) / (n ─ 1) =  (Aws).(Aws)T/ (n ─ 1) = (Aws sTAwT) / (n ─ 1) = I. This is because the variance-covariance matrix of s.sT = I as the components are independent to each other. Note that for any orthogonal matrix, Aw-1 = AwT , making the problem solving even simpler.

So now, our unmixing matrix is the inverse of the new orthogonal matrix Aw. The job now is to predict or estimate the inverse of Âw, such that ŝ = Âw-1 Yw is maximally independent. Note that ŝ and Âw are estimates since they are both unknown. The criteria of independence can be empirically found by using a certain algorithm, e.g. mutual information, negentropy, etc.

In dealing with fMRI data, specifically, FMRIB Oxford team proposed a probabilistic ICA model in the noisy dataset, a more robust independent component estimation where the noise properties follow is assumed to be Gaussian (Beckmann, et al., 2004). For a more detailed but readable explanation to write this post, I consult Chapter 10 in Ashby, 2011, textbook. I have been using FSL-MELODIC to perform a single-subject ICA and their website is also informative to read.

Theoretically speaking, when all ICs are added together (each one being a 4D signal formed by the outer product of the spatial map and timecourse) they equal the original data. Unlike PCA, ICA enforces independence between the components spatially, while PCA enforces orthogonality both spatially and temporally. Sometimes, the ICs may also share similar time courses. Once all the independent components have been identified, they are ordered according to the degree of importance, that is, % total variance explained.

The last step in the MELODIC is the thresholding stage that produces thresholded ICs overlayed on a background image of your choice. A threshold level of 0.5 (set by default) means that a voxel 'survives' as soon as the probability of being in the 'active' class exceeds the probability of being in the 'background' noise class. This 0.5 assumes we set an equal loss on false-positives and false-negatives. The user guide says that if instead we consider e.g. false-positives as being twice as bad as false-negatives you should change this value to 0.66.

Fig-2: ICA method to rs-fMRI data (TR = 450 msec) displays an independent component associated with the physiological signal. The thresholded spatial map illustrates the location of brain regions correlated with the heart rate (~1 Hz). The files are generated by FSL MELODIC with image cropping for practical reasons.


Fig-3: ICA method to rs-fMRI data (TR = 450 msec) displays an independent component associated with the visual and posterior parietal area (0.01 - 0.1 Hz). The files are generated by FSL MELODIC with image cropping for practical reasons.

By default, FSL MELODIC produces the following files: 

- HTML report with a logfile that collates all steps performed by melodic.

- The mask, the list of ICs, and other statistical information (smoothest : estimated smoothness, PPCA : estimated intrinsic dimensionality estimated from PPCA).

- melodic_mix: an ASCII text file that contains the estimated mixing matrix (in the noise-free case the ICA decomposition is typically written as X = A*C, where X is the original data, A is the mixing matrix, and C is the matrix containing the estimated independent components as its rows). melodic_mix contains #ICs time courses as its columns. Each time course is plotted in the IC report that melodic produces.

- melodic_FTmix: a matrix containing the power spectrum at different freq for the time courses contained in melodic_mix (plotted in the IC report under the time courses)

- Eigenvalues_adjusted : the set of eigenvalues from the initial PCA decomposition (after variance normalisation and adjusting for the dimensionality of X).
 

Group ICA Method
To assess resting-state network with ICA at higher-level or group application, group ICA (gICA) is introduced. One can apply gICA to a set of resting-state data through temporal concatenation gICA; or a set of task-based fMRI using tensor gICA. In temporal concatenation gICA, we are finding common spatial patterns across subjects with unpredictable or different time series. In tensor gICA, we assume that each subject carries a similar time-series pattern, e.g. in task-based fMRI.

Temporal concatenation is recommended for the group-level resting-state pipeline as an alternative method to seed-based analysis. Often, researchers are faced with ambiguity in seed selection and gICA can be useful, especially with the use of dual-regression or back-projection method. For a more detailed explanation, refer to the review paper by V. Calhoun's group (NeuroImage, 2011).

Smith et al., (PNAS, 2009) show that the functional connectivity networks obtained at rest from gICA method are similar to different task-based activation networks obtained from a separate analysis of thousands of subjects. This is quite an important study as the authors show that resting-state networks are not simply random fluctuation, but are recruited when subjects perform the actual task. Several functional networks are identified from resting-state data, e.g. the visual network, sensorimotor network, attention network, limbic network, cerebellar network, the default mode network (DMN), and auditory network. This so-called functional segregation shows that the brain is active even at rest.


Tuesday, July 1, 2014

fMRI Data Preprocessing

MRI has become one of the most popular non-invasive brain imaging tools in research. By default, MRI scanners save the data in DICOM format, a standard format for handling, managing, and distributing imaging data. Image data are basically matrices of numbers and can be analyzed further with any software (Matlab for example). We have to do something on the image data before any statistical inference. Why?
- Like many devices, there are sort of noises contaminating the recorded data.
- Participants are not able to stay still ... obviously.
- Our experiment requires them to perform something that may cause the head to move.
- For the hope to increase confidence in the presence of the BOLD signal changes.

There are basically three sets of data one has to acquire when one performs fMRI-related studies:
- The high-resolution anatomical image. This is a T1-weighted image.
- The functional images in 4D (3D spatial + 1D temporal). These are generally T2* EPI images.
- Local field information if one wants to take into account field inhomogeneity.

Data Preparation 
The very first step before any processing is to convert DICOM into the format that your software package is using. SPM8 has its own function to import DICOM > NIfTI format. There is one extremely useful conversion software here by Chris Rorden called dcm2nii. Note: People at the MNI prefers MINC format. We can check the post-conversion images by using any MRI viewer software, for example, FSLView. Once you get the files converted, you do some housekeeping. You have to carefully categorize the files into the appropriate folder. For example: T1 anatomical images can be put separately from the multiple runs of functional datasets. 

Once file housekeeping is done, let's preprocess the dataset. To me, the easiest one is to strip the skull from the brain. The skull is not of our interest and it should be removed to simplify the data crunching during the analysis. Apply this step to both the anatomical image and raw EPI images. In addition, you have to remove a few volumes from the raw EPI data. Why? Because the magnetic system requires ~1.5 - 2.0 seconds to reach a steady state.

Head motion detection and correction
Substantial and sudden head motion bring negative effects to image acquisition. First, it contaminates the signal of a particular voxel with the neighboring voxels, making the voxel time series inaccurate. Second, it disrupts the magnetic field homogeneity in the bore that has been adjusted prior to functional scans. Managing head motion should be done before any other signal processing and statistical processes and usually involves the detection and correction or realignment.

Motion detection can be done by first selecting the target or reference volume. The choice of this reference volume can vary. For example: FSL uses the middle volume, AFNI uses either the first or a preselected volume. Often, head motion occurs very abruptly giving rise to spikes or outliers. Such outliers can be identified using a motion outlier script that uses certain metrics, e.g. rms intensity difference between the test and reference volume, or between volume n and n + 1 (DVARS), etc. Once chosen, a threshold to determine an outlier can be found using a boxplot method. Another approach includes a data exploration technique called Independent Component Analysis (ICA) which is more popular for resting-state fMRI. AFNI, though, prefers a different approach called motion scrubbing which means you remove some particular volumes where spikes occur.

Knowing when and how much the movements occur is one thing, but correcting for the contamination is another step. In the "correction" step, all other volumes will go through a 6 DOF rigid body transformation to find the best possible alignment with respect to the reference. The six parameters (or DOF) include 3 XYZ translations and 3 roll/pitch/yaw rotations.

Fig-1: Example of how motion correction analysis estimates the amount of head movement (by FSL-MCFLIRT).

Traditionally, we would want to compute the sum of the squared difference between the target and reference images. We would then minimize this or some cost function through iterative procedures or optimization. Once motion parameters for getting the best alignment have been determined, the new and resampled volumes will be created through spatial interpolation. Why? Because we have to recalculate the new voxel values for each volume after the adjustment. Such "correction" parameters are especially useful as regressors in the next pipeline (statistical analysis using GLM). This is usually a list of 6 parameters according to the rigid body transformation. Often, motion outliers are also used as additional regressors, e.g in FSL.

Slice timing correction
To begin, one 3D image is also called a volume and one scan produces a full functional MRI dataset comprising multiple volumes. When analyzing one 3D image it is assumed that all slices are acquired simultaneously. In reality, however, this isn't the case. Rather, slices are obtained sequentially according to a certain slice order. Thus there is bound to shift in the individual time course across voxels of different slices. For example, imagine we did a scan with a TR = 2.5 sec, i.e. it takes 2.5 sec for a volume to be fully acquired. What this means is that the difference in time between the very first and last slice in a volume would be ~ 2.5 sec. The problem is worsened by the way we acquire the slice, e.g. using an interleaved slice acquisition.

The severity of slice time misalignment depends on the repetition time (TR) and paradigm involved. It is usually very important for an event-related design, which may not be that crucial for blocked-design. All software package that I know of has the correction feature. AFNI and FSL use an almost similar basic method of interpolating the time series. Such interpolation follows the way how the excitation happens in the magnet. It can be ascending or descending or interleaves. FSL provides options to do slice timing correction based on custom-made slice-timing data.

Although controversial, the slice timing correction is usually performed after the head motion correction because the overall effect of the first is less severe than the latter.

Note: Slice time information can be obtained directly from the DICOM header. This can be easily retrieved using a function in Matlab or SPM, but not FSL.

Spatial smoothing (blurring)
This step involves applying a Gaussian kernel to each voxel in a three-dimensional way, which is essentially averaging data points with the neighboring voxels. A Gaussian kernel is a mathematical function that is specified by its width σ (sigma) at half maximum or FWHM, whose center coincides with the voxel center. As a rule of thumb FWHM is selected to be 2x - 3x the voxel size of our functional data, although in practice any value 5-10 mm is very common. Smoothing is also similar to applying a low-pass filter to the original data.

Why do we do spatial smoothing? First, it improves the signal-to-noise ratio. This is related to the assumption that when a voxel is activated, the surrounding voxels are activated as well. Also, we assume the noise of a voxel is not correlated with the noise of the adjacent voxel. Second, it helps to maintain a statistical validity associated with Random Field Theory to solve a multiple comparison problem. Lastly, smoothing is useful during the group-analysis because it improves inter-subject registration by blurring any residual anatomical difference. Risk of overdoing it? The optimal spatial resolution and some desired frequency components are lost if we use too big a kernel. There is also the risk of incorrectly placing the location of the activity or even missing the whole activity itself if it averages out the signal too much.

High-pass filtering
This is basically a voxelwise temporal filter to remove unwanted drift and low-frequency components across time that obscure the actual BOLD changes. This drift is inherent in hardware design and the filtering is common ever since PET is used. High-pass filtering also removes the temporal dc value associated with the image. Common software tool such as SPM uses a discrete cosine function added to the General Linear Model (GLM) to model the drift. FSL uses a high-pass Gaussian filter.

What is the cutoff frequency? FSL gurus recommended that the filter period should be 2x or 3x the task duration per block. A cutoff of 100 sec or 0.01 Hz is therefore very common and also long enough to avoid removing meaningful signals. For an event-related design, there is no clear stimulation period. In order to assess what the cutoff should be, one has to analyze the frequency content of the expected activations.

Some practical examples
Okie, I managed to play around with my own data using different preprocessing steps. The data were obtained from Siemens Trio 3T scanner, with TR=1.69 sec EPI sequence, voxel size 2 × 2 × 2 mm. I stick with FSL5.0 as my software tool. Look how each preprocessing step has an impact on the raw EPI image, both for resting-state and task-based scans. The task is as followed. It is a blocked design between "Move" and "Rest". The subject puts her fist in front of her chest and makes repeated upward movements during "Move", and rests her hand on her chest during "Rest". It's clearly shown that a high-pass filter (HPF) removes low-frequency drift and spatial smoothing with σ = 5.0 mm attenuates the signal amplitude. These outcomes are clearly seen in resting-state data when the brain is not engaged in any special task.

Fig-2: The different outcomes of each preprocessing stage to resting-state data in FSL-FEAT.
Fig-3: Similar treatment but to task-based data, blocked design (Move & Rest, 30-sec each), right arm localizer task. Changes in BOLD signal outweigh the low-frequency drift. In practice, preprocessing is followed by a statistical analysis using the general linear modeling (GLM) framework.


What is the output of any fMRI data analyses? It is the statistical map that shows different activated brain regions associated with the task. What is the effect of different FWHM values on our map? Figure 5 below sums up everything by using FSL software package. Only five selected brain slices are shown through the bottom view:
(a) no preprocessing at all: no slice and motion correction, smoothing, and filtering;
(b) with preprocessing, σ = 0 mm;
(c) with preprocessing, σ = 2.0 mm; also refer to filtered EPI image on the right-hand panel.
(d) with preprocessing, σ = 4.5 mm, but without motion correction (MCFLIRT);
(e) with preprocessing, σ = 4.5 mm; the preferred FWHM is twice the voxel width, not too much!
(f) with preprocessing, σ = 15 mm; also refer to filtered EPI image on the right-hand panel.

Skipping motion correction gives a noisy map because it causes many false-positive clusters. Overblurring or oversmoothing the raw data yields to enlarged activation clusters that are most likely false positive.

Fig-4: Different brain activation maps due to different preprocessing steps. The task is the same as the one in the previous figures. The maps are rendered on the standard MNI 152 template. The M1, premotor, S1, S2, and vermis are shown to be activated during Task execution w.r.t Rest.

Fig-5: Activation maps with respect to rest in coronal, sagittal, and axial views (Z = 3.5, p < 0.05) overlayed on the standard MNI 152 template. Different spatial smoothing parameter is shown with different colors: 2.0 mm (green), 4.5 mm (blue), and 15.0 mm (red).  

  
References
[1]  Some course materials from: http://fsl.fmrib.ox.ac.uk/fslcourse/
[2]  http://support.brainvoyager.com/functional-analysis-preparation/
[3]  Lindquist M. (2008). The statistical analysis of fMRI data. Statistical Science 23: 439–464.
[4]  Ashby, F. Gregory. (2011). Statistical Analysis of fMRI Data, 1st ed. MIT Press.
*** Note: There are abundant piles of past literature dealing with each step of the preprocessing pipeline.