xtract 11.8 KB
Newer Older
Saad Jbabdi's avatar
Saad Jbabdi committed
1
2
3
4
5
6
7
8
9
#!/bin/bash

#   Copyright (C) 2019 University of Oxford
#
#   SHCOPYRIGHT

# Written by Saad Jbabdi & Stam Sotiropoulos (based on Marius de Groot autoPtx code)
# Protocols created by Rogier Mars et al.

10
# Location of probtrackx2_gpu binary
Saad Jbabdi's avatar
Saad Jbabdi committed
11
12
ptxbin_gpu=$FSLDIR/bin/probtrackx2_gpu

Saad Jbabdi's avatar
Saad Jbabdi committed
13
14
# Location of xtract data
datadir=$FSLDIR/etc/xtract_data
Saad Jbabdi's avatar
Saad Jbabdi committed
15
16
17
18

Usage() {
    cat << EOF

19
Usage:
20
21
    xtract -bpx <bedpostX_dir> -out <outputDir> -species <SPECIES> [options]
    xtract -list
Saad Jbabdi's avatar
Saad Jbabdi committed
22
23
24
25
26
27
28
29

    Compulsory arguments:

       -bpx <folder>                     Path to bedpostx folder
       -out <folder>                     Path to output folder
       -species <SPECIES>                One of HUMAN or MACAQUE

    Optional arguments:
30
31
32
       -list                             List the tract names used in XTRACT
       -str <file>                       Structures file (format: <tractName> per line OR format: <tractName> [samples=1], 1 means 1000, '#' to skip lines)
       -p   <folder>                     Protocols folder (all masks in same standard space) (Default=$FSLDIR/etc/xtract_data/<SPECIES>)
33
34
       -stdwarp <std2diff> <diff2std>    Standard2diff and Diff2standard transforms (Default=bedpostx_dir/xfms/{standard2diff,diff2standard})
       -gpu                              Use GPU version
Saad Jbabdi's avatar
Saad Jbabdi committed
35
36
       -native                           Run tractography in native (diffusion) space
       -res <mm>                         Output resolution (Default=same as in protocol folders unless '-native' used)
37
       -ptx_options <options.txt>	       Pass extra probtrackx2 options as a text file to override defaults, e.g. --steplength=0.2 --distthresh=10)
Saad Jbabdi's avatar
Saad Jbabdi committed
38
39
40
41
42

EOF
    exit 1
}

Saad Jbabdi's avatar
Saad Jbabdi committed
43
44
45
46
47

Splash (){

cat <<EOF

48
 __  _______ ____      _    ____ _____
Saad Jbabdi's avatar
Saad Jbabdi committed
49
 \ \/ /_   _|  _ \    / \  / ___|_   _|
50
51
52
53
  \  /  | | | |_) |  / _ \| |     | |
  /  \  | | |  _ <  / ___ \ |___  | |
 /_/\_\ |_| |_| \_\/_/   \_\____| |_|

Saad Jbabdi's avatar
Saad Jbabdi committed
54
55
56
57
58
59
EOF

}

Splash

Saad Jbabdi's avatar
Saad Jbabdi committed
60
61
62
63
64
65
66
67
68
[ "$1" = "" ] && Usage


# Set default options
bpx=""
out=""
str=""
p=""
std2diff=""
69
ptx_opts=""
70
stdref=""
Saad Jbabdi's avatar
Saad Jbabdi committed
71
72
73
74
gpu=0
nat=0
spec=""
res=-1
75
76
list=0
seedget=0 # if user specifies -str and format is <tractName>, try to get nseeds from default structureList (i.e. if using a subset of XTRACT tracts)
Saad Jbabdi's avatar
Saad Jbabdi committed
77
78
79
80
81
82
83
84

# Parse command-line arguments
while [ ! -z "$1" ];do
    case "$1" in
	-bpx) bpx=$2;shift;;
	-out) out=$2;shift;;
	-str) str=$2;shift;;
	-p)   p=$2;shift;;
Saad Jbabdi's avatar
Saad Jbabdi committed
85
	-species) spec=$2;shift;;  # converts to uppercase
Saad Jbabdi's avatar
Saad Jbabdi committed
86
87
88
89
	-stdwarp) std2diff=$2;diff2std=$3;shift;shift;;
	-gpu) gpu=1;;
	-native) nat=1;;
	-res) res=$2;shift;;
90
  -list) list=1;shift;;
91
92
	-ptx_options) ptx_opts=`cat $2`;shift;;

Saad Jbabdi's avatar
Saad Jbabdi committed
93
94
95
96
	*) echo "Unknown option '$1'";exit 1;;
    esac
    shift
done
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116

# list argument
if [ "$list" -eq 1 ];then
    tractNames=""
    echo ""
    echo "Tract names:"
    while read structstring; do
      struct=`echo $structstring | awk '{print $1}'`
      if [ "${struct:0:1}" == "#" ];then
    	foo=0
      elif [ "$struct" == "" ];then
    	foo=0
      else
      tractNames="${tractNames} ${struct},"
      fi
    done < $datadir/Human/structureList
    echo ${tractNames%,}
    exit 1
fi

Saad Jbabdi's avatar
Saad Jbabdi committed
117
118
119
120
121
122
123
124
125
126
127
128
129
130
# Default warps
if [ "$std2diff" == "" ];then
    std2diff=$bpx/xfms/standard2diff
    diff2std=$bpx/xfms/diff2standard
    if [ `$FSLDIR/bin/imtest $std2diff` -eq 0 ];then
	echo "Image $std2diff not found."
	exit 1
    fi
    if [ `$FSLDIR/bin/imtest $diff2std` -eq 0 ];then
	echo "Image $diff2std not found."
	exit 1
    fi
fi

Saad Jbabdi's avatar
Saad Jbabdi committed
131
echo SPECIES $spec
Saad Jbabdi's avatar
Saad Jbabdi committed
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154

# GPU stuff
if [ $gpu -eq 0 ];then
    ptxbin=$FSLDIR/bin/probtrackx2
else
    # Temp location of CUDA code
    ptxbin=${ptxbin_gpu}
fi

# Check compulsory arguments
errflag=0
if [ "$bpx" == "" ];then
    echo "Must set compulsory argument '-bpx'"
    errflag=1
elif [ ! -d $bpx ];then
    echo "Bedpostx folder $bpx not found"
    errflag=1
fi
if [ "$out" == "" ];then
    echo "Must set compulsory argument '-out'"
    errflag=1
fi

155
# Check which species and protocols to run
Saad Jbabdi's avatar
Saad Jbabdi committed
156
if [ "$spec" == "" ];then
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
  echo "Must set compulsory argument '-species'"
  errflag=1
elif [ "$spec" == "HUMAN" ];then
    stdref=$FSLDIR/data/standard/MNI152_T1_1mm
    strdef=$datadir/Human/structureList
    if [ "$p" == "" ];then
      p=$datadir/Human
      if [ "$str" == "" ];then
        str=$p/structureList
      fi
    elif [ ! "$p" == "" ];then
      if [ "$str" == "" ];then
        echo "If selecting a protocol folder, must set argument '-str'"
        errflag=1
      fi
Saad Jbabdi's avatar
Saad Jbabdi committed
172
    fi
173
174
175
elif [ "$spec" == "MACAQUE" ];then
    stdref=$datadir/standard/F99/mri/struct_brain
    strdef=$datadir/Macaque/structureList
Saad Jbabdi's avatar
Saad Jbabdi committed
176
    if [ "$p" == "" ];then
177
178
179
180
181
182
183
184
185
      p=$datadir/Macaque
      if [ "$str" == "" ];then
        str=$p/structureList
      fi
    elif [ ! "$p" == "" ];then
      if [ "$str" == "" ];then
        echo "If selecting a protocol folder, must set argument '-str'"
        errflag=1
      fi
Saad Jbabdi's avatar
Saad Jbabdi committed
186
187
    fi
else
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
    echo "Species must be one of HUMAN or MACAQUE"
    errflag=1
fi

# Check that -str and -p exist
if [ ! -d $p ];then
  echo "Protocol folder $p not found"
  errflag=1
fi
if [ ! -f $str ];then
  echo "Structures files $str not found"
  errflag=1
fi

# Check -str file format
tchk=()
if [ ! "$str" == "$strdef" ];then
  while read structstring; do
    struct=`echo $structstring | awk '{print $1}'`
    if [ "${struct:0:1}" == "#" ];then
  	foo=0
    elif [ "$struct" == "" ];then
  	foo=0
Saad Jbabdi's avatar
Saad Jbabdi committed
211
    else
212
213
214
215
216
      if [ "`echo $structstring | awk '{print $2}'`" == "" ];then
      tchk+=("1") # if empty, then 1 - get nseeds from default -str file
      else
      tchk+=("0")
      fi
Saad Jbabdi's avatar
Saad Jbabdi committed
217
    fi
218
219
220
221
222
223
224
225
226
227
228
229
230
231
  done < $str
  IFS=$'\n'; tchk=($(sort <<<"${tchk[*]}")); unset IFS
  if [ ! "${tchk[0]}" -eq "${tchk[${#tchk[@]}-1]}" ];then
    echo ""
    echo "-str file format is inconsistent. Format should be either:"
    echo "<tractName> per line"
    echo "OR"
    echo "<tractName> [samples=1] per line"
    echo "samples=1, 1 means 1000. Use '#' to skip lines"
    errflag=1
  elif [ ${tchk[0]} -eq 1 ];then
    echo "Using subset of XTRACT protocols and getting 'nsamples' from default structure file."
    seedget=1
  fi
Saad Jbabdi's avatar
Saad Jbabdi committed
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
fi

if [ "$errflag" -eq 1 ];then
    echo ""
    echo "Exit without doing anything.."
    exit 1
fi

# Create output folders
mkdir -p $out
mkdir -p $out/logs
mkdir -p $out/tracts

# Set common ptx options
opts=" -s $bpx/merged -m $bpx/nodif_brain_mask -V 1"
opts=" $opts --loopcheck --forcedir --opd --ompl --seedref=$stdref --sampvox=1 --randfib=1 "

249
250
251
# Add any user-defined ptx options
opts=" $opts $ptx_opts"

Saad Jbabdi's avatar
Saad Jbabdi committed
252
253
254
255
256
257
258
259
260
if [ "$nat" -eq 0 ];then
    opts="$opts --xfm=$std2diff --invxfm=$diff2std "
fi

# Loop over structures
commands=$out/commands.txt
rm -rf $commands
echo "Preparing submission script..."
while read structstring; do
261
262
263
  struct=`echo $structstring | awk '{print $1}'`
  # skip empty lines and lines that start with '#'
  if [ "${struct:0:1}" == "#" ];then
Saad Jbabdi's avatar
Saad Jbabdi committed
264
265
266
	# do nothing
	foo=0
	#echo "----- Skip line $structstring -----"
267
  elif [ "$struct" == "" ];then
Saad Jbabdi's avatar
Saad Jbabdi committed
268
269
270
	# do nothing
	foo=0
	#echo "----- Skip empty line -----"
271
  else
Saad Jbabdi's avatar
Saad Jbabdi committed
272
273
274
	#echo "autoTrack $struct"
	mkdir -p $out/tracts/$struct

275
276
277
278
279
280
281
282
283
284
285
  # if running a subset of tracts and -str format is <tractName>, check for nseed in str file
  if [ "$seedget" -eq 1 ]; then
      nseed=`grep -w ${struct} $strdef | awk '{print $2}'`
      if [ "$nseed" == "" ];then
        echo "Couldn't find number of samples for '$struct'. Exiting now."
        exit 1
      fi
  else
	    nseed=`echo $structstring | awk '{print $2}'`
	    if [ "$nseed" == "" ];then nseed=1; fi
  fi
Saad Jbabdi's avatar
Saad Jbabdi committed
286
	nseed=$(echo "scale=0; 1000 * ${nseed} / 1"|bc)
287

Saad Jbabdi's avatar
Saad Jbabdi committed
288
289
	maskdir=$p/$struct

290
291
	#  DEALING WITH RESAMPLING --
	# Pick space to run tractography in (diffusion or standard)
Saad Jbabdi's avatar
Saad Jbabdi committed
292
	if [ "$nat" -eq 1 ];then
293
	    echo " -- transforming masks into native space"
Saad Jbabdi's avatar
Saad Jbabdi committed
294
295
296
297
298
299
300
	    mkdir -p $out/masks/$struct
	    for m in seed stop exclude;do
		if [ `$FSLDIR/bin/imtest $maskdir/$m` -eq 1 ];then
		    $FSLDIR/bin/applywarp -i $maskdir/$m -o $out/masks/$struct/$m -w $std2diff -r $bpx/nodif_brain_mask -d float
		    $FSLDIR/bin/fslmaths $out/masks/$struct/$m -thr 0.1 -bin $out/masks/$struct/$m -odt char
		fi
		eval "${m}=$out/masks/$struct/$m"
301
	    done
Saad Jbabdi's avatar
Saad Jbabdi committed
302
	else
303
	    for m in seed stop exclude;do
Saad Jbabdi's avatar
Saad Jbabdi committed
304
305
306
307
		if [ $res -gt 0 ];then
		    # Resample at a different resolution
		    mkdir -p $out/masks/$struct
		    if [ `$FSLDIR/bin/imtest $maskdir/$m` -eq 1 ];then
308
			$FSLDIR/bin/flirt -in $maskdir/$m -out $out/masks/$struct/$m -applyisoxfm $res -ref $maskdir/$m
Saad Jbabdi's avatar
Saad Jbabdi committed
309
			$FSLDIR/bin/fslmaths $out/masks/$struct/$m -thr 0.1 -bin $out/masks/$struct/$m -odt char
310
311
		    fi
		    eval "${m}=$out/masks/$struct/$m"
Saad Jbabdi's avatar
Saad Jbabdi committed
312
313
314
315
		else
		    eval "${m}=$maskdir/$m"
		fi
	    done
316
	fi
Saad Jbabdi's avatar
Saad Jbabdi committed
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343

	# Deal with targets (in cases where there may be more than one)
	targets=`imglob $maskdir/target*`
	targetfile=$out/tracts/$struct/targets.txt
	if [ "$nat" -eq 1 ];then
	    for tfile in $targets;do
		t=`basename $tfile`
		$FSLDIR/bin/applywarp -i $tfile -o $out/masks/$struct/$t -w $std2diff -r $bpx/nodif_brain_mask -d float
		$FSLDIR/bin/fslmaths $out/masks/$struct/$t -thr 0.1 -bin $out/masks/$struct/$t -odt char
	    done
	    echo $out/masks/$struct/target* > $targetfile
	else
	    if [ $res -gt 0 ];then
		# Resample at a different resolution
		for tfile in $targets;do
		    t=`basename $tfile`
		    $FSLDIR/bin/flirt -in $tfile -out $out/masks/$struct/$t -applyisoxfm $res -ref $tfile
		    $FSLDIR/bin/fslmaths $out/masks/$struct/$t -thr 0.1 -bin $out/masks/$struct/$t -odt char
		done
		echo $out/masks/$struct/target* > $targetfile
	    else
		echo $targets > $targetfile
	    fi
	fi

	# Get generic options
	o=$opts
344

Saad Jbabdi's avatar
Saad Jbabdi committed
345
346
347
348
349
350
351
352
353
354
        # Add inclusion/exclusion masks
	if [ `$FSLDIR/bin/imtest $stop` -eq 1 ];then
	    o="$o --stop=$stop"
	fi
	if [ `$FSLDIR/bin/imtest $exclude` -eq 1 ];then
	    o="$o --avoid=$exclude"
	fi

	# Add seed/target
	o1="$o --nsamples=$nseed -x $seed "
355

Saad Jbabdi's avatar
Saad Jbabdi committed
356
357
358
359
360
361
362
	if [ "x${targets}" != "x" ];then #Add waypoints if there are any
       	    o1=" $o1 --waypoints=$targetfile "
	fi

        # Outputs
	o1=" $o1 -o density --dir=$out/tracts/$struct"

363
364

	# Does the protocol define a second run with inverted seed / target masks?
Saad Jbabdi's avatar
Saad Jbabdi committed
365
366
367
368
369
370
371
372
373
374
375
376
377
378
	if [ -e $maskdir/invert ]; then  #Invert-mode
	    if 	[ `$FSLDIR/bin/imtest $maskdir/target.nii.gz` -eq 1 ];then 	# Check if a target.nii.gz image exists when invert option has been selected.
		mkdir -p $out/tracts/$struct/tractsInv
		if [ `$FSLDIR/bin/imtest $out/masks/$struct/target.nii.gz` -eq 1 ]; then
		    target=$out/masks/$struct/target
		else
		    target=$maskdir/target
		fi
		o2="$o --nsamples=$nseed -x ${target} --waypoints=$seed -o density --dir=$out/tracts/$struct/tractsInv"

		# merge runs for forward and inverted tractography runs and then normalise (create commands but don't execute)
		mergecmd="$FSLDIR/bin/fslmaths $out/tracts/$struct/density -add $out/tracts/$struct/tractsInv/density $out/tracts/$struct/sum_density"

		#Add waypoints (create command but don't execute)
379
380
		addcmd="echo \"scale=5; \`cat $out/tracts/$struct/waytotal\` + \`cat $out/tracts/$struct/tractsInv/waytotal\` \"|bc > $out/tracts/$struct/sum_waytotal"

Saad Jbabdi's avatar
Saad Jbabdi committed
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
		# Waypoint normalisation (create command but don't execute)
		normcmd="$FSLDIR/bin/fslmaths $out/tracts/$struct/sum_density -div \`cat $out/tracts/$struct/sum_waytotal\` $out/tracts/$struct/densityNorm"

		# Append to command list
		echo "$ptxbin $o1; $ptxbin $o2; $mergecmd; $addcmd; $normcmd" >> $commands
	    else
		echo "Invert Option selected, but more than one target defined! A 'target.nii.gz' is expected. Exiting now"
		exit 1
	    fi
	else  	#No invert-mode
	    # Waypoint normalisation (create command but don't execute)
	    normcmd="$FSLDIR/bin/fslmaths $out/tracts/$struct/density -div \`cat $out/tracts/$struct/waytotal\` $out/tracts/$struct/densityNorm"

	    # Append to command list
	    echo "$ptxbin $o1; $normcmd" >> $commands
	fi
    fi

done < $str

chmod +x $commands

if [ "x$SGE_ROOT" != "x" ]; then  # Submit all commands to run in parallel on the cluster
    # One job per tract for a CPU cluster, one job for all tracts for a GPU cluster.
    if [ $gpu -eq 0 ];then
406
	fsl_sub -T 2160 -q long.q -l $out/logs -N xtract -t $commands
Saad Jbabdi's avatar
Saad Jbabdi committed
407
    else
408
	fsl_sub -T 300 -q $FSLGECUDAQ -l $out/logs -N xtract $commands
Saad Jbabdi's avatar
Saad Jbabdi committed
409
410
411
412
413
414
    fi
else   # If no SGE, run locally
    sh $commands
fi

#EOF