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