detrending and normalizing timecourses: afni 3dDetrend and 3dDeconvolve in R

This post introduces an expanded and updated demo of how detrending and normalizing ("scaling") of individual voxel (or vertex) timecourses works with afni 3dDeconvolve and 3dDetrend commands. 

In my original post I showed how to duplicate what 3dDetrend -normalize -polort 2 does with R code, not to replace the afni function, but to understand what it is doing. This expanded version simplifies the old code a bit, but more importantly, adds sections explaining another common method of preparing timecourses for analysis: using 3dDeconvolve with -num_stimts 0 -polort A -errts (plus censoring, motion regressors, etc.; see below) to create the residual error timeseries. 

The compiled demo is afni3dDetrend3dDeconvolve_R.pdf  and is a knitr file; its source code (with many comments) and files required for compilation are in a section of my osf site, DOI https://doi.org/10.17605/OSF.IO/NU324 (please include that DOI in any citation). 

starting point: a voxel's timecourse

The examples use a left motor grey matter voxel from a preprocessed (with fmriprep 1.3.2) fMRI task run from the DMCC55B dataset; I chose it arbitrarily. The plot below is directly from the preprocessed nifti; this voxel has values around 6900, and the timecourse vector is length 540. DMCC55B's TR was 1.2 s, so this is a 10.8 minute-long run. The grey vertical lines are at one-minute intervals, with TR (frame) number along the x-axis and BOLD amplitude along the y-axis. (See the knitr .rnw for details, plotting code, etc.; this blog post just has a few highlights.)

normalizing and detrending

Scaling alone isn't usually sufficient to prepare fMRI timeseries for analysis; we also need at least a bit of detrending. There's no universally correct degree or type of detrending to use. I generally recommend a modest amount of detrending before parcel-averaging types of analyses, specifically 3dDetrend -normalize -polort 2 . 

In the plot below, the same voxel timecourse as above is plotted after normalizing only (tc.np0, black, from 3dDetrend -normalize -polort 0), or with detrending at polort 2 (blue, tc.np2), and the more aggressive polort 5 (green, tc.np5). 

Notice that the spiky parts of the timecourses are pretty much the same in all three versions, but the slower changes vary more; e.g., in the first minute without no-detrending line (tc.np0) is furthest from zero, the tc.np2 line is closer, and tc.np5 line closest. 

It's sensible that larger polort numbers have more of an effect on the timecourse's shape, since, as explained in the afni help for 3dDetrend, -polort ppp gives "the Legendre polynomials of order up to and including 'ppp' in the list of vectors to remove", so larger -polort numbers means removing more complex trends. The R code below shows how to do this type of normalizing and detrending; Legendre() is from Gregor Gorjanc, and requires the orthopolynom and polynom R packages.

 # R commands for 3dDetrend -normalize -polort 2  
 lm.out <- lm(tc.raw ~ Legendre(x=seq(tc.raw), n=2)); # lm with two Legendre polynomials (polort 2)  
 tmp <- residuals(lm.out);  # extract the residuals  
 tc.Rnp2 <- (tmp-mean(tmp))/sqrt(sum((tmp-mean(tmp))^2));  # normalize the residuals, afni-style   

residual error via 3dDeconvolve

It's common to have the realignment parameters as nuisance regressors and censor high-motion frames prior to fMRI timecourse analyses, which can be done with 3dDeconvolve. Skipping rather a lot of explanations from the full demo, the afni command is:

 # errts.fname is the file made by 3dDeconvolve, from which the single voxel timecourse was extracted:  
 system2(paste0(afni.path, "3dDeconvolve"),   
     args=paste0("-input '", scale.fname, "' -polort A -float ",  
           "-censor '", c.fname, "' -num_stimts 0 ",  
           "-ortvec '", mot.fname, "' moveregs ",  
           "-nobucket -errts ", errts.fname), stdout=TRUE);  

where -num_stimts 0 means not to include any events in the model,  -errts that we want afni to write the residual error time series from the "full model fit to the input data" into file errts.fname (in this case, a .nii.gz), and -polort A that afni should set the polort level according to the run length (here, that gives 5).

Below is the 3dDetrend -normalize -polort 5 (tc.np5) timecourse again in green, with the new (-errts) version from the 3dDeconvolve command in pink: 

The errts timecourse is highly correlated with the np5 version, which makes sense, since both included polort 5 detrending. They're not perfectly correlated, though: the 3dDeconvolve command also did censoring and included the motion regressors. I don't have a simple way to describe the differences in the lines; they're clearly very similar, but not identical; sometimes one is more extreme or spiky, sometimes the other.

The errts image has 0 in the censored frames. This is obvious in a 4d nifti (entire frame filled with 0s), but ambiguous in a single voxel (or vertex) timecourse like this (in the plot the censored frames are circled on the tc.errts timecourse, squared on the np5 version). For some analyses in R (e.g., averaging frames after an event for temporal compression) it'd be sensible to use NA for the censored frames.

This R code matches the 3dDeconvolve calculations:

 # made in startup code chunk; download at https://osf.io/nu324/files/c6nax  
 mot.fname <- paste0(demo.path, "sub-f1027ao_ses-wave1bas_task-Stroop_6regressors_demean.txt");   
 mot.tbl <- read.delim(mot.fname, sep=" ", header=FALSE); # 540 x 6  
   
 # made in startup code chunk; download at https://osf.io/nu324/files/bjrk8  
 c.fname <- paste0(demo.path, "sub-f1027ao_ses-wave1bas_task-Stroop_FD_mask0.5.txt");  # 0 1 censor file  
 censor.vec <- read.table(c.fname, header=FALSE)[,1];  
 censor.TRs <- which(censor.vec == 0) # [1] 300 316 431 448 497  
   
 # first, remove censored frames from the motion regressors and input timecourse  
 # The input timecourse tc.scale is the voxel from file scale.fname: the input bold.fname after  
 # scaling with 3dcalc -expr 'min(200, a/b*100)*step(a)*step(b)'   
 scale.vec <- tc.scale[-censor.TRs];  
 col1.vec <- mot.tbl$V1[-censor.TRs]; # demeaned trans_x  
 col2.vec <- mot.tbl$V2[-censor.TRs];   
 col3.vec <- mot.tbl$V3[-censor.TRs];   
 col4.vec <- mot.tbl$V4[-censor.TRs];   
 col5.vec <- mot.tbl$V5[-censor.TRs];   
 col6.vec <- mot.tbl$V6[-censor.TRs];   
   
 # fit the lm, polort 5, on the censored scale.vec and including 6 censored motion regressors  
 lm.out <- lm(scale.vec ~ Legendre(x=seq(scale.vec), n=5) + col1.vec + col2.vec + col3.vec + col4.vec + col5.vec + col6.vec);   
 tc.Rerrts <- residuals(lm.out);  # extract the residuals  
   
 # put 0s back in where the censored frames were taken out.  
 for (i in 1:length(censor.TRs)) { tc.Rerrts <- append(tc.Rerrts, 0, (censor.TRs[i]-1)); }  
   
 # the R version matches the 3dDeconvolve errts version   
 cor(tc.errts, tc.Rerrts); # almost perfect  

musings

Working out the R code which matched the 3dDeconvolve -errts demystified it for me; the 3dDetrend -normalize -polort 2 detrending and normalizing ("np2") I've treated as a default (for timecourse-averaging type analyses) is closer to these 3dDeconvolve residuals ("errts") than I'd thought. 

Am I going to change my default detrending method? Is the 3dDeconvolve errts "better" than the 3dDetrend np2 for an analysis like this? I think incorporating the censoring (to NA, not 0) and polort-picking calculation (which gave 5 in this demo) from 3dDeconvolve (instead of using 2 regardless of run length) would be sensible, modest improvements.

I'm less confident about whether to add in the motion regressors. My sense was that including these would somehow "account for" or "clean up" any motion effects not "corrected" by preprocessing. And including the six realignment parameter columns does change the timecourse produced by the model a bit ... but it's not much different if the actual realignment parameters or random ones are included, making me suspect the change is more due to the change in model degrees of freedom than actually "fixing" the head motion.

Here's code for the random-motion-regressor models:

 # permute numbers in each motion regressor separately  
 lm.out <- lm(scale.vec ~ Legendre(x=seq(scale.vec), n=5) + sample(col1.vec) + sample(col2.vec) + sample(col3.vec) + sample(col4.vec) + sample(col5.vec) + sample(col6.vec));   
 tc.test1 <- residuals(lm.out);  # extract the residuals  
 for (i in 1:length(censor.TRs)) { tc.test1 <- append(tc.test1, 0, (censor.TRs[i]-1)); } # put 0s back in  
   
 # random numbers for the motion regressors  
 ct <- length(scale.vec);  # how long to make each fake motion regressor column  
 lm.out <- lm(scale.vec ~ Legendre(x=seq(scale.vec), n=5) + rnorm(ct) + rnorm(ct) + rnorm(ct) + rnorm(ct) + rnorm(ct) + rnorm(ct));   
 tc.test2 <- residuals(lm.out);  # extract the residuals  
 for (i in 1:length(censor.TRs)) { tc.test2 <- append(tc.test2, 0, (censor.TRs[i]-1)); } # put 0s back in  

DMCC GLMs: afni TENTzero, knots, and HDR curves, #2

This post is the second in a series, which begins with this post. In that first post I included examples of the estimated HDR curves for a positive control analysis, one of which is copied here for reference. 

What are the axis labels for the estimated HDR curves? How did we generate them? What made us think about changing them? 

At the first level the (curve) axis label question has a short answer: y-axes are average (across subjects) beta coefficient estimates, x-axes are knots. But this answer leads to another question: how do knots translate into time? Answering this question is more involved, because we can't translate knots to time without knowing both the TR and the afni 3dDeconvolve command. The TR is easy: 1.2 seconds. Listing the relevant afni commands is also easy:

3dDeconvolve \  
 -local_times \  
 -x1D_stop \  
 -GOFORIT 5 \  
 -input 'lpi_scale_blur4_tfMRI_AxcptBas1_AP.nii.gz lpi_scale_blur4_tfMRI_AxcptBas2_PA.nii.gz' \  
 -polort A \  
 -float \  
 -censor movregs_FD_mask.txt \  
 -num_stimts 3 \  
 -stim_times_AM1 1 f1031ax_Axcpt_baseline_block.txt 'dmBLOCK(1)' -stim_label 1 ON_BLOCKS \  
 -stim_times 2 f1031ax_Axcpt_baseline_blockONandOFF.txt 'TENTzero(0,16.8,8)' -stim_label 2 ON_blockONandOFF \  
 -stim_times 3 f1031ax_Axcpt_baseline_allTrials.txt 'TENTzero(0,21.6,10)' -stim_label 3 ON_TRIALS \  
 -ortvec motion_demean_baseline.1D movregs \  
 -x1D X.xmat.1D \  
 -xjpeg X.jpg \  
 -nobucket  
 
3dREMLfit \
 -matrix X.xmat.1D \
 -GOFORIT 5 \
 -input 'lpi_scale_blur4_tfMRI_AxcptBas1_AP.nii.gz lpi_scale_blur4_tfMRI_AxcptBas2_PA.nii.gz' \
 -Rvar stats_var_f1031ax_REML.nii.gz \
 -Rbuck STATS_f1031ax_REML.nii.gz \
 -fout \
 -tout \
 -nobout \
 -verb

Unpacking these commands is not so easy, but I will try to do so here for the most relevant parts; please comment if you spot anything not quite correct in my explanation.

First, the above commands are for the Axcpt task, baseline session, subject f1031ax (aside: their data, including the event timing text files named above, is in DMCC13benchmark). The 3dDeconvolve command has the part describing the knots and HDR estimation; I included the 3dREMLfit call for completeness, because the plotted estimates are from the Coef sub-bricks of the STATS image (produced by -Rbuck).

Since this is a control GLM in which all trials are given the same label there are only three stimulus time series: BLOCKS, blockONandOFF, and TRIALS (this is a mixed design; see introduction). The TRIALS part is what generates the estimated HDR curves plotted above, as defined by TENTzero(0,21.6,10).

We chose TENTzero instead of TENT for the HDR estimation because we do not expect anticipatory activity (the trial responses should start at zero) in this task. (see afni message board, e.g., here and here) To include the full response we decided to model the trial duration plus at least 14 seconds (not "exactly" 14 seconds because we want the durations to be a multiple of the TENT duration). My understanding is that it's not a problem if more time than needed is included in the duration; if too long the last few knots should just approximate zero. I doubt you'd often want to use a duration much shorter than that of the canonical HRF (there's probably some case when that would be useful, but for the DMCC we want to model the entire response).

From the AFNI 3dDeconvolve help: 

'TENT(b,c,n)': n parameter tent function expansion from times b..c after stimulus time [piecewise linear] [n must be at least 2; time step is (c-b)/(n-1)]

 

You can also use 'TENTzero' and 'CSPLINzero',which means to eliminate the first and last basis functions from each set. The effect of these omissions is to force the deconvolved HRF to be zero at t=b and t=c (to start and and end at zero response). With these 'zero' response models, there are n-2 parameters (thus for 'TENTzero', n must be at least 3).

The first TENTzero parameter is straightforward: b=0 to start at the event onset.

For the trial duration (c), we need to do some calculations. Here, the example is the DMCC Axcpt task, which has a trial duration of 5.9 seconds, so the modeled duration should be at least 14 + 5.9 = 19.9 seconds. That's not the middle value in the TENTzero command, though.

For these first analyses we decided to have the TENTs span two TRs (1.2 * 2 = 2.4 seconds), in the (mostly) mistaken hope it would improve our signal/noise, and also to make the file sizes and GLM estimation speeds more manageable (which it does). Thus, c=21.6, the shortest multiple of our desired TENT duration greater than the modeled duration (19.9/2.4 = 8.29, rounded up to 9; 9*2.4 = 21.6). 

Figuring out n requires a bit of mind bending but less calculation; I think the slide below (#5 in afni06_decon.pdf) helps: in tent function deconvolution, the n basis functions are divided into n-1 intervals. In the above calculations I rounded 8.29 up to 9; we want to model 9 intervals (of 2.4 s each). Thus,  n = 10 basis functions gives us the needed n-1 = 9 intervals.

Recall that since we're using TENTzero, the first and last tent functions are eliminated. Thus, having n=10 means that we will get 10-2=8 beta coefficients ("knots") out of the GLM. These, finally, are the x-axes in the estimated GLM curves above: the knots at which the HDR estimates are made, of which there are 8 for Axcpt. 

(Aside: 0-based afni can make some confusion for those of us more accustomed to 1-based systems, especially with TENTzero. The top plots put the value labeled as #0_Coef at 1; a point is added at 0=0 since TENTzero sets the curves to start at zero. I could also have added a zero for the last knot (9, for this Axcpt example), but did not.)

These GLMs seem to be working fine, producing sensible HDR estimates. We've been using these results in many analyses (multiple lab members with papers in preparation!), including the fMRIprep and HCP pipeline comparisons described in previous posts. We became curious, though, if we could set the knot duration to the TR, rather than twice the TR (as done here): would the estimated HDR then follow the trial timing even more precisely? In the next posts(s) I'll describe how we changed the GLMs to have the knot duration match the TR, and a (hopefully interesting) hiccup we had along the way.

UPDATE 4 January 2021: Corrected DMCC13benchmark openneuro links.

DMCC GLMs: afni TENTzero, knots, and HDR curves, #1

One of the things we've been busy with during the pandemic-induced data collection pause is reviewing  the DMCC (Dual Mechanisms of Cognitive Control; more descriptions here and here) project GLMs: did we set them up and are we interpreting them optimally? 

We decided that the models were not quite right, and are rerunning all of them. This post (and its successors) will describe how the original GLMs were specified, why we thought a change was needed, what the change is, and how the new results look. The posts are intended to summarize and clarify the logic and process for myself, but also for others, as I believe many people find modeling BOLD responses (here, with afni 3dDeconvolve TENTs) confusing.

 

I won't just start explaining the DMCC GLMs however, because they're complex: the DMCC has a mixed design (task trials embedded in long blocks), and we're estimating both event-related ("transient") and block-related ("sustained") responses. If you're not familiar with them, Petersen and Dubis (2012) review mixed block/event fMRI designs, and Figure 1 from Visscher, et. al (2003) (left) illustrates the idea. In addition to the event and block regressors we have "blockONandOFF" regressors in the DMCC, which are intended to capture the transition between the control and task periods (see Dosenbach et. al (2006)).

We fit many separate GLMs for the DMCC, some as controls and others to capture effects of interest. For these posts I'll describe one of the control analyses (the "ONs") to make it a bit more tractable. The ON GLMs include all three effect types (event (transient), block (sustained), blockONandOFF (block start and stop)), but do not distinguish between the trial types; all trials (events) are given the same label ("on"). The ON GLMs are a positive control; they should show areas with changes in activation between task and rest (30 seconds of fixation cross before and after task blocks). For example, there should be positive, HRF-type responses in visual areas, because all of our task trials include visual stimuli.

Below is an example of the estimated responses from our original ONs GLMs (generated by R knitr code similar to this). The GLMs were fit for every voxel individually, then averaged within the Schaefer, et al. (2018) parcels (400x7 version); these are results for four visual parcels and volumetric preprocessing. These are (robust) mean and SEM over 13 participants (DMCC13benchmark). The results have multiple columns both because there are four DMCC tasks (Axcpt (AX-CPT), Cuedts (Cued task-switching), Stern (Sternberg working memory), and Stroop (color-word variant)), and because of the mixed design, so we generate estimates for both event/transient (lines) and sustained (blue bars) effects (blockONandOFF are not included here).

The line graphs show the modeled response at each time point ("knot"). Note that the duration (and so number of estimated knots) for each task varies; e.g., Stroop has the shortest tasks and Stern the longest. (We set the modeled response to trial length + 14 seconds.)

The double-peak response shape for Axcpt, Cuedts, and Stern is expected, as the trials have two stimuli slides separated by a 4 second ITI; the Stroop response should resemble a canonical HRF since each trial has a single stimulus and is short. The task timing and some sample trials are shown below (this is  Figure 1 in a preprint (accepted version) which has much more task detail; see also the DMCC55B description). 

So far, so good: the estimated responses for each task have a sensible shape, reflecting both the HRF and task timing. But what exactly is along the x-axis for the curves? And how did we fit these GLMs, and then decide to change them a bit? ... stay tuned.

later posts in this series: #2

UPDATE 4 January 2021: Corrected DMCC13benchmark openneuro links.

detrending and normalizing timecourses: afni 3dDetrend in R

UPDATE 27 February 2026: An expanded version of this demo is summarized in this post, full version at https://doi.org/10.17605/OSF.IO/NU324.

I've been using the afni 3dDetrend command on datasets lately, but realized that I wasn't sure of the exact way the normalization and detrending was calculated. To help me understand (and reproduce the results on files that aren't easily read by afni), this post gives R code to match 3dDetrend output, plus shows some examples of how various normalizing and detrending schemes change timecourses (there's more in the knitr and pdf than shown in this post; files below the jump).

First, we need some data. I picked a run from a handy preprocessed dataset (it doesn't really matter for the demo, but it's from multibandACQtests), then ran the image through afni's 3dDetrend command twice: once with the options -normalize -polort 0, then with -normalize -polort 2. This gives three images, and I picked a voxel at random and extracted its timecourse from each image. Code for this step starts at line 22 of the knitr file (below the jump), and text files with the values are included in the archive for this post.

Here is the voxel's "raw" timecourse: after preprocessing (the HCP pipelines, in this case), but before running 3dDetrend:

And here is the same voxel's timecourse, from the image after 3dDetrend -normalize -polort 2: centered on zero, with the shift of baseline around TR 50 reduced a bit.

From the afni documentation, 3dDetrend -normalize "Normalize[s] each output voxel time series; that is, make the sum-of-squares equal to 1." Note that this is not what R scale() does by default (which is to mean 0 and standard deviation 1; see the full demo for a comparison). Line 113 of the demo code does the afni-style normalization: (tc.raw - mean(tc.raw)) / sqrt(sum((tc.raw - mean(tc.raw))^2)), where tc.raw is a vector with the voxel's timecourse.

For the -polort 2 part, the afni documentation says: "Add Legendre polynomials of order up to and including 'ppp' in the list of vectors to remove." I used the functions from Gregor Gorjanc's blog post (thanks!) and the R orthopolynom package for this; the key part is line 140: residuals(lm(tc.raw ~ seq(tc.raw) + Legendre(x=seq(tc.raw), n=1) + Legendre(x=seq(tc.raw), n=2))). The two Legendre function terms (n=1 and n=2) are to match -polort 2; more (or fewer) terms can be included as desired.

afni: 3dTstat with not-censored timepoints only

In a few recent posts I've shown images of the mean and standard deviation (calculated across time for each voxel), for QC tests. These are easy to calculate in afni (example here), but the 3dTstat command I used includes all timepoints (TRs), unless you specify otherwise. As described previously, we've been using a threshold of FD > 0.9 for censoring high-motion frames before doing GLMs. Thus, I wanted to calculate the mean and standard deviation images only including frames that were not marked for censoring (i.e., restrict the frames used by 3dTstat). This was a bit of a headache to code up, so R and afni code are after the jump, in the hopes it will be useful for others.

voxelwise standard deviation at different MB/SMS acquisitions

I've previously posted about using standard deviation and mean images to evaluate fMRI data quality, and practiCAL fMRI has a very useful series of posts explaining how various artifacts look in these images. In this post I use HCP (MB/SMS 8, RL/LR, customized scanner), and some images from one of our studies (MB4 and MB8, AP/PA, Siemens Prisma scanner); see links for acquisition details, and afni commands to generate these images is after the jump.

These tests were inspired by the work of Benjamin Risk, particularly his illustrations of banding and poor middle-of-the-brain sensitivity in the residuals (e.g., slides 28 and 29 of his OHBM talk). He mentioned that you can see some of the same patterns in simple standard deviation (stdev) images, which are in this post.

Here is a typical example of what I've been seeing with raw (not preprocessed) images. The HCP and MB8 scans are of the same person. The HCP scan is of the WM_RL task; the other two are of one of our tasks and the runs are about 12 minutes in duration (much longer than the HCP task). These MB8 and MB4 runs have low overall motion.

 Looking closely, horizontal bands are visible in the HCP and MB8 images (marked with green arrows). None of the MB4 images I've checked have such banding (though all the HCP ones I've looked at do); sensible anatomy (big vessels, brain edges, etc.) is brightest at MB4, as it should be.

Here are the same runs again, after preprocessing (HCP pipelines minimal preprocessing pipelines), first with all three at a low color threshold (800), second, with all three at a high color threshold (2000).

The HCP run is "washed out" at the 800 threshold with much higher standard deviation in the middle of the brain. Increasing the color scaling makes some anatomy visible in the HCP stdev map, but not as much as the others, and with a hint of the banding (marked with green arrows; easier to see in 3d). The MB4 and MB8 runs don't have as dramatic a "wash out" at any color scaling, with more anatomic structure visible at MB8 and especially MB4.The horizontal banding is still somewhat visible in the MB8 run (marked with an arrow), and the MB4 run has much lower stdev at the tip of the frontal lobe (marked with a green arc). tSNR versions of these images are below the jump.

These patterns are in line with Todd et. al (2017), as well as Benjamin Risk's residual maps. I'm encouraged that it looks like we're getting better signal-to-noise with our MB4 scans (though will be investigating the frontal dropout). Other metrics (innovation variance? GLM residuals?) may be even more useful for exploring these patterns. I suspect that some of the difference between the HCP and MB8 runs above may be due to the longer MB8 runs, but haven't confirmed.

getting afni sub-bricks by name in knitr

We've been using afni for GLMs, and using knitr to generate various summary documents (my demo for plotting NIfTI images in knitr is here).

A very nice feature of afni is that the 4th dimension in the statistics output files (even NIfTIs) can be assigned labels. The labels are viewed by the afni command 3dinfo -verb filename.nii.gz, which prints a list including the "sub-brick" (afni-speak for the 4th dimension) names and (0-based) indices, such as:

  -- At sub-brick #102 'Bng#5_Coef' datum type is float:   -13.6568 to    14.6825  
  -- At sub-brick #103 'Bng#5_Tstat' datum type is float:   -5.06214 to    4.55812 statcode = fitt; statpar = 1153  
  -- At sub-brick #104 'Bng#6_Coef' datum type is float:   -9.25095 to    13.3367  
  -- At sub-brick #105 'Bng#6_Tstat' datum type is float:   -4.56816 to    3.55517 statcode = fitt; statpar = 1153  
  -- At sub-brick #106 'Bng_Fstat' datum type is float:      0 to    5.56191 statcode = fift; statpar = 7 1153  

What if I want to display the Bng_Fstat image? It's possible to hard-code the corresponding number (107, in this case: R is 1-based) in the knitr R code, but that's dangerous, as the indexing could change. But we can run 3dinfo -label2index label dataset from within knitr to look up the indexes by name:

 img <- readNIfTI(fname, reorient=FALSE);  
 brick.num <- as.numeric(system2("/usr/local/pkg/linux_openmp_64/3dinfo", args=paste("-label2index Bng_Fstat", fname), stdout=TRUE));  
 if (!is.na(brick.num)) { make.plots(in.img=img[,,,1,brick.num+1]); }  

That the make.plots function is from my demo, and I provide the full path to 3dinfo in the system2 call to avoid command-not-found issues. Also note the indexing: afni-created images are often read by oro.nifti with an extra dimension (the sub-bricks are in the 5th dimension rather than the 4th). 
dim(img)
[1]  75  90  75   1 152
Also, note again that it's brick.num+1: afni has a sub-brick #0, but R does not have an index 0!

resampling images with afni 3dresample

This is the third post showing how to resample an image: I first covered SPM, then wb_command. afni's 3dresample command is the shortest version yet. I highly suggest you verify that the resampling worked regardless of the method you choose; the best way I know of is simply visually checking landmarks, as described at the end of the SPM post.

The setup is the same as the previous examples:

• inImage.nii.gz is the image you want to resample (for example, the 1x1x1 mm anatomical image)
• matchImage.nii.gz is the image with the dimensions you want the output image to have - what inImage should be transformed to match (for example, the 3x3x3 mm functional image)
• outImage.nii.gz is the new image that will be written: inImage resampled to match matchImage.

The command just lists those three files, with the proper flags:

 3dresample -master matchImage.nii.gz -prefix outImage.nii.gz -inset inImage.nii.gz  

Not being especially linux-savvy (afni does not play nice with Windows, so I run it in NeuroDebian), I sometimes get stuck with path issues when running afni commands. When in doubt, navigate to the directory in which the afni command programs (3dresample, in this case) are located, and execute the commands from that directory, typing./ at the start of the command. You can specify the full paths to images if needed; the example command above assumes 3dresample is on the path and the images are in the directory from which the command is typed. If you're using NeuroDebian, you can add afni to the session's path by typing . /etc/afni/afni.sh at the command prompt before trying any afni commands.

demo: clustering in afni with Clusterize

My go-to program for separating clusters out of an image is the Clusterize routine in AFNI. This little tutorial steps you through getting a NIfTI image into AFNI, using Clusterize, then getting a NIfTI out again. A word of warning: be sure to check laterality in the post-Clusterize NIfTI; sometimes things get flipped when you use multiple analysis programs. Also, I have a Windows box, so run AFNI within NeuroDebian (you should, too, especially if you run Windows), as the screenshots and notes below reflect. 

First, you need to get your NIfTI image into AFNI. Since I use NeuroDebian I started by putting the NIfTI I want to open into the for_afni subdirectory of the shared folder. Then you need to tell AFNI which directory to find the images in, which you do by clicking the Read button in the DataDir window (top red arrow). The Read Session window appears (right side of the screenshot, and, since I'm in NeuroDebian, I find my for_afni subdirectory under /home/brain/host/. Clicking the Set button (bottom red arrow) makes AFNI look for images in the directory.

Now we need to display the image that we want to clusterize. The image needs to be loaded as an OverLay, but AFNI is happiest if it has both an UnderLay and the OverLay, which are loaded via the circled buttons. Clicking the UnderLay button brings up a list of images, from both the for_afni subdirectory (since it was Read in the previous step) and standard anatomies (in my installation). In the screenshot I picked a standard anatomy; it also works if you use the overlay for the underlay (but you need something for the underlay). Then click the OverLay button and select the image you want to cluster. After setting both images you should see colored blobs on top of a greyscale background image: the colored image (the OverLay) will be the one clustered. Then click the Define OverLay button (arrow) to bring up the display shown in the upper right corner of the screenshot.

Next, set the threshold so that only the voxels you want to cluster are shown. Here, my overlay image consists of integers, and I want to identify clusters of at least 10 voxels with the value of 6 or higher. The screenshot shows how to set the threshold of 6: I put the little ** dropdown menu to 1 so that the values in the color slider are the actual numerical values (rather than a statistic). Next, I uncheck the autoRange box, also since I want the slider to be the actual numerical values. Finally, I move the slider (top arrow) to be exactly at 6 (use the up and down arrows for fine-tuning). The overlay changed as I moved the slider: now only voxels with values of 6 or larger are shown, and the overlay color scaling shifted.

Now we can do the clustering. First, click the Clear button (circled), in case any previous clustering is still in memory. Then click the Clusterize button (also circled), which brings up the menu dialog box (shown at left). Adjust the NN level and Voxels boxes to match your clustering parameters; the screenshot is set to find clusters with at least 10 voxels, and voxels must share a side to be in the same cluster. Click the Set button to close the menu dialog box. The main display won't change, except that the Rpt button (circled) will be enabled.

Clicking the Rpt button brings up the AFNI Cluster Results dialog box, as shown here. The display shows that AFNI found 8 clusters in my mask, ranging in size from 277 to 10 voxels. The coordinates are shown for the peak voxel in each cluster (since the XYZ dropdown is set to Peak), and clicking the Jump button in each row changes the coordinates in the display accordingly. To save the clustered version of the image, type a name into the box to the left of the SaveMsk button (marked with an arrow), then click the SaveMsk button. It doesn't look like anything happened, but there should now be a pair of images in the AFNI output directory (/brain/ by default in NeuroDebian) starting with the name specified.

Last, we need to convert the clustered mask back in NIfTI. I do this at the command line with 3dcopy. Not liking to mess about with configuration files, when I first open up the terminal I run . /etc/afni/afni.sh so that it can find 3dcopy. In the screenshot the terminal window opened up in the /brain/ directory, which is where the AFNI files are, so running 3dcopy outImage_mask+tlrc outImage.nii.gz writes the NIfTI file in /brain/ as well. Then I copy-paste outImage.nii.gz into /host/for_afni/ so that I can get to the file in Windows.

None of these steps are particularly difficult, but navigating back and forth can be a bit tricky, and the steps need to be done in the proper order. Good luck!


UPDATE 19 November 2015: You can also use afni at the command prompt for clustering; try the different options of 3dmerge and 3dmask_tool.