#!/bin/tcsh

# makes AFNI compatible atlases from the Keuken MNI05 probabilistic datases

# where is the original Keuken data (assume current directory)
set origdir = `pwd`
set outatlasdir = ${origdir}/KeukenAtlasDist/

rm -rf $outatlasdir
mkdir $outatlasdir
@AfniEnv -set AFNI_SUPP_ATLAS_DIR $outatlasdir


foreach group ( young middle_age elderly )
   mkdir mpm_${group}
   cd mpm_${group}

   rm labellist.txt
   set labeli = 1
   foreach roiset ( ${origdir}/*_non_linear_${group}.nii.gz )
      set label = `basename -s _prob_mni_non_linear_elderly.nii.gz $roiset`
      echo $labeli $label >> labellist.txt
      @ labeli ++
   end

   # create map of maximum probabilities
   3dMean -max -overwrite -prefix probmax_nonlinear_${group}.nii.gz \
      ${origdir}/*_non_linear_${group}.nii.gz
   # each group's "probability map" is in terms of number of subject with differing numbers of subjects overall
   # so get the maximum number here 
   set maxvalue = `3dBrickStat -max -slow probmax_nonlinear_${group}.nii.gz`
   # scale and threshold map to at least 20%
   3dcalc -a probmax_nonlinear_${group}.nii.gz -exp "step(a/$maxvalue-0.2)" -overwrite \
      -prefix probmax_thrmask_nonlinear_${group}.nii.gz
   set dsetlist = ( ${origdir}/*_non_linear_${group}.nii.gz )
   # map maximum index as maximum probability
   3dTstat -argmax1 -prefix mpm_nonlinear_${group}.nii.gz -mask probmax_thrmask_nonlinear_${group}.nii.gz \
      -overwrite "$dsetlist"

   # a little modal smoothing is nicer here to fix alignment, drawing and probability thresholding issues
   3dLocalstat -stat mode -nbhd 'SPHERE(2)' -prefix mpm_nonlinear_${group}_mode.nii.gz \
       mpm_nonlinear_${group}.nii.gz

   # make formal atlas
   @Atlasize -space MNI -dset mpm_nonlinear_${group}.nii.gz -lab_file labellist.txt 1 0 -atlas_type  G \
       -atlas_name keuken_mpm_${group} -atlas_description  \
       "sub-cortical regions, Keuken, et al, https://doi.org/10.1016/j.neuroimage.2014.03.032"

   # make formal atlas - also modal smoothed atlas
   @Atlasize -space MNI -dset mpm_nonlinear_${group}_mode.nii.gz -lab_file labellist.txt 1 0 -atlas_type  G \
       -atlas_name keuken_mpm_mode_${group} -atlas_description    \
       "sub-cortical regions (modally smoothed 2mm),Keuken, et al, https://doi.org/10.1016/j.neuroimage.2014.03.032"

   # create surfaces for each region
   IsoSurface -isorois+dsets -o keuken_${group}.gii -input mpm_nonlinear_${group}.nii.gz \
      -Tsmooth 0.1 100 -overwrite
   cd ..

   #show how to use these surfaces
   echo "cd mpm_${group}; suma -i *.gii -onestate"
   mv mpm_${group} $outatlasdir
end
