1#!/bin/tcsh
2
3## help the pitiful luser?
4
5set dohelp = 0
6if( $#argv == 0 ) set dohelp = 1
7if( dohelp == 0 )then
8  if( $argv[1] == "-help" ) set dohelp = 1
9endif
10
11if( $dohelp )then
12  echo " "
13  echo "Usage: @banded_clustsim list_of_directories"
14  echo " "
15  echo "The input to this script is a list of afni_proc.py results directories."
16  echo "The steps followed are:"
17  echo " a) determine the collective passband for the stimuli from the"
18  echo "    X.nocensor.xmat.1D matrix file in each directory (program stimband)."
19  echo " b) bandpass each directory's pre-whitened residual file, and determine"
20  echo "    its spatial ACF = AutoCorrelation Function (program 3dFWHMx)"
21  echo " c) average the spatial ACF parameters across directories"
22  echo " d) use the resulting smoothness to estimate the cluster-size"
23  echo "    thresholds (program 3dClustSim)"
24  echo "The goal is to provide cluster significance levels for group tests."
25  echo " "
26  echo "** This script is currently in development and is experimental -- RWCox **"
27  exit 0
28endif
29
30## scan inputs and assemble the lists of files to process
31
32set nbad    = 0
33set matlist = ( )
34set errlist = ( ) ; set nerrWH = 0 ; set nerrOT = 0
35set msklist = ( ) ; set nmskGM = 0 ; set nmskGR = 0
36
37## We prefer whitened residuals (from 3dREMLfit -Rwherr), but
38## will process non-whitened residuals if that's all there is
39
40set errts_template_WH1   = 'errts.\*.WH+tlrc.HEAD'
41set errts_template_WH2   = 'errts_whitened.\*+tlrc.HEAD'
42set errts_template_OTHER = 'errts.\*+tlrc.HEAD'
43
44echo "========== Scanning for input files =========="
45
46foreach aaa ( $argv )
47
48  if( ! -d $aaa )then
49    echo "** ERROR: directory $aaa does not exist" ; @ nbad++ ; continue
50  endif
51
52  if( ! -f $aaa/X.nocensor.xmat.1D )then
53    echo "** ERROR: matrix $aaa/X.nocensor.xmat.1D does not exist" ; @ nbad++
54  else
55    set matlist = ( $matlist $aaa/X.nocensor.xmat.1D )
56  endif
57
58  unset mmm
59  if( -f $aaa/mask_GM_resam+tlrc.HEAD )then
60    set mmm = $aaa/mask_GM_resam+tlrc.HEAD ; @ nmskGM++
61  else if( -f $aaa/mask_group+tlrc.HEAD )then
62    set mmm = $aaa/mask_group+tlrc.HEAD    ; @ nmskGR++
63  endif
64  if( $?mmm == 0 )then
65    echo "** ERROR: neither mask $aaa/mask_GM_resam+tlrc.HEAD nor $aaa/mask_group+tlrc.HEAD is found"
66    @ nbad++
67  else
68    set msklist = ( $msklist $mmm )
69  endif
70
71  unset fff
72  if( $?fff == 0 )then
73    set eee = ( `find $aaa -name errts\*WH\*+tlrc.HEAD` )
74    if( $#eee > 0 )then
75      set fff = $eee[1]
76      @ nerrWH++
77    endif
78  endif
79
80  if( $?fff == 0 )then
81    set eee = ( `find $aaa -name errts\*whitened\*+tlrc.HEAD` )
82    if( $#eee > 0 )then
83      set fff = $eee[1]
84      @ nerrWH++
85    endif
86  endif
87
88  if( $?fff == 0 )then
89    set eee = ( `find $aaa -name errts\*+tlrc.HEAD` )
90    if( $#eee > 0 )
91      set fff = $eee[1]
92      @ nerrOT++
93    endif
94  endif
95
96  if( $?fff == 0 )then
97    echo "** ERROR: no file of form $errts'*'+tlrc.HEAD in directory $aaa"
98    @ nbad++
99  else
100    set errlist = ( $errlist $fff )
101  endif
102
103end
104
105if( ! -f $argv[1]/mask_group+tlrc.HEAD )then
106  echo "** ERROR: mask file $argv[1]/mask_group+tlrc.HEAD does not exist" ; nbad++
107endif
108
109if( $nbad > 0 )then
110  echo "** FATAL ERROR: Can't continue after such problems :-("
111  exit 1
112endif
113
114echo "-- Using $nerrWH prewhitened errts datasets and $nerrOT other errts datasets"
115echo "-- Using $nmskGM gray-matter mask datasets  and $nmskGR group mask datasets"
116
117## echo "-- list of errts datasets:"
118## echo "   $errlist"
119## echo "-- list of mask datasets:"
120## echo "   $msklist"
121
122## process all the matrices at once to get the collective passband
123
124echo "========== Computing stimulus passband =========="
125
126set bbb = ( `stimband $matlist` )
127
128if( $#bbb < 2 || "$bbb[1]" == "$bbb[2]" )then
129  echo "program stimband failed on the matrices for some reason :-("
130  exit 1
131endif
132
133echo "="
134echo "Stimulus passband = $bbb[1] to $bbb[2] Hz"
135echo "="
136
137## bandpass and 3dFWHMx each noise file separately
138
139set rrr = `3dnewid -fun`
140
141echo "========== Bandpassing and computing ACF smoothness =========="
142
143echo "# 3dFWHMx output" > BPtemp.$rrr.1D
144
145foreach nnn ( `count -dig 1 1 $#errlist` )
146
147  set eee = $errlist[$nnn]
148  set mmm = $msklist[$nnn]
149
150  3dTproject -input $eee -mask $mmm -prefix BPtemp.$rrr.nii -passband $bbb[1] $bbb[2]
151
152# the 'tail -1' is to keep only the 2nd (ACF) line from 3dFWHMx
153
154  3dFWHMx -acf BPtemp.$nnn.$rrr.1D -mask $mmm BPtemp.$rrr.nii | tail -1 >> BPtemp.$rrr.1D
155
156  \rm BPtemp.$rrr.nii
157
158end
159
160## compute the average of the ACF outputs
161
162set bbb = ( `1dsum -mean BPtemp.$rrr.1D` )
163
164cp BPtemp.$rrr.1D BPacf.1D
165
166echo "Average ACF FWHM = $bbb[4] mm"
167
168## run 3dClustSim
169
170echo "========== Running 3dClustSim =========="
171
1723dClustSim -acf $bbb[1] $bbb[2] $bbb[3]        \
173           -mask $argv[1]/mask_group+tlrc.HEAD \
174           -both -MEGA -sumup -sumup -prefix BPclust
175
176\rm BPtemp*${rrr}*
177
178echo "========== e Finito =========="
179exit 0
180