Python fsl.findimage() function

bin/dwibiascorrect

@@ -51,10 +51,6 @@ if app.args.fsl:
       app.error('Could not find FSL program fast; please verify FSL install')
 
   fsl_suffix = fsl.suffix()
-  if fast_cmd == 'fast':
-    fast_suffix = fsl_suffix
-  else:
-    fast_suffix = '.nii.gz'
 
 elif app.args.ants:
 
@@ -101,10 +97,12 @@ else:
 run.command('dwiextract in.mif - -bzero | mrmath - mean mean_bzero.mif -axis 3')
 
 if app.args.fsl:
+
   # FAST doesn't accept a mask input; therefore need to explicitly mask the input image
   run.command('mrcalc mean_bzero.mif mask.mif -mult - | mrconvert - mean_bzero_masked.nii -stride -1,+2,+3')
   run.command(fast_cmd + ' -t 2 -o fast -n 3 -b mean_bzero_masked.nii')
-  bias_path = 'fast_bias' + fast_suffix
+  bias_path = fsl.findImage('fast_bias')
+
 elif app.args.ants:
 
   # If the mask image was provided manually, and doesn't match the input image perfectly

bin/dwipreproc

@@ -473,14 +473,19 @@ if do_topup:
       run.command('dwiextract dwi.mif' + import_dwi_pe_table_option + ' -pe ' + ','.join(line.split()) + ' - | mrconvert - ' + input_path + ' -json_export ' + json_path)
       run.command(applytopup_cmd + ' --imain=' + input_path + ' --datain=applytopup_config.txt --inindex=' + str(index) + ' --topup=field --out=' + temp_path + ' --method=jac')
       file.delTempFile(input_path)
+      temp_path = fsl.findImage(temp_path)
       run.command('mrconvert ' + temp_path + ' ' + output_path + ' -json_import ' + json_path)
       file.delTempFile(json_path)
       file.delTempFile(temp_path)
       applytopup_image_list.append(output_path)
       index += 1
 
   # Use the initial corrected volumes to derive a brain mask for eddy
-  run.command('mrcat ' + ' '.join(applytopup_image_list) + ' - | dwi2mask - - | maskfilter - dilate - | mrconvert - mask.nii -datatype float32 -stride -1,+2,+3')
+  if len(applytopup_image_list) == 1:
+    run.command('dwi2mask ' + applytopup_image_list[0] + ' - | maskfilter - dilate - | mrconvert - mask.nii -datatype float32 -stride -1,+2,+3')
+  else:
+    run.command('mrcat ' + ' '.join(applytopup_image_list) + ' - | dwi2mask - - | maskfilter - dilate - | mrconvert - mask.nii -datatype float32 -stride -1,+2,+3')
+
   for entry in applytopup_image_list:
     file.delTempFile(entry)
 
@@ -539,7 +544,7 @@ for index1 in range(num_volumes):
 if not len(volume_pairs) == int(num_volumes/2):
 
   # Convert the resulting volume to the output image, and re-insert the diffusion encoding
-  run.command('mrconvert dwi_post_eddy' + fsl_suffix + ' result.mif' + stride_option + ' -fslgrad ' + bvecs_path + ' bvals')
+  run.command('mrconvert ' + fsl.findImage('dwi_post_eddy') + ' result.mif' + stride_option + ' -fslgrad ' + bvecs_path + ' bvals')
 
 else:
   app.console('Detected matching DWI volumes with opposing phase encoding; performing explicit volume recombination')
@@ -569,7 +574,7 @@ else:
     # If one diffusion sensitisation gradient direction is reversed with respect to
     #   the other, still want to enable their recombination; but need to explicitly
     #   account for this when averaging the two directions
-    if norm2 < 0.999:
+    if norm2 < 0.0:
       bvec_sum = [ float(bvecs[0][pair[0]]) - float(bvecs[0][pair[1]]),
                    float(bvecs[1][pair[0]]) - float(bvecs[1][pair[1]]),
                    float(bvecs[2][pair[0]]) - float(bvecs[2][pair[1]]) ]
@@ -600,12 +605,12 @@ else:
   # Detect this, and manually replace the transform if necessary
   #   (even if this doesn't cause an issue with the subsequent mrcalc command, it may in the future, it's better for
   #   visualising the script temporary files, and it gives the user a warning about an out-of-date FSL)
-  field_map_image = 'field_map' + fsl_suffix
+  field_map_image = fsl.findImage('field_map')
   if not image.match('topup_in.nii', field_map_image):
     app.warn('topup output field image has erroneous header; recommend updating FSL to version 5.0.8 or later')
-    run.command('mrtransform ' + field_map_image + ' -replace topup_in.nii field_map_fix' + fsl_suffix)
+    run.command('mrtransform ' + field_map_image + ' -replace topup_in.nii field_map_fix.mif')
     file.delTempFile(field_map_image)
-    field_map_image = 'field_map_fix' + fsl_suffix
+    field_map_image = 'field_map_fix.mif'
 
 
 
@@ -638,8 +643,8 @@ else:
   # If eddy provides its main image output in a compressed format, the code block below will need to
   #   uncompress that image independently for every volume pair. Instead, if this is the case, let's
   #   convert it to an uncompressed format before we do anything with it.
-  eddy_output = 'dwi_post_eddy' + fsl_suffix
-  if fsl_suffix.endswith('.gz'):
+  eddy_output = fsl.findImage('dwi_post_eddy')
+  if eddy_output.endswith('.gz'):
     run.command('mrconvert ' + eddy_output + ' dwi_post_eddy.nii')
     file.delTempFile(eddy_output)
     eddy_output = 'dwi_post_eddy.nii'

bin/labelsgmfix

@@ -83,7 +83,7 @@ if app.args.premasked:
   first_input_is_brain_extracted = ' -b'
 run.command(first_cmd + ' -s ' + ','.join(structure_map.keys()) + ' -i T1.nii' + first_input_is_brain_extracted + ' -o first')
 
-# Generate an empty image that will be used to contruct the new SGM nodes
+# Generate an empty image that will be used to construct the new SGM nodes
 run.command('mrcalc parc.mif 0 -min sgm.mif')
 
 # Read the local connectome LUT file

lib/mrtrix3/_5ttgen/fsl.py

@@ -135,14 +135,11 @@ def execute():
       run.command(ssroi_cmd + ' T1.nii T1_preBET' + fsl_suffix + ' -maskMASK mni_mask.nii' + ssroi_roi_option, False)
     else:
       run.command(ssroi_cmd + ' T1.nii T1_preBET' + fsl_suffix + ' -b', False)
-
-    # For whatever reason, the output file from standard_space_roi may not be
-    #   completed before BET is run
-    file.waitFor('T1_preBET' + fsl_suffix)
+    pre_bet_image = fsl.findImage('T1_preBET')
 
     # BET
-    fast_t1_input = 'T1_BET' + fsl_suffix
-    run.command(bet_cmd + ' T1_preBET' + fsl_suffix + ' ' + fast_t1_input + ' -f 0.15 -R')
+    run.command(bet_cmd + ' ' + pre_bet_image + ' T1_BET' + fsl_suffix + ' -f 0.15 -R')
+    fast_t1_input = fsl.findImage('T1_BET' + fsl_suffix)
 
     if os.path.exists('T2.nii'):
       if app.args.nocrop:
@@ -159,7 +156,6 @@ def execute():
     run.command(fast_cmd + ' -S 2 ' + fast_t2_input + ' ' + fast_t1_input)
   else:
     run.command(fast_cmd + ' ' + fast_t1_input)
-  fast_output_prefix = fast_t1_input.split('.')[0]
 
   # FIRST
   first_input_is_brain_extracted = ''
@@ -178,7 +174,10 @@ def execute():
       app.console('(note however that FIRST may fail silently, and hence this script may hang indefinitely)')
       file.waitFor(combined_image_path)
     else:
-      app.error('FSL FIRST has failed; not all structures were segmented successfully (check ' + path.toTemp('first.logs', False) + ')')
+      combined_image_path = fsl.findImage('first_all_none_firstseg')
+      if not os.path.isfile(combined_image_path):
+        app.error('FSL FIRST has failed; not all structures were segmented successfully (check ' + 
+path.toTemp('first.logs', False) + ')')
 
   # Convert FIRST meshes to partial volume images
   pve_image_list = [ ]
@@ -192,23 +191,23 @@ def execute():
   pve_cat = ' '.join(pve_image_list)
   run.command('mrmath ' + pve_cat + ' sum - | mrcalc - 1.0 -min all_sgms.mif')
 
-  # Looks like FAST in 5.0 ignores FSLOUTPUTTYPE when writing the PVE images
-  # Will have to wait and see whether this changes, and update the script accordingly
-  if fast_cmd == 'fast':
-    fast_suffix = fsl_suffix
-  else:
-    fast_suffix = '.nii.gz'
-
   # Combine the tissue images into the 5TT format within the script itself
+  fast_output_prefix = fast_t1_input.split('.')[0]
+  fast_csf_output = fsl.findImage(fast_output_prefix + '_pve_0')
+  fast_gm_output = fsl.findImage(fast_output_prefix + '_pve_1')
+  fast_wm_output = fsl.findImage(fast_output_prefix + '_pve_2')
   # Step 1: Run LCC on the WM image
-  run.command('mrthreshold ' + fast_output_prefix + '_pve_2' + fast_suffix + ' - -abs 0.001 | maskfilter - connect - -connectivity | mrcalc 1 - 1 -gt -sub remove_unconnected_wm_mask.mif -datatype bit')
-  # Step 2: Generate the images in the same fashion as the 5ttgen command
-  run.command('mrcalc ' + fast_output_prefix + '_pve_0' + fast_suffix + ' remove_unconnected_wm_mask.mif -mult csf.mif')
+  run.command('mrthreshold ' + fast_wm_output + ' - -abs 0.001 | maskfilter - connect - -connectivity | mrcalc 1 - 1 -gt -sub remove_unconnected_wm_mask.mif -datatype bit')
+  # Step 2: Generate the images in the same fashion as the old 5ttgen binary used to:
+  #   - Preserve CSF as-is
+  #   - Preserve SGM, unless it results in a sum of volume fractions greater than 1, in which case clamp
+  #   - Multiply the FAST volume fractions of GM and CSF, so that the sum of CSF, SGM, CGM and WM is 1.0
+  run.command('mrcalc ' + fast_csf_output + ' remove_unconnected_wm_mask.mif -mult csf.mif')
   run.command('mrcalc 1.0 csf.mif -sub all_sgms.mif -min sgm.mif')
-  run.command('mrcalc 1.0 csf.mif sgm.mif -add -sub ' + fast_output_prefix + '_pve_1' + fast_suffix + ' ' + fast_output_prefix + '_pve_2' + fast_suffix + ' -add -div multiplier.mif')
+  run.command('mrcalc 1.0 csf.mif sgm.mif -add -sub ' + fast_gm_output + ' ' + fast_wm_output + ' -add -div multiplier.mif')
   run.command('mrcalc multiplier.mif -finite multiplier.mif 0.0 -if multiplier_noNAN.mif')
-  run.command('mrcalc ' + fast_output_prefix + '_pve_1' + fast_suffix + ' multiplier_noNAN.mif -mult remove_unconnected_wm_mask.mif -mult cgm.mif')
-  run.command('mrcalc ' + fast_output_prefix + '_pve_2' + fast_suffix + ' multiplier_noNAN.mif -mult remove_unconnected_wm_mask.mif -mult wm.mif')
+  run.command('mrcalc ' + fast_gm_output + ' multiplier_noNAN.mif -mult remove_unconnected_wm_mask.mif -mult cgm.mif')
+  run.command('mrcalc ' + fast_wm_output + ' multiplier_noNAN.mif -mult remove_unconnected_wm_mask.mif -mult wm.mif')
   run.command('mrcalc 0 wm.mif -min path.mif')
   run.command('mrcat cgm.mif sgm.mif wm.mif csf.mif path.mif - -axis 3 | mrconvert - combined_precrop.mif -stride +2,+3,+4,+1')
 

lib/mrtrix3/fsl.py

@@ -1,3 +1,5 @@
+_suffix = ''
+
 # Functions that may be useful for scripts that interface with FMRIB FSL tools
 
 # Get the name of the binary file that should be invoked to run eddy;
@@ -27,25 +29,52 @@ def eddyBinary(cuda):
 
 
 
+# In some versions of FSL, even though we try to predict the names of image files that
+#   FSL commands will generate based on the suffix() function, the FSL binaries themselves
+#   ignore the FSLOUTPUTTYPE environment variable. Therefore, the safest approach is:
+# Whenever receiving an output image from an FSL command, explicitly search for the path
+def findImage(name):
+  import os
+  from mrtrix3 import app
+  basename = name.split('.')[0]
+  if os.path.isfile(basename + suffix()):
+    app.debug('Image at expected location: \"' + basename + suffix() + '\"')
+    return basename + suffix()
+  for suf in ['.nii', '.nii.gz', '.img']:
+    if os.path.isfile(basename + suf):
+      app.debug('Expected image at \"' + basename + suffix() + '\", but found at \"' + basename + suf + '\"')
+      return basename + suf
+  app.error('Unable to find FSL output file for path \"' + name + '\"')
+
+
+
 # For many FSL commands, the format of any output images will depend on the string
 #   stored in 'FSLOUTPUTTYPE'. This may even override a filename extension provided
 #   to the relevant command. Therefore use this function to 'guess' what the names
 #   of images provided by FSL commands will be.
 def suffix():
   import os
   from mrtrix3 import app
+  global _suffix
+  if _suffix:
+    return _suffix
   fsl_output_type = os.environ.get('FSLOUTPUTTYPE', '')
   if fsl_output_type == 'NIFTI':
     app.debug('NIFTI -> .nii')
-    return '.nii'
-  if fsl_output_type == 'NIFTI_GZ':
+    _suffix = '.nii'
+  elif fsl_output_type == 'NIFTI_GZ':
     app.debug('NIFTI_GZ -> .nii.gz')
-    return '.nii.gz'
-  if fsl_output_type == 'NIFTI_PAIR':
+    _suffix = '.nii.gz'
+  elif fsl_output_type == 'NIFTI_PAIR':
     app.debug('NIFTI_PAIR -> .img')
-    return '.img'
-  if fsl_output_type == 'NIFTI_PAIR_GZ':
+    _suffix = '.img'
+  elif fsl_output_type == 'NIFTI_PAIR_GZ':
     app.error('MRtrix3 does not support compressed NIFTI pairs; please change FSLOUTPUTTYPE environment variable')
-  app.warn('Environment variable FSLOUTPUTTYPE not set; FSL commands may fail, or script may fail to locate FSL command outputs')
-  return '.nii.gz'
+  elif fsl_output_type:
+    app.warn('Unrecognised value for environment variable FSLOUTPUTTYPE (\"' + fsl_output_type + '\"): Expecting compressed NIfTIs, but FSL commands may fail')
+    _suffix = '.nii.gz'
+  else:
+    app.warn('Environment variable FSLOUTPUTTYPE not set; FSL commands may fail, or script may fail to locate FSL command outputs')
+    _suffix = '.nii.gz'
+  return _suffix