#!/bin/tcsh
# rank_betas.csh
# ranks the betas in an impulse response function dataset

# modify the values here for your data
set roidset = mask+orig
set irespdset = iresp+orig
# only look at the first nvals sub-bricks
set nvals = 15
# top percentile rank of voxels
set topvoxels = 25

# make masked, integer version of iresp dset with integer values only
3dcalc -a $irespdset -b $roidset -expr 'int(1000 * a*step(b))' \
-prefix iresp_mask_scaled -overwrite -datum short -nscale

# for every sub-brick, compute ranking and extract top voxels
foreach sbi (`count -digits 2 0 $nvals `)
# extract sub-brick
3dTcat -prefix tempsb -overwrite iresp_mask_scaled+orig."[${sbi}]"

# compute of rank of voxels from low to high at each voxel
3dRank -prefix irank2_${sbi} -overwrite -input tempsb+orig

# get value that corresponds to "top" cutoff from rank mapping
set uvals = `cat irank2_${sbi}.rankmap.1D | wc -l`
@ uvals = $uvals - 3
set minval = `tail -${topvoxels} irank2_${sbi}.rankmap.1D |head -1`
set mv = $minval[1]

# get ranked mask of just the top values
3dcalc -overwrite -a irank2_${sbi}+orig -expr "a*step(a-${mv}+1)" \
-prefix top${topvoxels}_${sbi}
end

# put the rankings and the top voxel masks back together across sub-bricks
3dTcat -overwrite -prefix irank2_all irank2_??+orig.HEAD
3dTcat -overwrite -prefix top${topvoxels}_all top${topvoxels}_??+orig.HEAD