#!/bin/tcsh -f

if ("$1" == "-clean") then
   \rm -f AVG.ts.*
   \rm -f lh.OccROIs.FULL.niml.dset
   \rm -f v2s.lh.TS.FULL.niml.dset
   \rm ALLv2*.ts.1D
   \rm lh.OccROIs.labeled.niml.dset
   \rm toylut.niml.cmap
   goto END
endif 

if ( ! -f v2s.lh.TS.niml.dset || \
     ! -f lh.OccROIs.niml.dset ) then
    echo "Missing v2s.lh.TS.niml.dset and/or lh.OccROIs.niml.dset" 
    echo "Make sure you first run:"
    echo " tcsh run_3dVol2Surf"
    echo " tcsh run_ROI2dataset"
    echo ""
    exit 1 
endif

#To make use of the maskdumping tools in AFNI , we should
#make sure we have a value for each node in both datasets
#and mask.
if ( ! -f v2s.lh.TS.FULL.niml.dset) then
   ConvertDset -i v2s.lh.TS.niml.dset  \
              -o_NIML -pad_to_node ld120 \
               -prefix v2s.lh.TS.FULL   
else
   echo "Reusing v2s.lh.TS.FULL.niml.dset"
endif              

if ( ! -f lh.OccROIs.FULL.niml.dset) then
   ConvertDset -i lh.OccROIs.niml.dset  \
              -o_NIML -pad_to_node ld120 \
               -prefix lh.OccROIs.FULL   
else
   echo "Reusing lh.OccROIs.FULL.niml.dset"
endif

#To get the average time series from each mask value:              
#in readable form
3dROIstats  -nzmean -nomeanout  \
            -mask lh.OccROIs.FULL.niml.dset  \
            v2s.lh.TS.FULL.niml.dset > AVG.ts.text

#in less readable but easily plottable form
3dROIstats  -nzmean -nomeanout -1Dformat \
            -mask lh.OccROIs.FULL.niml.dset  \
            v2s.lh.TS.FULL.niml.dset > AVG.ts.1D      

1dplot -title "Normalized average TS from each ROI" -one -norm2 AVG.ts.1D      &

#If you want a dump of all the values in ROI valued 2
3dmaskdump  -cmask "-a lh.OccROIs.FULL.niml.dset -expr equals(a,2)"   \
            -index -noijk -o ALLv2a.ts.1D \
            v2s.lh.TS.FULL.niml.dset

         #or
3dmaskdump  -mask lh.OccROIs.FULL.niml.dset   \
            -mrange 2 2 \
            -index -noijk -o ALLv2b.ts.1D \
            v2s.lh.TS.FULL.niml.dset
            
#Each row in ALLv2a.ts.1D contains the index of a node, followed by the 
#the timeseries at that node. To plot them all in one swoop
#with 1dplot we skip column 0, then transpose (\' thingy) the file because
#1dplot plots columns, not rows
1dplot -title "All time series from region 2" \
       -one ALLv2a.ts.1D'[1-$]'\' 

#Getting data from ROIs defined in the annotation files
set zone1 = G_cuneus
set zone2 = G_occipital_middle
3dcalc -a ../SurfData/SUMA/std.lh.aparc.a2005s.annot.niml.dset"<$zone1>" \
       -b ../SurfData/SUMA/std.lh.aparc.a2005s.annot.niml.dset"<$zone2>" \
       -expr 'step(a)+2*step(b)' \
       -overwrite -prefix auto.lh.from_annot.gii.dset 
3dROIstats  -nzmean -nomeanout -1Dformat \
            -mask auto.lh.from_annot.gii.dset   \
            v2s.lh.TS.FULL.niml.dset > AVG.ts.from_annot.1D  

#Turning a bunch of blobs into a labeled dataset
#  Say you want to label the 3 rois (values 1, 2, and 3) in lh.OccROIs.niml.dset
#  First we decide on the name and color for each ROI and write them into a
#  text file. Here I will use echo commands to create the file colut.txt in this
#  script. Notice that the first line is not necessary. See MakeColorMap -help
#  and its example 5 for more information."
echo "     #integer label    String Label      R    G    B    A" > colut.txt
echo "     1                Big_House         0.3  0.1  1    1" >> colut.txt
echo "     2                Small_Face        1    0.2  0.4  1" >> colut.txt
echo "     3                Electric          1    1    0    1" >> colut.txt
echo "     4                Atomic            0.1  1    0.3  1" >> colut.txt

#  Now turn the color lookup file into a label table
MakeColorMap -usercolutfile colut.txt -suma_cmap toylut

#  Create a labeled dataset from lh.OccROIs.niml.dset and this colorscheme
ConvertDset -i lh.OccROIs.niml.dset -o lh.OccROIs.labeled.niml.dset \
            -labelize toylut.niml.cmap
#  And that's it, now load both of lh.OccROIs.niml.dset
#and lh.OccROIs.labeled.niml.dset into SUMA, switch between them and notice
#the differences. A labeled dataset no longer allows you to set the colormapping,
#However, once loaded into SUMA and you click on a node in the labeled region
#you get the label of that region displayed in the viewer.
#Also you can extract regions by name from a labeled dataset as shown in the 
#examples using the annot files above

END: