(All the hard work and the write up, compliments of Jen Johnson and Pamela Perschler)

You need to average both anatomy and functional data.
In both cases you want to use +tlrc (talariached data).

For ease, make a directory for your averaged files.
They may all need to be in the same directory to average them.

Anatomy
Copy 3D+tlrc.BRIK & .HEAD into your "average" directory
>3dmerge -doall 3d*+tlrc*
(where "3d*+tlrc* is using enough wild cards to get all the 3D anatomy files you want to average)

The output name is mrg+tlrc.BRIK and mrg+tlrc.HEAD by default, but you can use the -prefix flag to specify an output prefix (other than mrg) that you prefer.

Functionals
You are making a binary functional image file from your clustered, talairached data: Every active voxel (member of any cluster) will be called "1" and all nonactive voxels are called "0." 3dcalc can actually do a number of different algebraic manipulations on the data, see -expr below.

So really "average" is a misnomer here, because what you really do is to sum all of these "one/zero" files. Then every voxel that is in common across subjects will be given the value of the number of subjects in common.

How To Create Individual Binary "onezero" Files

3dcalc -prefix onezero_subj# -a filename -expr 'astep (a,0)'

Example:
>3dcalc -prefix onezero1234_5 -a deconblur1234_5+tlrc -expr 'astep (a,0)'

-prefix precedes the name of your new output file, onezero1234_5+tlrc.BRIK &HEAD

-a dname = Read dataset 'dname' and call the voxel values 'a' in the -expr that follows. 'a' may be any single letter from 'a' to 'z'.

If some letter name is used in the expression, but not present in one of the dataset options here, then that variable is set to 0.

If the letter is followed by a number, then that number is used to select the sub-brick of the dataset which will be used in the calculations. For example, '-b3 dname' specifies that the variable 'b' refers to sub-brick 3 of the dataset (indexes start at 0). N.B.: Another way to achieve the effect of '-b3' is described in the 'INPUT DATASET SPECIFICATION' section of the 3dcalc help.

-expr "expression"
Apply the expression within quotes to the input datasets,
one voxel at a time, to produce the output dataset.
("sqrt(a*b)" to compute the geometric mean, for example).

The following functions are designed to help implement logical
functions, such as masking of 3D volumes against some criterion:
step(x) = {1 if x>0 , 0 if x<=0},
astep(x,y) = {1 if abs(x) > y , 0 otherwise} = step(abs(x)-y)
rect(x) = {1 if abs(x)<=0.5, 0 if abs(x)>0.5},
bool(x) = {1 if x != 0.0 , 0 if x == 0.0},
and(a,b,...,c) = {1 if all arguments are nonzero, 0 if any are zero}
or(a,b,...,c) = {1 if any arguments are nonzero, 0 if all are zero}
mofn(m,a,...,c) = {1 if at least 'm' arguments are nonzero, 0 otherwise}
[These last 3 functions take a variable number of arguments.]

How to Sum the Binary Files

3dcalc -prefix name -a onezero_subj#1 -b onezero_subj#2 ... (list each) -expr 'a+b+c...'

Example:
>3dcalc -prefix sum1234post -a onezero1234_9+tlrc -b onezero1234_11+tlrc -c onezero1234_14+tlrc -d onezero1234_21+tlrc -e onezero1234_24+tlrc -f onezero1234_29+tlrc -expr 'a+b+c+d+e+f'

IN this command, you first specify an output name (-prefix sum1234post), then define the files you want to sum, by assigning a letter to each. These letters are used in the expression part of the command. The letters do not have to be the same as those specified when you created the individual binary files. There you can run that "onezero" command for each file, and keep calling each one "a," but here you have to call each file a different letter, it doesn't matter which letter, or which order, as long as the letters you use for each file match the letters you use in the -expr (expression) part of the command. For example if you want to sum 5 files, and you name them a,b,c,d,and e, then your expression would be "-expr 'a+b+c+d+e'"