Ding, dong.
There are several ways to accomplish this. Here's one that uses 3dcalc's 1D input applied slicewise. You may use 3dAutomask instead of 3dSkullStrip. Also 3dcalc on the single slice volumes and then 3dZcat'ing them all back together is one of the other ways.

#!/bin/tcsh

set anatdset = anat+orig
set lastslice = 174
rm -f temp_anat_ns+orig.*
3dSkullStrip -prefix temp_anat_ns -orig_vol -input $anatdset

rm -f medslices.1D
touch medslices.1D

echo computing median for each slice
foreach slice (`count 0 $lastslice `)
echo "Slice: $slice"
3dZcutup -keep $slice $slice -prefix tt_slice -overwrite temp_anat_ns+orig
set sl_med = `3dBrickStat -median -non-zero tt_slice+orig`
echo $sl_med[2] >> medslices.1D
end

3dcalc -a $anatdset -b K:medslices.1D'[0]' -ok_1D_text -expr 'a*1000/b' \
-prefix slicemedian2 -overwrite