** As of Dec 18, 2012, is (essentially) no longer being
     supported (by its author, Rayus Kuplicki of the University of Tulsa).

     Consider using (written by Rick Reynolds of the NIH).

This script takes care of preprocessing commonly done prior to resting state
analyses. The process is similar to ANATICOR and RETROICOR and consists of
removing a number of regressors of noninterest. Some regressors need to be
computed prior to running this script, others do not. These regressors include
all or a subset of (with requirements in parentheses):

        motion parameters (-align on, default)
        local or global WM signal (supply appropriate mask or segmentation file)
        average ventricle signal (supply appropriate mask or segmentation file)
        average whole brain signal (-includebrain)
        RVT signal (compute first using RetroTS.m)
        arbitrary ROI average regressors (supply masks)
        other arbitrary regressors supplied by the user (compute first)

Additionally, this script takes care of alignment, registration, and slice
timing correction via  It will also talairach the data if
desired-either before or after processing, your choice. Additionally, this
script implements two censoring methods.  One is based on outliers and the
other is identical to the method used by Power et al. (Neuroimage 2012)

Required Options:
        -anat aa        :aa is the high resolution anatomy file
        -epi ee         :ee is the epi timeseries to process

Optional Options:

        Informational Options:
                -help           :Display this help message
                -changelog      :Display a log of changes to this script

        Output Options:
                -prefix p       :Prepend p to each of the final output files.
                                Default is rest_proc.
                -dest d         :Put the results in directory d. Default is to
                                match p.
                -tsnr           :Compute the tSNR of the EPI as described below.
                -snr n c        :Compute the SNR of the EPI. n should be a scan
                                collected with the RF turned off. This will look
                                like static on an old TV and is used to estimate
                                the variance of the noise inherent in the
                                system. c should be the number of channels in
                                the coil you used and will determine a
                                correction factor.
                -corrmap        :Run 3dTcorrMap on the clean data.
                -corrmapt t     :Use t as the threshold when computing average
                                correlation strengths. Default is 0.3. The idea
                                here is that you may be interested in the
                                average correlation between each voxel and all
                                other voxels it is connected to, but below
                                a certain threshold two voxels could be
                                considered disconnected, so discard those weak
                                correlations. Regardless of what t is, mean
                                correlations without thresholding are also
                -script sc      :Write all commands to a script named sc. This
                                script can be modified and run later, similar
                                to the output from

        Alignment Options:
                -tlrc           :Do the preprocessing in talairach space.
                                Default is to stay in orig space.
                -tlrclast       :Do all preprocessing in orig space, but then
                                talairach the results. Pick either -tlrc or
                                -tlrclast (or neither).
                -episize mm     :Cubic voxel size of all datasets (other than
                                the anatomy) after transforming them to
                                talairach space. Default is 2. This only works
                                with -tlrclast or -tlrc.
                -align [on]/off :Do the Alignment etc. step. Turn this option
                                off if all of the datasets are already aligned.
                -alignbase b    :Align both the epi and anat datasets to b.
                                b should be an epi dataset and the first image
                                will be used for alignment. This option makes
                                sense if you have several epi runs and you want
                                all of them aligned to the same base.  Also,
                                this option only makes sense when not using
                                -tlrc or -epi2anat (-tlrclast is ok).
                -epi2anat       :Align epi to anat instead of anat to epi.
                                This only makes sense when not using -tlrc.
                -uniformize     :Uniformize anat before alignment.  Sometimes
                                This helps with skull stripping problems.
                -anat_has_skull [yes]/no
                                :Set this option to no if the anatomy has
                                already been skull stripped (useful when default
                                skull stripping doesn't work right).

        Regressor Options:
                -aseg a         :a is the aseg segmentation file from
                                freesurfer. It should be aligned with the
                                anatomy supplied as -anat and can be in .mgz,
                                .nii, or .BRIK format.
                -wmsize w       :Radius (in mm) of the sphere to use when
                                computing the local white matter regressors.
                                Default is 15mm
                -globalwm       :Use the global wm average as a single
                                regressor instead of computing local wm
                -venterode v    :Number of nonmask neighbors required to cause
                                erosion in the ventricles. Default is 2
                -wmerode we     :Number of nonmask neighbors required to cause
                                erosion in the WM mask. Default is 1.
                -rvt r          :r is the RVT file produced by RetroTS.m
                -includebrain   :Include the whole brain average regressor.
                -dreg           :Add the derivatives of all regressors as
                -regressor re   :Use re as a regressor. re will be processed
                                in the same way the other regressors are
                                (detrended, catenated). If you do alignment
                                and registration outside of this script, it may
                                be a good idea to provide the motion parameters
                                as a regressor. re can be either a 3d+t volume
                                (specifying a different regressor for each
                                voxel) or a .1D file (specifying a single global
                -globalregmask g:Use the average signal extracted from the mask
                                g as a global regressor of noninterest. This
                                will produce one regressor used for all voxels.
                -localregmask rm rr
                                :Use the local average signal extracted from
                                rm as a regressor of noninterest. This will
                                produce a different regressor for each voxel.
                                rm should be a mask defining the ROI to use
                                and rr is the radius in mm to use when computing
                                local average signals.
                                -regressor, -globalregmask and -localregmask
                                can be used multiple times to supply an
                                arbitrary number of regressors.

        Censoring Options:
                -outcensor      :Censor timepoints based on their number of
                                outliers and head motion magnitude. Censored
                                time points are cut out.
                -fraclimit f    :When using -outcensor, fraction of voxels
                                identified as outliers needed to censor a time
                                point. Default is 0.05.
                -motlimit m     :When using -outcensor, limit on rms motion to
                                censor a point. Default is 0.3.
                -dvarscensor    :Create a censor file based on FD (framewise
                                displacement) and DVARS as from Power et. al,
                                Neuroimage 2012.
                -fdlimit ff     :Set the FD limit to be ff
                -dvarslimit dd  :Set the DVARS limit to be dd
                -censorleft s   :Censor s steps to the left of bad time points.
                                Default is 1.
                -censorright ss :Censor ss steps to the right of bad time
                                points. Default is 2.
                -censorunion    :Censor the union of fraclimit and motlimit or
                                FD and DVARS, instead of the intersection.
                -keepuncensored :Keep a copy of the uncensored timeseries. It
                                will be called

        Normalization Options:
                -localnorm      :Normalize based on voxelwise mean
                -globalnorm     :Normalize based on global mean
                -modenorm       :Normalize based on global mode using 100 bins
                -normval n      :Scale the selected attribute to be n

        Smoothing Options:
                -smooth [on]/off:Smooth the clean timeseries data.
                -smoothrad s    :FWHM size of smoothing to apply after cleaning
                                the data.  Default is 4mm. Smoothing is done
                                using a grey/nongrey matter mask by default.
                -smoothtogether :Smooth everything inside a brain mask together,
                                rather than smoothing the grey/nongrey matter
                -smoothfirst    :Smooth the data before doing regression,
                                instead of after.

        Misc. Processing Options:
                -despike[on]/off:Despike the timeseries as the fist
                                preprocessing step.
                -trcut t        :Number of TRs to throw away.  Default is 4.
                -polort p       :Polynomial to detrend from the regressors and
                                the timeseries.  Similar to 3dDeconvolve
                                -polort A, default is floor(1 + TR*nVOLS / 150).
                -bandpass       :Do bandpass filtering with LHz < f < HHz.
                                Default is 0.009 and 0.08.
                -setbands L H   :Set L and H for bandpass filtering
                -bpassregs      :Also bandpass filter the regressors.
                -exec [on]/off  :Execute the commands. Turn this off and use
                                -script to get things setup without running

        Other Options:
                -apply_censor e c p
                                :This option is used to apply a censor file to
                                remove timepoints from a timeseries. If it is
                                given, this is the only option that will be
                                processed. e is the timeseries to censor. c is
                                a 1D file consisting of a single column of
                                0's and 1's which must be the same length as e.
                                Time points with 1's in c will be kept, 0's will
                                be discarded. p is the prefix to use for the
                                output timeseries.

The following steps should be done before running this script:
        Create anatomical regressor masks:
                If you want to remove anatomical regressors or noninterest,
                the average ventricle signal, for example, you will need to
                provide masks used to extract these signals. This can be done
                in two different ways. Either supply the aseg file produced by
                freesurfer, which is used to extract the ventricle and white
                matter ROIs, or supply your own arbitrary masks with the
                -globalregmask or -localregmask options.
        Align the segmentation file with the experimental anatomy:
                If given, The aseg file from freesurfer (it can be in mgz, nii,
                or BRIK format) is assumed to be in alignment with the
                experimental anatomy.  The aseg and anat files will already be
                aligned if the anatomy is the one used by freesurfer.  If it is
                not, you may need to use something like @SUMA_AlignToExperiment.
        Create the RVT file:
                This is done by processing the experiment's cardiac and
                respiratory files using RetroTS.m, available in the AFNI Matlab
                library. While it is probably beneficial to remove the
                estimated cardiac and respiratory signals, this step is not
                necessary and this script will run fine without them.

Processing is done in the following steps:
        Copy Files:
                The output directory (specified by -dest) is created and input
                files are copied to dest/tmp
                This step is done first, if at all, so that spikes are not
                'smeared around' by registration and slice timing correction
        Alignment etc.:
                This step aligns the epi and anat datasets while also taking
                care of slice timing correction and the talairach
                transformation, if requested. These steps are combined using
       to minimize the number of interpolations
                required. If processing is done in +orig space, the anat is
                aligned to the epi by default. Using -epi2anat will cause the
                epi to be aligned to the anat. The appropriate transformation
                is also applied to the aseg file and any masks provided by
                -localregmask and -globalregmask to keep them aligned with the
                This is where the first few time points are thrown out.  This
                step is delayed until after alignment so that the first high
                contrast epi image can be used in the alignment process.
                At this stage in processing the tSNR of the EPI data is
                computed if requested.  It is taken to be
                mean(pEPI) / stdev(det(pEPI)) where pEPI is the processed EPI
                (despiked, aligned, tshifted, catenated) and det represents
                detrending with polynomial order polort.
                If you have collected a scan with the RF turned off (used to
                estimate the varience of the noise) the SNR can be computed
                for each voxel. It is taken to be S/(sigma*corr) where S is
                the signal in the first frame of the epi timeseries, sigma is
                the standard deviation of the values in the noise scan, and
                corr is a correction factor based on the number of channels in
                the coil. After computing the SNR of the original EPI volume,
                it is transformed to be aligned with and have the same voxel
                size as the final EPI. Correction factors are:
                        corr1 = 1.5263997
                        corr8 = 1.4257312
                        corr16 = 1.4198559
                        corr32 = 1.4170053
        Normalize EPI:
                If a method of normalization is chosen, it is applied here.
                Normalization is done to scale the selected attribute to be
                -normval inside a brain mask. Specifically:
                -globalnorm computes the global mean signal across time and
                space in the brain and scales accordingly.
                normval * (voxel intensity)/(global mean)
                -modenorm computes the global mode intensity across time and
                space in the brain and scales accordingly.
                normval * (voxel intensity)/(global mode)
                -localnorm computes the temporal average for each voxel and
                scales accordingly.
                normval * (voxel intensity)/(voxel mean)
        Prep WM Mask:
                The WM mask is taken from the aseg file from freesurfer.  It is
                first taken to be all voxels with labels: 2,7,16,41,46,251,252,
                253,254,255. This mask is resampled to match the resolution of
                the epi dataset. After resampling, the mask is eroded so that
                it is less likely to contain any grey matter. By default,
                typical erosion is done which removes voxels which have a
                single non-mask neighbor from the mask.
        Prep Ventricle Mask:
                The ventricle mask is taken from the aseg file using labels 4
                and 43. This mask is resampled to match the epi and then eroded.
                Erosion of the ventricle mask is by default less conservative
                than typical erosion. Voxels require two neighbors to be
                non-mask voxels in order to be eroded. This is done because a
                large number of subjects end up without any voxels in the
                ventricle mask using standard erosion and a 64x64 matrix. If
                your data are higher resolution, you may want to use -venterode
                1. It is a good idea to check both the WM and ventricle masks
                to make sure they look good.
        Prep Blurring Mask:
                The last step (after regression) is to apply gaussian smoothing.
                By default (if -aseg was specified), this smoothing is done in
                the grey and nongrey matter seperately via 3dBlurInMask. The
                blurring mask is created so that grey matter voxels are labeled 1
                and nongrey voxels (inside the brain) are labeld 2. The labeling
                is simply (automask + WM mask + Vent mask). If -aseg was not
                specified, or -smoothtogether was given, the smoothing is done
                using the whole brain as one region.
                If -smoothfirst was selected, this is where smoothing takes
                place. It is done as described below.
        Extract Regressors from masks:
                Regressor timeseries are extracted from the WM and ventricle
                masks as well as any masks supplied by the user with the
                -globalregmask and -localregmask options. A single regressor
                timeseries is computed as the average value for each supplied
                global masks and the ventricle mask. For the WM mask and any
                masks supplied with -localregmask, different regressors are
                computed for each voxel. For a single voxel, the regressor is
                defined as the average timeseries in voxels which are both
                within the supplied radius and included in the mask.
                -globalregmask and -localregmask are useful if you want to use
                software other than Freesurfer for the segmentation step. For
                example, using whatever method you like, you can create a
                ventricle mask and supply it as a global mask. Likewise, you can
                create a whitematter mask and supply it as a local regressor.
                Multiple local and global masks can be supplied.
        Differentiate Regressors:
                If desired, the temporal derivatives of each regressor are
                computed and added to the list of regressors.
        Detrend Regressors:
                The polynomial of order -polort is removed from each of the
                regressors (RVT, WM, Vent, Motion). This is done so there are
                no competing polynomial terms during the regression step.
        Bandpass Filtering:
                If bandpass filtering is selected, it is applied to the EPI
                data here after regression. This is also where the regressors
                are bandpass filtered if -bpassregs was selected.
                The regressors of noninterest (RVT,WM,Vent,Motion,other
                arbitrary regressors) are taken out of the epi timeseries using
                3dTfitter, which also removes the polynomial selected using
        Create Censor File:
                If a censoring method was chosen, the offending time points are
                identified here.
                        Make a temporal mask marking frames with more outliers
                        than the threshold specified by -fraclimit.
                        Make a temporal mask marking frames with more RMS
                        motion than specified by -motlimit.
                        Mark frames -censorleft and -censorright steps to the
                        left and right of time points flagged in the two masks.
                        Take the intersection of the two masks created above
                        (or the union, if -censorunion was specified)
                        Make a temporal mask marking frames with FD greater
                        than the threshold specified by -fdlimit.
                        Make a temporal mask marking frames with DVARS greater
                        than the threshold specified by -dvarslimit.
                        Mark frames -censorleft and -censorright steps to the
                        left and right of time points flagged in the two masks.
                        Take the intersection of the two masks created above
                        (or the union, if -censorunion was specified)
                By default, 3dBlurInMask is used to smooth the timeseries in
                the grey and nongrey matter separately. Grey matter voxels are
                likely the interesting ones, but it can't hurt to apply the
                same process to nongrey voxels to see what they look like. The
                -smoothtogether flag can be used to apply uniform smoothing to
                all voxels in the brain instead. -smooth off skips this step.
                If a censoring method was chosen, the censored time points are
                removed here.
                At this point, the data have been preprocessed to remove
                uninteresting signals. It is now appropriate to do resting
                state functional connectivity analysis on the clean data. One
                thing to examine is the result of running 3dTcorrMap. For a
                full description of what this does, see the help from
                3dTcorrMap. This script uses it as follows:
                        3dTcorrMap -input cleanEPI -mask automask -polort -1
                                -mean prefix.MeanCorr -Hist 400 prefix.CorHist
                                -Cexpr 'step(r-t)*r' prefix.MeanCorrGT
                where t can be specified using -corrmapt and is 0.3 by default.

Things to check after running this script:
                Make sure the various mask datasets are in good alignment with
                the anatomical dataset.
        Mask Coverage:
                Make sure the ventricle and white matter masks cover what you
                think are appropriate voxels.

Example Usage:

        #Basic usage:
        #Remove RVT, motion parameters, WM and ventricle signals from
        #Store the results in a directory named preproc
        #Prefix each result file with subjX
        #Processing is done in orig space -anat mprage+orig. \
                -epi epi+orig. \
                -rvt RVT.slibase.1D \
                -aseg aseg.mgz \
                -dest preproc \
                -prefix subjX

        #Produce a tsnr map and results from 3dTcorrMap using a threshold of .15
        #Write a script called proc.tcsh but don't execute it yet
        #This script can be modified and executed at your leisure -anat mprage+orig. \
                -epi epi+orig. \
                -rvt RVT.slibase.1D \
                -aseg aseg.mgz \
                -dest preproc \
                -prefix subjX \
                -corrmap \
                -corrmapt .15 \
                -tsnr \
                -script proc.tcsh \
                -exec off

        #Alignment and talairaching were done already, so skip those steps
        #Use the provided motion parameter file as a regressor \
                -epi epi+tlrc \
                -rvt RVT.slibase.1D \
                -anat mprage+tlrc \
                -aseg aseg+tlrc \
                -regressor epi_tsh_vr_motion.1D \
                -dest prealigned \
                -prefix subjX \
                -align off

        #Do processing like it was done in Power et al. Neuroimage 2012 \
                -despike off \
                -aseg aseg.mgz \
                -anat mprage+orig \
                -epi rest+orig \
                -script power_method.tcsh \
                -dest power_method_subjx \
                -prefix pm \
                -dvarscensor \
                -tlrc \
                -episize 3 \
                -dreg \
                -smoothfirst \
                -smoothrad 6 \
                -smoothtogether \
                -bandpass \
                -includebrain \
                -polort 0 \
                -globalwm \
                -censorleft 1 \
                -censorright 2 \
                -fdlimit 0.5 \
                -dvarslimit 5 \

        #Apply a censor file to a timeseries.  This will output a file called
        #epi.censored+orig that has TRs cut out wherever censor.1D is 0. \
                -apply_censor \
                epi+orig \
                censor.1D \

Original version by Rayus Kuplicki.
University of Tulsa
Laureate Institute for Brain Research
Report problems or feature requests to