Merge pull request #5 from SPMIC-UoN/xpy

Bug fixes and add xtract_divergence

Makefile

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

README.md

@@ -7,27 +7,22 @@ define a set of masks in standard space (e.g. MNI152).
 The script was written by Shaun Warrington, Saad Jbabdi & Stamatios Sotiropoulos
 (based on the autoPtx tool by Marius de Groot - see https://fsl.fmrib.ox.ac.uk/fsl/fslwiki/AutoPtx)
 
-The HUMAN and MACAQUE tractography protocols were created by:
+The HUMAN and MACAQUE tractography protocols were created by: Rogier Mars, Stamatios Sotiropoulos, Saad Jbabdi, Kathryn Bryant, Shaun Warrington, Marina Charquero-Ballester, Gwenaelle Douaud
 
-Rogier Mars, Stamatios Sotiropoulos, Saad Jbabdi, Kathryn Bryant, Shaun Warrington, Marina Charquero-Ballester, Gwenaelle Douaud
-
-The BABY tractography protocols were created by:
-
-Elinor Thompson, Shaun Warrington, Matteo Bastiani, Jessica Dubois, Stamatios Sotiropoulos
-
-The MACAQUE tractography protocols were generalised across templates by:
-
-Stephania Assimopoulos, Shaun Warrington, Stamatios Sotiropoulos
+The BABY tractography protocols were created by: Elinor Thompson, Shaun Warrington, Matteo Bastiani, Jessica Dubois, Stamatios Sotiropoulos
 
+The MACAQUE tractography protocols were generalised across templates by: Stephania Assimopoulos, Shaun Warrington, Stamatios Sotiropoulos
 
 The XTRACT viewer helper script was written by Shaun Warrington
 
 The XTRACT stats helper script was written by Shaun Warrington
 
-The XTRACT blueprint helper script was written by Shaun Warrington and Saad Jbabdi
+The XTRACT blueprint script was written by Shaun Warrington and Saad Jbabdi
 
 The XTRACT QC helper script was written by Shaun Warrington
 
+The XTRACT divergence script was written by Shaun Warrington and Stephania Assimopoulos
+
 ---------------------------------------------------------------------
 
 ## Citations:
@@ -42,8 +37,8 @@ NeuroImage, 76(1), 400-411. DOI: 10.1016/j.neuroimage.2013.03.015
 
 If you use the neonatal protocols, please cite:
 
-Warrington, S.*, Thompson, E.*, Bastiani, M., Dubois, J., Baxter, L., Slater, R., Jbabdi, S., Mars, R.B., 
-and Sotiropoulos, S.N. (2022) “Concurrent mapping of brain ontogeny and phylogeny within a common space: 
+Warrington, S.*, Thompson, E.*, Bastiani, M., Dubois, J., Baxter, L., Slater, R., Jbabdi, S., Mars, R.B.,
+and Sotiropoulos, S.N. (2022) “Concurrent mapping of brain ontogeny and phylogeny within a common space:
 Standardized tractography and applications”, Science Advances, vol. 8(42), doi.org/10.1101/2022.03.03.482776
 
 If you use the macaque protocols, please cite:
@@ -54,4 +49,3 @@ Assimopoulos, S., Warrington, S., Bryant, K.L., Pszczolkowski, S., Jbabdi S., Ma
 ---------------------------------------------------------------------
 
 For full details and usage, see the FSL wiki page: https://fsl.fmrib.ox.ac.uk/fsl/docs/#/diffusion/xtract
-

xtract

@@ -23,7 +23,7 @@ def errchk(errflag):
 def imgtest(fname):
     r = subprocess.run([f'{os.path.join(FSLbin, "imtest")} {fname}'], capture_output=True, text=True, shell=True)
     return int(r.stdout)
-    
+
 def applywarp(fin, out, ref, warp, dtype, interp):
     if interp == 'tri':
         r = subprocess.run([os.path.join(FSLbin, "applywarp"), '-i', fin, '-o', out, '-r', ref, '-w', warp, '-d', dtype, '--interp=trilinear'])
@@ -103,36 +103,39 @@ parser = MyParser(prog='XTRACT',
                   epilog=textwrap.dedent('''Example usage:
                         xtract -bpx /data/Diffusion.bedpostX -out /data/xtract -species HUMAN -stdwarp std2diff.nii.gz diff2std.nii.gz -gpu
                   '''))
-
+
+required = parser.add_argument_group('Required arguments')
+optional = parser.add_argument_group('Optional arguments')
+
 # Compulsory arguments:
-parser.add_argument("-bpx", metavar='<folder>', help="Path to bedpostx folder")
-parser.add_argument("-out", metavar='<folder>', help="Path to output folder")
-parser.add_argument("-species", metavar='<SPECIES>', choices=['HUMAN', 'MACAQUE', 'HUMAN_BABY', 'CUSTOM'], help="One of HUMAN or MACAQUE or HUMAN_BABY or CUSTOM.")
-# parser.add_argument("-species", metavar='<SPECIES>', choices=['HUMAN', 'MACAQUE', 'MACAQUE_F99', 'MACAQUE_D99', 'MACAQUE_INIA', 'MACAQUE_NMT', 'MACAQUE_YRK', 'HUMAN_BABY', 'CUSTOM'], help="One of HUMAN or MACAQUE or HUMAN_BABY or CUSTOM. Can also specify the macaque template space to be used by appending _[D99,INIA,NMT,YRK] (default is F99).")
+required.add_argument("-bpx", metavar='<folder>', help="Path to bedpostx folder")
+required.add_argument("-out", metavar='<folder>', help="Path to output folder")
+required.add_argument("-species", metavar='<SPECIES>', choices=['HUMAN', 'MACAQUE', 'HUMAN_BABY', 'CUSTOM'], help="One of HUMAN or MACAQUE or HUMAN_BABY or CUSTOM.")
 
 # Optional protocol arguments:
-parser.add_argument("-tract_list", metavar='<list>', dest='tract_list', help="Comma separated tract list defining the subset of tracts to run (must follow xtract naming convention)")
-parser.add_argument("-str", metavar='<file>', dest='structureList', help="Structures file (format: <tractName> [samples=1], 1 means 1000, '#' to skip lines)")
-parser.add_argument("-p", metavar='<folder>', help=f"Protocols folder (all masks in same standard space) (Default={os.path.join(FSLDIR, 'data', 'xtract_data', '<SPECIES>')})")
-parser.add_argument("-stdref", metavar='<reference>', help=f"Standard space reference image (Default = {os.path.join(FSLDIR, 'data', 'standard', 'MNI152_T1_1mm [HUMAN]')}, {os.path.join(datadir, 'standard', 'F99', 'mri', 'struct_brain [MACAQUE]')}, {os.path.join(datadir, 'standard', 'neonate', 'mri', 'schuh_template [HUMAN_BABY]')}")
-parser.add_argument("-stdwarp", metavar='<path>', nargs=2, help=f"<std2diff> <diff2std> Non-linear Standard2diff and Diff2standard transforms (Default={os.path.join('bedpostx_dir', 'xfms', '{standard2diff_warp.nii.gz,diff2standard_warp.nii.gz}')})")
+optional.add_argument("-tract_list", metavar='<list>', dest='tract_list', help="Comma separated tract list defining the subset of tracts to run (must follow xtract naming convention)")
+optional.add_argument("-str", metavar='<file>', dest='structureList', help="Structures file (format: <tractName> [samples=1], 1 means 1000, '#' to skip lines)")
+optional.add_argument("-p", metavar='<folder>', help=f"Protocols folder (all masks in same standard space) (Default={os.path.join(FSLDIR, 'data', 'xtract_data', '<SPECIES>')})")
+optional.add_argument("-stdref", metavar='<reference>', help=f"Standard space reference image (Default = {os.path.join(FSLDIR, 'data', 'standard', 'MNI152_T1_1mm [HUMAN]')}, {os.path.join(datadir, 'standard', 'F99', 'mri', 'struct_brain [MACAQUE]')}, {os.path.join(datadir, 'standard', 'neonate', 'mri', 'schuh_template [HUMAN_BABY]')}")
+optional.add_argument("-stdwarp", metavar='<path>', nargs=2, help=f"<std2diff> <diff2std> Non-linear Standard2diff and Diff2standard transforms (Default={os.path.join('bedpostx_dir', 'xfms', '{standard2diff_warp.nii.gz,diff2standard_warp.nii.gz}')})")
 
 # Optional space arguments:
-parser.add_argument("-native", action="store_true", default=False, help="Run tractography in native (diffusion) space")
-parser.add_argument("-ref", metavar='<path>', default=False, nargs=3, help="<refimage> <diff2ref> <ref2diff> Reference image for running tractography in reference space, Diff2Reference and Reference2Diff transforms")
+optional.add_argument("-native", action="store_true", default=False, help="Run tractography in native (diffusion) space")
+optional.add_argument("-ref", metavar='<path>', default=False, nargs=3, help="<refimage> <diff2ref> <ref2diff> Reference image for running tractography in reference space, Diff2Reference and Reference2Diff transforms")
 
 # Optional tracking/computing arguments:
-parser.add_argument("-res", metavar='<mm>', default=0, help="Output resolution (Default=same as in protocol folders unless '-native' used)")
-parser.add_argument("-ptx_options", metavar='<file>', help="Pass extra probtrackx2 options as a text file to override defaults, e.g. --steplength=0.2 --distthresh=10)")
-parser.add_argument("-interp", metavar='<str>', choices=['nn', 'tri'], default='nn', help="If native/ref: default interpolation of protocol ROIs is nearest neighbour ('nn', since vXXX). For backwards compatability, trilinear plus thresholding ('tri') is available")
-parser.add_argument("-gpu", action="store_true", default=False, help="Use GPU version")
-parser.add_argument("-par", action="store_true", default=False, help="If cluster, run in parallel: submit 1 job per tract")
-parser.add_argument("-print_list", dest='print_list', action="store_true", default=False, help="List the tract names used in XTRACT")
-parser.add_argument("-version", dest='version', action="store_true", default=False, help="List the package versions")
+optional.add_argument("-res", metavar='<mm>', default=0, help="Output resolution (Default=same as in protocol folders unless '-native' used)")
+optional.add_argument("-ptx_options", metavar='<file>', help="Pass extra probtrackx2 options as a text file to override defaults, e.g. --steplength=0.2 --distthresh=10)")
+optional.add_argument("-interp", metavar='<str>', choices=['nn', 'tri'], default='nn', help="If native/ref: default interpolation of protocol ROIs is nearest neighbour ('nn', since vXXX). For backwards compatability, trilinear plus thresholding ('tri') is available")
+optional.add_argument("-gpu", action="store_true", default=False, help="Use GPU version")
+optional.add_argument("-par", action="store_true", default=False, help="If cluster, run in parallel: submit 1 job per tract")
+optional.add_argument("-print_list", dest='print_list', action="store_true", default=False, help="List the tract names used in XTRACT")
+optional.add_argument("-version", dest='version', action="store_true", default=False, help="List the package versions")
 
 if FSLINFMRIB is not None:
-    parser.add_argument("-queue", metavar='<str>', default=None, help="For job schedulers, specify the queue fsl_sub submits to")
-    parser.add_argument("-mac_ext", metavar='<str>', choices=['_F99', '_D99', '_INIA', '_NMT', '_YRK'], default=None, help="Macaque template extension")
+    hidden = parser.add_argument_group('Hidden arguments')
+    hidden.add_argument("-queue", metavar='<str>', default=None, help="For job schedulers, specify the queue fsl_sub submits to")
+    hidden.add_argument("-mac_ext", metavar='<str>', choices=['_F99', '_D99', '_INIA', '_NMT', '_YRK'], default=None, help="Macaque template extension")
 
 argsa = parser.parse_args()
 
@@ -177,6 +180,7 @@ if FSLINFMRIB is not None and argsa.mac_ext is not None:
 if not os.path.isdir(bpx):
     print(f'Bedpostx folder {bpx} not found')
     errflag = 1
+
 if os.path.isdir(out):
     print('Warning: output directory already exists. XTRACT may overwrite existing content.')
 
@@ -226,9 +230,11 @@ errchk(errflag)
 if imgtest(stdref) == 0:
     print(f"Standard space reference image '-stdref' {stdref} not found")
     errflag = 1
+
 if not os.path.isdir(p):
     print(f"Protocol folder {p} not found")
     errflag = 1
+
 if not os.path.isfile(structureList):
     print(f"Structures files {structureList} not found")
     errflag = 1
@@ -250,6 +256,7 @@ else:
 if imgtest(std2diff) == 0:
     print(f'std2diff {std2diff} not found.')
     errflag = 1
+
 if imgtest(diff2std) == 0:
     print(f'diff2std {diff2std} not found.')
     errflag = 1
@@ -299,11 +306,17 @@ errchk(errflag)
 # Set common ptx options
 nodif_bm = os.path.join(bpx, 'nodif_brain_mask')
 opts = f" -s {os.path.join(bpx, 'merged')} -m {nodif_bm} -V 1"
-opts += f" --loopcheck --forcedir --opd --ompl --sampvox=1 --randfib=1 "
+opts += " --loopcheck --forcedir --opd --ompl"
 
 # Add any user-defined ptx options
 opts += f" {ptx_opts}"
 
+if 'sampvox' not in ptx_opts: 
+    opts += " --sampvox=1"
+
+if 'randfib' not in ptx_opts: 
+    opts += " --randfib=1"
+
 # space options
 if not argsa.native and argsa.ref is False:
     opts += f" --seedref={stdref} --xfm={std2diff} --invxfm={diff2std} "
@@ -348,9 +361,7 @@ print('Preparing submission script...')
 for struct, nseed in structure_arr:    
     structdir = os.path.join(out, 'tracts', struct)
     os.makedirs(structdir, exist_ok=True)
-
     nseed = int(float(nseed)*1000)
-
     maskdir = os.path.join(p, struct)
     maskdict = {'seed': '', 'target': '', 'exclude': '', 'stop': ''} # mask paths stored here
     #  DEALING WITH RESAMPLING --
@@ -410,11 +421,12 @@ for struct, nseed in structure_arr:
             with open(targetfile, 'w') as f:
                 print(*targets, file=f)
             maskdict['target'] = targetfile
-                
+
     # Get generic options
     struct_opts = opts
     if imgtest(maskdict['stop']):
         struct_opts += f' --stop={maskdict["stop"]}'
+
     if imgtest(maskdict['exclude']):
         struct_opts += f' --avoid={maskdict["exclude"]}'
 
@@ -436,7 +448,11 @@ for struct, nseed in structure_arr:
         os.makedirs(structdir_inv, exist_ok=True)
 
         # invert targets/seed
-        targets_inv = targets.copy()
+        if argsa.native:
+            targets_inv = open(targetfile, "r").read().rstrip().split('\n')
+        else:
+            targets_inv = targets.copy()
+
         targets_inv.reverse()
         targets_inv.append(maskdict["seed"])
         seed_inv = targets_inv.pop(0)