GSReg_demo

#!/bin/tcsh -f
set sdir = $PWD
set stat = 0
set pname = `basename $0`

SETUP:
set gngain = 4.0  #base nz gain factor
set gnr = (1.0 2.0 1.0 2.0) #nz gain factor per group, 30 subjects per group
set gbr = (1.0 1.0 2.0 2.0) #background gain factor for each group
set n_grps = $#gbr
set imx = `ccalc -i "$n_grps * 30 -1"`  #-i以交互方式运行

set sext = .Rnd
set no_corr_in_bkg = 0  #if 1 don't compute correlations in background region
set regen = 0     #if 1 then recreate existing data files
#if 0 then reuse when possible. Set regen to 1 when
#Changing parameters under SETUP section or you'll
#end up looking at previously generated data
set killall = 1   #if 1 (recommended) then close then reopen all
#exisiting AFNI/3dGroupInCorr when script
#is rerun
set all_auto_seeds = 0  #if 1 then automatically go through all seed regions
#and save results. Otherwise, just stop at 3st seed
set wdir = /Volumes/afni/pub/dist/edu
goto PARSE
RET_PARSE:

SKIP:
#goto SHOWME2
#goto SHOWME_GRP_DIFF

CREATE_LATTICE:
foreach sub (`count -digits 2 0 $imx`)
set ig = `ccalc -fint $sub/30+1`
set sout = s$sub$sext

if (-f $sout/msk+tlrc.HEAD && $regen == 0) then
echo "Reusing lattice for $sout"
goto NEXT_LATTICE
endif

rm -rf ./$sout
if (! -d ${sout}) mkdir ${sout}
cd ${sout}

echo "Creating subject $sub"
echo "10 10 1 2 10" > udm.1D
echo "10 50 1 3 10" >> udm.1D
echo "40 30 1 4 20" >> udm.1D

3dUndump -dimen 64 64 2 \
-dval 1  \
-ijk  \
-prefix tmplt  \
-overwrite  \
udm.1D
3drefit -view +tlrc tmplt+orig

3dcalc -a tmplt+tlrc    \
-expr '(1-step(a))+a'  \
-overwrite  \
-prefix msk
cd -
NEXT_LATTICE:
end

CREATE_RAND_DATA:
foreach sub (`count -digits 2 0 $imx`)
set ig = `ccalc -fint $sub/30+1`
set sout = s$sub$sext

if (-f $sout/rsort+tlrc.HEAD && $regen == 0) then
echo "Reusing data for $sout"
goto NEXT_DATA
endif

cd ${sout}

echo "Creating subject $sub"
1deval -num 300 -dt 2.0 -expr 'gran(0,1)+(i-i)+0.7*z' > ARsig1.$sub.1D
1deval -num 300 -dt 2.0 -expr 'gran(0,1)+(i-i)+0.7*z' > ARsig2.$sub.1D
1deval -num 300 -dt 2.0 -expr 'gran(0,1)+(i-i)+0.7*z' > ARsig3.$sub.1D
1deval -num 300 -dt 2.0 -expr 'gran(0,1)+(i-i)+0.7*z' > ARsig4.$sub.1D

#Background sharing
set fn = `ccalc "$gngain * $gnr[$ig]"`   #Noise gain
set fb = $gbr[$ig]   #Background gain
3dcalc   -dt 2.0 \
-a msk+tlrc  \
-b ARsig1.$sub.1D \
-c ARsig2.$sub.1D \
-d ARsig3.$sub.1D \
-e ARsig4.$sub.1D \
-expr " step(a)     * (b*$fb) +\
($fn*gran(0,1)) +\
equals(a,2) * (c ) +\
equals(a,3) * (d ) +\
equals(a,4) * (e ) "\
-overwrite \
-prefix rs
3dROIstats -quiet -mask msk+tlrc. rs+tlrc.BRIK > rall.1D
3dmaskave -quiet rs+tlrc. > gs.1D
3dmaskave -quiet -mask msk+tlrc -mrange 1 1 rs+tlrc. > r1.1D
@ROI_Corr_Mat -ts rs+tlrc -roi msk+tlrc -prefix xcm.rs

#Background + Some cross talk
3dcalc   -dt 2.0 \
-a msk+tlrc  \
-b ARsig1.$sub.1D \
-c ARsig2.$sub.1D \
-d ARsig3.$sub.1D \
-e ARsig4.$sub.1D \
-expr " step(a)     * (b*$fb) +\
($fn*gran(0,1)) +\
equals(a,2) * (c ) +\
equals(a,3) * ((d+c)/2.0) +\
equals(a,4) * (e ) "\
-overwrite \
-prefix rs_23
3dROIstats -quiet -mask msk+tlrc. rs_23+tlrc.BRIK > rall_23.1D
3dmaskave -quiet rs_23+tlrc. > gs_23.1D
3dmaskave -quiet -mask msk+tlrc -mrange 1 1 rs_23+tlrc. > r1_23.1D
@ROI_Corr_Mat -ts rs_23+tlrc -roi msk+tlrc -prefix xcm.rs_23

#No background sharing
3dcalc   -dt 2.0 \
-a msk+tlrc  \
-b ARsig1.$sub.1D \
-c ARsig2.$sub.1D \
-d ARsig3.$sub.1D \
-e ARsig4.$sub.1D \
-expr " equals(a,1) * (b*$fb) +\
($fn*gran(0,1)) +\
equals(a,2) * (c) +\
equals(a,3) * (d) +\
equals(a,4) * (e) "\
-overwrite \
-prefix rsort
3dROIstats -quiet -mask msk+tlrc. rsort+tlrc.BRIK > rallort.1D
3dmaskave -quiet rsort+tlrc. > gsort.1D
3dmaskave -quiet -mask msk+tlrc -mrange 1 1 rsort+tlrc. > r1ort.1D
@ROI_Corr_Mat -ts rsort+tlrc -roi msk+tlrc -prefix xcm.rsort

cd $sdir
NEXT_DATA:
end
goto GINCOR

GINCOR:
BANDPASS:
foreach sub (`count -digits 2 0 $imx`)
set ig = `ccalc -fint $sub/30+1`
set sout = s$sub$sext

if ( -f $sout/rs_bp_23.r1+tlrc.HEAD && $regen == 0) then
echo "Reusing filtered for $sout"
goto NEXT_PASS
endif

cd ${sout}
echo "Bandpassing subject $sout ... "
3dBandpass  -overwrite -blur 2.0 -band 0.01 0.10 -prefix rs_bp \
-nodetrend -notrans -quiet   -input rs+tlrc
3dBandpass  -overwrite -blur 2.0 -band 0.01 0.10 -prefix rs_bp.gs \
-nodetrend -notrans -quiet   -input rs+tlrc -ort gs.1D
3dBandpass  -overwrite -blur 2.0 -band 0.01 0.10 -prefix rs_bp.r1 \
-nodetrend -notrans -quiet   -input rs+tlrc -ort r1.1D

3dBandpass  -overwrite -blur 2.0 -band 0.01 0.10 -prefix rs_bport \
-overwrite -nodetrend -notrans -quiet   -input rsort+tlrc
3dBandpass  -overwrite -blur 2.0 -band 0.01 0.10 -prefix rs_bport.gs \
-nodetrend -notrans -quiet   -input rsort+tlrc -ort gsort.1D
3dBandpass  -overwrite -blur 2.0 -band 0.01 0.10 -prefix rs_bport.r1 \
-nodetrend -notrans -quiet   -input rsort+tlrc -ort r1ort.1D

3dBandpass  -overwrite -blur 2.0 -band 0.01 0.10 -prefix rs_bp_23 \
-overwrite -nodetrend -notrans -quiet   -input rs_23+tlrc
3dBandpass  -overwrite -blur 2.0 -band 0.01 0.10 -prefix rs_bp_23.gs \
-nodetrend -notrans -quiet   -input rs_23+tlrc -ort gs_23.1D
3dBandpass  -overwrite -blur 2.0 -band 0.01 0.10 -prefix rs_bp_23.r1 \
-nodetrend -notrans -quiet   -input rs_23+tlrc -ort r1_23.1D

cd $sdir

NEXT_PASS:
end

ASSEMBLE:
set mopt = ''
if ($no_corr_in_bkg == 1) then
if (!(-f s00$sext/msk_nobkg+tlrc.HEAD) || $regen == 1) then
3dcalc -a s00$sext/msk+tlrc -expr 'a*step(a-1)' -overwrite \
-prefix s00$sext/msk_nobkg
endif
set mopt = "-mask s00$sext/msk_nobkg+tlrc"
endif

foreach ig (`count -digits 1 1 $n_grps`)
if ($ig == 1) then
set gflg = ''
set sel = ("s[012]?$sext/rs_bp+tlrc.HEAD")
else if ($ig == 2) then
set gflg = 'G2.'
set sel = ("s[345]?$sext/rs_bp+tlrc.HEAD")
else if ($ig == 3) then
set gflg = 'G3.'
set sel = ("s[678]?$sext/rs_bp+tlrc.HEAD")
else if ($ig == 4) then
set gflg = 'G4.'
set sel = ("s[9]?$sext/rs_bp+tlrc.HEAD" "s1[01]?$sext/rs_bp+tlrc.HEAD")
else
echo "ig = $ig is no good"
goto BEND
endif

if ( -f ${gflg}bp_23$sext.r1.grpincorr.data && $regen == 0) then
echo "Reusing GIC dataset ${gflg}bp_23$sext.*.grpincorr.data"
goto NEXT_SET
endif

set seln = (`echo "$sel"`)
3dSetupGroupInCorr -overwrite $mopt -prefix ${gflg}bp$sext \
$seln
set seln = (`echo "$sel" | sed 's/_bp/_bp.gs/'`)
3dSetupGroupInCorr -overwrite $mopt -prefix ${gflg}bp$sext.gs \
$seln
set seln = (`echo "$sel" | sed 's/_bp/_bp.r1/'`)
3dSetupGroupInCorr -overwrite $mopt -prefix ${gflg}bp$sext.r1 \
$seln

set seln = (`echo "$sel" | sed 's/_bp/_bport/'`)
3dSetupGroupInCorr -overwrite $mopt -prefix ${gflg}bport$sext \
$seln
set seln = (`echo "$sel" | sed 's/_bp/_bport.gs/'`)
3dSetupGroupInCorr -overwrite $mopt -prefix ${gflg}bport$sext.gs \
$seln
set seln = (`echo "$sel" | sed 's/_bp/_bport.r1/'`)
3dSetupGroupInCorr -overwrite $mopt -prefix ${gflg}bport$sext.r1 \
$seln

set seln = (`echo "$sel" | sed 's/_bp/_bp_23/'`)
3dSetupGroupInCorr -overwrite $mopt -prefix ${gflg}bp_23$sext \
$seln
set seln = (`echo "$sel" | sed 's/_bp/_bp_23.gs/'`)
3dSetupGroupInCorr -overwrite $mopt -prefix ${gflg}bp_23$sext.gs \
$seln
set seln = (`echo "$sel" | sed 's/_bp/_bp_23.r1/'`)
3dSetupGroupInCorr -overwrite $mopt -prefix ${gflg}bp_23$sext.r1 \
$seln

NEXT_SET:
end
goto SHOWME2

SHOWME2:
set vols = (bp$sext     bp$sext.gs     bp$sext.r1 \
bport$sext  bport$sext.gs  bport$sext.r1  \
bp_23$sext  bp_23$sext.gs  bp_23$sext.r1 )
set ccl = ( Mix$sext    Mix$sext.gs    Mix$sext.r1 \
Ort$sext    Ort$sext.gs    Ort$sext.r1 \
C23$sext    C23$sext.gs    C23$sext.r1)

if ($killall) then
@Quiet_talkers -npb_range 1 9
endif

set cnt = 1
foreach cc ($ccl)
set pid = `ps -a | \grep afni | \grep niml | \grep s00 | \grep "npb $cnt"`
if ($status) then
afni -npb $cnt -niml -yesplugouts \
-tbar $cc   \
-layout layout.$cnt s00$sext/ >& /dev/null &
sleep 2
plugout_drive  -npb $cnt -maxwait 20 \
-com 'SET_PBAR_NUMBER A.10' \
-com 'SET_PBAR_SIGN A.-' \
-com 'SET_THRESHNEW A.0.0001 *p' \
-com 'SET_FUNC_RANGE A.3' \
-com 'SET_IJK 40 30 1'  \
-quit >& /dev/null &
endif

set pid = `ps -a | \grep 3dGroupInCorr | \grep setA | \
\grep grpincorr.niml | \
\grep "npb $cnt"`
if ($status) then
echo "Plan $cc"
echo "   3dGroupInCorr -npb $cnt -setA $vols[$cnt].grpincorr.niml >& /dev/null"
3dGroupInCorr -npb $cnt -setA $vols[$cnt].grpincorr.niml >& /dev/null &
endif
@ cnt ++
end

INSTA_SETS:
set seedl = ()
foreach loc (0 1 2 3)
if ($loc < 3) then
set IJK = (`1dcat s00$sext/udm.1D"[0..2]{$loc}"`)
set seedloc = "Reg_`ccalc -i $loc + 2`"
else
set IJK = ( 10 30 1 ) #in background
set seedloc = "Backg"
endif
set seedl = ($seedl $seedloc)
set cnt = 1
foreach cc ($ccl)
plugout_drive  -npb $cnt \
-com "SET_IJK $IJK"  \
-com 'INSTACORR SET' \
-com 'SET_THRESHNEW A.0.0001 *p' \
-com 'SET_FUNC_RANGE A.3' \
-quit
plugout_drive  -npb $cnt \
-com "SET_XHAIRS A.OFF" \
-com "SAVE_PNG A.axialimage $cc.$seedloc.png" \
-com "SET_XHAIRS A.MULTI" \
-quit
@ cnt ++
end
if ($all_auto_seeds == 0) goto INFO
end

ASSEMBLE_IMAGES:
foreach seedloc ($seedl)
\rm $seedloc$sext.ppm
imcat    -matrix 3 3 -prefix $seedloc$sext.ppm -gap 1 -gap_col 0 0 0  \
Mix$sext.*$seedloc.png \
Ort$sext.*$seedloc.png \
C23$sext.*$seedloc.png
end

INFO:
set l = `prompt_user -pause \
"Each window shows one-sample ttest on correlations with the seed location.\n\
To create a new map:\n\
Move to a seed location\n\
Right-click and select 'InstaCorr set'\n\
or\n\
Shift+ctrl+right-click (and drag if you like).\n\
(To control whether or not the background region is masked, \n\
use option -mask_background 1 with $pname.)\n\
\n\
The windows are arranged to show 9 cases. Each row has a slightly different \n\
model, and each column process the data differently.\n\
Row 1: Uncorrelated regions + shared background\n\
Row 2: Uncorrelated regions + No shared background\n\
Row 3: Regions 2&3 correlated + shared background\n\
Col 1: No Global Signal Regression\n\
Col 2: Global Signal Regression\n\
Col 3: Background Signal Regression\n\
\n\
Each of these windows is a separate afni session and communicates with a \n\
separate 3dGroupInCorr program. To bring up the AFNI control window for one\n\
of the images, right-click on the 'Disp' button in the lower left corner\n\
of the image window. The title bars can also be used to match image window \n\
to AFNI controller.\n\
\n\
Press OK to close all and move to two sample results with different SNR"`

if ($l == 0) goto BEND
if ($killall == 0) then
#Must kill when going to SHOWME_GRP_DIFF from here
@Quiet_talkers -npb_range 1 9
endif

SHOWME_GRP_DIFF:
if ($killall) then
@Quiet_talkers -npb_range 1 9
endif
set cnt = 1
set g1 = ( bp$sext     bp$sext.gs  \
bp$sext     bp$sext.gs  \
bp$sext     bp$sext.gs )
set ccl = ( G1-G2.+Noise.Mix.$sext    G1-G2.+Noise.Mix$sext.gs  \
G1-G3.+Back.Mix$sext    G1-G3.+Back.Mix$sext.gs  \
G1-G4.+Noise+Back.Mix$sext    G1-G4.+Noise+Back.Mix$sext.gs)
set lays = ( 1 2 4 5 7 8)
foreach cc ($ccl)
if ($cnt < 3) then
set gg = G2
else if ($cnt < 5) then
set gg = G3
else if ($cnt < 7) then
set gg = G4
endif
set pid = `ps -a | \grep afni | \grep niml | \grep s00 | \grep "npb $cnt"`
if ($status) then
afni  -npb $cnt -niml -yesplugouts -tbar $cc   \
-layout layout.$lays[$cnt] s00$sext/ \
>& /dev/null &
sleep 2
plugout_drive  -npb $cnt -maxwait 20 \
-com 'SET_PBAR_NUMBER A.10' \
-com 'SET_PBAR_SIGN A.-' \
-com 'SET_THRESHNEW A.0.0001 *p' \
-com 'SET_FUNC_RANGE A.3' \
-com 'SET_IJK 40 30 1'  \
-quit >& /dev/null
endif

set pid = `ps -a | \grep 3dGroupInCorr | \grep setA | \
\grep grpincorr.niml | \
\grep "npb $cnt"`
if ($status) then
echo "Plan $cc"
echo "   3dGroupInCorr -npb $cnt \
-setA $g1[$cnt].grpincorr.niml \
-setB $gg.$g1[$cnt].grpincorr.niml \
>& /dev/null"
3dGroupInCorr -npb $cnt    \
-setA $g1[$cnt].grpincorr.niml \
-setB $gg.$g1[$cnt].grpincorr.niml \
>& /dev/null &
endif
set loc = 2
set IJK = (`1dcat s00$sext/udm.1D"[0..2]{$loc}"`)
set seedloc = "Reg_`ccalc -i $loc + 2`"
plugout_drive  -npb $cnt \
-com "SET_IJK $IJK"  \
-com 'INSTACORR SET' \
-com 'SET_THRESHNEW A.0.0001 *p' \
-com 'SET_FUNC_RANGE A.3' \
-quit >& /dev/null
sleep 1
plugout_drive  -npb $cnt \
-com 'SET_FUNC_RANGE A.3' \
-com 'SET_THRESHNEW A.0.0001 *p' \
-quit >& /dev/null

if ($all_auto_seeds == 1) then
set icon = 1
set con = (g1-g2 g1 g2)
foreach sb ( 0 2 4 )
set sb1 = `ccalc -i $sb + 1`
plugout_drive  -npb $cnt \
-com "SET_SUBBRICKS A -1 $sb $sb1" \
-quit >& /dev/null
plugout_drive  -npb $cnt \
-com 'SET_THRESHNEW A.0.0001 *p' \
-quit >& /dev/null
plugout_drive  -npb $cnt \
-com "SET_XHAIRS A.OFF" \
-com "SAVE_PNG A.axialimage $cc.$con[$icon].$seedloc.png" \
-com "SET_XHAIRS A.MULTI" \
-quit >& /dev/null
@ icon ++
end
endif

@ cnt ++
end

INFO2:
set l = `prompt_user -pause \
"Each window shows two-sample ttest results on correlations from two groups\n\
at the seed location. There are no differences between the group pairs, except\n\
for levels of noise and/or shared background fluctuations.\n\
To create a new map:\n\
Move to a seed location\n\
Right-click and select 'InstaCorr set'\n\
or\n\
Shift+ctrl+right-click (and drag if you like).\n\
\n\
The windows are arranged to show 6 cases. Each row is a contrast between two\n\ groups. The 2nd column is processed with Global Signal Regression prior to \n\
the two-sample ttest.\n\
Row 1: G1 vs. G2, G2 has twice the noise as G1\n\
Row 2: G1 vs. G3, G3 has twice the background fluctuations as G1\n\
Row 3: G1 vs. G4, G4 has twice the background fluctuations and noise as G1\n\
Col 1: No Global Signal Regression\n\
Col 2: Global Signal Regression\n\
\n\
Press OK to close all and end demo. Cancel to keep everything up."`

if ($l == 0) goto BEND
@Quiet_talkers -npb_range 1 9

goto END

PARSE:
set Narg = $#
set cnt = 1
while ($cnt <= $Narg)

if ( "$argv[$cnt]" == "-echo") then
set echo
goto NEXT
endif

if ( ("$argv[$cnt]" == "-h" || "$argv[$cnt]" == "-help") ) then
goto HELP
endif

if ("$argv[$cnt]" == "-recreate") then
set regen = 1
goto NEXT
endif

if ("$argv[$cnt]" == "-webit") then
if (-f ./$pname) then
echo ""
echo "Enter c if you want to copy ./$pname to: $wdir"
echo ""
set jnk = $<
if ($jnk == 'c' || $jnk == 'C') then
cp -p ./$pname $wdir
if ($status) then
echo "Failed to copy ./$pname to: $wdir"
goto BEND
endif
goto END
endif
else
echo ""
echo "$pname is not in the current directory."
echo "This option should be run from the directory."
echo "where $pname resides."
echo "Perhaps: `which $pname` "
echo ""
goto BEND
endif
goto NEXT
endif

if ("$argv[$cnt]" == "-mask_background") then
set regen = 1
@ cnt ++
if ($cnt > $Narg) then
echo "Need 0 or 1 after -mask_background"
goto BEND
endif
set no_corr_in_bkg = `printf %d $argv[$cnt]`
goto NEXT
endif

if ("$argv[$cnt]" == "-noise_gain") then
set regen = 1
@ cnt ++
if ($cnt > $Narg) then
echo "Need a value after -noise_gain"
goto BEND
endif
set gngain = `printf %f $argv[$cnt]`
goto NEXT
endif

if (1) then
echo "Error: Option or parameter '$argv[$cnt]' not understood"
goto END
endif

NEXT:
@ cnt ++
end

goto RET_PARSE

HELP:
echo ""
echo "Usage: $pname "
echo "A script to illustrate the effects of Global Signal Regression"
echo "on a toy model of a brain made."
echo "The group correlation maps are generated interactively with"
echo "AFNI's 3dGroupInCorr"
echo "The script takes time to setup the first time it is run to generate the"
echo "data. Once the data are created the script launches AFNI and "
echo "3dGroupInCorr sessions for interactive amusement."
echo ""
echo "You can change default parameters by directly editing the section"
echo "SETUP in the beginning of the script, or by using the options below"
echo ""
echo "Options:"
echo " -mask_background MSK: MSK=0 => Do include background region in "
echo "                       corelation computations. "
echo "                       MSK=1 => Do not include background."
echo "                       Default is $no_corr_in_bkg"
echo " -noise_gain GN: Noise gain factor common to all groups."
echo "                 Default GN is $gngain"
echo ""
echo " -recreate: Recreate simulated data. Default is to reuse"
echo "            whatever is available, unless you use options"
echo "            -mask_background or -noise_gain on the command line"
echo ""
echo " -help: This message"
echo " -webit: Put copy of this script on afni website"
echo "         This option only works from Ziad's computer."
echo ""
echo "Questions, comments should be addressed to:"
echo "         saadz@mail.nih.gov"
echo "This script can be obtained from:"
echo "    http://afni.nimh.nih.gov/pub/dist/edu/@GSReg_demo" echo "and requires AFNI versions compiled after July 07/2011"
echo "If you have AFNI installed already, you can update with"
echo "    @update_binaries -defaults"
echo ""
goto END

BEND:
echo "$pname : Error, or User Interruption"
set stat = 1
goto END

END:
exit $stat