#!/bin/tcsh -f @global_parse `basename $0` "$*" ; if ($status) exit 0 ######################################################################## # # Auto grad-flip tester # by PA Taylor (NIH, UCT, AIMS) # # Sep 2015, v2.1: # + linear DTI fits, for speed # # Sep 2015, v2.2: # + echo edu to print commands # # May 2016, v2.3: # + use '-out_grad_cols_bwt' to allow for non-constant DW factors. # # Jan 2017, v2.4: # + command line optioning # + a real help output! # #set version = "2.5"; set rev_dat = "May 5, 2017" # + changed the way outdir works: can create a dir # + also, make sure output text file goes into output dir # #set version = "2.6"; set rev_dat = "June 6, 2017" # + better results summary, no auto-prending of './' on output path # names (-> can be bad/inappropriate) # #set version = "2.7"; set rev_dat = "June 6, 2017" # + in 3dAutomask, take abs val of DWI[0], bc TORTOISE can output # huge neg values that apparently can make a big difference. # #set version = "2.8"; set rev_dat = "July 31, 2017" # + echo the useful recommendation stuff into a textfile whose name # is based on the grad-recording file. # + add opt for "-wdir *" name # + also echo "@GradFlipTest ..." command used # #set version = "2.9"; set rev_dat = "Aug 24, 2017" # + output text file with single grad will overwrite previous, not # append to # #set version = "2.91"; set rev_dat = "Sep 20, 2017" # + 'more' -> 'cat', to not wait for pause in terminal! # ---> Thanks, Tiffany Nash for finding this. # #set version = "2.92"; set rev_dat = "Oct 04, 2017" # + one case of needing to: $out_ff -> $odir/$out_ff # + update using '-prefix ...' better, no more separate 'outdir' # #set version = "2.93"; set rev_dat = "Oct 12, 2017" # + update rules of: $out_ff -> $odir/$out_ff # + better line alignment of text dumping # set version = "2.94"; set rev_dat = "Feb 1 2018" # + internal change in format of input/vars to allow subbrick # selection in bval/bvecs # ######################################################################## set this_prog = "@GradFlipTest" set here = $PWD # ----------------- find AFNI and set viewer --------------------- # find AFNI binaries directory and viewer location set adir = "" set my_viewer = "" which afni >& /dev/null if ( $status ) then echo "** Cannot find 'afni' (??!?!)." goto BAD_EXIT else set aa = `which afni` set adir = $aa:h endif # default location of viewer: user could modify! set my_viewer = "$adir/@chauffeur_afni" # ----------------------- set defaults -------------------------- set in_DSET = "" # DWI set #set in_CMD = "" # inputting grads/matr #set in_BVAL = "" # inputting bval, a la 1dDW*++ set in_MASK = "" # inputting bval, a la 1dDW*++ set odir = "." set invecmat = ( "" "" "" "" ) # switches+names for bvecs and bvals set out_ff = "" set min_tr_len = 30 # default min tract length set min_fa = 0.2 # default min FA set SOT = "" # unnec ability to scale DT by 1000 set DO_CLEAN = 0 # don't delete tmpdir, by def set FLIPS = ( "-no_flip" '-flip_x' "-flip_y" "-flip_z" ) set FLIP_sp = ( "" " " " " " " ) set FNOMS = ( "F0" "FX" "FY" "FZ" ) set tmpdir = "_tmp_TESTFLIP" set ALL_ARR = ( ) set SUMA_CALL = "" # ------------------- process options, a la rr ---------------------- if ( $#argv == 0 ) goto SHOW_HELP set ac = 1 while ( $ac <= $#argv ) # terminal options if ( ("$argv[$ac]" == "-h" ) || ("$argv[$ac]" == "-help" )) then goto SHOW_HELP endif if ( "$argv[$ac]" == "-ver" ) then goto SHOW_VERSION endif # ----- input file of DWIs ------------- if ( "$argv[$ac]" == "-in_dwi" ) then if ( $ac >= $#argv ) goto FAIL_MISSING_ARG @ ac += 1 set in_DSET = "$argv[$ac]" else if ( "$argv[$ac]" == "-mask" ) then if ( $ac >= $#argv ) goto FAIL_MISSING_ARG @ ac += 1 set in_MASK = "$argv[$ac]" # ------------- input vecmat and bval -------------- else if ( "$argv[$ac]" == "-in_col_matA" ) then if ( $ac >= $#argv ) goto FAIL_MISSING_ARG set invecmat[1] = $argv[$ac] @ ac += 1 set invecmat[2] = "$argv[$ac]" else if ( "$argv[$ac]" == "-in_col_matT" ) then if ( $ac >= $#argv ) goto FAIL_MISSING_ARG set invecmat[1] = $argv[$ac] @ ac += 1 set invecmat[2] = "$argv[$ac]" else if ( "$argv[$ac]" == "-in_col_vec" ) then if ( $ac >= $#argv ) goto FAIL_MISSING_ARG set invecmat[1] = $argv[$ac] @ ac += 1 set invecmat[2] = "$argv[$ac]" else if ( "$argv[$ac]" == "-in_row_vec" ) then if ( $ac >= $#argv ) goto FAIL_MISSING_ARG set invecmat[1] = $argv[$ac] @ ac += 1 set invecmat[2] = "$argv[$ac]" # not necessary; default is just empty else if ( "$argv[$ac]" == "-in_bvals" ) then if ( $ac >= $#argv ) goto FAIL_MISSING_ARG set invecmat[3] = $argv[$ac] @ ac += 1 set invecmat[4] = "$argv[$ac]" # --- control tracking, i.e., for infants, non-humans, etc. --- else if ( "$argv[$ac]" == "-alg_Thresh_FA" ) then if ( $ac >= $#argv ) goto FAIL_MISSING_ARG @ ac += 1 set min_FA = $argv[$ac] else if ( "$argv[$ac]" == "-alg_Thresh_Len" ) then if ( $ac >= $#argv ) goto FAIL_MISSING_ARG @ ac += 1 set min_tr_len = $argv[$ac] # ---------- other opts ----------------- # [PT: Oct 04, 2017] Don't need separate, just get from '-prefix ...' part # else if ( "$argv[$ac]" == "-outdir" ) then # if ( $ac >= $#argv ) goto FAIL_MISSING_ARG # @ ac += 1 # set odir = "$argv[$ac]" else if ( "$argv[$ac]" == "-outdir" ) then if ( $ac >= $#argv ) goto FAIL_MISSING_ARG echo "** ERROR: this option no longer exists-- just specify output" echo " directory from path part of '-prefix ...'" goto BAD_EXIT # [PT: July 31, 2017] useful if running a lot else if ( "$argv[$ac]" == "-wdir" ) then if ( $ac >= $#argv ) goto FAIL_MISSING_ARG @ ac += 1 set tmpdir = "$argv[$ac]" # [PT: Oct 04, 2017] Get $odir solely from this else if ( "$argv[$ac]" == "-prefix" ) then if ( $ac >= $#argv ) goto FAIL_MISSING_ARG @ ac += 1 set out_ff = "$argv[$ac]" set odir = `dirname "$out_ff"` else if ( "$argv[$ac]" == "-scale_out_1000" ) then set SOT = "-scale_out_1000" else if ( "$argv[$ac]" == "-do_clean" ) then set DO_CLEAN = 1 else echo "\n\n** ERROR: unexpected option #$ac = '$argv[$ac]'\n\n" goto BAD_EXIT endif @ ac += 1 end # ======================================================================= # ============================ ** SETUP ** ============================== # ======================================================================= # ============================ input files ============================== echo "++ Start script version: $version" # check for DWI set check = `3dinfo "$in_DSET"` if ( "$#check" == "0" ) then echo "** ERROR: can't find inset file $in_DSET !" goto BAD_EXIT endif # check mask, if one has been input if ( "$in_MASK" != "" ) then set check = `3dinfo "$in_MASK"` if ( "$#check" == "0" ) then echo "** ERROR: can't find inset file $in_MASK !" goto BAD_EXIT endif endif # [PT: Feb 1, 2018] New format of this if ( "$invecmat[2]" == "" ) then echo "** ERROR: no gradient/matrix file input?!" goto BAD_EXIT endif # default output dir, if nothing input. if ( ! -e "$odir" ) then echo "+* Output directory didn't exist. Trying to make '$odir' now." mkdir "$odir" endif # in case user wants to change this loc set tmpdir = "$odir/$tmpdir" if ( "$out_ff" == "" ) then set out_ff = "$odir/GradFlipTest_rec.txt" endif # make a temporary output directory if( -d "$tmpdir" ) then echo "++ Directory '$tmpdir' exists already." else echo "++ Making the temporary directory '$tmpdir' for the files." mkdir $tmpdir endif if ( ! -e "$tmpdir" ) then echo "** ERROR: Couldn't make working dir '$tmpdir' for some reason." goto BAD_EXIT endif # [PT: July 31, 2017]: echo the useful parts dumped to the terminal # into a text file set eee = ${out_ff:t} set out_echo = `basename $eee .dat` set out_echo = `basename $out_echo .txt` set out_echo = `basename $out_echo .log` set out_echo = `basename $out_echo .flip` set out_echo = `basename $out_echo .tex` set out_echo = "$odir/${out_echo}_echo.txt" echo "OUTECHO: $out_echo" # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - # # Start the looping # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - # printf "\n\n\n++ Starting to test...\n" set i = '1' # Setup names set out_MATA = "$tmpdir/$FNOMS[$i]_matA.txt" set out_DSET = "$tmpdir/$FNOMS[$i]_DSET.nii.gz" set out_DT = "$tmpdir/$FNOMS[$i]_DT" set out_TRK = "$tmpdir/$FNOMS[$i]_WB" set out_GRID = "$tmpdir/$FNOMS[$i]_WB_000.grid" set out_TRACT = "$tmpdir/$FNOMS[$i]_WB_000.niml.tract" set out_mask = "$tmpdir/mask.nii.gz" printf "\n\n++ Grad_o_Matting (just needed once).\n\n\n" # -------------- copy over grads and data sets ----------------------- 1dDW_Grad_o_Mat++ \ -overwrite \ -echo_edu \ -out_col_matA $out_MATA \ "$invecmat[1]" "$invecmat[2]" \ "$invecmat[3]" "$invecmat[4]" 3dcalc \ -echo_edu \ -a "$in_DSET" \ -expr "a" \ -prefix "$out_DSET" \ -overwrite # check that grad and vol numbers match set Nvm = `wc $out_MATA` set Nvol = `3dinfo -nv $out_DSET` if ( $Nvm[1] != $Nvol ) then echo "** ERROR! Number of grads ($Nvm[1]) != number of vols ($Nvol)!" goto BAD_EXIT else echo "++ Good: found the same number of grads ($Nvm[1]) and vols ($Nvol)" endif # Only do Grad_o_Mat with no flipping-- only fit tensors once, and # then flip eigenvectors individually if ( "$in_MASK" != "" ) then echo "++ Copying mask" 3dcalc \ -echo_edu \ -a $in_MASK \ -expr "a" \ -prefix $out_mask \ -overwrite else echo "++ Automasking (-> $out_mask)" # mask based on 0th brick of 1st set # [June 6, 2017]: updated because TORTOISE can have large # neg values in DWIs, so take abs value 3dAutomask \ -echo_edu \ -overwrite \ -prefix $out_mask \ "3dcalc( -a ${out_DSET}[0] -expr (a*step(a)) )" #$out_DSET'[0]' endif printf "\n\n++ Quick tensor fitting (just needed once):" # Tensor fitting of non-flipped 3dDWItoDT \ -echo_edu \ -eigs -linear -sep_dsets \ -mask $out_mask \ $SOT \ -prefix $out_DT \ -bmatrix_FULL $out_MATA \ $out_DSET \ -overwrite set ext = `3dinfo -av_space ${out_DT}_FA*HEAD` printf "\n\n++ Make the new eigenvectors for each flip.\n\n\n" # Make other flips foreach j ( `seq 1 1 3` ) # flip_x (i.e., 0th brick in V? files) set out_DTx = "$tmpdir/$FNOMS[2]_DT_V${j}" 3dcalc \ -a ${out_DT}_V${j}$ext \ -expr 'a*(-equals(l,0)+equals(l,1)+equals(l,2))' \ -prefix $out_DTx \ -overwrite set out_DTy = "$tmpdir/$FNOMS[3]_DT_V${j}" 3dcalc \ -a ${out_DT}_V${j}$ext \ -expr 'a*(equals(l,0)-equals(l,1)+equals(l,2))' \ -prefix $out_DTy \ -overwrite set out_DTz = "$tmpdir/$FNOMS[4]_DT_V${j}" 3dcalc \ -a ${out_DT}_V${j}$ext \ -expr 'a*(equals(l,0)+equals(l,1)-equals(l,2))' \ -prefix $out_DTz \ -overwrite end printf "\n\n++ Now for some tracking with each set of eigenvectors.\n\n" foreach j ( `seq 1 1 4` ) set out_DTx_niml = "$tmpdir/$FNOMS[$j].niml.opts" # Clear file first, then fill printf '' > $out_DTx_niml printf "> $out_DTx_niml printf "dti_V1='$tmpdir/$FNOMS[$j]_DT_V1$ext'\n" >> $out_DTx_niml printf "dti_V2='$tmpdir/$FNOMS[$j]_DT_V2$ext'\n" >> $out_DTx_niml printf "dti_V3='$tmpdir/$FNOMS[$j]_DT_V3$ext'\n" >> $out_DTx_niml printf "dti_FA='${out_DT}_FA$ext'\n" >> $out_DTx_niml printf "dti_MD='${out_DT}_MD$ext'\n" >> $out_DTx_niml printf "dti_L1='${out_DT}_L1$ext'\n" >> $out_DTx_niml printf "dti_RD='${out_DT}_RD$ext' />\n" >> $out_DTx_niml # long tracts set out_TRK = "$tmpdir/$FNOMS[$j]_WB" 3dTrackID -mode DET \ -echo_edu \ -no_indipair_out \ -logic OR \ -alg_Nseed_X 1 \ -alg_Nseed_Y 1 \ -alg_Nseed_Z 1 \ -alg_Thresh_Len $min_tr_len \ -alg_Thresh_FA $min_fa \ -dti_list $out_DTx_niml \ -netrois $out_mask \ -mask $out_mask \ -prefix $out_TRK \ -overwrite end ######################################################################## # Summary stuff ######################################################################## foreach i ( `seq 1 1 4` ) # always underlay single FA, which wouldn't change set out_DT = "$tmpdir/$FNOMS[1]_DT" set out_TRK = "$tmpdir/$FNOMS[$i]_WB" set out_GRID = "$tmpdir/$FNOMS[$i]_WB_000.grid" set out_TRACT = "$tmpdir/$FNOMS[$i]_WB_000.niml.tract" # grep for the numbers of tracts in each case set HasNT = `grep -A 1 "\# NT" $out_GRID` printf "++ Checking GRID file: $out_GRID\n" set ALL_ARR = ( $ALL_ARR "$HasNT[3]" ) # for building output command set SUMA_CALL = "${SUMA_CALL}# For '$FLIPS[$i]':\n" set SUMA_CALL = "${SUMA_CALL}suma -vol ${out_DT}_FA*HEAD -tract $out_TRACT &\n" end # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - # # Find the max of the tracts # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - # printf "++ Calculate max of tracts.\n" set L = `printf $#ALL_ARR` set MaxI = "1" set MaxV = "$ALL_ARR[$MaxI]" # prepare the echo file, clear out anything in their if this is being # run again printf "" > $out_echo echo "# This command:" >> $out_echo echo "$this_prog $argv\n\n" >> $out_echo foreach i ( `seq 1 1 $L` ) if ( "$ALL_ARR[$i]" > "$MaxV" ) then set MaxI = "$i" set MaxV = "$ALL_ARR[$MaxI]" endif printf "\n The number of tracks for '$FLIPS[$i]' is: $FLIP_sp[$i]:q %12s" \ "$ALL_ARR[$i]" >> $out_echo end # [PT: Aug 24, 2017] overwrite any previous run, if it exists printf "" > $out_ff echo "$FLIPS[$MaxI]" >> $out_ff printf "\n\n -->" >> $out_echo printf "Therefore, I *guess* that the best flip is: " >> $out_echo printf "'$FLIPS[$MaxI]'\n\n" >> $out_echo printf "\n" >> $out_echo printf " NB: If that isn't a *very* clear winner, then check \n" >> $out_echo printf " for any problems, such as DICOM conversion, etc.\n" >> $out_echo printf "\n" >> $out_echo printf "That recommendation is stored in '$out_ff'\n\n" >> $out_echo printf "Volume + tracking results are in directory '$tmpdir'\n" >> $out_echo echo "" >> $out_echo printf "You may (should!) check the tracking results to verify:\n" >> $out_echo printf "$SUMA_CALL\n" >> $out_echo printf "DONE.\n\n" >> $out_echo echo "\n" echo "# --------------------- RESULTS ---------------------------\n" cat $out_echo echo "\n" echo "# ---------------------------------------------------------" if ( $DO_CLEAN ) then echo "++ ... and cleaning up by removing the temp dir: $tmpdir" \rm -rf $tmpdir endif goto GOOD_EXIT # ======================================================================== # ======================================================================== SHOW_HELP: cat << EOF #------------------------------------------------------------------------ Simple script to test what 'flip', if any, should likely be performed for a data set when using 1dDW_Grad_o_Mat++. **Majorly updated in Jan, 2017-- otherwise you wouldn't even be reading this help file description!** When using this function and looking at the number of tracts per flip, there should be a *very* clear winner. If there isn't, then probably something is not correct in the data (something inconsistent in bvals or bvecs, large noise, etc.). Please make sure to look at the results in SUMA when prompted at the end, to make sure that everything makes sense! ver $version; revision date $rev_dat. Written by PA Taylor (NIH). # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - OUTPUT On a good day, this function will: + Recommend using either '-no_flip', or one of the {-flip_x|-flip_y|-flip_z} options for 1dDW_Grad_o_Mat++. + It will store this snippet of code in a file called $out_ff (default name), which the User could be used in scripting later. + It will produce a temporary working directory called '_tmp_TESTFLIP/' to store intermediate files, of which there are many (could be wiped away with '-do_clean'). + It will also prompt you, O User, to visually check the tract results with some simple example scripts (some day it might automatically make snapshots!). # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - COMMAND $this_prog \ -in_dwi DWI \ { -in_row_vec | -in_col_vec | -in_col_matA | -in_col_matT } FF \ { -mask MASK } \ { -in_bvals BB } \ { -alg_Thresh_FA X } \ { -alg_Thresh_Len L } \ { -prefix PPP } \ { -scale_out_1000 } \ { -wdir WWW } \ { -do_clean } USAGE (*must* input 1 set of DWIs *and* 1 set of grads-- choice of format): -in_dwi DWI :set of DWIs (N total volumes) -in_row_vec FF :set of row-wise gradient vectors -in_col_vec FF :set of column-wise gradient vectors -in_col_matA FF :set of column-wise g- or b-matrix elements ("AFNI"-style format, "diagonal-first") -in_col_matT FF :set of column-wise g- or b-matrix elements ("TORTOISE"-style format, "row-first") -mask MASK :option mask (probably whole brain); otherwise, automasking is performed -in_bvals BB :can input bvals, as in 1dDW_Grad_o_Mat++, if necessary (but shouldn't be necessary?) -alg_Thresh_FA X :set minimum FA value for tracking (default X=0.2 as for adult, healthy WM parenchyma) -alg_Thresh_Len L :set minimum tract length to require to keep a tract when propagating (default L=30mm ; probably want it to be a bit on the longside for clear counting and comparison) -prefix PPP :output name of text file that stores recommended flip opt (default is $out_ff). This option is now also used to determine the directory for all outputs of this program, via the path of PPP. NB: The previous, separate option for specifying output directory was '-outdir OUT', but this no longer is used; the path of an output directory is specified by taking the path-part of the '-prefix PPP' input. -scale_out_1000 :as in 3dDWItoDT. Probably not necessary, since we are just checking out trackability -wdir WWW :rename working directory output; useful if running multiple iterations. Default: ${tmpdir}. NB: WWW should *only* be the name of the directory, not contain path info-- the location of WWW is just determined by the path for output, which comes from the path part of PPP/ -do_clean :remove temporary directory # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - EXAMPLES $this_prog \ -in_dwi DWI.nii.gz \ -in_col_matA BMTXT_AFNI.txt or (perhaps if scanning infants, who have less developed myelin) $this_prog \ -in_dwi DWI.nii.gz \ -in_col_vec GRADS.txt \ -mask mask_DWI.nii.gz \ -alg_Thresh_FA 0.1 # ----------------------------------------------------------------------- EOF goto GOOD_EXIT SHOW_VERSION: echo "version $version (${rev_dat})" goto GOOD_EXIT FAIL_MISSING_ARG: echo "** ERROR! Missing an argument after option flag: '$argv[$ac]'" goto BAD_EXIT BAD_EXIT: exit 1 GOOD_EXIT: exit 0