Merge pull request #506 from MRtrix3/fix_strides_for_FSL

Changing strides in scripts for input to FSL

scripts/5ttgen

@@ -42,7 +42,7 @@ algorithm = getattr(_5ttgen, lib.app.args.algorithm)
 lib.app.checkOutputFile(lib.app.args.output)
 algorithm.checkOutputFiles()
 
-runCommand('mrconvert ' + lib.app.args.input + ' ' + os.path.join(lib.app.tempDir, 'input.mif') + ' -stride +1,+2,+3')
+runCommand('mrconvert ' + lib.app.args.input + ' ' + os.path.join(lib.app.tempDir, 'input.mif'))
 algorithm.getInputFiles()
 
 lib.app.gotoTempDir()

scripts/dwibiascorrect

@@ -69,7 +69,7 @@ lib.app.checkOutputFile(lib.app.args.bias)
 
 runCommand('mrconvert ' + lib.app.args.input + ' ' + os.path.join(lib.app.tempDir, 'in.mif') + grad_import_option)
 if lib.app.args.mask:
-  runCommand('mrconvert ' + lib.app.args.mask + ' ' + os.path.join(lib.app.tempDir, 'mask.nii') + ' -stride +1,+2,+3')
+  runCommand('mrconvert ' + lib.app.args.mask + ' ' + os.path.join(lib.app.tempDir, 'mask.mif'))
 
 lib.app.gotoTempDir()
 
@@ -83,28 +83,30 @@ if len(DW_scheme) != int(dwi_sizes[3]):
 
 # Generate a brain mask if required, or check the mask if provided
 if lib.app.args.mask:
-  mask_sizes = getHeaderInfo('mask.nii', 'size').split()
+  mask_sizes = getHeaderInfo('mask.mif', 'size').split()
   if not mask_sizes[:3] == dwi_sizes[:3]:
     errorMessage('Provided mask image does not match input DWI')
 else:
-  runCommand('dwi2mask in.mif - | mrconvert - mask.nii -stride +1,+2,+3')
+  runCommand('dwi2mask in.mif mask.mif')
 
 # Make a NiFTI image since this is what the external softwares need
 runCommand('dwiextract in.mif bzeros.mif -bzero')
 bzero_sizes = getHeaderInfo('bzeros.mif', 'size')
 if len(bzero_sizes.split()) == 4:
-  runCommand('mrmath bzeros.mif mean - -axis 3 | mrconvert - mean_bzero.nii -stride +1,+2,+3')
+  runCommand('mrmath bzeros.mif mean - -axis 3 | mrconvert - mean_bzero.mif -stride +1,+2,+3')
 else:
-  runCommand('mrconvert bzeros.mif mean_bzero.nii -stride +1,+2,+3')
+  runCommand('mrconvert bzeros.mif mean_bzero.mif')
 
 if lib.app.args.fsl:
   # FAST doesn't accept a mask input; therefore need to explicitly mask the input image
-  runCommand('mrcalc mean_bzero.nii mask.nii -mult mean_bzero_masked.nii')
+  runCommand('mrcalc mean_bzero.mif mask.mif -mult - | mrconvert - mean_bzero_masked.nii -stride -1,+2,+3')
   runCommand(fast_cmd + ' -t 2 -o fast -n 3 -b mean_bzero_masked.nii')
   bias_path = 'fast_bias' + fast_suffix
 elif lib.app.args.ants:
   # Use the brain mask as a weights image rather than a mask; means that voxels at the edge of the mask
   #   will have a smoothly-varying bias field correction applied, rather than multiplying by 1.0 outside the mask
+  runCommand('mrconvert mean_bzero.mif mean_bzero.nii -stride +1,+2,+3')
+  runCommand('mrconvert mask.mif mask.nii -stride +1,+2,+3')
   bias_path = 'bias.nii'
   runCommand('N4BiasFieldCorrection -d 3 -i mean_bzero.nii -w mask.nii -o [corrected.nii,' + bias_path + '] -s 2 -b [150] -c [200x200,0.0]')
 

scripts/dwipreproc

@@ -103,7 +103,6 @@ lib.app.gotoTempDir()
 series_size = getHeaderInfo('series.mif', 'size').split()
 grad = getHeaderInfo('series.mif', 'dwgrad').split('\n')
 stride = getHeaderInfo('series.mif', 'stride')
-transform = getHeaderInfo('series.mif', 'transform')
 if PE_design == 'Pair':
   Pair1_size  = getHeaderInfo('pair1.mif',  'size').split()
   Pair2_size  = getHeaderInfo('pair2.mif',  'size').split()
@@ -183,17 +182,17 @@ if not PE_design == 'None':
 
 # Convert the input files as necessary for FSL tools
 if PE_design == 'None':
-  runCommand('mrconvert ' + series_path + ' dwi_pre_topup.nii -stride +1,+2,+3,+4')
+  runCommand('mrconvert ' + series_path + ' dwi_pre_topup.nii -stride -1,+2,+3,+4')
 if PE_design == 'Pair':
-  runCommand('mrconvert ' + series_path + ' dwi_pre_topup.nii -stride +1,+2,+3,+4')
-  runCommand('mrcat ' + pair1_path + ' ' + pair2_path + ' - -axis 3 | mrconvert - topup_in.nii -stride +1,+2,+3,+4')
+  runCommand('mrconvert ' + series_path + ' dwi_pre_topup.nii -stride -1,+2,+3,+4')
+  runCommand('mrcat ' + pair1_path + ' ' + pair2_path + ' - -axis 3 | mrconvert - topup_in.nii -stride -1,+2,+3,+4')
 elif PE_design == 'All':
   runCommand('dwiextract ' + series_path + ' -bzero pair1.mif')
   runCommand('dwiextract ' + series2_path + ' -bzero pair2.mif')
-  runCommand('mrcat pair1.mif pair2.mif - -axis 3 | mrconvert - topup_in.nii -stride +1,+2,+3,+4')
-  runCommand('mrconvert ' + series_path + ' dwi1_pre_topup.nii -stride +1,+2,+3,+4')
+  runCommand('mrcat pair1.mif pair2.mif - -axis 3 | mrconvert - topup_in.nii -stride -1,+2,+3,+4')
+  runCommand('mrconvert ' + series_path + ' dwi1_pre_topup.nii -stride -1,+2,+3,+4')
   delFile(series_path)
-  runCommand('mrconvert ' + series2_path + ' dwi2_pre_topup.nii -stride +1,+2,+3,+4')
+  runCommand('mrconvert ' + series2_path + ' dwi2_pre_topup.nii -stride -1,+2,+3,+4')
   delFile(series2_path)
   runCommand('mrcat dwi1_pre_topup.nii dwi2_pre_topup.nii dwi_pre_eddy.nii -axis 3')
 
@@ -249,8 +248,8 @@ eddy_in_topup = ''
 if PE_design == 'None':
 
   # Generate a processing mask for eddy based on the input series
-  runCommand('dwi2mask ' + series_path + ' - | maskfilter - dilate - | mrconvert - mask.nii -datatype float32 -stride +1,+2,+3')
-  runCommand('mrconvert ' + series_path + ' - -stride +1,+2,+3,+4 | mrinfo - -export_grad_fsl bvecs bvals')
+  runCommand('dwi2mask ' + series_path + ' - | maskfilter - dilate - | mrconvert - mask.nii -datatype float32 -stride -1,+2,+3')
+  runCommand('mrconvert ' + series_path + ' - -stride -1,+2,+3,+4 | mrinfo - -export_grad_fsl bvecs bvals')
 
 else:
 
@@ -276,7 +275,7 @@ else:
   # Create the diffusion gradient table in FSL format
   # Make sure the strides are identical to the image actually being passed to eddy before exporting the gradient table
   if PE_design == 'Pair':
-    runCommand('mrconvert ' + series_path + ' - -stride +1,+2,+3,+4 | mrinfo - -export_grad_fsl bvecs bvals')
+    runCommand('mrconvert ' + series_path + ' - -stride -1,+2,+3,+4 | mrinfo - -export_grad_fsl bvecs bvals')
   else:
     # Concatenate the diffusion gradient table twice
     with open('grad.b', 'w') as outfile:
@@ -288,7 +287,7 @@ else:
 
 
   # Use the initial corrected image series from applytopup to derive a processing mask for eddy
-  runCommand('mrconvert dwi_post_topup' + fsl_suffix + ' -fslgrad bvecs bvals - | dwi2mask - - | maskfilter - dilate - | mrconvert - mask.nii -datatype float32 -stride +1,+2,+3')
+  runCommand('mrconvert dwi_post_topup' + fsl_suffix + ' -fslgrad bvecs bvals - | dwi2mask - - | maskfilter - dilate - | mrconvert - mask.nii -datatype float32 -stride -1,+2,+3')
 
   eddy_in_topup = ' --topup=field'
 
@@ -356,27 +355,36 @@ else: # 'All'
   bvals = bvals[0:num_volumes]
   with open('bvals_combined', 'w') as f:
     f.write(' '.join(bvals))
-
 
-  # Derive the weight images
-  # To properly get a voxel displacement field, need to scale the field map (which is output in Hz) appropriately
-  # Scaling appears to be correct; however FSL is a bit ambiguous when it comes to image transforms & strides
-  # Looks like topup flips the x-axis on output (both corrected images and field), and erases image transform
-  # Eddy seems to expect this and operates accordingly; but since we're doing things external to FSL and
-  #   voxel-by-voxel, need to do the flip explicitly; also import the transform from the first input image
-  #   for the sake of visualisation
-  with open( 'transform.txt', 'w' ) as f:
-    for line in transform:
-      f.write (line)
-  runCommand('mrtransform field_map' + fsl_suffix + ' - -flip 0 -linear transform.txt -replace | mrconvert - field_map_flip.mif' + stride_option)
-  delFile('field_map' + fsl_suffix)
+  # Prior to 5.0.8, a bug resulted in the output field map image from topup having an identity transform,
+  #   regardless of the transform of the input image
+  # 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)
+  input_transform_text = getHeaderInfo('topup_in.nii', 'transform')
+  input_transform = [ float(f) for f in input_transform_text.replace('\n', ' ').replace(',', ' ').split() ]
+  topup_transform = [ float(f) for f in getHeaderInfo('field_map' + fsl_suffix, 'transform').replace('\n', ' ').replace(',', ' ').split() ]
+  transform_error = False
+  field_map_image = 'field_map' + fsl_suffix
+  for i, t in zip(input_transform, topup_transform):
+    if (abs(i-t) / (0.5*(i+t))) > 0.001:
+      transform_error = True
+  if transform_error:
+    warnMessage('topup output field image has incorrect transform; recommend updating FSL to version 5.0.8 or later')
+    with open( 'transform.txt', 'w' ) as f:
+      for line in input_transform_text:
+        f.write (line)
+    field_map_image = 'field_map_fix' + fsl_suffix
+    runCommand('mrtransform field_map' + fsl_suffix + ' -linear transform.txt -replace ' + field_map_image)
+    delFile('field_map' + fsl_suffix)
 
 
+  # Derive the weight images
   # Scaling term for field map is identical to the bandwidth provided in the topup config file
   #   (converts Hz to pixel count; that way a simple image gradient can be used to get the Jacobians)
   # Let mrfilter apply the default 1 voxel size gaussian smoothing filter before calculating the field gradient
-  runCommand('mrcalc field_map_flip.mif 0.1 -mult - | mrfilter - gradient - | mrconvert - field_deriv_pe.mif -coord 3 1 -axes 0,1,2')
-  delFile('field_map_flip.mif')
+  runCommand('mrcalc ' + field_map_image + ' 0.1 -mult - | mrfilter - gradient - | mrconvert - field_deriv_pe.mif -coord 3 1 -axes 0,1,2')
+  delFile(field_map_image)
   runCommand('mrcalc 1.0 field_deriv_pe.mif -add 0.0 -max jacobian1.mif')
   runCommand('mrcalc 1.0 field_deriv_pe.mif -sub 0.0 -max jacobian2.mif')
   delFile('field_deriv_pe.mif')

scripts/labelsgmfix

@@ -76,7 +76,7 @@ with open(lib.app.args.config) as f:
 
 # Get the parcellation and T1 images into the temporary directory, with conversion of the T1 into the correct format for FSL
 runCommand('mrconvert ' + lib.app.args.parc + ' ' + os.path.join(lib.app.tempDir, 'parc.mif'))
-runCommand('mrconvert ' + lib.app.args.t1   + ' ' + os.path.join(lib.app.tempDir, 'T1.nii')   + ' -stride +1,+2,+3')
+runCommand('mrconvert ' + lib.app.args.t1   + ' ' + os.path.join(lib.app.tempDir, 'T1.nii') + ' -stride -1,+2,+3')
 
 lib.app.gotoTempDir()
 
@@ -90,7 +90,7 @@ runCommand(first_cmd + ' -s ' + ','.join(structure_map.keys()) + ' -i T1.nii' +
 # In this use case, don't want the PVE images; want to threshold at 0.5
 mask_list = [ ]
 for struct in structure_map.keys():
-  image_path = struct + '_mask.nii'
+  image_path = struct + '_mask.mif'
   mask_list.append (image_path)
   vtk_in_path = 'first-' + struct + '_first.vtk'
   if not os.path.exists(vtk_in_path):
@@ -111,8 +111,8 @@ for struct in structure_map.keys():
 # Detect any overlapping voxels between these masks
 # These will be set to 0 in the final parcellated image
 result_path = 'result' + os.path.splitext(lib.app.args.output)[1]
-runCommand('mrmath ' + ' '.join(mask_list) + ' sum - | mrcalc - 1 -gt overlap_mask.nii')
-runCommand('mrcalc overlap_mask.nii 0 parc.mif -if result.mif')
+runCommand('mrmath ' + ' '.join(mask_list) + ' sum - | mrcalc - 1 -gt overlap_mask.mif')
+runCommand('mrcalc overlap_mask.mif 0 parc.mif -if result.mif')
 runCommand('mrconvert result.mif ' + os.path.join(lib.app.workingDir, lib.app.args.output) + lib.app.mrtrixForce)
 
 lib.app.complete()

scripts/src/_5ttgen/fsl.py

@@ -19,7 +19,7 @@ def getInputFiles():
   import lib.app
   from lib.runCommand   import runCommand
   if hasattr(lib.app.args, 'mask') and lib.app.args.mask is not None:
-    runCommand('mrconvert ' + lib.app.args.mask + ' ' + os.path.join(lib.app.tempDir, 'mask.mif') + ' -datatype bit -stride +1,+2,+3')
+    runCommand('mrconvert ' + lib.app.args.mask + ' ' + os.path.join(lib.app.tempDir, 'mask.mif') + ' -datatype bit -stride -1,+2,+3')
 
 
 
@@ -76,23 +76,23 @@ def execute():
   if lib.app.args.sgm_amyg_hipp:
     sgm_structures.extend([ 'L_Amyg', 'R_Amyg', 'L_Hipp', 'R_Hipp' ])
 
-  runCommand('mrconvert input.mif T1.nii')  
+  runCommand('mrconvert input.mif T1.nii -stride -1,+2,+3')
 
   # Decide whether or not we're going to do any brain masking
   if os.path.exists('mask.mif'):
 
     # Check to see if the dimensions match the T1 image
-    T1_size = getHeaderInfo('input.mif', 'size')
+    T1_size = getHeaderInfo('T1.nii', 'size')
     mask_size = getHeaderInfo('mask.mif', 'size')
     if mask_size == T1_size:
-      runCommand('mrcalc input.mif mask.mif -mult T1_masked.nii')
+      runCommand('mrcalc input.mif mask.mif -mult - | mrconvert - T1_masked' + fsl_suffix + ' -stride -1,+2,+3')
     else:
-      runCommand('mrtransform mask.mif mask_regrid.mif -template input.mif')
-      runCommand('mrcalc input.mif mask_regrid.mif -mult T1_masked' + fsl_suffix)
+      runCommand('mrtransform mask.mif mask_regrid.mif -template T1.nii')
+      runCommand('mrcalc input.mif mask_regrid.mif -mult - | mrconvert - T1_masked' + fsl_suffix)
 
   elif lib.app.args.premasked:
 
-    runCommand('mrconvert input.mif T1_masked' + fsl_suffix)
+    runCommand('mrconvert input.mif T1_masked' + fsl_suffix + ' -stride -1,+2,+3')
 
   else:
 
@@ -121,7 +121,7 @@ def execute():
 
     if not os.path.isfile('T1_preBET' + fsl_suffix):
       warnMessage('FSL command ' + ssroi_cmd + ' appears to have failed; passing T1 directly to BET')
-      runCommand('mrconvert input.mif T1_preBET' + fsl_suffix)
+      runCommand('mrconvert input.mif T1_preBET' + fsl_suffix + ' -stride -1,+2,+3')
 
     # BET
     runCommand(bet_cmd + ' T1_preBET' + fsl_suffix + ' T1_masked' + fsl_suffix + ' -f 0.15 -R')
@@ -140,16 +140,16 @@ def execute():
   # Convert FIRST meshes to partial volume images
   pve_image_list = [ ]
   for struct in sgm_structures:
-    pve_image_path = 'mesh2pve_' + struct + '.nii'
+    pve_image_path = 'mesh2pve_' + struct + '.mif'
     vtk_in_path = 'first-' + struct + '_first.vtk'
     vtk_temp_path = struct + '.vtk'
     if not os.path.exists(vtk_in_path):
       errorMessage('Missing .vtk file for structure ' + struct + '; run_first_all must have failed')
-    runCommand('meshconvert ' + vtk_in_path + ' ' + vtk_temp_path + ' -transform_first2real input.mif')
+    runCommand('meshconvert ' + vtk_in_path + ' ' + vtk_temp_path + ' -transform_first2real T1.nii')
     runCommand('mesh2pve ' + vtk_temp_path + ' T1_masked' + fsl_suffix + ' ' + pve_image_path)
     pve_image_list.append(pve_image_path)
   pve_cat = ' '.join(pve_image_list)
-  runCommand('mrmath ' + pve_cat + ' sum - | mrcalc - 1.0 -min all_sgms.nii')
+  runCommand('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
@@ -163,7 +163,7 @@ def execute():
   runCommand('mrthreshold T1_masked_pve_2' + fast_suffix + ' - -abs 0.001 | maskfilter - connect wm_mask.mif -largest')
   # Step 2: Generate the images in the same fashion as the 5ttgen command
   runCommand('mrconvert T1_masked_pve_0' + fast_suffix + ' csf.mif')
-  runCommand('mrcalc 1 csf.mif -sub all_sgms.nii -min sgm.mif')
+  runCommand('mrcalc 1 csf.mif -sub all_sgms.mif -min sgm.mif')
   runCommand('mrcalc 1.0 csf.mif sgm.mif -add -sub T1_masked_pve_1' + fast_suffix + ' T1_masked_pve_2' + fast_suffix + ' -add -div multiplier.mif')
   runCommand('mrcalc multiplier.mif -finite multiplier.mif 0.0 -if multiplier_noNAN.mif')
   runCommand('mrcalc T1_masked_pve_1' + fast_suffix + ' multiplier_noNAN.mif -mult cgm.mif')

src/dwi/gradient.cpp

@@ -12,7 +12,9 @@
  * For more details, see www.mrtrix.org
  * 
  */
+
 #include "dwi/gradient.h"
+#include "file/nifti1_utils.h"
 
 namespace MR
 {
@@ -80,10 +82,17 @@ namespace MR
       if (bvals.cols() != bvecs.cols() || bvals.cols() != header.size (3))
         throw Exception ("bvals and bvecs files must have same number of diffusion directions as DW-image");
 
+      // bvecs format actually assumes a LHS coordinate system even if image is
+      // stored using RHS - x axis is flipped to make linear 3x3 part of
+      // transform have negative determinant:
+      std::vector<size_t> order;
+      auto adjusted_transform = File::NIfTI::adjust_transform (header, order);
+      if (adjusted_transform.linear().determinant() > 0.0) 
+        bvecs.row(0) = -bvecs.row(0);
+
       // account for the fact that bvecs are specified wrt original image axes,
       // which may have been re-ordered and/or inverted by MRtrix to match the
       // expected anatomical frame of reference:
-      std::vector<size_t> order = Stride::order (header, 0, 3);
       MatrixXd G (bvecs.cols(), 3);
       for (ssize_t n = 0; n < G.rows(); ++n) {
         G(n,order[0]) = header.stride(order[0]) > 0 ? bvecs(0,n) : -bvecs(0,n);
@@ -116,7 +125,8 @@ namespace MR
 
       // deal with FSL requiring gradient directions to coincide with data strides
       // also transpose matrices in preparation for file output
-      auto order = Stride::order (header, 0, 3);
+      std::vector<size_t> order;
+      auto adjusted_transform = File::NIfTI::adjust_transform (header, order);
       MatrixXd bvecs (3, grad.rows());
       MatrixXd bvals (1, grad.rows());
       for (ssize_t n = 0; n < G.rows(); ++n) {
@@ -126,6 +136,12 @@ namespace MR
         bvals(0,n) = grad(n,3);
       }
 
+      // bvecs format actually assumes a LHS coordinate system even if image is
+      // stored using RHS - x axis is flipped to make linear 3x3 part of
+      // transform have negative determinant:
+      if (adjusted_transform.linear().determinant() > 0.0) 
+        bvecs.row(0) = -bvecs.row(0);
+
       save_matrix (bvecs, bvecs_path);
       save_matrix (bvals, bvals_path);
     }