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

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

7
# Written by Saad Jbabdi, Stam Sotiropoulos & Shaun Warrington (based on Marius de Groot autoPtx code)
Saad Jbabdi's avatar
Saad Jbabdi committed
8
9
# 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
# Location of xtract data
14
datadir=$FSLDIR/data/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

    Compulsory arguments:

25
26
27
       -bpx <folder>                          Path to bedpostx folder
       -out <folder>                          Path to output folder
       -species <SPECIES>                     One of HUMAN or MACAQUE
Saad Jbabdi's avatar
Saad Jbabdi committed
28
29

    Optional arguments:
30
31
       -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)
32
       -p <folder>                            Protocols folder (all masks in same standard space) (Default=$FSLDIR/data/xtract_data/<SPECIES>)
33
34
35
36
37
38
39
40
41
42
       -stdwarp <std2diff> <diff2std>         Standard2diff and Diff2standard transforms (Default=bedpostx_dir/xfms/{standard2diff,diff2standard})
       -gpu                                   Use GPU version
       -res <mm>                              Output resolution (Default=same as in protocol folders unless '-native' used)
       -ptx_options <options.txt>	            Pass extra probtrackx2 options as a text file to override defaults, e.g. --steplength=0.2 --distthresh=10)

       And EITHER:
       -native                                Run tractography in native (diffusion) space

       OR:
       -ref <refimage> <diff2ref> <ref2diff>  Reference image for running tractography in reference space, Diff2Reference and Reference2Diff transforms
Saad Jbabdi's avatar
Saad Jbabdi committed
43
44
45
46
47

EOF
    exit 1
}

Saad Jbabdi's avatar
Saad Jbabdi committed
48
49
50
51
52

Splash (){

cat <<EOF

53
 __  _______ ____      _    ____ _____
Saad Jbabdi's avatar
Saad Jbabdi committed
54
 \ \/ /_   _|  _ \    / \  / ___|_   _|
55
56
57
58
  \  /  | | | |_) |  / _ \| |     | |
  /  \  | | |  _ <  / ___ \ |___  | |
 /_/\_\ |_| |_| \_\/_/   \_\____| |_|

Saad Jbabdi's avatar
Saad Jbabdi committed
59
60
61
62
63
64
EOF

}

Splash

Saad Jbabdi's avatar
Saad Jbabdi committed
65
66
67
68
69
70
71
72
73
[ "$1" = "" ] && Usage


# Set default options
bpx=""
out=""
str=""
p=""
std2diff=""
74
ptx_opts=""
75
stdref=""
Saad Jbabdi's avatar
Saad Jbabdi committed
76
77
gpu=0
nat=0
78
79
80
ref=""
diff2ref=""
ref2diff=""
Saad Jbabdi's avatar
Saad Jbabdi committed
81
82
spec=""
res=-1
83
84
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
85
86
87
88
89
90
91
92

# 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
93
	-species) spec=$2;shift;;  # converts to uppercase
Saad Jbabdi's avatar
Saad Jbabdi committed
94
95
96
	-stdwarp) std2diff=$2;diff2std=$3;shift;shift;;
	-gpu) gpu=1;;
	-native) nat=1;;
97
  -ref) ref=$2;diff2ref=$3;ref2diff=$4;shift;shift;shift;;
Saad Jbabdi's avatar
Saad Jbabdi committed
98
	-res) res=$2;shift;;
99
  -list) list=1;shift;;
100
101
	-ptx_options) ptx_opts=`cat $2`;shift;;

Saad Jbabdi's avatar
Saad Jbabdi committed
102
103
104
105
	*) echo "Unknown option '$1'";exit 1;;
    esac
    shift
done
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125

# 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
126
127
128
129
130
131
132
133
134
135
136
137
# 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
138
139
140
141
142
143
144
145
146
147
148
149
150
elif [ ! "$ref" == "" ];then
    if [ "$diff2ref" == "" ] || [ "$ref2diff" == "" ];then
  echo "If running in ref space, you must specify '-ref <refimage> <diff2ref> <ref2diff>'"
  exit 1
    fi
    if [ `$FSLDIR/bin/imtest $diff2ref` -eq 0 ];then
  echo "Image $diff2ref not found."
  exit 1
    fi
    if [ `$FSLDIR/bin/imtest $ref2diff` -eq 0 ];then
  echo "Image $ref2diff not found."
  exit 1
    fi
Saad Jbabdi's avatar
Saad Jbabdi committed
151
152
fi

Saad Jbabdi's avatar
Saad Jbabdi committed
153
echo SPECIES $spec
Saad Jbabdi's avatar
Saad Jbabdi committed
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176

# 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

177
# Check which species and protocols to run
Saad Jbabdi's avatar
Saad Jbabdi committed
178
if [ "$spec" == "" ];then
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
  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
194
    fi
195
196
197
elif [ "$spec" == "MACAQUE" ];then
    stdref=$datadir/standard/F99/mri/struct_brain
    strdef=$datadir/Macaque/structureList
Saad Jbabdi's avatar
Saad Jbabdi committed
198
    if [ "$p" == "" ];then
199
200
201
202
203
204
205
206
207
      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
208
209
    fi
else
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
    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
233
    else
234
235
236
237
238
      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
239
    fi
240
241
242
243
244
245
246
247
248
249
250
  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
251
    echo " -- getting 'nsamples' from default structure file"
252
253
    seedget=1
  fi
Saad Jbabdi's avatar
Saad Jbabdi committed
254
255
fi

256
257
258
259
260
261
262
# Check space option
if [ $nat -eq 1 ] && [ ! "$ref" == "" ];then
    echo "You have selected the native space and ref space options"
    echo "Must select EITHER '-native', '-ref <refimage> <diff2ref> <ref2diff>', OR use the default standard space"
    errflag=1
fi

Saad Jbabdi's avatar
Saad Jbabdi committed
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
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 "

278
279
280
# Add any user-defined ptx options
opts=" $opts $ptx_opts"

281
282
283
284
285
286
287
288
289
290
291
if [ "$nat" -eq 0 ] && [ "$ref" == "" ];then
    opts="$opts --seedref=$stdref --xfm=$std2diff --invxfm=$diff2std "
elif [ ! "$ref" == "" ]; then
    opts="$opts --seedref=$ref --xfm=$ref2diff --invxfm=$diff2ref "
fi

# If running in reference space, combine std2diff and diff2ref for std2ref
if [ ! "$ref" == "" ];then
    echo " -- combining standard-to-diffusion and diffusion-to-reference transforms"
    std2ref=$out/standard2ref
    $FSLDIR/bin/convertwarp -o $std2ref -r $ref --warp1=$std2diff --warp2=$diff2ref
Saad Jbabdi's avatar
Saad Jbabdi committed
292
293
294
295
296
297
298
fi

# Loop over structures
commands=$out/commands.txt
rm -rf $commands
echo "Preparing submission script..."
while read structstring; do
299
300
301
  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
302
303
304
	# do nothing
	foo=0
	#echo "----- Skip line $structstring -----"
305
  elif [ "$struct" == "" ];then
Saad Jbabdi's avatar
Saad Jbabdi committed
306
307
308
	# do nothing
	foo=0
	#echo "----- Skip empty line -----"
309
  else
Saad Jbabdi's avatar
Saad Jbabdi committed
310
311
312
	#echo "autoTrack $struct"
	mkdir -p $out/tracts/$struct

313
314
315
316
317
318
319
320
321
322
323
  # 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
324
	nseed=$(echo "scale=0; 1000 * ${nseed} / 1"|bc)
325

Saad Jbabdi's avatar
Saad Jbabdi committed
326
327
	maskdir=$p/$struct

328
329
	#  DEALING WITH RESAMPLING --
	# Pick space to run tractography in (diffusion or standard)
Saad Jbabdi's avatar
Saad Jbabdi committed
330
	if [ "$nat" -eq 1 ];then
331
	    echo "${struct} -- transforming masks into native space"
Saad Jbabdi's avatar
Saad Jbabdi committed
332
333
334
335
336
337
338
	    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"
339
	    done
340
341
342
343
344
345
346
347
348
349
  elif [ ! "$ref" == "" ];then
      echo "${struct} -- transforming masks into ref space"
      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 $std2ref -r $ref -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"
      done
Saad Jbabdi's avatar
Saad Jbabdi committed
350
	else
351
	    for m in seed stop exclude;do
Saad Jbabdi's avatar
Saad Jbabdi committed
352
353
354
355
		if [ $res -gt 0 ];then
		    # Resample at a different resolution
		    mkdir -p $out/masks/$struct
		    if [ `$FSLDIR/bin/imtest $maskdir/$m` -eq 1 ];then
356
			$FSLDIR/bin/flirt -in $maskdir/$m -out $out/masks/$struct/$m -applyisoxfm $res -ref $maskdir/$m
Saad Jbabdi's avatar
Saad Jbabdi committed
357
			$FSLDIR/bin/fslmaths $out/masks/$struct/$m -thr 0.1 -bin $out/masks/$struct/$m -odt char
358
359
		    fi
		    eval "${m}=$out/masks/$struct/$m"
Saad Jbabdi's avatar
Saad Jbabdi committed
360
361
362
363
		else
		    eval "${m}=$maskdir/$m"
		fi
	    done
364
	fi
Saad Jbabdi's avatar
Saad Jbabdi committed
365
366
367
368
369
370
371
372
373
374
375

	# 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
376
377
378
379
380
381
382
  elif [ ! "$ref" == "" ];then
      for tfile in $targets;do
    t=`basename $tfile`
    $FSLDIR/bin/applywarp -i $tfile -o $out/masks/$struct/$t -w $std2ref -r $ref -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
Saad Jbabdi's avatar
Saad Jbabdi committed
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
	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
399

Saad Jbabdi's avatar
Saad Jbabdi committed
400
401
402
403
404
405
406
407
408
409
        # 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 "
410

Saad Jbabdi's avatar
Saad Jbabdi committed
411
412
413
414
415
416
417
	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"

418
419

	# Does the protocol define a second run with inverted seed / target masks?
Saad Jbabdi's avatar
Saad Jbabdi committed
420
421
422
423
424
425
426
427
428
429
430
431
432
433
	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)
434
435
		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
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
		# 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
461
	fsl_sub -T 2160 -q long.q -l $out/logs -N xtract -t $commands
Saad Jbabdi's avatar
Saad Jbabdi committed
462
    else
463
	fsl_sub -T 300 -q $FSLGECUDAQ -l $out/logs -N xtract $commands
Saad Jbabdi's avatar
Saad Jbabdi committed
464
465
466
467
468
469
    fi
else   # If no SGE, run locally
    sh $commands
fi

#EOF