Hi Soo Hyun,

Here's a short script where almost all the action happens in just two steps. If you are careful, you can actually do this mostly in a single 3dcalc command using relative indexing, but I think this is a little clearer and more flexible.

#!/bin/tcsh
# @edgy_rois
# make ROIs edgy by creating dataset of just the edges of the ROI
# usage:
# @edgy_rois mydset.nii
# output is edgy_rois.nii.gz

set dset = $1
set prefix = "edgy_rois"

# see if voxels are on edge by looking for any differences in a voxel neighborhood
3dLocalstat -stat stdev -nbhd 'SPHERE(-1.42)' -prefix temp_sd.nii.gz -overwrite \
$dset

# alternatively, use a bigger neighborhood, particularly for anisotropic voxels - thicker edges
#3dLocalstat -stat stdev -nbhd 'RECT(-1,-1,-1)' -prefix temp_sd.nii.gz -overwrite \
# $dset

# color the edges only by the ROI value
3dcalc -a $dset -b temp_sd.nii.gz -prefix $prefix.nii.gz -overwrite \
-expr 'a*step(b)'
3drefit -cmap INT_CMAP $prefix.nii.gz