# make MPM atlases for the Desai MPM datasets
set dkmpm = "desikan_killiany.desai_prob.at_subcortical_TT_N27.bl6.max.nii.gz"
set dkmpmout = "TT_desai_dk_mpm"
set ddmpm = "destrieux.desai_prob.at_subcortical_TT_N27.bl6.max.nii.gz"
set ddmpmout = "TT_desai_dd_mpm"

set nimloutmpm = dkmpm.niml

# build maximum probability map atlases
rm $ddmpmout* $dkmpmout*
rm temp_*

# make these datasets as small as possible
setenv AFNI_COMPRESSOR GZIP

# make concatentated copy of dd and fs maps and dk with fs maps
3dTcat -prefix temp_tcat_dd_fs ~/abin/TT_desai_ddpmaps+tlrc. ~/abin/TT_desai_fspmaps+tlrc.
3dTcat -prefix temp_tcat_dk_fs ~/abin/TT_desai_dkpmaps+tlrc. ~/abin/TT_desai_fspmaps+tlrc.

# find all values that have a value greater than 0 to use for a mask
3dTstat -max -prefix temp_ts_dd_fs_max temp_tcat_dd_fs+tlrc
3dTstat -max -prefix temp_ts_dk_fs_max temp_tcat_dk_fs+tlrc

# apply the mask to the MPM atlas and put output in byte integers
# for data with larger than 255 ROIs, would need integers
# also threshold to a minimum of 5% probability
3dcalc -byte -nscale \
       -a $ddmpm \
       -b temp_ts_dd_fs_max+tlrc -overwrite -expr 'a*step(b-0.04999)' -prefix $ddmpmout
3dcalc -byte -nscale \
       -a $dkmpm \
       -b temp_ts_dk_fs_max+tlrc -overwrite -expr 'a*step(b-0.04999)' -prefix $dkmpmout



# use atlas_label_table attributes from probability maps to build MPM table
set dd_structlist = `3dAttribute -name ATLAS_LABEL_TABLE ~/abin/TT_desai_ddpmaps+tlrc. | \
                 grep STRUCT | awk -F= '{print $2}'`
set dk_structlist = `3dAttribute -name ATLAS_LABEL_TABLE ~/abin/TT_desai_dkpmaps+tlrc. | \
                 grep STRUCT | awk -F= '{print $2}'`
set fs_structlist = `3dAttribute -name ATLAS_LABEL_TABLE ~/abin/TT_desai_fspmaps+tlrc. | \
                 grep STRUCT | awk -F= '{print $2}'`

set mpmdset = $ddmpmout+tlrc
set nimloutmpm = ddmpm.niml
set si = 1
echo "<atlas_point_list"  > $nimloutmpm
echo ' ni_form="ni_group" >' >> $nimloutmpm

foreach struct ( $dd_structlist $fs_structlist )
  # MPM indices run from 1 to nv while PM indices run from 0 to nv-1
  # write out NIML point structure
  # the STRUCT field contains the structure name
  # the SB_LABEL field contains the sub-brick label that is associated with
  #  this structure
  set cm = `3dCM $mpmdset"<$si..$si>"`
  echo '<ATLAS_POINT' >> $nimloutmpm
  echo 'data_type="atlas_point"' >> $nimloutmpm
  echo 'STRUCT='$struct >> $nimloutmpm
  echo 'VAL='\"$si\" >> $nimloutmpm
  echo 'OKEY='\"$si\" >> $nimloutmpm
  echo 'GYoAR="0"' >> $nimloutmpm
  echo 'COG='\"$cm\" >> $nimloutmpm
  echo 'SB_LABEL='$struct' />' >> $nimloutmpm
  set si = `ccalc -int -expr "$si + 1"`
end
echo '</atlas_point_list>' >> $nimloutmpm

# 3dcopy $ddmpm $ddmpmout
3drefit -atrstring ATLAS_LABEL_TABLE file:$nimloutmpm \
  ${ddmpmout}+tlrc
3drefit -cmap INT_CMAP ${ddmpmout}+tlrc
3drefit -space TT_N27 ${ddmpmout}+tlrc


# now do the same for the dessikan-killiany MPM atlas
set mpmdset = $dkmpmout+tlrc
set nimloutmpm = dkmpm.niml
set si = 1
echo "<atlas_point_list" > $nimloutmpm
echo ' ni_form="ni_group" >' >> $nimloutmpm

foreach struct ( $dk_structlist $fs_structlist )
  # MPM indices run from 1 to nv while PM indices run from 0 to nv-1
  # write out NIML point structure
  # the STRUCT field contains the structure name
  # the SB_LABEL field contains the sub-brick label that is associated with
  #  this structure
  # compute the center of mass for each structure (note left/right are separate)
  set cm = `3dCM $mpmdset"<$si..$si>"`
  echo '<ATLAS_POINT' >> $nimloutmpm
  echo 'data_type="atlas_point"' >> $nimloutmpm
  echo 'STRUCT='$struct >> $nimloutmpm
  echo 'VAL='\"$si\" >> $nimloutmpm
  echo 'OKEY='\"$si\" >> $nimloutmpm
  echo 'GYoAR="0"' >> $nimloutmpm
  echo 'COG='\"$cm\" >> $nimloutmpm
  echo 'SB_LABEL='$struct' />' >> $nimloutmpm
  set si = `ccalc -int -expr "$si + 1"`
end
echo '</atlas_point_list>' >> $nimloutmpm

# 3dcopy $dkmpm $dkmpmout
3drefit -atrstring ATLAS_LABEL_TABLE file:$nimloutmpm \
  ${dkmpmout}+tlrc
3drefit -cmap INT_CMAP ${dkmpmout}+tlrc
3drefit -space TT_N27 ${dkmpmout}+tlrc
