1#!/bin/tcsh -f
2#last edited June 03 2008
3
4set cleanup = 0
5set resample = 1
6set fast_grouping = 1
7
8#count the total number of command line parameters
9set al=$#argv
10#echo $al
11
12#if you want to overwrite pre-existing files, you can do
13setenv AFNI_DECONFLICT OVERWRITE
14   #or just add a -overwrite to most AFNI's commandline programs
15
16if (${al} < 1 ) then
17   echo ""
18   echo "**************************************"
19   echo "give me mprage volume for segmentation"
20   echo "**************************************"
21   echo ""
22   goto END
23endif
24
25if (${al} != 1 ) then
26   goto BAD_INPUT
27endif
28
29PREP:
30   set starttime=`date`
31
32   set ipref=`@GetAfniPrefix $1`
33   echo ${ipref}
34
35   echo "##########################################"
36   echo "Gathering input volume information ..."
37
38   set olddelta = `3dAttribute DELTA ${ipref}+orig \
39                    | sed 's/-//g'`
40
41   set olddim = `3dAttribute DATASET_DIMENSIONS ${ipref}+orig \
42                   | sed 's/-//g'`
43
44   ### 30 becouse 30mm is the biggest square!!!
45
46   set newx = `ccalc -form int -eval "${olddim[1]}+30/${olddelta[1]}"`
47   set newy = `ccalc -form int -eval "${olddim[2]}+30/${olddelta[2]}"`
48   set newz = `ccalc -form int -eval "${olddim[3]}+30/${olddelta[3]}"`
49
50   echo ${newx}
51   echo ${newy}
52   echo ${newz}
53
54#Good place for a jumping point
55#goto GROUPING
56
57ZEROPADDING:
58   if (0) then
59      echo "Zeropadding ..."
60      #I am not sure what this is for here,
61      #but is fails at times (anat0149_ss+orig) and it makes dset huge
62      set iprefzp = ${ipref}_ZP
63      3dZeropad   -RL $newx -AP $newy -IS $newz \
64                  -mm -prefix ${iprefzp}+orig \
65                  ${ipref}+orig
66   else
67      echo "NO Zeropadding ..."
68      set iprefzp = ${ipref}
69   endif
70LOCAL_STATS:
71   echo "RESAMPLING TO 1mm ..."
72   3dresample  -rmode Linear -dxyz 1.0 1.0 1.0 \
73               -inset ${iprefzp}+orig \
74               -prefix ${ipref}_1mm
75
76   echo "CALCULATING LOCALSTATS (1-5mm) ..."
77   set in_vol = "${ipref}_1mm+orig"
78    set voxl_size = 1mm
79
80   foreach i (1 2 3 4 5)
81       3dLocalstat   -nbhd "RECT($i,$i,$i)" \
82                     -stat median -stat MAD -stat mean \
83                     -datum short \
84                     -prefix "${ipref}_stat0${i}_${voxl_size}" \
85                     "${in_vol}"
86
87       ####chewing gum statistics
88       3dLocalstat   -nbhd "RECT($i,0,0)" \
89                     -stat median -stat MAD \
90                     -datum short \
91                     -prefix "${ipref}_chewX0${i}_${voxl_size}" \
92                     ${in_vol}
93       3dLocalstat   -nbhd "RECT(0,$i,0)" \
94                     -stat  median -stat MAD \
95                     -datum short \
96                     -prefix "${ipref}_chewY0${i}_${voxl_size}" \
97                     ${in_vol}
98       3dLocalstat   -nbhd "RECT(0,0,$i)" \
99                     -stat median -stat MAD \
100                     -datum short \
101                     -prefix "${ipref}_chewZ0${i}_${voxl_size}" \
102                     ${in_vol}
103
104      #### skewness
105      3dcalc   -a "${ipref}_stat0${i}_${voxl_size}+orig"'[$]' \
106               -b "${ipref}_stat0${i}_${voxl_size}+orig"'[0]' \
107               -datum short \
108               -expr 'a-b' \
109               -prefix "${ipref}_skew0${i}_${voxl_size}"
110
111   end
112
113    if ($resample == 0) then
114	goto LOCALSTATS6
115    endif
116
117   echo "RESAMPLING to 2mm and CALCULATING LOCALSTATS (6 and 8mm) ..."
118
119    set voxl_size = 2mm
120   3dresample  -rmode Linear \
121               -dxyz 2.0 2.0 2.0 \
122               -inset ${iprefzp}+orig \
123               -prefix ${ipref}_${voxl_size}
124
125   set in_vol="${ipref}_${voxl_size}+orig"
126
127
128LOCALSTATS6:
129   foreach i (6 8)
130      3dLocalstat    -nbhd "RECT($i,$i,$i)" \
131                     -stat median -stat MAD -stat mean \
132                     -datum short \
133                     -prefix "${ipref}_stat0${i}_${voxl_size}" \
134                     "${in_vol}"
135
136      ####chewing gum statistics
137      3dLocalstat    -nbhd "RECT($i,0,0)" \
138                     -stat median -stat MAD \
139                     -datum short \
140                     -prefix "${ipref}_chewX0${i}_${voxl_size}" \
141                     ${in_vol}
142      3dLocalstat    -nbhd "RECT(0,$i,0)" \
143                     -stat  median -stat MAD \
144                     -datum short \
145                     -prefix "${ipref}_chewY0${i}_${voxl_size}" \
146                     ${in_vol}
147      3dLocalstat    -nbhd "RECT(0,0,$i)" \
148                     -stat median -stat MAD \
149                     -datum short \
150                     -prefix "${ipref}_chewZ0${i}_${voxl_size}" \
151                     ${in_vol}
152
153      #### skewness
154      3dcalc   -a "${ipref}_stat0${i}_${voxl_size}+orig"'[$]' \
155               -b "${ipref}_stat0${i}_${voxl_size}+orig"'[0]' \
156               -datum short \
157               -expr 'a-b' \
158               -prefix "${ipref}_skew0${i}_${voxl_size}"
159
160   end
161
162    if ($resample == 0) then
163	goto LOCALSTATS10
164    endif
165
166   echo "RESAMPLING to 3mm and CALCULATING LOCALSTATS (10,13 and 16mm) ..."
167
168    set voxl_size = 3mm
169   3dresample  -rmode Linear \
170               -dxyz 3.0 3.0 3.0 \
171               -inset ${iprefzp}+orig \
172               -prefix ${ipref}_${voxl_size}
173
174   set in_vol="${ipref}_${voxl_size}+orig"
175
176LOCALSTATS10:
177   foreach i (10 13 16)
178      3dLocalstat    -nbhd "RECT($i,$i,$i)" \
179                     -stat median -stat MAD -stat mean \
180                     -datum short \
181                     -prefix "${ipref}_stat${i}_${voxl_size}" \
182                     "${in_vol}"
183
184      ####chewing gum statistics
185      3dLocalstat    -nbhd "RECT($i,0,0)" \
186                     -stat median -stat MAD \
187                     -datum short \
188                     -prefix "${ipref}_chewX${i}_${voxl_size}" \
189                     ${in_vol}
190      3dLocalstat    -nbhd "RECT(0,$i,0)" \
191                     -stat  median -stat MAD \
192                     -datum short \
193                     -prefix "${ipref}_chewY${i}_${voxl_size}" \
194                     ${in_vol}
195      3dLocalstat    -nbhd "RECT(0,0,$i)" \
196                     -stat median -stat MAD \
197                     -datum short \
198                     -prefix "${ipref}_chewZ${i}_${voxl_size}" \
199                     ${in_vol}
200
201      #### skewness
202      3dcalc   -a "${ipref}_stat${i}_${voxl_size}+orig"'[$]' \
203               -b "${ipref}_stat${i}_${voxl_size}+orig"'[0]' \
204               -datum short \
205               -expr 'a-b' \
206               -prefix "${ipref}_skew${i}_${voxl_size}"
207
208   end
209
210    if ($resample == 0) then
211	goto LOCALSTATS20
212    endif
213
214   echo "RESAMPLING to 4mm and CALCULATING LOCALSTATS (20,25 and 30mm) ..."
215
216    set voxl_size = 4mm
217   3dresample  -rmode Linear \
218               -dxyz 4.0 4.0 4.0 \
219               -inset ${iprefzp}+orig \
220               -prefix ${ipref}_4mm
221   set in_vol="${ipref}_${voxl_size}+orig"
222
223LOCALSTATS20:
224   foreach i (20 25 30)
225      3dLocalstat    -nbhd "RECT($i,$i,$i)" \
226                     -stat median -stat MAD -stat mean \
227                     -datum short \
228                     -prefix "${ipref}_stat${i}_${voxl_size}" \
229                     "${in_vol}"
230
231      ####chewing gum statistics
232      3dLocalstat    -nbhd "RECT($i,0,0)" \
233                     -stat median -stat MAD \
234                     -datum short \
235                     -prefix "${ipref}_chewX${i}_${voxl_size}" \
236                     ${in_vol}
237      3dLocalstat    -nbhd "RECT(0,$i,0)" \
238                     -stat  median -stat MAD \
239                     -datum short \
240                     -prefix "${ipref}_chewY${i}_${voxl_size}" \
241                     ${in_vol}
242      3dLocalstat    -nbhd "RECT(0,0,$i)" \
243                     -stat median -stat MAD \
244                     -datum short \
245                     -prefix "${ipref}_chewZ${i}_${voxl_size}" \
246                     ${in_vol}
247      #### skewness
248      3dcalc   -a "${ipref}_stat${i}_${voxl_size}+orig"'[$]' \
249               -b "${ipref}_stat${i}_${voxl_size}+orig"'[0]' \
250               -datum short \
251               -expr 'a-b' \
252               -prefix "${ipref}_skew${i}_${voxl_size}"
253
254   end
255
256    if ($resample == 0) then
257	goto GROUPING
258    endif
259
260RESAMPLING:
261   #########################################################################
262   ###!!! I'm grouping medians, MADs... together...  first resample all to 1mm###
263   echo "RESAMPLING ALL VOLUMES TO 1mm ..."
264
265   foreach j (06 08 10 13 16 20 25 30)
266	foreach k (2 3 4)
267	    if ( -f ${ipref}_stat${j}_${k}mm+orig.HEAD) then
268		3dresample  -rmode Linear \
269			    -dxyz 1.0 1.0 1.0 \
270			    -master ${ipref}_1mm+orig \
271			    -inset ${ipref}_stat${j}_${k}mm+orig \
272			    -prefix ${ipref}_stat${j}_1mm
273	    endif
274
275	    if ( -f ${ipref}_skew${j}_${k}mm+orig.HEAD) then
276		3dresample  -rmode Linear \
277			    -dxyz 1.0 1.0 1.0 \
278			    -master ${ipref}_1mm+orig \
279			    -inset ${ipref}_skew${j}_${k}mm+orig \
280			    -prefix ${ipref}_skew${j}_1mm
281	    endif
282
283	    if ( -f ${ipref}_chewX${j}_${k}mm+orig.HEAD) then
284		3dresample  -rmode Linear \
285			    -dxyz 1.0 1.0 1.0 \
286			    -master ${ipref}_1mm+orig \
287			    -inset ${ipref}_chewX${j}_${k}mm+orig \
288			    -prefix ${ipref}_chewX${j}_1mm
289	    endif
290
291	    if ( -f ${ipref}_chewY${j}_${k}mm+orig.HEAD) then
292		3dresample  -rmode Linear \
293			    -dxyz 1.0 1.0 1.0 \
294			    -master ${ipref}_1mm+orig \
295			    -inset ${ipref}_chewY${j}_${k}mm+orig \
296			    -prefix ${ipref}_chewY${j}_1mm
297	    endif
298
299	   if ( -f ${ipref}_chewZ${j}_${k}mm+orig.HEAD) then
300		3dresample  -rmode Linear \
301			    -dxyz 1.0 1.0 1.0 \
302			    -master ${ipref}_1mm+orig \
303			    -inset ${ipref}_chewZ${j}_${k}mm+orig \
304			    -prefix ${ipref}_chewZ${j}_1mm
305	    endif
306	end
307   end
308
309
310GROUPING:
311   echo "GROUPING TOGETHER STATS OF THE SAME TYPE"
312
313   3dcalc   -a "${ipref}_stat01_1mm+orig[0]" \
314            -expr 'a' -datum short \
315            -prefix "${ipref}_ALL_med+orig"
316   3dcalc   -a "${ipref}_stat01_1mm+orig[1]" \
317            -expr 'a' -datum short \
318            -prefix "${ipref}_ALL_MAD+orig"
319   3dcalc   -a "${ipref}_skew01_1mm+orig" \
320            -expr 'a' -datum short \
321            -prefix "${ipref}_ALL_skew+orig"
322
323   if ($fast_grouping == 1) then
324      #potentially faster but untested
325      set li1 = ()
326      set li2 = ()
327      set li3 = ()
328      foreach j ( 02 03 04 05 06 08 10 13 16 20 25 30 )
329         set li1 = ("$li1" "${ipref}_stat${j}_1mm+orig[0]")
330         set li2 = ("$li2" "${ipref}_stat${j}_1mm+orig[1]")
331         set li3 = ("$li3" "${ipref}_skew${j}_1mm+orig")
332      end
333      3dTcat -glueto "${ipref}_ALL_med+orig"    `echo "$li1"`
334      3dTcat -glueto "${ipref}_ALL_MAD+orig"    `echo "$li2"`
335      3dTcat -glueto "${ipref}_ALL_skew+orig"   `echo "$li3"`
336   else
337      #This is most inefficient but tested
338      foreach j ( 02 03 04 05 06 08 10 13 16 20 25 30 )
339		   set fname=${ipref}_stat${j}+orig
340		   echo $fname
341		   3dTcat -glueto "${ipref}_ALL_med+orig"   ${ipref}_stat${j}_1mm+orig'[0]'
342		   3dTcat -glueto "${ipref}_ALL_MAD+orig"   ${ipref}_stat${j}_1mm+orig'[1]'
343		   3dTcat -glueto "${ipref}_ALL_skew+orig"  ${ipref}_skew${j}_1mm+orig
344	   end
345   endif
346
347
348
349ANISO:
350   #############################################################################
351   ######### group chewing gum stat like median XYZ, MAD XYZ,... ###############
352   ## did median max -min along "timeline" (between X Y Z) and named it aniso...
353
354
355   echo "DOING ANISOMETRIC STATISTICS ..."
356
357   foreach k (01 02 03 04 05 06 08 10 13 16 20 25 30)
358      3dcalc   -a "${ipref}_chewX${k}_1mm+orig[0]" \
359               -b "${ipref}_chewY${k}_1mm+orig[0]"\
360               -expr 'max(a,b)' -datum short\
361               -prefix "${ipref}_chew_mmxy"
362      3dcalc   -a "${ipref}_chew_mmxy+orig"\
363               -b "${ipref}_chewZ${k}_1mm+orig[0]"\
364               -expr 'max(a,b)' -datum short\
365               -prefix "${ipref}_chew_mmax${k}"
366      rm ${ipref}_chew_mmxy+orig*
367      3dcalc   -a "${ipref}_chewX${k}_1mm+orig[0]"\
368               -b "${ipref}_chewY${k}_1mm+orig[0]"\
369               -expr 'min(a,b)' -datum short\
370               -prefix "${ipref}_chew_mmxy"
371      3dcalc   -a "${ipref}_chew_mmxy+orig"\
372               -b "${ipref}_chewZ${k}_1mm+orig[0]"\
373               -expr 'min(a,b)' -datum short\
374               -prefix "${ipref}_chew_mmin${k}"
375      rm ${ipref}_chew_mmxy+orig*
376      3dcalc   -a "${ipref}_chew_mmax${k}+orig"\
377               -b "${ipref}_chew_mmin${k}+orig"\
378               -expr 'a-b' -datum short\
379               -prefix "${ipref}_chew_medaniso${k}"
380
381      ###### ZIAD told me to include MAD
382      3dcalc   -a "${ipref}_chewX${k}_1mm+orig[1]"\
383               -b "${ipref}_chewY${k}_1mm+orig[1]"\
384               -expr 'max(a,b)' -datum short\
385               -prefix "${ipref}_chew_mmxy"
386         3dcalc   -a "${ipref}_chew_mmxy+orig"\
387               -b "${ipref}_chewZ${k}_1mm+orig[1]"\
388               -expr 'max(a,b)' -datum short\
389               -prefix "${ipref}_chew_MADmax${k}"
390         rm ${ipref}_chew_mmxy+orig*
391         3dcalc   -a "${ipref}_chewX${k}_1mm+orig[1]"\
392               -b "${ipref}_chewY${k}_1mm+orig[1]"\
393               -expr 'min(a,b)' -datum short\
394               -prefix "${ipref}_chew_mmxy"
395         3dcalc   -a "${ipref}_chew_mmxy+orig"\
396               -b "${ipref}_chewZ${k}_1mm+orig[1]"\
397               -expr 'min(a,b)' -datum short\
398               -prefix "${ipref}_chew_MADmin${k}"
399         rm ${ipref}_chew_mmxy+orig*
400         3dcalc   -a "${ipref}_chew_MADmax${k}+orig"\
401               -b "${ipref}_chew_MADmin${k}+orig"\
402               -expr 'a-b' -datum short\
403               -prefix "${ipref}_chew_MADaniso${k}"
404
405   end
406
407GROUP_ANISO:
408   echo "GROUPING ANISOMETRIC STATISTICS ..."
409   #DO NOT USE something like "${ipref}_chew_medaniso"??+orig.BRIK
410   #because a .BRIK.gz would not get read! Better stick to .HEAD
411   3dTcat -prefix "${ipref}_chew_medaniso" "${ipref}_chew_medaniso"??+orig.HEAD
412
413   3dTcat -prefix "${ipref}_chew_MADaniso" "${ipref}_chew_MADaniso"??+orig.HEAD
414
415
416GROUP_AVM:
417   echo "Creating AVM"
418   3dcalc   -a ${ipref}_1mm+orig \
419            -b ${ipref}_ALL_med+orig   \
420            -expr 'a-b' \
421            -prefix "${ipref}_AvMed"
422
423if ($cleanup == 0) then
424   goto TIME_REPORT
425endif
426
427CLEANUP:
428   echo "Cleaning up ..."
429
430   foreach j (01 02 03 04 05 06 08 10 13 16 20 25 30)
431      rm ${ipref}_stat${j}+orig.*
432      rm ${ipref}_stat${j}_1mm+orig.*
433      rm ${ipref}_skew${j}+orig.*
434      rm ${ipref}_skew${j}_1mm+orig.*
435   end
436
437   rm ${ipref}_2mm*
438   rm ${ipref}_3mm*
439   rm ${ipref}_4mm*
440
441   #se rm ...
442   foreach k ( 01 02 03 04 05 06 08 10 13 16 20 25 30 )
443     rm ${ipref}_chewX${k}*
444     rm ${ipref}_chewY${k}*
445     rm ${ipref}_chewZ${k}*
446     rm ${ipref}_chew_medaniso${k}*
447     rm ${ipref}_chew_MADaniso${k}*
448     rm ${ipref}_chew_MADmax${k}+orig*
449     rm ${ipref}_chew_mmax${k}+orig*
450     rm ${ipref}_chew_MADmin${k}+orig*
451     rm ${ipref}_chew_mmin${k}+orig*
452     #rm ${ipref}_chew_ISOx${k}*
453     #rm ${ipref}_chew_ISOy${k}*
454     #rm ${ipref}_chew_ISOz${k}*
455   end
456   goto TIME_REPORT
457
458TIME_REPORT:
459set endtime=`date`
460echo ${starttime}
461echo ${endtime}
462
463goto END
464
465BAD_INPUT:
466   echo "bad input"
467   goto END
468
469END:
470