MRtrix 3.0 (release candidate 2)

MRtrix 3.0 (release candidate 2)

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
@@ -118,7 +116,7 @@ elif app.args.ants:
   run.command('mrconvert mean_bzero.mif mean_bzero.nii -stride +1,+2,+3')
   run.command('mrconvert mask.mif mask.nii -stride +1,+2,+3')
   bias_path = 'bias.nii'
-  run.command('N4BiasFieldCorrection -d 3 -i mean_bzero.nii -w mask.nii -o [corrected.nii,' + bias_path + '] -s 2 -b [150] -c [200x200,0.0]')
+  run.command('N4BiasFieldCorrection -d 3 -i mean_bzero.nii -w mask.nii -o [corrected.nii,' + bias_path + '] -b [150,3] -c [1000x1000,0.0]')
 
 run.command('mrcalc in.mif ' + bias_path + ' -div result.mif')
 run.command('mrconvert result.mif ' + path.fromUser(app.args.output, True) + (' -force' if app.force else ''))

bin/dwipreproc

@@ -2,7 +2,7 @@
 
 # Script for performing DWI pre-processing using FSL 5.0 tools eddy / topup / applytopup
 
-# This script is generally the first operation that will be applied to diffusion image data (with the possible exception of dwidenoise). The precise details of how this image pre-processing takes place depends heavily on the DWI acquisition; specifically, the presence or absence of reversed phase-encoding data for the purposes of EPI susceptibility distortion correction.
+# This script is generally one of the first operations that will be applied to diffusion image data. The precise details of how this image pre-processing takes place depends heavily on the DWI acquisition; specifically, the presence or absence of reversed phase-encoding data for the purposes of EPI susceptibility distortion correction.
 
 # The script is capable of handling a wide range of DWI acquisitions with respect to the design of phase encoding directions. This is dependent upon information regarding the phase encoding being embedded within theimage headers. The relevant information should be captured by MRtrix when importing DICOM images; it should also be the case for BIDS-compatible datasets. If the user expects this information to be present within the image headers, the -rpe_header option must be specified.
 
@@ -24,7 +24,7 @@ if not os.path.isdir(lib_folder):
 sys.path.insert(0, lib_folder)
 
 
-import math
+import math, shutil
 from distutils.spawn import find_executable
 from mrtrix3 import app, file, fsl, image, path, phaseEncoding, run
 
@@ -182,6 +182,27 @@ if not len(grad) == num_volumes:
 do_topup = (not PE_design == 'None')
 
 
+def grads_match(one, two):
+  # Dot product between gradient directions
+  # First, need to check for zero-norm vectors:
+  # - If both are zero, skip this check
+  # - If one is zero and the other is not, volumes don't match
+  # - If neither is zero, test the dot product
+  if any([val for val in one[0:3]]):
+    if not any([val for val in two[0:3]]):
+      return False
+    dot_product = one[0]*two[0] + one[1]*two[1] + one[2]*two[2]
+    if abs(dot_product) < 0.999:
+      return False
+  elif any([val for val in two[0:3]]):
+    return False
+  # b-value
+  if abs(one[3]-two[3]) > 10.0:
+    return False
+  return True
+
+
+
 # Manually generate a phase-encoding table for the input DWI based on user input
 dwi_manual_pe_scheme = None
 topup_manual_pe_scheme = None
@@ -220,26 +241,39 @@ if manual_pe_dir:
   # If -rpe_all, need to scan through grad and figure out the pairings
   # This will be required if relying on user-specified phase encode direction
   # It will also be required at the end of the script for the manual recombination
+  # Update: The possible permutations of volume-matched acquisition is limited within the
+  #   context of the -rpe_all option. In particular, the potential for having more
+  #   than one b=0 volume within each half means that it is not possible to permit
+  #   arbitrary ordering of those pairs, since b=0 volumes would then be matched
+  #   despite having the same phase-encoding direction. Instead, explicitly enforce
+  #   that volumes must be matched between the first and second halves of the DWI data.
   elif PE_design == 'All':
-    grad_matchings = [ num_volumes ] * num_volumes
+    if num_volumes%2:
+      app.error('If using -rpe_all option, input image must contain an even number of volumes')
+    grads_matched = [ num_volumes ] * num_volumes
     grad_pairs = [ ]
-    for index1 in range(num_volumes):
-      if grad_matchings[index1] == num_volumes: # As yet unpaired
-        for index2 in range(index1+1, num_volumes):
-          if grad_matchings[index2] == num_volumes: # Also as yet unpaired
-            if grad[index1] == grad[index2]:
-              grad_matchings[index1] = index2;
-              grad_matchings[index2] = index1;
+    app.debug('Commencing gradient direction matching; ' + str(num_volumes) + ' volumes')
+    for index1 in range(int(num_volumes/2)):
+      if grads_matched[index1] == num_volumes: # As yet unpaired
+        for index2 in range(int(num_volumes/2), num_volumes):
+          if grads_matched[index2] == num_volumes: # Also as yet unpaired
+            if grads_match(grad[index1], grad[index2]):
+              grads_matched[index1] = index2;
+              grads_matched[index2] = index1;
               grad_pairs.append([index1, index2])
+              app.debug('Matched volume ' + str(index1) + ' with ' + str(index2) + ': ' + str(grad[index1]) + ' ' + str(grad[index2]))
+              break
+        else:
+          app.error('Unable to determine matching reversed phase-encode direction volume for DWI volume ' + str(index1))
     if not len(grad_pairs) == num_volumes/2:
-      app.error('Unable to determine matching DWI volume pairs for reversed phase-encode combination')
+      app.error('Unable to determine complete matching DWI volume pairs for reversed phase-encode combination')
     # Construct manual PE scheme here:
     #   Regardless of whether or not there's a scheme in the header, need to have it:
     #   if there's one in the header, want to compare to the manually-generated one
     dwi_manual_pe_scheme = [ ]
     for index in range(0, num_volumes):
       line = list(manual_pe_dir)
-      if grad_matchings[index] < index:
+      if index >= int(num_volumes/2):
         line = [ (-i if i else 0.0) for i in line ]
       line.append(trt)
       dwi_manual_pe_scheme.append(line)
@@ -435,22 +469,36 @@ if do_topup:
   else:
     run.command('mrinfo dwi.mif -export_pe_eddy applytopup_config.txt applytopup_indices.txt')
 
-  applytopup_index_list = [ ]
+  # Update: Call applytopup separately for each unique phase-encoding
+  # This should be the most compatible option with more complex phase-encoding acquisition designs,
+  #   since we don't need to worry about applytopup performing volume recombination
+  # Plus, recombination doesn't need to be optimal; we're only using this to derive a brain mask
   applytopup_image_list = [ ]
   index = 1
   with open('applytopup_config.txt', 'r') as f:
     for line in f:
-      image_path = 'dwi_pe_' + str(index) + '.nii'
-      run.command('dwiextract dwi.mif' + import_dwi_pe_table_option + ' -pe ' + ','.join(line.split()) + ' ' + image_path)
-      applytopup_index_list.append(str(index))
-      applytopup_image_list.append(image_path)
+      input_path = 'dwi_pe_' + str(index) + '.nii'
+      json_path = 'dwi_pe_' + str(index) + '.json'
+      temp_path = 'dwi_pe_' + str(index) + '_topup' + fsl_suffix
+      output_path = 'dwi_pe_' + str(index) + '_topup.mif'
+      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
 
-  # Finally ready to run applytopup
-  run.command(applytopup_cmd + ' --imain=' + ','.join(applytopup_image_list) + ' --datain=applytopup_config.txt --inindex=' + ','.join(applytopup_index_list) + ' --topup=field --out=dwi_post_topup' + fsl_suffix + ' --method=jac')
-
   # Use the initial corrected volumes to derive a brain mask for eddy
-  run.command('mrconvert dwi_post_topup' + fsl_suffix + ' -grad grad.b - | 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)
 
   eddy_in_topup_option = ' --topup=field'
 
@@ -491,22 +539,26 @@ if not os.path.isfile(bvecs_path):
 # The phase-encoding scheme needs to be checked also
 volume_matchings = [ num_volumes ] * num_volumes
 volume_pairs = [ ]
+app.debug('Commencing gradient direction matching; ' + str(num_volumes) + ' volumes')
 for index1 in range(num_volumes):
   if volume_matchings[index1] == num_volumes: # As yet unpaired
     for index2 in range(index1+1, num_volumes):
       if volume_matchings[index2] == num_volumes: # Also as yet unpaired
         # Here, need to check both gradient matching and reversed phase-encode direction
-        if not any(dwi_pe_scheme[index1][i] + dwi_pe_scheme[index2][i] for i in range(0,3)) and grad[index1] == grad[index2]:
+        if not any(dwi_pe_scheme[index1][i] + dwi_pe_scheme[index2][i] for i in range(0,3)) and grads_match(grad[index1], grad[index2]):
           volume_matchings[index1] = index2;
           volume_matchings[index2] = index1;
           volume_pairs.append([index1, index2])
+          app.debug('Matched volume ' + str(index1) + ' with ' + str(index2) + '\n' +
+                    'Phase encoding: ' + str(dwi_pe_scheme[index1]) + ' ' + str(dwi_pe_scheme[index2]) + '\n' +
+                    'Gradients: ' + str(grad[index1]) + ' ' + str(grad[index2]))
+          break
 
 
-
-if not len(volume_pairs) == num_volumes/2:
+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')
@@ -533,6 +585,14 @@ else:
                  float(bvecs[1][pair[0]]) + float(bvecs[1][pair[1]]),
                  float(bvecs[2][pair[0]]) + float(bvecs[2][pair[1]]) ]
     norm2 = bvec_sum[0]*bvec_sum[0] + bvec_sum[1]*bvec_sum[1] + bvec_sum[2]*bvec_sum[2]
+    # 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.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]]) ]
+      norm2 = bvec_sum[0]*bvec_sum[0] + bvec_sum[1]*bvec_sum[1] + bvec_sum[2]*bvec_sum[2]
     # Occasionally a bzero volume can have a zero vector
     if norm2:
       factor = 1.0 / math.sqrt(norm2)
@@ -559,12 +619,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'
 
 
 
@@ -594,27 +654,37 @@ else:
     run.command('mrcalc ' + jacobian_path + ' ' + jacobian_path + ' -mult weight' + str(index+1) + '.mif')
     file.delTempFile(jacobian_path)
 
+  # 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 = 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'
+
   # This section extracts the two volumes corresponding to each reversed phase-encoded volume pair, and
   #   derives a single image volume based on the recombination equation
   combined_image_list = [ ]
   for index, volumes in enumerate(volume_pairs):
     pe_indices = [ eddy_indices[i] for i in volumes ]
-    run.command('mrconvert dwi_post_eddy' + fsl_suffix + ' volume0.mif -coord 3 ' + str(volumes[0]))
-    run.command('mrconvert dwi_post_eddy' + fsl_suffix + ' volume1.mif -coord 3 ' + str(volumes[1]))
+    run.command('mrconvert ' + eddy_output + ' volume0.mif -coord 3 ' + str(volumes[0]))
+    run.command('mrconvert ' + eddy_output + ' volume1.mif -coord 3 ' + str(volumes[1]))
     # Volume recombination equation described in Skare and Bammer 2010
     run.command('mrcalc volume0.mif weight' + str(pe_indices[0]) + '.mif -mult volume1.mif weight' + str(pe_indices[1]) + '.mif -mult -add weight' + str(pe_indices[0]) + '.mif weight' + str(pe_indices[1]) + '.mif -add -divide 0.0 -max combined' + str(index) + '.mif')
     combined_image_list.append('combined' + str(index) + '.mif')
-    file.delTempFile('volume0.mif')
-    file.delTempFile('volume1.mif')
+    os.remove('volume0.mif')
+    os.remove('volume1.mif')
+  file.delTempFile(eddy_output)
 
   for index in range(0, len(eddy_config)):
     file.delTempFile('weight' + str(index+1) + '.mif')
 
   # Finally the recombined volumes must be concatenated to produce the resulting image series
   run.command('mrcat ' + ' '.join(combined_image_list) + ' - -axis 3 | mrconvert - result.mif -fslgrad bvecs_combined bvals_combined' + stride_option)
 
-  for path in combined_image_list:
-    file.delTempFile(path)
+  for entry in combined_image_list:
+    file.delTempFile(entry)
 
 
 

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

bin/population_template

@@ -188,6 +188,7 @@ nloptions.add_argument('-nl_grad_step', default='0.5', help='The gradient step s
 
 options = app.cmdline.add_argument_group('Input, output and general options')
 options.add_argument('-type', help='Specifiy the types of registration stages to perform. Options are "rigid" (perform rigid registration only which might be useful for intra-subject registration in longitudinal analysis), "affine" (perform affine registration) and "nonlinear" as well as cominations of registration types: %s. Default: rigid_affine_nonlinear' % ', '.join('"'+x+'"' for x in registration_modes if "_" in x), default='rigid_affine_nonlinear')
+options.add_argument('-voxel_size', help='Define the template voxel size in mm. Use either a single value for isotropic voxels or 3 comma separated values.')
 options.add_argument('-initial_alignment', default='mass', help='Method of alignment to form the initial template. Options are "mass" (default), "geometric" and "none".')
 options.add_argument('-mask_dir', help='Optionally input a set of masks inside a single directory, one per input image (with the same file name prefix). Using masks will speed up registration significantly')
 options.add_argument('-warp_dir', help='Output a directory containing warps from each input to the template. If the folder does not exist it will be created')
@@ -216,6 +217,18 @@ if len(inFiles) <= 1:
 else:
   app.console('Generating a population-average template from ' + str(len(inFiles)) + ' input images')
 
+voxel_size = None
+if app.args.voxel_size:
+  voxel_size = app.args.voxel_size.split()
+  if len(voxel_size) == 1:
+    voxel_size = voxel_size * 3
+  try:
+    if len(voxel_size) != 3:
+      raise
+    [float(v) for v in voxel_size]
+  except:
+    app.error('voxel size needs to be a single or three comma separated floating point numbers, received: '+str(app.args.voxel_size))
+
 initial_alignment = app.args.initial_alignment
 if initial_alignment not in ["mass", "geometric", "none"]:
   message.error('initial_alignment must be one of ' + " ".join(["mass", "geometric", "none"]));
@@ -428,7 +441,10 @@ app.console('Generating initial template')
 input_filenames = []
 for i in input:
   input_filenames.append (abspath(i.directory, i.filename));
-run.command('mraverageheader ' + ' '.join(input_filenames) + ' average_header.mif' ' -fill')
+if voxel_size is None:
+  run.command('mraverageheader ' + ' '.join(input_filenames) + ' average_header.mif -fill')
+else:
+  run.command('mraverageheader -fill ' + ' '.join(input_filenames) + ' - | mrresize - -voxel '+','.join(voxel_size)+' average_header.mif')
 
 # crop average space to extent defined by original masks
 if useMasks:
@@ -480,7 +496,10 @@ else:
              ' ' + os.path.join('masks_transformed', i.prefix + '_translated.mif'))
   # update average space to new extent
   run.command('mraverageheader ' + ' '.join([os.path.join('input_transformed', i.prefix + '_translated.mif') for i in input]) + ' average_header_tight.mif')
-  run.command('mrpad -uniform 10 average_header_tight.mif average_header.mif -force')
+  if voxel_size is None:
+    run.command('mrpad -uniform 10 average_header_tight.mif average_header.mif -force')
+  else:
+    run.command('mrpad -uniform 10 average_header_tight.mif - | mrresize - -voxel '+','.join(voxel_size)+' average_header.mif -force')
   run.function(remove, 'average_header_tight.mif')
   if useMasks:
     # reslice masks