Ideally, the following process would be quite excellent, not sure how hard it would be to build the code though.

STEPS:
1) use 3dmerge with the -1thresh and -1clust to correct for multiple comparisons (ie 3dmerge -thresh 3.00 -1clust 4 400) whatever parameters are required to do the correction

2) find all local extrema (peaks) greater than a specified value (z>3.00) upwards of a 1000 voxels

3) for peaks closer than specified radius (ie 13.5 mm) collapse the two regions based on the center of mass

4) finally, after all peaks are found that are further apart than the specified radius; make spheres of a specified radius around them. In the case that spheres overlap, allocate the voxels in the overlapping region to the nearest peak -- half will go to each region.

5) The output is regions that do not overlap and whose peak/center of mass are more than the specified radius apart. Then one could run 3dROIStats on those regions (would be nice to number them based on their ranking of the ir maximal intensity). Alternatively, one might want the option to mask the regions with the original data -- combine the original dataset and the new dataset only when both datasets are nonzero.

I know this is done with other analysis programs. Not sure if it has been or can be implemented into AFNI. Let me know if its possible or not. More details on the process can be provided.