Merge branch 'xtract_blueprint' into 'master'

Xtract blueprint

See merge request !19

Makefile

@@ -2,5 +2,4 @@ include ${FSLCONFDIR}/default.mk
 
 PROJNAME = xtract
 
-SCRIPTS = xtract xtract_viewer xtract_stats
-
+SCRIPTS = xtract xtract_viewer xtract_stats xtract_blueprint create_blueprint.py

README.md

@@ -262,7 +262,7 @@ data of interest, the directory containing the XTRACT output, the warp field (or
 if tracts are already in diffusion space). If tracts are not in diffusion space, you must also
 provide a reference image in diffusion space (e.g. FA map).
 
-e.g. call: xtract_stats -d /home/DTI/dti_ -xtract /home/xtract -w /home/warp/standard2diff -r /home/DTI/dti_FA
+e.g. call: `xtract_stats -d /home/DTI/dti_ -xtract /home/xtract -w /home/warp/standard2diff -r /home/DTI/dti_FA`
 
 The output (a .csv file) by default contains the tract volume (mm3) and the mean, median and
 standard deviation of the probability, length, FA and MD for each tract.
@@ -297,3 +297,155 @@ Usage:
        -keepfiles                             Keep temporary files
 
 ```
+
+---------------------------------------------------------------------
+
+## Generating Connectivity Blueprints with xtract_blueprint
+
+xtract_blueprint calculates the "connectivity blueprint" that represents how different grey matter areas are connected to a set of white matter bundles (Mars et al. 2018 eLife, Warrington et al. 2020 NeuroImage). Formally, this is a seeds x tracts matrix where each row describes how each (cortical) seed location is connected to major white matter fibre bundles, and each column describes the cortical terminations of each of the white matter fibre bundles. xtract_blueprint can, in principle, build a connectivity blueprint for any brain (e.g. different species), at any resolution (i.e. number of surface vertices) and for any region (whole whole or a given ROI, a lobe for example). It inherently supports surface input files (for instance a GIFTI representing the WM/GM boundary interface).
+
+The script was written by Shaun Warrington (University of Nottingham), Saad Jbabdi (Oxford University) and Stamatios Sotiropoulos (University of Nottingham).
+
+---------------------------------------------------------------------
+
+**Citations:**
+
+Mars R, Sotiropoulos SN, Passingham RE, Sallet J, Verhagen L, Khrapitchev AA, Sibson N, Jbabdi S (2018) Whole brain comparative anatomy using connectivity blueprints. eLife. DOI: 10.7554/eLife.35237
+
+Warrington S, Bryant K, Khrapitchev A, Sallet J, Charquero-Ballester M, Douaud G, Jbabdi S*, Mars R*, Sotiropoulos SN* (2020) XTRACT - Standardised protocols for automated tractography in the human and macaque brain. NeuroImage. DOI: 10.1016/j.neuroimage.2020.116923
+
+---------------------------------------------------------------------
+**Usage:**
+```
+
+__  _______ ____      _    ____ _____ _     _                       _       _
+\ \/ /_   _|  _ \    / \  / ___|_   _| |__ | |_   _  ___ _ __  _ __(_)_ __ | |_
+ \  /  | | | |_) |  / _ \| |     | | | '_ \| | | | |/ _ \ '_ \| '__| | '_ \| __|
+ /  \  | | |  _ <  / ___ \ |___  | | | |_) | | |_| |  __/ |_) | |  | | | | | |_
+/_/\_\ |_| |_| \_\/_/   \_\____| |_| |_.__/|_|\__,_|\___| .__/|_|  |_|_| |_|\__|
+                                                        |_|
+
+
+Usage:
+    xtract_blueprint -bpx <folder> -xtract <folder> -seeds <list> -warps <ref> <xtract2diff> <diff2xtract> -out <folder> [options]
+
+    Compulsory arguments:
+
+       -bpx       <folder>                          Path to bedpostx folder
+       -out       <folder>                          Path to output folder
+       -xtract    <folder>                          Path to xtract folder
+       -seeds     <list>                            Comma separated list of seeds for which a blueprint is requested (e.g. left and right cortex in standard space)
+       -warps     <ref> <xtract2diff> <diff2xtract> Standard space reference image and transforms between xtract space and diffusion space
+
+    Optional arguments:
+       -stage                                       What to run. 1:matrix2, 2:blueprint, all:everythng (default)
+       -gpu                                         Use GPU version
+       -savetxt                                     Save blueprint as txt file (nseed by ntracts) instead of CIFTI
+       -prefix    <string>                          Specify a prefix for the final blueprint filename (e.g. <prefix>_BP.LR.dscalar.nii)
+
+       -rois  <list>                                Comma separated list (1 per seed): ROIs (gifti) to restrict seeding (e.g. medial wall masks)
+       -stops     <stop.txt>                        Text file containing line separated list
+       -wtstops   <wtstop.txt>                      Text file containing line separated list
+       -tract_list                                  Comma separated list of tracts to include (default = all found under -xtract <folder>)
+
+       -thr                                         Threshold applied to XTRACT tracts prior to blueprint calculation (default = 0.001, i.e. 0.1% probability).
+       -nsamples                                    Number of samples per seed used in tractography (default = 1000)
+       -res       <mm>                              Resolution of matrix2 output (Default = 3 mm)
+       -ptx_options <options.txt>                   Pass extra probtrackx2 options as a text file to override defaults
+
+   Example (recommended) usage:
+      xtract_blueprint -bpx sub001/dMRI.bedpostx -out sub001/blueprint -xtract sub001/xtract -seeds sub001/l.white.surf.gii,sub001/r.white.surf.gii \
+           -rois sub001/l.medwall.shape.gii,sub001/r.medwall.shape.gii -warps MNI152_brain.nii.gz sub001/xtract2diff.nii.gz sub001/diff2xtract.nii.gz -gpu \
+
+
+```
+
+---------------------------------------------------------------------
+
+## Running xtract_blueprint:
+
+In order to use xtract_blueprint, you need to have run xtract first. To construct the connectivity blueprints, xtract_blueprint expects the same warp fields as used in the running of xtract.
+
+xtract_blueprint currently supports the construction of connectivity blueprints using surface (GIFTI) files only. As such, you must provide GIFTI surface files containing the white-grey matter boundary surface data. (We plan to extend this to support volume seeding in the future.)
+
+**Required input:**
+-	bedpostx folder - crossing fibre modelled diffusion data (expects to find nodif_brain_mask)
+-	xtract folder   - xtract tract folder (the output from xtract)
+-	seed 			      - the comma separated seed masks to use in tractography (e.g. the white-grey matter boundary surfaces), e.g. `L.white.surf.gii,R.white.surf.gii`
+-	warps			      - a reference standard space image and warps to and from the native diffusion and standard spaces, e.g. `MNI152.nii.gz standard2diff.nii.gz diff2standard.nii.gz`
+
+Note: if running whole-brain (recommended), you must provide the left seed first, as exampled.
+xtract_blueprint uses the seed GIFTI header information to build the CIFTI blueprint. Ensure that your structure is set in the seed GIFTI: either CortexLeft or CortexRight.
+
+**Running modes:**
+- Stage 1 only		- only run seed-based tractography
+- Stage 2 only		- only run blueprint processing (requires xtract output and tractography from stage 1)
+- All			        - runs both stage 1 and stage 2 processing
+
+xtract_blueprint is capable of GPU acceleration (`'-gpu'` flag) and detects $SGE_ROOT to work with fsl_sub. If using the CPU version, expect tractography to take many hours.
+
+
+**Tractography details and options:**
+
+Tractography is performed for each seed separately. The resultant connectivity matrices (fdt_matrix) are stacked in order to construct a whole-brain connectivity blueprint. You may also provide a single hemisphere if you wish.
+
+Optionally, you may also provide stop (stop tracking at locations given by this mask file) and wtstop (allow propagation within mask but terminate on exit) masks. Stop is typically the pial surface. wtstop is typically subcortical structures and/or the white surface. These should be specified as line separated text files. e.g. `-seeds <l.white.surf.gii,r.white.surf.gii> -stop stop.txt -wtstop wtstop.txt`
+
+Spatial resolution: by default tractography will be ran and stored using a resolution of 3 mm. This may be adjusted using the `'-res'` argument. Note: if required, xtract_blueprint will resample the xtract tracts. Warning: connectivity matrices are very large and require a lot of memory to handle - 3 mm is usually sufficient for the adult human brain.
+
+Additional probtrackx options may also be supplied. Simply add the probtrackx arguments to a text file and direct xtract_blueprint to this using the `'-ptx_options'` argument.
+
+Connectivity blueprints are primarily concerned with the connectivity of the cortex to white matter tracts. As such, we offer the option the mask out the medial wall. To do so, provide a single medial wall mask per supplied seed: e.g. `-seeds l.white.surf.gii,r.white.surf.gii -rois l.roi,r.roi`. By default, the medial wall is included in the calculation of the connectivity blueprint: we recommend the use of the medial wall mask to prevent this.
+
+The '-roi' argument may be used to restrict the blueprint to any region of interest, not just to exclude the medial wall. For example, you may provide an ROI restricting the blueprint to the temporal or frontal lobe.
+
+If you wish to use a stop/wtstop surface mask, you must ensure that the number of vertices matches the seed mask. This means that, if you are providing a seed mask and medial wall mask to xtract_blueprint, and wish to provide a surface stop mask, you must convert the stop mask to asc, restricting the data points to the medial wall mask, e.g.:
+
+	${FSLDIR}/bin/surf2surf -i stop.L.surf.gii -o stop.L.asc --outputtype=ASCII --values=l.roi.shape.gii
+	${FSLDIR}/bin/surf2surf -i stop.R.surf.gii -o stop.R.asc --outputtype=ASCII --values=r.roi.shape.gii
+
+Then supply "stop.L.asc" and "stop.R.asc" in a text file to xtract_blueprint using the `'-stop'` argument. This conversion is automatically performed for the seed mask in xtract_blueprint if a medial wall mask is supplied.
+
+
+**Which tracts are included?**
+
+Connectivity blueprints may be constructed using the provided XTRACT tracts or using your own. By default, xtract_blueprint will use all tracts it finds under the xtract folder. You can specify a subset, or you own tracts, by providing a comma separated list of tracts using the `'-tracts <str,str,str>'` argument.
+
+Certain tracts, e.g. the Middle Cerebellar Peduncle (MCP), do not project to the cortex. As such, they should be disregarded when interpreting the connectivity blueprint, or excluded from its construction.
+
+
+**Only interested in the connectivity of a specific area?**
+
+In many cases, the connectivity to a particular lobe, e.g. temporal or frontal, is of interest. You can use xtract_blueprint to obtain a connectivity blueprint for such a region:
+
+1. Define a binary mask which contains the region of interest as a shape.gii or func.gii file
+2. Select the tracts of interest: in all likelihood, only a subset of XTRACT's tracts will project to the ROI
+3. Supply the whole white matter surface file along with the ROI to xtract_blueprint, e.g. for the temporal lobe
+
+  xtract_blueprint -bpx sub001/dMRI.bedpostx -out sub001/blueprint -xtract sub001/xtract \
+  -warps MNI152_brain.nii.gz sub001/xtract2diff.nii.gz sub001/diff2xtract.nii.gz -gpu \
+  -seeds sub001/l.white.surf.gii -rois sub001/l.temporal_lobe.shape.gii -tract_list af_l,ilf_l,ifo_l,mdlf_l
+
+
+**Outputs of xtract_blueprint**
+
+xtract_blueprint will create an output directory specified by the `'-out'` argument. This will contain any log and command files along with a sub-directory per seed. Each sub-directory contains the resultant connectivity matrix from stage 1. The connectivity blueprint will be saved in the parent output directory (a CIFTI dscalar.nii file).
+
+Under outputDir:
+- Stage 1 (matrix2 tractography) output
+    - omatrix2:
+        - "ptx_commands.txt" - the probtrackx commands for tractography
+        - "`<seed>`" - sub-directory containing tractography results for each seed supplied, each containing:
+            - "coords_fdt_matrix2" - the coordinates of the seed mask
+            - "lookup_tractspace_fdt_matrix.nii.gz" - the target lookup space
+            - "tract_space_coords_for_fdt_matrix2" - the coordinates of the target mask
+            - "probtrackx.log" - the probtrackx log file
+            - "fdt_paths.nii.gz" - a NIFTI file containing the generated streamlines
+            - "fdt_matrix.dot" - the sparse format connectivity matrix (used to calculated the blueprint)
+            - "waytotal" - txt file continaing the number of valid streamlines
+- Stage 2 (blueprint calculation) output
+    - "bp_commands.txt" - the blueprint calculation commands
+    - "BP.`<L,R,LR>`.dscalar.nii" - CIFTI file containing the whole-brain connectivity blueprint - if running both hemispheres
+- "logs" - sub-directory containing job scheduler log files for both stages
+
+Alternatively, the `'-savetxt'` option may be used to override this. In this case, two txt files will be saved: the first (BP.<L,R,LR>.txt) will be an n_seed by n_tracts array containing the blueprint; the second (tract_order.`<L,R,LR>`.txt) is an n_tracts by 1 array containing the tract order in which the blueprint is structured. Note: no CIFTI file will be generated.

create_blueprint.py

+#!/usr/bin/env fslpython
+
+# Multiply a bunch of fdt_paths with a matrix2 to form a blueprint
+#
+# Author: Saad Jbabdi <saad@fmrib.ox.ac.uk>
+#
+# Copyright (C) 2020 University of Oxford
+# SHBASECOPYRIGHT
+
+import numpy as np
+import sys,os,glob
+import scipy.sparse as sps
+
+# Image stuff
+import nibabel as nib
+from nibabel import cifti2
+from fsl.data.image import Image
+from fsl.data.cifti import cifti2_axes
+from fsl.data.cifti import Cifti
+
+xtract_folder = sys.argv[1]
+ptx_folder    = sys.argv[2] # the ptx_folder(s)- if 2, comma separated
+seed_path      = sys.argv[3] # the seed mask(s) - if 2, comma separated
+roi_path  = sys.argv[4] # using the medial wall as a mask(s)? - if 2, comma separated
+tracts        = sys.argv[5]
+distnorm      = int(sys.argv[6])
+savetxt       = int(sys.argv[7]) # cii (0) or txt (1)?
+prefix        = sys.argv[8]
+
+# split the arguments here
+ptx_folder=ptx_folder.split(",")
+seed_path=seed_path.split(",")
+tracts=tracts.split(",")
+if prefix != "x":
+    prefix=f'{prefix}_'
+else:
+    prefix=''
+
+if len(ptx_folder) == 2:
+    print('Building whole-brain connectivity blueprint')
+else:
+    print('Building single hemisphere connectivity blueprint')
+
+# load the seed ROIs
+seed = []
+for p in seed_path:
+    temp = nib.load(p).darrays[0].data != 0
+    seed.append(temp[:,0])
+
+# if using medial wall, load the masks
+if roi_path != "x":
+    roi_path=roi_path.split(",")
+    roi = []
+    for p in roi_path:
+        temp = nib.load(p).darrays[0].data
+        if np.unique(temp).shape[0] > 2:
+            print('Warning!! Medial wall mask is not binary.')
+            print('Binarising...')
+            temp = (temp > 0)*1
+        roi.append(temp)
+    # and stack
+    if len(ptx_folder) == 2:
+        roi = np.concatenate((roi[0],roi[1]))
+    else:
+        roi = roi[0]
+
+print('')
+print('')
+
+# mask and lut are equal across hemispheres, so only need 1
+maskfile = os.path.join(ptx_folder[0],'lookup_tractspace_fdt_matrix2')
+mask     = Image(maskfile)
+lut      = Image(os.path.join(ptx_folder[0],'lookup_tractspace_fdt_matrix2.nii.gz'))
+lut      = lut.data[mask.data>0]-1
+
+
+# Get num voxels in mask and num tracts
+nvoxels = np.sum(mask.data>0)
+ntracts = len(tracts)
+
+print(f'Generating blueprint with {ntracts} tracts')
+# Collect tracts
+tracts_mat = np.zeros( (nvoxels, ntracts) )
+print('Reading tracts...')
+trm=[]
+for idx,t in enumerate(tracts):
+    t = os.path.join(xtract_folder,f'{t}.nii.gz')
+    if os.path.exists(t):
+        tract = Image(t)
+        tracts_mat[lut,idx] = tract[mask.data>0]
+    else:
+        print(f'Could not find {tracts[idx]}. Skipping... (will remove from final output)')
+        trm.append(int(idx))
+
+# remove missing tract columns
+tracts_mat = np.delete(tracts_mat, [trm], 1)
+for i in trm:
+    tracts.pop(i)
+
+def normalise(M):
+    D = np.sum(M,axis=1)
+    D[D==0] = 1
+    M = M / D[:,None]
+    return M
+
+def load_fdt(fdt_path):
+    print('Reading tractography matrix. This may take a few minutes...')
+    mat         = np.loadtxt(fdt_path)
+    data        = mat[:-1, -1]
+    rows        = np.array(mat[:-1, 0]-1, dtype=int)
+    cols        = np.array(mat[:-1, 1]-1, dtype=int)
+    nrows,ncols = int(mat[-1, 0]), int(mat[-1,1])
+    M           = sps.csc_matrix((data, (rows,cols)), shape=(nrows,ncols)).toarray()
+    M           = normalise(M)
+    return M
+
+# Open matrix2 file
+M = []
+for p in ptx_folder:
+    M.append(load_fdt(os.path.join(p, "fdt_matrix2.dot")))
+
+if len(ptx_folder) == 2:
+    print('Stacking hemispheres...')
+    M = np.concatenate((M[0], M[1]))
+else:
+    M = M[0]
+
+print(f'Tractography matrix dimensions: {M.shape[0]} vertices, {M.shape[1]} voxels')
+
+# Create blueprint
+print('Calculating blueprint...')
+BP = M@tracts_mat
+BP = normalise(BP)
+
+# Convert to full cortex structure here (i.e. add in empty medial wall)
+if roi_path != "x":
+    # and stack
+    if len(ptx_folder) == 2:
+        seed_temp = np.concatenate((seed[0],seed[1]))
+    else:
+        seed_temp = seed[0]
+    full_BP = np.zeros([np.shape(seed_temp)[0], np.shape(BP)[1]])
+    full_BP[roi == 1, :] = BP
+    BP = full_BP
+
+# function to find which hemisphere/cifti structure we're working with
+def get_structure(p):
+    f = open(p, 'r')
+    text = f.read()
+    text = text.split('Cortex',1)[1]
+    cstruct = text.split(']]></Value>',1)[0]
+    return cstruct
+
+out_folder = os.path.dirname(os.path.dirname(ptx_folder[0]))
+if savetxt == 1:
+    if len(ptx_folder) == 1:
+        cstruct = get_structure(seed_path[0])
+        if cstruct == 'Left':
+            side='L'
+        else:
+            side='R'
+
+        new_fname = os.path.join(out_folder, f'{prefix}BP.{side}.txt')
+        np.savetxt(new_fname, BP)
+        new_fname_tord = os.path.join(out_folder, f'{prefix}tract_order.{side}.txt')
+        np.savetxt(new_fname_tord, tracts, fmt="%s")
+    elif len(ptx_folder) == 2:
+        new_fname = os.path.join(out_folder, f'{prefix}BP.LR.txt')
+        np.savetxt(new_fname, BP)
+        new_fname_tord = os.path.join(out_folder, f'{prefix}tract_order.LR.txt')
+        np.savetxt(new_fname_tord, tracts, fmt="%s")
+else:
+    # set up cifti brain model axes
+    if len(ptx_folder) == 1:
+        cstruct = get_structure(seed_path[0])
+        if cstruct == 'Left':
+            side='L'
+        else:
+            side='R'
+        bm        = cifti2_axes.BrainModelAxis.from_mask(seed[0], name=f'Cortex{cstruct}')
+        new_fname = os.path.join(out_folder, f'{prefix}BP.{side}.dscalar.nii')
+    elif len(ptx_folder) == 2:
+        bm_l      = cifti2_axes.BrainModelAxis.from_mask(seed[0], name=f'CortexLeft')
+        bm_r      = cifti2_axes.BrainModelAxis.from_mask(seed[1], name=f'CortexRight')
+        bm        = bm_l + bm_r
+        new_fname = os.path.join(out_folder, f'{prefix}BP.LR.dscalar.nii')
+    # save cifti
+    sc        = cifti2_axes.ScalarAxis(tracts)
+    hdr       = cifti2.Cifti2Header.from_axes((sc, bm))
+    img       = cifti2.Cifti2Image(BP.T, hdr)
+    nib.save(img, new_fname)
+
+print(f'Saved: {new_fname}')
+print('')

xtract_blueprint

+#!/bin/bash
+
+#   Copyright (C) 2020 University of Oxford
+#
+#   SHCOPYRIGHT
+
+# Written by Saad Jbabdi, Stam Sotiropoulos & Shaun Warrington
+
+shopt -s extglob
+
+code_folder="$FSLDIR/bin"
+
+# Location of probtrackx2_gpu binary
+ptxbin_gpu=$FSLDIR/bin/probtrackx2_gpu
+
+Splash (){
+cat <<EOF
+
+__  _______ ____      _    ____ _____ _     _                       _       _
+\ \/ /_   _|  _ \    / \  / ___|_   _| |__ | |_   _  ___ _ __  _ __(_)_ __ | |_
+ \  /  | | | |_) |  / _ \| |     | | | '_ \| | | | |/ _ \ '_ \| '__| | '_ \| __|
+ /  \  | | |  _ <  / ___ \ |___  | | | |_) | | |_| |  __/ |_) | |  | | | | | |_
+/_/\_\ |_| |_| \_\/_/   \_\____| |_| |_.__/|_|\__,_|\___| .__/|_|  |_|_| |_|\__|
+                                                        |_|
+
+EOF
+}
+Splash
+
+
+Usage() {
+    cat << EOF
+
+Usage:
+    xtract_blueprint -bpx <folder> -out <folder> -xtract <folder> -seeds <list> [options]
+
+    Compulsory arguments:
+
+       -bpx       <folder>                          Path to bedpostx folder
+       -out       <folder>                          Path to output folder
+       -xtract    <folder>                          Path to xtract folder
+       -seeds     <list>                            Comma separated list of seeds for which a blueprint is requested (e.g. left and right cortex in standard space)
+       -warps     <ref> <xtract2diff> <diff2xtract> Standard space reference image, and transforms between xtract space and diffusion space
+
+    Optional arguments:
+       -stage                                       What to run. 1:matrix2, 2:blueprint, all:everythng (default)
+       -gpu                                         Use GPU version
+       -savetxt                                     Save blueprint as txt file (nseed by ntracts) instead of CIFTI
+       -prefix    <string>                          Specify a prefix for the final blueprint filename (e.g. <prefix>_BP.LR.dscalar.nii)
+
+       -rois      <list>                            Comma separated list (1 per seed): ROIs (gifti) to restrict seeding (e.g. medial wall masks)
+       -stops     <stop.txt>                        Text file containing line separated list
+       -wtstops   <wtstop.txt>                      Text file containing line separated list
+       -tract_list                                  Comma separated list of tracts to include (default = all found under -xtract <folder>)
+
+       -thr                                         Threshold applied to XTRACT tracts prior to blueprint calculation (default = 0.001, i.e. 0.1% probability).
+       -nsamples                                    Number of samples per seed used in tractography (default = 1000)
+       -res       <mm>                              Resolution of matrix2 output (Default = 3 mm)
+       -ptx_options <options.txt>                   Pass extra probtrackx2 options as a text file to override defaults
+
+   Example (recommended) usage:
+      xtract_blueprint -bpx <bpxdir> -out <outdir> -xtract <xtractdir> -seeds <l.white.surf.gii,r.white.surf.gii> \
+           -rois <l.medwall.shape.gii,r.medwall.shape.gii> -warps <ref> <xtract2diff> <diff2xtract> -gpu \
+
+EOF
+    exit 1
+}
+#       -distnorm                                    Normalise tracts by distance
+
+[ "$1" = "" ] && Usage
+
+# Set default options
+bpx=""
+out=""
+xtract=""
+seeds=""
+stops=""
+wtstops=""
+rois=""
+distnorm=0
+ptx_opts=""
+stdref=""
+gpu=0
+ref=""
+diff2xtract=""
+xtract2diff=""
+diff2seed=""
+seed2diff=""
+spec=""
+res=3
+stage=all
+savetxt=0
+thr=0.001
+nsamples=1000
+tract_list=""
+prefix="x" # if prefix not specified, pass 'x' to create_blueprint which is nulled
+
+# Parse command-line arguments
+while [ ! -z "$1" ];do
+    case "$1" in
+	-bpx) bpx=$2;shift;;
+	-out) out=$2;shift;;
+	-seeds) seeds=(`echo $2 | sed s@","@" "@g`);shift;;
+  -stops) stops=(`echo $2 | sed s@","@" "@g`);shift;;
+  -wtstops) wtstops=(`echo $2 | sed s@","@" "@g`);shift;;
+  -rois) rois=(`echo $2 | sed s@","@" "@g`);shift;;
+	-xtract) xtract=$2;shift;;
+	-warps) ref=$2;xtract2diff=$3;diff2xtract=$4;shift;shift;shift;;
+	-gpu) gpu=1;;
+	-res) res=$2;shift;;
+	-stage) stage=$2;shift;;
+	-ptx_options) ptx_options=`cat $2`;shift;;
+  -savetxt) savetxt=1;;
+  -thr) thr=$2;shift;;
+  -tract_list) tract_list=$2;shift;;
+  -prefix) prefix=$2;shift;;
+  -nsamples) nsamples=$2;shift;;
+	*) echo "Unknown option '$1'";exit 1;;
+    esac
+    shift
+done
+
+#	-distnorm) distnorm=1;;
+
+# Step 0 : Prepare
+# 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 [ "$xtract" == "" ];then
+    echo "Must set compulsory argument '-xtract'"
+    errflag=1
+elif [ ! -d $xtract ];then
+    echo "XTRACT folder $xtract not found"
+    errflag=1
+fi
+if [ "$out" == "" ];then
+    echo "Must set compulsory argument '-out'"
+    errflag=1
+fi
+if [ "$seeds" == "" ];then
+    echo "Must set compulsory argument '-seeds'"
+    errflag=1
+fi
+
+# Check seeds and stops (check they exist and 1 stop per seed if using)
+for seed in ${seeds[@]};do
+  if [ ! -f $seed ];then
+    echo "Seed file $seed not found"
+    errflag=1
+  fi
+done
+
+nseeds=${#seeds[@]}
+nrois=${#rois[@]}
+if [ ! "$stops" == "" ]; then
+    if [ ! -f $stop ];then
+      echo "Stop file $stop not found"
+      errflag=1
+    fi
+fi
+if [ ! "$wtstops" == "" ]; then
+    if [ ! -f $stop ];then
+      echo "wtstop file $wtstop not found"
+      errflag=1
+    fi
+fi
+if [ ! "$rois" == "" ]; then
+    if [ ! $nrois -eq $nseeds ]; then
+      echo "Supplied $nseeds seeds but $nrois ROIs"
+      echo "If using ROIs, must supply 1 per seed"
+      errflag=1
+    fi
+    for roi in ${rois[@]};do
+      if [ ! -f $roi ];then
+        echo "ROI file $roi not found"
+        errflag=1
+      fi
+    done
+fi
+
+# check ref and warps
+if [ ! -f $ref ];then
+  echo "Reference file $ref not found"
+  errflag=1
+fi
+if [ ! -f $xtract2diff ];then
+  echo "Warp file $xtract2diff not found"
+  errflag=1
+fi
+if [ ! -f $diff2xtract ];then
+  echo "Warp file $diff2xtract not found"
+  errflag=1
+fi
+
+if [ "$errflag" -eq 1 ];then
+    echo ""
+    echo "Exit without doing anything..."
+    exit 1
+fi
+
+# GPU stuff
+if [ "$gpu" == 0 ];then
+    ptxbin=$FSLDIR/bin/probtrackx2
+    echo "Warning: not using GPU mode - this may take a while. Consider downsampling the seed mask."
+else
+    ptxbin=${ptxbin_gpu}
+fi
+
+# what are we running?
+if [ $stage == 1 ];then
+    echo "Running tractography only"
+elif [ $stage == 2 ];then
+    echo "Running blueprint processing only"
+elif [ $stage == "all" ];then
+    echo "Running tractography and blueprint processing"
+else
+    echo "Unknown stage option. Either: '1' for matrix2 tractography, '2' for blueprint or 'all' for both"
+fi
+
+if [ -d $out ];then
+    echo "Warning: Output folder already exists. Some of the files may be overwritten"
+fi
+mkdir -p $out
+mkdir -p $out/omatrix2
+
+# Get tract list here
+# if not provided by user, get list from under xtract folder
+function join_str { local IFS="$1"; shift; echo "$*"; }
+if [ "$tract_list" == "" ]; then
+    # if not supplied, get tracts from xtract folder
+    struct=(`ls ${xtract}/tracts/*/densityNorm.nii.gz`)
+    unset tracts
+    for t in ${struct[@]}; do
+        t=$( echo ${t%/*} )
+        t=$( echo ${t##*/} )
+        tracts+=($t)
+    done
+    # also build comma separated tract_list for python usage later
+    tract_list=`join_str , ${tracts[@]}`
+else
+    tracts=(`echo $tract_list | sed s@","@" "@g`)
+fi
+
+# function for image resmapling
+resample(){
+    # resample <src_img> <mm_resolution> <out_img> <interp>
+    i=$1
+    r=$2
+    o=$3