Actually, I've figured out another way to do roughly the same thing. Try this in your script.

# 3dclust with whatever parameters you typically use
3dclust -nosum ... stats+orig > clust.1D
# count the number of rows actually containing the clusters
@ endrow = `wc clust.1D | awk '{print $1}'` - 11
# create a simple series of numbers in a row
count 1 $endrow > serialrow.1D
# transpose the row to column
1dtranspose serialrow.1D > serialcol.1D
# extract just the coordinates for the centers of mass for the clusters
# and add the serial column as an index for each cluster
# (only columns 1 to 3 starting in row 9 up to the end row)
1dcat 'clust.1D[1..3]{9..$}' serialcol.1D > clustcoords.1D
# create new dataset with spheres with a radius of 5 mm here around
# each coordinate
3dUndump -prefix clustspheres -master stats+orig -srad 5 -xyz clustcoords.1D
# Now compute the means for the spheres applied to the original dataset
3dROIstats -mask clustspheres+orig stats+orig > clustsphere_means.1D