xtract 9.7 KB
Newer Older
Saad Jbabdi's avatar
Saad Jbabdi committed
1
2
3
4
5
6
7
8
9
10
11
12
#!/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.


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
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41

Usage() {
    cat << EOF

Usage: 
    xtract -bpx <bedpostX_dir> -out <outputDir> -str <structuresFile> -p <protocolsFolder> [options]
    xtract -bpx <bedpostX_dir> -out <outputDir> -species HUMAN [options]
    xtract -bpx <bedpostX_dir> -out <outputDir> -species MACAQUE [options]

    Compulsory arguments:

       -bpx <folder>                     Path to bedpostx folder
       -out <folder>                     Path to output folder
       
       And EITHER:
       -str <file>                       Structures file (format: <tractName> [samples=1], 1 means 1000, '#' to skip lines)
       -p   <folder>                     Protocols folder (all masks in same standard space)

       Or:
       -species <SPECIES>                One of HUMAN or MACAQUE

    Optional arguments:

       -stdwarp <std2diff> <diff2std>    Standard2diff and Diff2standard transforms (Default=bedpostx_dir/xfms/{standard2diff,diff2standard}) 
       -gpu                              Use GPU version 
       -native                           Run tractography in native (diffusion) space
       -res <mm>                         Output resolution (Default=same as in protocol folders unless '-native' used)
42
       -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
43
44
45
46
47

EOF
    exit 1
}

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

Splash (){

cat <<EOF

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

}

63
64
65
66
67
68
69
70
Warning (){
cat <<EOF
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
  WARNING!!!! MACAQUE TRACTS ARE A WORK IN PROGRESS STILL....
  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
EOF
}

Saad Jbabdi's avatar
Saad Jbabdi committed
71
72
Splash

Saad Jbabdi's avatar
Saad Jbabdi committed
73
74
75
76
77
78
79
80
81
[ "$1" = "" ] && Usage


# Set default options
bpx=""
out=""
str=""
p=""
std2diff=""
82
ptx_opts=""
Saad Jbabdi's avatar
Saad Jbabdi committed
83
84
85
86
87
88
89
90
91
92
93
94
95
stdref="$FSLDIR/data/standard/MNI152_T1_1mm"
gpu=0
nat=0
spec=""
res=-1

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

Saad Jbabdi's avatar
Saad Jbabdi committed
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
	*) echo "Unknown option '$1'";exit 1;;
    esac
    shift
done
# 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
121
echo SPECIES $spec
Saad Jbabdi's avatar
Saad Jbabdi committed
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146


# 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

Saad Jbabdi's avatar
Saad Jbabdi committed
147
if [ "$spec" == "" ];then
Saad Jbabdi's avatar
Saad Jbabdi committed
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
    if [ "$str" == "" ];then
	echo "Must set compulsory argument '-str'"
	errflag=1
    elif [ ! -f $str ];then
	echo "Structure file $str not found"
	errflag=1
    fi
    if [ "$p" == "" ];then
	echo "Must set compulsory argument '-p'"
	errflag=1
    elif [ ! -d $p ];then
	echo "Protocol folder $p not found"
	errflag=1
    fi
else
    if [ "$spec" == "HUMAN" ];then
Saad Jbabdi's avatar
Saad Jbabdi committed
164
165
	p=$datadir/Human
	str=$p/structureList
Saad Jbabdi's avatar
Saad Jbabdi committed
166
    elif [ "$spec" == "MACAQUE" ];then
Saad Jbabdi's avatar
Saad Jbabdi committed
167
168
	p=$datadir/Macaque
	str=$p/structureList
169
	Warning
Saad Jbabdi's avatar
Saad Jbabdi committed
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
    else
	echo "Species must be one of HUMAN or MACAQUE"
	errflag=1
    fi
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 "

191
192
193
# Add any user-defined ptx options
opts=" $opts $ptx_opts"

Saad Jbabdi's avatar
Saad Jbabdi committed
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
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
    struct=`echo $structstring | awk '{print $1}'`
    # skip empty lines and lines that start with '#'
    if [ "${struct:0:1}" == "#" ];then
	# do nothing
	foo=0
	#echo "----- Skip line $structstring -----"
    elif [ "$struct" == "" ];then
	# do nothing
	foo=0
	#echo "----- Skip empty line -----"
    else
	#echo "autoTrack $struct"
	mkdir -p $out/tracts/$struct

	nseed=`echo $structstring | awk '{print $2}'`
	if [ "$nseed" == "" ];then
	    nseed=1
	fi
	nseed=$(echo "scale=0; 1000 * ${nseed} / 1"|bc)
	
	maskdir=$p/$struct

	#  DEALING WITH RESAMPLING -- 
	# Pick space to run tractography in (diffusion or standard)	
	if [ "$nat" -eq 1 ];then
	    echo " -- transforming masks into native 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 $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"
	    done	    
	else
	    for m in seed stop exclude;do		
		if [ $res -gt 0 ];then
		    # Resample at a different resolution
		    mkdir -p $out/masks/$struct
		    if [ `$FSLDIR/bin/imtest $maskdir/$m` -eq 1 ];then
			$FSLDIR/bin/flirt -in $maskdir/$m -out $out/masks/$struct/$m -applyisoxfm $res -ref $maskdir/$m 
			$FSLDIR/bin/fslmaths $out/masks/$struct/$m -thr 0.1 -bin $out/masks/$struct/$m -odt char
		    fi		
		    eval "${m}=$out/masks/$struct/$m"		    
		else
		    eval "${m}=$maskdir/$m"
		fi
	    done
	fi	

	# 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
	
        # 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 "
	
	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"

	
	# Does the protocol define a second run with inverted seed / target masks? 
	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)
		addcmd="echo \"scale=5; \`cat $out/tracts/$struct/waytotal\` + \`cat $out/tracts/$struct/tractsInv/waytotal\` \"|bc > $out/tracts/$struct/sum_waytotal" 
		
		# 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
Saad Jbabdi's avatar
Saad Jbabdi committed
341
	fsl_sub -q long.q -l $out/logs -N xtract -t $commands	
Saad Jbabdi's avatar
Saad Jbabdi committed
342
    else
343
	fsl_sub -q $FSLGECUDAQ -T 300 -l $out/logs -N xtract $commands
Saad Jbabdi's avatar
Saad Jbabdi committed
344
345
346
347
348
349
    fi
else   # If no SGE, run locally
    sh $commands
fi

#EOF