1#!/bin/tcsh 2 3@global_parse `basename $0` "$*" ; if ($status) exit 0 4 5# --------------------- revision history ------------------------- 6# Jan, 2017 7# + rename 8# 9# Jan 27, 2017 10# + update opts 11# 12#set version = "1.3" 13#set rev_dat = "Feb, 2017" 14# + drive viewer 15# 16#set version = "1.4" 17#set rev_dat = "Feb 20, 2017" 18# + supplementary anat input 19# 20#set version = "1.5" 21#set rev_dat = "Feb 26, 2017" 22# + more with anat input 23# + wdir 24# + pre-alignment/matching stuff 25# 26#set version = "1.6"; set rev_dat = "Feb 27, 2017" 27# + pre-alignment/matching stuff 28# 29#set version = "1.7"; set rev_dat = "March 10, 2017" 30# + fix when res struct and dwi don't have same origin to start 31# + add QC snapshot of b0 edges overlayed on orig anat 32# 33#set version = "1.8"; set rev_dat = "Apr 13, 2017" 34# + alter I/O to be consistent across progs 35# + wdir update 36# 37#set version = "1.9"; set rev_dat = "Apr 26, 2017" 38# + wdir update 39# 40#set version = "2.0"; set rev_dat = "May 1, 2017" 41# + add in ability to mask struc set 42# 43#set version = "2.1"; set rev_dat = "May 12, 2017" 44# + fixed error if no anat entered 45# 46#set version = "2.2"; set rev_dat = "May 17, 2017" 47# + now '3dDWItoDT -debug_briks ...' by default 48# + Has some fit-type parameters there. 49# 50#set version = "2.3"; set rev_dat = "Aug 1, 2017" 51# + now '-reweight' on by default 52# + now '-cumulative_wts' on by default 53# 54#set version = "2.4"; set rev_dat = "Sep 04, 2017" 55# + work with new @chauffeur, -prefix only 56# + updated format of help file, so that user can be prompted for 57# more opts 58# + output QC images for uncert (V12 and FA), and also for showing 59# anat lines on the dwi vol 60# 61#set version = "2.5b"; set rev_dat = "Sep 06, 2017" 62# + for showing anat lines on DWI, allow anat lines to be clearer 63# --> quickchange in how this is done: don't use -master. 64# 65#set version = "2.6"; set rev_dat = "Feb 08, 2018" 66# + fixed bug: wasn't resampling mask if input when the DWI was, 67# which created an erroneous mismatch between them. Sigh. Fixed 68# now. Thanks, Sandra H., for finding this! 69# 70#set version = "2.7"; set rev_dat = "Feb 15, 2018" 71# + make QC dir for images 72# + rename output QC files more succinctly 73# 74#set version = "2.8"; set rev_dat = "Mar 14, 2018" 75# + fix minor bug/crash if no "ref" set input 76# --> mask error about recentering mask. 77# 78#set version = "2.9"; set rev_dat = "July 25, 2018" 79# + match func_range_nz to updated @chauffeur 80# 81#set version = "3.0"; set rev_dat = "July 31, 2018" 82# + add in 1dDW_Grad_o_Mat++'s "-check_abs_min .." opt 83# 84#set version = "3.1"; set rev_dat = "Feb 12, 2019" 85# + [PT] change "checks" to use '3dinfo -prefix ...' as a better 86# methodology 87# 88#set version = "3.2"; set rev_dat = "Jan 29, 2020" 89# + [PT] make output images with olay=b0 have 95%ile (in volume) 90# value for top of colorbar: easier to see structure. 91# 92set version = "3.21"; set rev_dat = "Sep 27, 2021" 93# + [PT] chauffeur label_size 3 -> 4, bc imseq.c shifted all sizes 94# down one level 95# 96# ---------------------------------------------------------------- 97 98set this_prog = "fat_proc_dwi_to_dt" 99set tpname = "${this_prog:gas/fat_proc_//}" 100set here = $PWD 101 102# ----------------- find AFNI and set viewer --------------------- 103 104# find AFNI binaries directory and viewer location 105set adir = "" 106set my_viewer = "" 107which afni >& /dev/null 108if ( $status ) then 109 echo "** Cannot find 'afni' (?!)." 110 goto BAD_EXIT 111else 112 set aa = `which afni` 113 set adir = $aa:h 114endif 115 116# default location of viewer: user could modify! 117set my_viewer = "$adir/@chauffeur_afni" 118 119# ----------------------- set defaults -------------------------- 120 121set idwi = "" # necessary input: DWI dset 122set invecmat = ( "" "" "" "" ) # switches+names for bvecs and bvals 123set imask = "" # optional mask; otherwise, automask 124set istrc = "" # optional struc vol (like in TORT) 125set iref = "" # option, for getting orient + orig 126 127set iflip = "" # for 1dDW_*++: can flip (yaaawn) 128set sca = "-scale_out_1000" # for 3dDWItoDT: -scale_out_1000 129set do_rwt = "-reweight" # on by def 130set do_cwt = "-cumulative_wts" # on by def 131 132set odir = "" 133set opref = "" 134set dtipref = "dt" # NIFTI prefix for NII output 135set ori_new = "RPI" # default file output orientation 136 137set DO_VIEWER = 1 # def: do viewing 138set qc_prefix = "" # def: autoname; user can enter 139set qc_fa_thr = 0.2 # def: FA olay thr for QC 140set qc_fa_max = 0.9 # def: FA olay max for QC 141set qc_v12unc_max = 0.349 # def: V12 uncert stdev max, ~20deg 142set qc_faunc_max = 0.05 # def: FA uncert stdev max 143set output_cmd = 1 # def: output copy of this command 144 # -> should expand to have all cmds 145set cmd_file = "" # def: same name as viewer 146 147set DO_UNCERT = 1 148set Niters = 300 149set extra_unc_cmds = "" 150 151set DO_CLEAN = "1" 152set wdir = "__WORKING_$tpname" 153set WDIR_EX = "1" # put opref on wdir (unless user names) 154 155set DO_STRUC_MASK = "0" 156set sss = "" # just null val; later is CoM 157 158set ca_min_arg = ( ) # for '-check_abs_min ..' 159 160# ------------------- process options, a la rr ---------------------- 161 162if ( $#argv == 0 ) goto SHOW_HELP 163 164set ac = 1 165while ( $ac <= $#argv ) 166 # terminal options 167 if ( ("$argv[$ac]" == "-h" ) || ("$argv[$ac]" == "-help" )) then 168 goto SHOW_HELP 169 endif 170 if ( "$argv[$ac]" == "-ver" ) then 171 goto SHOW_VERSION 172 endif 173 174 # --------------- input dset(s) ---------------- 175 if ( "$argv[$ac]" == "-in_dwi" ) then 176 if ( $ac >= $#argv ) goto FAIL_MISSING_ARG 177 @ ac += 1 178 set idwi = "$argv[$ac]" 179 180 else if ( "$argv[$ac]" == "-mask" ) then 181 if ( $ac >= $#argv ) goto FAIL_MISSING_ARG 182 @ ac += 1 183 set imask = "$argv[$ac]" 184 185 else if ( "$argv[$ac]" == "-in_struc_res" ) then 186 if ( $ac >= $#argv ) goto FAIL_MISSING_ARG 187 @ ac += 1 188 set istrc = "$argv[$ac]" 189 190 else if ( "$argv[$ac]" == "-in_ref_orig" ) then 191 if ( $ac >= $#argv ) goto FAIL_MISSING_ARG 192 @ ac += 1 193 set iref = "$argv[$ac]" 194 195 # ------------- make mask from struc -------------- 196 else if ( "$argv[$ac]" == "-mask_from_struc" ) then 197 set DO_STRUC_MASK = "1" 198 199 # ------------- input vecmat and bval -------------- 200 else if ( "$argv[$ac]" == "-in_col_matA" ) then 201 if ( $ac >= $#argv ) goto FAIL_MISSING_ARG 202 set invecmat[1] = $argv[$ac] 203 @ ac += 1 204 set invecmat[2] = "$argv[$ac]" 205 206 else if ( "$argv[$ac]" == "-in_col_matT" ) then 207 if ( $ac >= $#argv ) goto FAIL_MISSING_ARG 208 set invecmat[1] = $argv[$ac] 209 @ ac += 1 210 set invecmat[2] = "$argv[$ac]" 211 212 else if ( "$argv[$ac]" == "-in_col_vec" ) then 213 if ( $ac >= $#argv ) goto FAIL_MISSING_ARG 214 set invecmat[1] = $argv[$ac] 215 @ ac += 1 216 set invecmat[2] = "$argv[$ac]" 217 218 else if ( "$argv[$ac]" == "-in_row_vec" ) then 219 if ( $ac >= $#argv ) goto FAIL_MISSING_ARG 220 set invecmat[1] = $argv[$ac] 221 @ ac += 1 222 set invecmat[2] = "$argv[$ac]" 223 224 # not necessary; default is just empty 225 else if ( "$argv[$ac]" == "-in_bvals" ) then 226 if ( $ac >= $#argv ) goto FAIL_MISSING_ARG 227 set invecmat[3] = $argv[$ac] 228 @ ac += 1 229 set invecmat[4] = "$argv[$ac]" 230 231 # ----------------- outputs --------------------- 232 # *THIS* is the main prefix that can determine 233 # output file location-- required! 234 else if ( "$argv[$ac]" == "-prefix" ) then 235 if ( $ac >= $#argv ) goto FAIL_MISSING_ARG 236 @ ac += 1 237 set opref = "$argv[$ac]" 238 239 # will just take the basename of this-- optional 240 else if ( "$argv[$ac]" == "-prefix_dti" ) then 241 if ( $ac >= $#argv ) goto FAIL_MISSING_ARG 242 @ ac += 1 243 set dtipref = "$argv[$ac]" 244 245 else if ( "$argv[$ac]" == "-workdir" ) then 246 if ( $ac >= $#argv ) goto FAIL_MISSING_ARG 247 @ ac += 1 248 set wdir = "$argv[$ac]" 249 set WDIR_EX = "0" 250 251 # ----------------- other opts --------------------- 252 else if ( ("$argv[$ac]" == "-flip_x") || \ 253 ("$argv[$ac]" == "-flip_y") || \ 254 ("$argv[$ac]" == "-flip_z") || \ 255 ("$argv[$ac]" == "-no_flip") ) then 256 set iflip = "$argv[$ac]" 257 258 # [PT: July 31, 2018] add in for 1dDW_Grad*++; pass both opt flag 259 # and argvalue along. 260 else if ( "$argv[$ac]" == "-check_abs_min" ) then 261 if ( $ac >= $#argv ) goto FAIL_MISSING_ARG 262 @ ac += 1 263 set ca_min_arg = ( "-check_abs_min" $argv[$ac] ) 264 265 else if ( "$argv[$ac]" == "-no_scale_out_1000" ) then 266 set sca = "" 267 268 else if ( "$argv[$ac]" == "-no_reweight" ) then 269 set do_rwt = "" 270 271 else if ( "$argv[$ac]" == "-no_cumulative_wts" ) then 272 set do_cwt = "" 273 274 else if ( "$argv[$ac]" == "-qc_fa_thr" ) then 275 if ( $ac >= $#argv ) goto FAIL_MISSING_ARG 276 @ ac += 1 277 set qc_fa_thr = "$argv[$ac]" 278 279 else if ( "$argv[$ac]" == "-qc_fa_max" ) then 280 if ( $ac >= $#argv ) goto FAIL_MISSING_ARG 281 @ ac += 1 282 set qc_fa_max = "$argv[$ac]" 283 284 else if ( "$argv[$ac]" == "-qc_fa_unc_max" ) then 285 if ( $ac >= $#argv ) goto FAIL_MISSING_ARG 286 @ ac += 1 287 set qc_faunc_max = "$argv[$ac]" 288 289 else if ( "$argv[$ac]" == "-qc_v12_unc_max" ) then 290 if ( $ac >= $#argv ) goto FAIL_MISSING_ARG 291 @ ac += 1 292 set qc_v12unc_max = "$argv[$ac]" 293 294 else if ( "$argv[$ac]" == "-qc_prefix" ) then 295 if ( $ac >= $#argv ) goto FAIL_MISSING_ARG 296 @ ac += 1 297 set qc_prefix = "$argv[$ac]" 298 299 else if ( "$argv[$ac]" == "-no_qc_view" ) then 300 set DO_VIEWER = 0 301 302 else if ( "$argv[$ac]" == "-no_cmd_out" ) then 303 set output_cmd = 0 304 305 else if ( "$argv[$ac]" == "-no_clean" ) then 306 set DO_CLEAN = "0" 307 308 # ------ dti par uncert -------------- 309 else if ( "$argv[$ac]" == "-uncert_off" ) then 310 set DO_UNCERT = 0 311 312 else if ( "$argv[$ac]" == "-uncert_iters" ) then 313 if ( $ac >= $#argv ) goto FAIL_MISSING_ARG 314 @ ac += 1 315 set Niters = "$argv[$ac]" 316 317 else if ( "$argv[$ac]" == "-uncert_extra_cmds" ) then 318 if ( $ac >= $#argv ) goto FAIL_MISSING_ARG 319 @ ac += 1 320 set extra_unc_cmds = "$argv[$ac]" 321 322 else 323 echo "** unexpected option #$ac = '$argv[$ac]'" 324 goto BAD_EXIT 325 326 endif 327 @ ac += 1 328end 329 330# ======================================================================= 331# ============================ ** SETUP ** ============================== 332# ======================================================================= 333 334# ============================ input files ============================== 335 336echo "++ Start script version: $version" 337 338# NEED these two inputs 339if ( "$idwi" == "" ) then 340 echo "** ERROR: no DWI file input?" 341 goto BAD_EXIT 342else if ( "$invecmat[2]" == "" ) then 343 echo "** ERROR: no gradient/matrix file input?" 344 goto BAD_EXIT 345endif 346 347# make sure we can read DWI OK 348set check = `3dinfo -prefix "$idwi"` 349if ( "$check" == "NO-DSET" ) then 350 echo "** ERROR: can't find inset file: $idwi" 351 goto BAD_EXIT 352else 353 echo "++ Found inset DWI file: $idwi" 354endif 355 356# check for vec/mat file ???????????????/ 357if ( 0 ) then 358if ( -f "$invecmat[2]" ) then 359 echo "++ Found input vec/mat file: $invecmat[2]" 360else 361 echo "** ERROR: can't find entered vec/mat file $invecmat[2]" 362 goto BAD_EXIT 363endif 364endif 365 366# ------------------------------- 367 368# check if MASK was input, then if can be read 369if ( "$imask" != "" ) then 370 set check = `3dinfo -prefix "$imask"` 371 if ( "$check" == "NO-DSET" ) then 372 echo "** ERROR: can't find input mask file: $imask" 373 goto BAD_EXIT 374 else 375 echo "++ Found inset input file: $imask" 376 endif 377endif 378 379# check if ANAT was input, then if can be read 380if ( "$istrc" != "" ) then 381 set check = `3dinfo -prefix "$istrc"` 382 if ( "$check" == "NO-DSET" ) then 383 echo "** ERROR: can't find input anatomical file: $istrc" 384 goto BAD_EXIT 385 else 386 echo "++ Found inset input file: $istrc" 387 endif 388endif 389 390# check if REF was input, then if can be read 391if ( "$iref" != "" ) then 392 set check = `3dinfo -prefix "$iref"` 393 if ( "$check" == "NO-DSET" ) then 394 echo "** ERROR: can't find input ref file: $iref" 395 goto BAD_EXIT 396 else 397 echo "++ Found input ref file: $iref" 398 set ori_new = `3dinfo -orient "$iref"` 399 endif 400endif 401 402# can't have struc mask opt used if no struc vol, and can't input 403# another mask 404if ( "$DO_STRUC_MASK" == "1" ) then 405 if ( "$istrc" == "" ) then 406 echo "** ERROR: you asked for masking a struc file, but no such" 407 echo "** struc file was input using '-in_struc_res ...'" 408 goto BAD_EXIT 409 endif 410 if ( "$imask" != "" ) then 411 echo "** ERROR: you asked for masking a struc file, but you also" 412 echo "** input a mask file using '-mask ...'." 413 echo "** Sorry, mate, can't have both!" 414 goto BAD_EXIT 415 endif 416endif 417 418# check for bval file, if entered 419if ( 0 ) then 420if ( "$invecmat[4]" != "" ) then 421 if ( -f "$invecmat[4]" ) then 422 echo "++ Found entered bval file: $invecmat[4]" 423 else 424 echo "** ERROR: can't find entered bval file: $invecmat[4] !" 425 goto BAD_EXIT 426 endif 427endif 428endif 429 430# ========================= output/working dir ========================== 431 432if ( "$opref" == "" ) then 433 echo "** ERROR: need '-prefix ...' option provided!" 434 echo " See the helpfile for more information." 435 goto BAD_EXIT 436else 437 set odir = `dirname $opref` 438 set opref = `basename $opref` 439 echo "" 440 echo "++ Based on prefix, the output directory will be:" 441 echo " $odir" 442 echo "++ Based on prefix, the output prefix will be:" 443 echo " $opref" 444 echo "" 445endif 446 447# check output directory, use input one if nothing given 448 449# default output dir, if nothing input. 450if ( ! -e "$odir" ) then 451 echo "+* Output directory didn't exist. Trying to make '$odir' now." 452 mkdir "$odir" 453endif 454 455set odirqc = "$odir/QC" 456if ( "$DO_VIEWER" == "1" ) then 457 if ( ! -e "$odirqc" ) then 458 mkdir "$odirqc" 459 endif 460endif 461 462# and put working directory as subdirectory. 463if ( "$WDIR_EX" == "1" ) then 464 set wdir = $odir/${wdir}_$opref 465else 466 set wdir = $odir/$wdir 467endif 468 469# make the working directory 470if ( ! -e $wdir ) then 471 echo "++ Making working directory: $wdir" 472 mkdir $wdir 473else 474 echo "+* WARNING: Somehow found a premade working directory (?):" 475 echo " $wdir" 476 477 # don't clean preexisting directories-- could be user mistake. 478 echo " NB: will *not* clean it afterwards." 479 set DO_CLEAN = "0" 480endif 481 482# file names for lots of outputs 483 484set ocmd = "${opref}_cmd.txt" # name for output command 485set omata = "${opref}_bmatA.txt" # name for full afni bmatrix 486set obvec = "${opref}_bvec.txt" # name for full afni grads 487set obval = "${opref}_bval.txt" # name for full afni bvals 488set odwi = "${opref}_dwi.nii.gz" # name for dwis 489set omask = "${opref}_mask.nii.gz" # name for mask 490set oanat = "${opref}_anat.nii.gz" # name for anat in output 491 492#set dti = "DTI" 493set dti_par = "${dtipref}.nii.gz" # NIFTI prefix for NII output 494set dti_unc = "${dtipref}_UNC.nii.gz" # name for mask 495 496# ======================================================================= 497# =========================== ** PROCESS ** ============================= 498# ======================================================================= 499 500echo "\n-----> STARTING $this_prog ---->" 501 502# ---------------------------- CMD --------------------------------- 503 504echo "\n\nThis command:" 505echo "$this_prog $argv\n\n" 506 507if ( "$cmd_file" == "" ) then 508 set cmd_file = "$odir/$ocmd" 509endif 510 511# copy original command: 512# dump copy of command into workdir/.. 513if ( $output_cmd == 1 ) then 514 echo "++ Echoing the command to: $cmd_file" 515 516 set rec_afni_ver = `afni -ver` 517 echo "### AFNI version:" > $cmd_file 518 echo "# $rec_afni_ver\n" >> $cmd_file 519 520 echo "### Executed from the directory location:" >> $cmd_file 521 echo "# $here\n" >> $cmd_file 522 echo "### The command was:" >> $cmd_file 523 echo "# $this_prog $argv" >> $cmd_file 524 echo "\n" >> $cmd_file 525endif 526 527# ----------------- grads/bmats --------------------- 528# Make grad file: bmatrix, bvector and bvalues 529 530# make grads ... 5311dDW_Grad_o_Mat++ \ 532 -overwrite \ 533 -echo_edu \ 534 $iflip \ 535 -out_col_vec $odir/$obvec \ 536 -out_col_bval_sep $odir/$obval \ 537 $ca_min_arg \ 538 "$invecmat[1]" "$invecmat[2]" \ 539 "$invecmat[3]" "$invecmat[4]" 540# NB: the above works somewhat cheatingly by having the empty strings 541# at the end of the function call. 542 543# ... and matching matA 5441dDW_Grad_o_Mat++ \ 545 -echo_edu \ 546 -in_col_vec $odir/$obvec \ 547 -out_col_matA $odir/$omata \ 548 -overwrite 549 550# -------------- copy DWI/reset ori^2 -------------- 551 552# [PT: Feb 8, 2018] Check that mask matches grid of DWI set! 553if ( "$imask" != "") then 554 set checkgrid = `3dinfo -same_grid $imask $idwi` 555 if ( "$checkgrid[1]" == "1" ) then 556 echo "++ OK, the entered mask's grid matches that of the DWI set." 557 else 558 echo "** ERROR: need mask needs to match the grid of the DWI set!" 559 goto BAD_EXIT 560 endif 561endif 562 563 564if ( "$iref" != "" ) then 565 566 # make sure this is in same orientation with ref 567 set ori_idwi = `3dinfo -orient "$idwi"` 568 if ( "$ori_idwi" != "$ori_new" ) then 569 570 echo "++ Resampling/fitting input DWI vols." 571 set tmp_idwi = $wdir/DWI_res.nii 572 3dresample \ 573 -echo_edu \ 574 -overwrite \ 575 -inset "$idwi" \ 576 -orient "$ori_new" \ 577 -prefix $tmp_idwi 578 579 # [PT: Feb 8, 2018] ... and apply the same to the mask-- 580 # thanks, Sandra H! 581 if ( "$imask" != "") then 582 583 echo "++ Resampling the input mask." 584 set tmp_imask = $wdir/dwi_mask_res.nii 585 3dresample \ 586 -echo_edu \ 587 -overwrite \ 588 -inset "$imask" \ 589 -orient "$ori_new" \ 590 -prefix $tmp_imask 591 592 set imask = $tmp_imask 593 594 endif 595 else 596 # ... or just use input file as is 597 set tmp_idwi = "$idwi" 598 endif 599 600 if ( "$istrc" != "" ) then 601 602 # resample anat, if nec; then use it to find shift to get to 603 # refspace 604 set tmp_istrc = $wdir/anat_res.nii 605 606 set ori_anat = `3dinfo -orient "$istrc"` 607 608 if ( "$ori_anat" != "$ori_new" ) then 609 echo "++ Resampling/fitting input anat vol." 610 3dresample \ 611 -echo_edu \ 612 -overwrite \ 613 -inset "$istrc" \ 614 -orient "$ori_new" \ 615 -prefix $tmp_istrc 616 else 617 # ... or just use input file as is 618 echo "++ Copying input anat vol." 619 3dcalc \ 620 -echo_edu \ 621 -overwrite \ 622 -a "$istrc" \ 623 -expr 'a' \ 624 -prefix $tmp_istrc 625 endif 626 627 # check IJK dimensions of anat and dwi. [March 10, 2017] 628 set dimA = `3dinfo -n4 $tmp_istrc` 629 set dimD = `3dinfo -n4 $tmp_idwi` 630 foreach i ( `seq 1 1 3` ) 631 if ( $dimA[$i] != $dimD[$i] ) then 632 echo "** ERROR! IJK dimensions of $istrc and $idwi don't match!" 633 goto BAD_EXIT 634 endif 635 end 636 637 # [March 10, 2017] 638 echo "++ Making sure the anat and DWI have the same origin!" 639 3drefit \ 640 -echo_edu \ 641 -duporigin $tmp_idwi \ 642 $tmp_istrc 643 644 # [May 2, 2017] 645 if ( "$DO_STRUC_MASK" == "1" ) then 646 set ss_istrc = $wdir/anat_ss.nii 647 set imask = $wdir/anat_mask.nii 648 649 3dSkullStrip \ 650 -echo_edu \ 651 -overwrite \ 652 -mask_vol \ 653 -prefix $ss_istrc \ 654 -input $tmp_istrc 655 656 # smooth out a bit, maybe fill small holes, and make 657 # output be smoothed 658 3dmask_tool \ 659 -overwrite \ 660 -prefix $imask \ 661 -dilate_inputs 2 -2 \ 662 -inputs $ss_istrc 663 endif 664 665 # now do alignment to ref anat 666 set isetinref = $wdir/anat_in_ref.nii 667 set map_to_ref = $wdir/map_to_ref.param.1D 668 3dAllineate \ 669 -lpa \ 670 -source_automask \ 671 -autoweight \ 672 -cmass \ 673 -final wsinc5 \ 674 -1Dparam_save $map_to_ref \ 675 -warp shift_only \ 676 -source $tmp_istrc \ 677 -base "$iref" \ 678 -prefix $isetinref \ 679 -overwrite 680 else 681 # use dwi[0] for it. 682 set isetinref = $wdir/dwi0_in_ref.nii 683 set map_to_ref = $wdir/map_to_ref.param.1D 684 3dAllineate \ 685 -lpa \ 686 -source_automask \ 687 -autoweight \ 688 -cmass \ 689 -final wsinc5 \ 690 -1Dparam_save $map_to_ref \ 691 -warp shift_only \ 692 -source $tmp_idwi"[0]" \ 693 -base "$iref" \ 694 -prefix $isetinref \ 695 -overwrite 696 endif 697 698 if ( $DO_VIEWER ) then 699 700 echo "++ Initial QC images: how did we align to ref?" 701 702 if ( "$istrc" == "" ) then 703 set ffo = "b0" 704 else 705 set ffo = "struc" 706 endif 707 708 if ( $qc_prefix == "" ) then 709 set vpref0 = ${opref}_qc_ref_$ffo 710 else 711 set vpref0 = ${qc_prefix}_qc_ref_$ffo 712 endif 713 714 echo "\n\n" 715 echo "++ QC image 00 ($isetinref on $iref): $vpref0" 716 echo "\n\n" 717 # need to put '[0]' on $iref? 718 # [PT: Jan 29, 2020] make the perc range of olay to be 95, to 719 # be more useful in output img 720 $my_viewer \ 721 -ulay "$iref" \ 722 -ulay_range "2%" "98%" \ 723 -olay "$isetinref" \ 724 -func_range_perc 95 \ 725 -pbar_posonly \ 726 -opacity 4 \ 727 -prefix "$odirqc/$vpref0" \ 728 -montx 5 -monty 3 \ 729 -set_xhairs OFF \ 730 -label_mode 1 -label_size 4 \ 731 -do_clean 732 733 endif 734 735 # --> from either of the above maps, we just care about the 736 # --> $map_to_ref, to apply those shifts to the DWI set, so it 737 # --> isn't resampled/smoothed. 738 739 # les shifts; need to reverse signs to apply 740 set shsh = `grep -v -h "#" $map_to_ref` 741 set sss = ( 0 0 0 ) 742 foreach i ( `seq 1 1 3` ) 743 set sss[$i] = `echo "scale=2; -1.0 * $shsh[$i]" | bc` 744 end 745 746 # copy dwi, and apply shifts with refit 747 3dcalc \ 748 -echo_edu \ 749 -overwrite \ 750 -a "$tmp_idwi" \ 751 -expr 'a' \ 752 -prefix $odir/$odwi 753 754 # apply to both DWI and, if applicable, the anat 755 3drefit \ 756 -dxorigin "$sss[1]" \ 757 -dyorigin "$sss[2]" \ 758 -dzorigin "$sss[3]" \ 759 $odir/$odwi 760 761 if ( "$istrc" != "" ) then 762 3dcalc \ 763 -echo_edu \ 764 -overwrite \ 765 -a "$tmp_istrc" \ 766 -expr 'a' \ 767 -prefix $odir/$oanat 768 3drefit \ 769 -dxorigin "$sss[1]" \ 770 -dyorigin "$sss[2]" \ 771 -dzorigin "$sss[3]" \ 772 $odir/$oanat 773 endif 774 775else 776 # just copy 777 3dcalc \ 778 -echo_edu \ 779 -overwrite \ 780 -a "$idwi" \ 781 -expr 'a' \ 782 -prefix $odir/$odwi 783 784 if ( "$istrc" != "" ) then 785 3dcalc \ 786 -echo_edu \ 787 -overwrite \ 788 -a "$istrc" \ 789 -expr 'a' \ 790 -prefix $odir/$oanat 791 endif 792endif 793 794# ----------------- WB mask: make or copy ------------- 795 796if ( "$imask" == "" ) then 797 # make mask. 798 # [June 6, 2017]: updated because TORTOISE can have large 799 # neg values in DWIs, so take abs value 800 echo "++ Automasking to make: $odir/$omask" 801 3dAutomask \ 802 -echo_edu \ 803 -overwrite \ 804 -prefix $odir/$omask \ 805 "3dcalc( -a ${odir}/${odwi}[0] -expr (a*step(a)) )" 806 #$odir/$odwi'[0]' 807else 808 # just copy mask and apply transform [May 2, 2017] 809 3dcalc \ 810 -echo_edu \ 811 -overwrite \ 812 -a $imask \ 813 -expr 'a' \ 814 -prefix $odir/$omask 815 816 if ( "$sss" != "" ) then 817 818 3drefit \ 819 -dxorigin "$sss[1]" \ 820 -dyorigin "$sss[2]" \ 821 -dzorigin "$sss[3]" \ 822 $odir/$omask 823 endif 824endif 825 826# [March 10, 2017] 827if ( $DO_VIEWER ) then 828 if ( "$iref" != "" ) then 829 830 echo "++ More QC images: b0 on initial ref." 831 832 set b0mskd = "$wdir/b0mskd.nii" 833 set b0edge = "$wdir/b0edge.nii" 834 835 3dcalc \ 836 -a $odir/$omask \ 837 -b $odir/$odwi'[0]' \ 838 -expr 'step(a)*b' \ 839 -prefix $b0mskd \ 840 -float \ 841 -overwrite 842 843 3dedge3 \ 844 -overwrite \ 845 -prefix $b0edge \ 846 -input $b0mskd 847 848 if ( $qc_prefix == "" ) then 849 set vpref0 = ${opref}_qc_ref_b0e 850 else 851 set vpref0 = ${qc_prefix}_qc_ref_b0e 852 endif 853 854 echo "\n\n" 855 echo "++ QC image 01 ($odir/$odwi'[0]' on $iref): $vpref0" 856 echo "\n\n" 857 # need to put '[0]' on $iref? 858 $my_viewer \ 859 -ulay "$iref" \ 860 -ulay_range "2%" "98%" \ 861 -olay "$b0edge" \ 862 -func_range_perc_nz 50 \ 863 -pbar_posonly \ 864 -cbar "red_monochrome" \ 865 -opacity 6 \ 866 -prefix "$odirqc/$vpref0" \ 867 -montx 5 -monty 3 \ 868 -set_xhairs OFF \ 869 -label_mode 1 -label_size 4 \ 870 -do_clean 871 872 # ------------- 873 874 echo "++ More QC images: initial ref on b0." 875 876 set refedge = "$wdir/refedge.nii" 877 3dedge3 \ 878 -overwrite \ 879 -prefix $refedge \ 880 -input $iref 881 882 # we want to view real b0 at its resolution, but also want the 883 # olay anat edges to be clear; so make this resampled version 884 # that looks the same (b/c of NN interp) but has higher res 885 # (so we see the olay nicerly). 886 887 set newres = `3dinfo -ad3 $iref` 888 889 set b0resam = "$wdir/b0resam.nii" 890 3dresample \ 891 -echo_edu \ 892 -prefix $b0resam \ 893 -rmode NN \ 894 -dxyz $newres \ 895 -input $odir/$odwi'[0]' 896 897 898 if ( $qc_prefix == "" ) then 899 set vpref0b = ${opref}_qc_b0_refe 900 else 901 set vpref0b = ${qc_prefix}_qc_b0_refe 902 endif 903 904 echo "\n\n" 905 echo "++ QC image 01 (edgey $iref on $odir/$odwi'[0]'): $vpref0b" 906 echo "\n\n" 907 # need to put '[0]' on $iref? 908 $my_viewer \ 909 -ulay $b0resam \ 910 -ulay_range "2%" "98%" \ 911 -olay "$refedge" \ 912 -func_range_perc_nz 50 \ 913 -pbar_posonly \ 914 -cbar "red_monochrome" \ 915 -opacity 6 \ 916 -prefix "$odirqc/$vpref0b" \ 917 -montx 5 -monty 3 \ 918 -set_xhairs OFF \ 919 -label_mode 1 -label_size 4 \ 920 -do_clean 921 922 endif 923endif 924 925 926# ------------------- tensor fit ----------------- 927 928echo "++ Fitting tensor, estimating params and fit: $odir/${dti_par}*" 9293dDWItoDT \ 930 -echo_edu \ 931 -overwrite \ 932 -debug_briks \ 933 -eigs -nonlinear -sep_dsets \ 934 $do_rwt \ 935 $do_cwt \ 936 $sca \ 937 -prefix $odir/$dti_par \ 938 -mask $odir/$omask \ 939 -bmatrix_FULL $odir/$omata \ 940 $odir/$odwi 941 942# ---------------------------- VIEWER --------------------------------- 943 944if ( $DO_VIEWER ) then 945 946 echo "++ Making QC images." 947 948 if ( $qc_prefix == "" ) then 949 set vpref0 = ${opref}_qc_ref_b0 950 set vpref1 = ${opref}_qc_MD_FA${qc_fa_thr:gas/.//} 951 else 952 set vpref0 = ${qc_prefix}_qc_ref_b0 953 set vpref1 = ${qc_prefix}_qc_MD_FA${qc_fa_thr:gas/.//} 954 endif 955 956 if ( "$iref" != "" ) then 957 958 echo "\n\n" 959 echo "++ QC image 02 ($odir/${odwi}[0] on $iref): $vpref0" 960 echo "\n\n" 961 # need to put '[0]' on $iref? 962 $my_viewer \ 963 -ulay "$iref" \ 964 -ulay_range "2%" "98%" \ 965 -olay "$odir/${odwi}""[0]" \ 966 -pbar_posonly \ 967 -opacity 4 \ 968 -prefix "$odirqc/$vpref0" \ 969 -montx 5 -monty 3 \ 970 -set_xhairs OFF \ 971 -label_mode 1 -label_size 4 \ 972 -do_clean 973 endif 974 975 echo "\n\n" 976 echo "++ QC image 03 (final output FA>${qc_fa_thr} on MD): $vpref1" 977 echo " (olay colorbar is 'Plasma', with max value ${qc_fa_max})" 978 echo "\n\n" 979 980 $my_viewer \ 981 -ulay "$odir/${dtipref}_MD.nii.gz" \ 982 -ulay_range "2%" "98%" \ 983 -olay "$odir/${dtipref}_FA.nii.gz" \ 984 -func_range $qc_fa_max \ 985 -thr_olay $qc_fa_thr \ 986 -pbar_posonly \ 987 -opacity 7 \ 988 -prefix "$odirqc/$vpref1" \ 989 -montx 5 -monty 3 \ 990 -set_xhairs OFF \ 991 -label_mode 1 -label_size 4 \ 992 -do_clean 993 994endif 995 996echo "++ Tensor par. uncertainty with $Niters iterations: $odir/${dti_unc}*" 997 998if ( $DO_UNCERT ) then 999 3dDWUncert \ 1000 -echo_edu \ 1001 -overwrite \ 1002 -inset $odir/$odwi \ 1003 -input $odir/$dtipref \ 1004 -bmatrix_FULL $odir/$omata \ 1005 -mask $odir/$omask \ 1006 -prefix $odir/$dti_unc \ 1007 -iters $Niters \ 1008 $extra_unc_cmds 1009 1010 if ( $DO_VIEWER ) then 1011 1012 if ( $qc_prefix == "" ) then 1013 set vpref4 = ${opref}_qc_b0_uncV12 1014 set vpref5 = ${opref}_qc_b0_uncFA 1015 else 1016 set vpref4 = ${qc_prefix}_qc_b0_uncV12 1017 set vpref5 = ${qc_prefix}_qc_b0_uncFA 1018 endif 1019 1020 echo "\n\n" 1021 echo "++ QC image 04 (final output uncert-stdev of V1 on b0): $vpref4" 1022 echo " (olay colorbar is 'Viridis', with max value ${qc_v12unc_max})" 1023 echo "\n\n" 1024 1025 $my_viewer \ 1026 -ulay "$odir/$odwi""[0]" \ 1027 -ulay_range "2%" "98%" \ 1028 -olay "$odir/${dtipref}_UNC.nii.gz""[1]" \ 1029 -func_range $qc_v12unc_max \ 1030 -pbar_posonly \ 1031 -cbar "Viridis" \ 1032 -opacity 6 \ 1033 -prefix "$odirqc/$vpref4" \ 1034 -montx 5 -monty 3 \ 1035 -set_xhairs OFF \ 1036 -label_mode 1 -label_size 4 \ 1037 -do_clean 1038 1039 echo "\n\n" 1040 echo "++ QC image 05 (final output uncert-stdev of FA on b0): $vpref5" 1041 echo " (olay colorbar is 'Viridis', with max value ${qc_faunc_max})" 1042 echo "\n\n" 1043 1044 $my_viewer \ 1045 -ulay "$odir/$odwi""[0]" \ 1046 -ulay_range "2%" "98%" \ 1047 -olay "$odir/${dtipref}_UNC.nii.gz""[5]" \ 1048 -func_range $qc_faunc_max \ 1049 -pbar_posonly \ 1050 -cbar "Viridis" \ 1051 -opacity 6 \ 1052 -prefix "$odirqc/$vpref5" \ 1053 -montx 5 -monty 3 \ 1054 -set_xhairs OFF \ 1055 -label_mode 1 -label_size 4 \ 1056 -do_clean 1057 endif 1058endif 1059 1060# ------------------------- clean (y/n) ----------------------------- 1061 1062# clean, by default 1063if ( "$DO_CLEAN" == "1" ) then 1064 echo "\n++ Cleaning working directory!\n" 1065 \rm -rf $wdir 1066else 1067 echo "\n++ NOT removing working directory '$wdir'.\n" 1068endif 1069 1070# final messages 1071echo "" 1072echo "++ The final DWI data are here: $odir/${opref}*" 1073echo "++ The final DTI data are here: $odir/${dtipref}*" 1074if ( "$DO_VIEWER" == "1") then 1075 echo "++ The final QC images are here: ${odirqc}/${opref}*" 1076endif 1077echo "" 1078 1079# --------------------------------------------------------------------- 1080 1081goto GOOD_EXIT 1082 1083# ======================================================================== 1084# ======================================================================== 1085 1086SHOW_HELP: 1087cat << EOF 1088------------------------------------------------------------------------- 1089 1090 This program is for doing tensor and DT parameter fitting, as well as 1091 the uncertainty of DT parameters that are needed for tractography. 1092 1093 Ver. $version (PA Taylor, ${rev_dat}) 1094 1095------------------------------------------------------------------------- 1096 1097 RUNNING: 1098 1099 This script has two *required* arguments ('-in_dwi ...' and some 1100 kind of gradient/matrix file input. 1101 1102 The rest are optional, but it is highly recommended to input a 1103 reference data set ('-in_ref ...') if you have used a processing 1104 tool that resets origin+orientation (such as TORTOISE), as well as 1105 using '-scale_out_1000' to make the output units of the physical DT 1106 measures nicer. 1107 1108 $this_prog \ 1109 -in_dwi DWI \ 1110 {-in_col_matA | -in_col_matT | \ 1111 -in_col_vec | -in_row_vec} GRADMAT \ 1112 -prefix PPP \ 1113 {-in_bvals BVAL} \ 1114 {-mask MASK} \ 1115 {-mask_from_struc} \ 1116 {-in_struc_res STRUC} \ 1117 {-in_ref_orig REF} \ 1118 {-prefix_dti PREFIX_D} \ 1119 {-flip_x | -flip_y | -flip_z | -no_flip} \ 1120 {-no_scale_out_1000} \ 1121 {-no_reweight} \ 1122 {-no_cumulative_wts} \ 1123 {-qc_prefix QCPREF} \ 1124 {-qc_fa_thr TTT} \ 1125 {-qc_fa_max MMM} \ 1126 {-qc_fa_unc_max UM} \ 1127 {-qc_v12_unc_max V} \ 1128 {-no_qc_view} \ 1129 {-no_cmd_out} \ 1130 {-workdir WWW} \ 1131 {-no_clean} \ 1132 {-uncert_off} \ 1133 {-uncert_iters NN} \ 1134 {-uncert_extra_cmds STR} 1135 1136 where: 1137 1138 -in_dwi DWI :4D volume of N DWIs. Required. 1139 1140 -in_col_matA | 1141 -in_col_matT | 1142 -in_col_vec | 1143 -in_row_vec GRADMAT 1144 :input text file of N gradient vectors or 1145 bmatrices. By default, it is assumed that 1146 these still have physical units in them (or that 1147 there is an accompanying BVAL file input), so 1148 scaling physical values by 1000 is on by default; 1149 see turning this scaling off, if unnecessary, by 1150 using '-no_scale_out_1000', below. 1151 1152 -prefix PPP :set prefix for output DWI data; required. 1153 1154 -in_bvals BVAL :optional, if bvalue information is 1155 in a separate file from the b-vectors 1156 or matrices; should have same number N as 1157 volumes and vectors/matrices. 1158 -flip_x | 1159 -flip_y | 1160 -flip_z | 1161 -no_flip :can flip the DW grads, if needed; for example, 1162 based on the recommendation of @GradFlipTest. 1163 1164 -check_abs_min VVV :briefly, this can help the program push through 1165 finding tiny negative values (that miiiight be 1166 due to rounding errors or small numerical 1167 things) in columns that should only contain 1168 numbers >=0. 'VVV' is basically a tolerance for 1169 the magnitude of negative values you are willing 1170 to allow: anything between [-VVV, 0) gets zeroed 1171 for further calcs. See 1dDW_Grad_o_Mat++'s help 1172 for more information on this option (of the same 1173 name). 1174 1175 -mask MASK :optional whole brain mask can be input; 1176 otherwise, automasking is performed for the 1177 region to be tensor and parameter fit. 1178 -mask_from_struc :flag to make a mask using 3dSkullStrip+3dmask_tool 1179 from the STRUC file. 1180 NB ---> If no "-mask*" option is given, then 3dAutomask is run on 1181 the DWI set. This often ain't great, so if TORTOISE isn't 1182 producing a mask, 1) email Okan and ask him about that, and 1183 2) try '-mask_from_struc'. 1184 ALSO, if you want the whole volume to be estimated 1185 tensorially for some reason, then make a volume fully 1186 filled with 1s and pass that in as the MASK, et voila 1187 (but then calcs will likely be slooow). 1188 1189 -in_ref_orig REF :use another data set to adjust the DWI (and 1190 subsequent parameter) dsets' orientation and 1191 origin; for example, TORTOISE has default 1192 orientation and origin for all output DWIs-- it 1193 would be very advisable to use the anatomical 1194 volume that you had input into TORTOISE as REF, 1195 so that the DWIs should be viewable overlaying 1196 it afterwards; if an ANAT (below) that has been 1197 merely resampled is *not* used, then you really, 1198 really want REF to have the same contrast as the 1199 b=0 DWI volume. *Highly recommended to include!* 1200 -in_struc_res STRUC :accomplish the alignment of the output DWI to the 1201 REF data set via ANAT: a version of the anatomical 1202 that has been resampled to match the DWI set (in 1203 both orientation and origin); for example, in 1204 TORTOISE there is a 'structural.nii' file that should 1205 match this description. Both ANAT and DWI should 1206 then be well aligned to the original REF (and to 1207 each other). *Highly recommended to include!* 1208 1209 -prefix_dti PREFIX2 :set prefix for output DTI data; optional, 1210 default is '${dtipref}'. Do *not* include path 1211 information here-- that is only supplied using 1212 '-prefix ..'. 1213 1214 -no_scale_out_1000 :by default, for tensor fitting it is assumed 1215 that 1) the DW b-value information is included 1216 in the gradient vectors or grads, and 2) you are 1217 happy to have tiny numbers of physical 1218 diffusion, which in standard units are like 1219 MD~0.001 "mm^2/s", scaled by 1000 so that they 1220 are returned as MD~1 "10^{-3} mm^2/s". Isn't 1221 that nicer? I thought you'd agree-- therefore, 1222 such a kind of scaling is *on* by default. To 1223 turn that *off*, use this option flag. 1224 See the 3dDWItoDT help file for what this 1225 entails. Basically, you will likely have nicer 1226 numeric values (from scaling physical length 1227 units by 1000); otherwise, you might have small 1228 numerical values leading to issues with 1229 statistical modeling. 1230 1231 -no_reweight :by default, we *do* reweight+refit tensors during 1232 estimation; should improve fit. But what do I 1233 know? This option turns that functionality *off*. 1234 -no_cumulative_wts :by default, show overall weight factors for each 1235 gradient; may be useful as a quality control, but 1236 this option will turn that functionality *off*. 1237 1238 -qc_fa_thr TTT :set threshold for overlay FA volume in QC image 1239 (default: TTT=0.2, as for healthy adult human 1240 parenchyma). 1241 -qc_fa_max MMM :set cbar max for overlay FA volume in QC image 1242 (default: MMM=0.9, a very large value even for 1243 healthy adult human parenchyma). 1244 -qc_fa_unc_max UM :set cbar max for overlay uncert (stdev) of FA 1245 in QC image (default: UM=0.05). 1246 -qc_v12_unc_max V :set cbar max for overlay uncert (stdev) of V1 1247 towards the V2 direction for DTs, in QC image 1248 (default: UM=0.349 rads, which corresponds to 1249 20 deg). 1250 1251 -qc_prefix QCPREF :can set the prefix of the QC image files separately 1252 (default is '$opref'). 1253 1254 -no_qc_view :can turn off generating QC image files (why?) 1255 1256 -no_cmd_out :don't save the command line call of this program 1257 and the location where it was run (otherwise, it is 1258 saved by default in the ODIR/). 1259 1260 -no_clean :is an optional switch to NOT remove working 1261 directory: 1262 '$wdir' 1263 (default: remove working dir). 1264 -workdir WWW :specify a working directory, which can be removed; 1265 (default name = '$wdir'). 1266 1267 -uncert_off :don't do uncertainty calc (default is to do so); 1268 perhaps if it is slow or you want *very* different 1269 options. 1270 -uncert_iters NN :set the number of Monte Carlo iterations for the 1271 uncertainty calc (default NN=300). 1272-uncert_extra_cmds STR:put in extra commands for the uncertainty calcs 1273 (see the 3dDWUncert helpfile for more opts). 1274 1275# ----------------------------------------------------------------------- 1276 1277 EXAMPLE 1278 1279 $this_prog \ 1280 -in_dwi DWI.nii \ 1281 -in_col_matA BMTXT_AFNI.txt \ 1282 -in_struc_res ../structural.nii \ 1283 -in_ref_orig t2w.nii \ 1284 -mask mask_DWI.nii.gz \ 1285 -prefix OUTPUT/dwi 1286 1287 or 1288 1289 $this_prog \ 1290 -in_dwi ap_proc_DRBUDDI_final.nii \ 1291 -in_col_matT ap_proc_DRBUDDI_final.bmtxt \ 1292 -in_struc_res structural.nii \ 1293 -in_ref_orig t2w.nii \ 1294 -mask_from_struc \ 1295 -prefix dwi_03/dwi 1296 1297 1298------------------------------------------------------------------------- 1299 1300EOF 1301 1302 goto GOOD_EXIT 1303 1304SHOW_VERSION: 1305 echo "version $version (${rev_dat})" 1306 goto GOOD_EXIT 1307 1308FAIL_MISSING_ARG: 1309 echo "** ERROR! Missing an argument after option flag: '$argv[$ac]'" 1310 goto BAD_EXIT 1311 1312BAD_EXIT: 1313 exit 1 1314 1315# send everyone here, in case there is any cleanup to do 1316GOOD_EXIT: 1317 exit 0 1318