How about using 3dcalc in the way you say to create a copy of the original for all included voxels? Then create a summed image with 3dTstat -sum. Then use 3dcalc to divide by the total 1's in your censor file.
3dcalc -a data+orig -b 'templist.1D' -expr 'a*b' -prefix 'censordata'
3dTstat -sum -prefix testsum censordata+orig
set nvals = `awk 'BEGIN {sum=0}; {sum=sum+$1}; END {print sum};' templist.1D`
3dcalc -a testsum+orig -expr "a/$nvals" -prefix censored_avg
rm testsum+orig.*
rm censordata+orig.*
There's no doubt more than one way to do this, but there you go.
