Merge pull request #3146 from MRtrix3/fix_pehandling_307

Fix phase encoding handling on 3.0.7

Author: Lestropie

bin/dwifslpreproc

@@ -506,7 +506,7 @@ def execute(): #pylint: disable=unused-variable
   # This may be required by -rpe_all for extracting b=0 volumes while retaining phase-encoding information
   import_dwi_pe_table_option = ''
   if dwi_manual_pe_scheme:
-    phaseencoding.save('dwi_manual_pe_scheme.txt', dwi_manual_pe_scheme)
+    phaseencoding.save_table('dwi_manual_pe_scheme.txt', dwi_manual_pe_scheme)
     import_dwi_pe_table_option = ' -import_pe_table dwi_manual_pe_scheme.txt'
 
 
@@ -744,7 +744,7 @@ def execute(): #pylint: disable=unused-variable
   # This may be required when setting up the topup call
   se_epi_manual_pe_table_option = ''
   if se_epi_manual_pe_scheme:
-    phaseencoding.save('se_epi_manual_pe_scheme.txt', se_epi_manual_pe_scheme)
+    phaseencoding.save_table('se_epi_manual_pe_scheme.txt', se_epi_manual_pe_scheme)
     se_epi_manual_pe_table_option = ' -import_pe_table se_epi_manual_pe_scheme.txt'
 
 
@@ -825,8 +825,12 @@ def execute(): #pylint: disable=unused-variable
     # Apply the warp field to the input image series to get an initial corrected volume estimate
     # applytopup can't receive the complete DWI input and correct it as a whole, because the phase-encoding
     #   details may vary between volumes
-    run.command('mrconvert ' + dwi_path + import_dwi_pe_table_option + STRIDES_OPTION_FSL_COMPAT_4D + ' - | '
-                'mrinfo - -export_pe_eddy applytopup_config.txt applytopup_indices.txt')
+    run.command('mrconvert ' + dwi_path + ' applytopup_in.nii'
+                + import_dwi_pe_table_option
+                + STRIDES_OPTION_FSL_COMPAT_4D
+                + ' -export_pe_eddy applytopup_config.txt applytopup_indices.txt'
+                + ' -export_grad_fsl applytopup_in.bvec applytopup_in.bval'
+                + ' -json_export applytopup_in.json')
 
     applytopup_config = matrix.load_matrix('applytopup_config.txt')
     if any(row[2] for row in applytopup_config):
@@ -846,20 +850,36 @@ def execute(): #pylint: disable=unused-variable
     app.debug('applytopup_indices: ' + str(applytopup_indices))
     app.debug('applytopup_volumegroups: ' + str(applytopup_volumegroups))
     for index, group in enumerate(applytopup_volumegroups):
-      prefix = os.path.splitext(dwi_path)[0] + '_pe' + str(index+1)
+      prefix = 'applytopup_in_pe' + str(index+1)
       input_path = prefix + '.nii'
       json_path = prefix + '.json'
+      bvec_path = prefix + '.bvec'
+      bval_path = prefix + '.bval'
       temp_path = prefix + '_applytopup.nii'
       output_path = prefix + '_applytopup.mif'
-      run.command('mrconvert ' + dwi_path + ' ' + input_path + ' -coord 3 ' + ','.join(str(value) for value in group) + STRIDES_OPTION_FSL_COMPAT_4D + ' -json_export ' + json_path)
+      run.command('mrconvert applytopup_in.nii ' + input_path
+                  + ' -json_import applytopup_in.json'
+                  + ' -fslgrad applytopup_in.bvec applytopup_in.bval'
+                  + ' -coord 3 ' + ','.join(str(value) for value in group)
+                  + STRIDES_OPTION_FSL_COMPAT_4D
+                  + ' -json_export ' + json_path
+                  + ' -export_grad_fsl ' + bvec_path + ' ' + bval_path)
       run.command(applytopup_cmd + ' --imain=' + input_path + ' --datain=applytopup_config.txt --inindex=' + str(index+1) + ' --topup=field --out=' + temp_path + ' --method=jac')
       app.cleanup(input_path)
       temp_path = fsl.find_image(temp_path)
-      run.command('mrconvert ' + temp_path + ' ' + output_path + ' -json_import ' + json_path)
+      run.command('mrconvert ' + temp_path + ' ' + output_path
+                  + ' -json_import ' + json_path
+                  + ' -fslgrad ' + bvec_path + ' ' + bval_path)
       app.cleanup(json_path)
       app.cleanup(temp_path)
+      app.cleanup(bvec_path)
+      app.cleanup(bval_path)
       applytopup_image_list.append(output_path)
       index += 1
+    app.cleanup('applytopup_in.nii')
+    app.cleanup('applytopup_in.json')
+    app.cleanup('applytopup_in.bvec')
+    app.cleanup('applytopup_in.bval')
 
     # Use the initial corrected volumes to derive a brain mask for eddy
     if not app.ARGS.eddy_mask:

cmd/mrinfo.cpp

@@ -78,13 +78,6 @@ void usage ()
        "raw gradient table information as stored in the image header, i.e. without "
        "MRtrix3 back-end processing, use \"-property dw_scheme\"."
 
-    + "The -petable option exports the MRtrix3 internal representation of the phase encoding table, "
-      "regardless of whether the relevant metadata are stored in the BIDS fields \"PhaseEncodingDirection\" "
-      "and \"TotalReadoutTime\" or the MRtrix3-specific \"pe_scheme\". The contents of this query "
-      "should however *not* be provided to FSL tools: despite the contents being of the same format, "
-      "the phase encoding directions may be erroneously interpreted. If extracting phase encoding information "
-      "to interface with FSL tools, use the -export_pe_topup or -export_pe_eddy options."
-
     + DWI::bvalue_scaling_description;
 
   ARGUMENTS
@@ -115,7 +108,7 @@ void usage ()
     +   Option ("shell_indices", "list the image volumes attributed to each b-value shell")
 
     + Metadata::PhaseEncoding::ExportOptions
-    +   Option ("petable", "print the MRtrix internal representation of the phase encoding table (see Description)")
+    +   Option ("petable", "print the phase encoding table")
 
     + OptionGroup ("Handling of piped images")
     +   Option ("nodelete", "don't delete temporary images or images passed to mrinfo via Unix pipes");

core/header.cpp

@@ -108,7 +108,7 @@ namespace MR
               auto scheme = DWI::resolve_DW_scheme (parse_matrix (item.second), parse_matrix (it->second));
               DWI::set_DW_scheme (new_keyval, scheme);
             } catch (Exception& e) {
-              WARN("Unable to merge inconsistent DW gradient tables between headers");
+              INFO("Unable to merge inconsistent DW gradient tables between headers");
               new_keyval["dw_scheme"] = "variable";
             }
           }
@@ -342,7 +342,8 @@ namespace MR
       }
       parser.calculate_padding (Pdim);
 
-
+      // FIXME This fails to appropriately assign rows of these schemes to images
+      //   if splitting 4D image into 2D images
       const bool split_4d_schemes = (parser.ndim() == 1 && template_header.ndim() == 4);
       Eigen::MatrixXd dw_scheme, pe_scheme;
       try {

core/metadata/phase_encoding.cpp

@@ -27,6 +27,8 @@ namespace MR {
       using namespace App;
       const OptionGroup ImportOptions =
         OptionGroup("Options for importing phase-encode tables")
+        + Option("import_pe_table", "import a phase-encoding table from file")
+          + Argument("file").type_file_in()
         + Option("import_pe_topup", "import a phase-encoding table intended for FSL TOPUP from file")
           + Argument("file").type_file_in()
         + Option("import_pe_eddy", "import phase-encoding information from an EDDY-style config / index file pair")
@@ -44,6 +46,8 @@ namespace MR {
 
       const OptionGroup ExportOptions =
         OptionGroup("Options for exporting phase-encode tables")
+        + Option("export_pe_table", "export phase-encoding table to file")
+          + Argument("file").type_file_out()
         + Option("export_pe_topup", "export phase-encoding table to a file intended for FSL topup")
           + Argument("file").type_file_out()
         + Option("export_pe_eddy", "export phase-encoding information to an EDDY-style config / index file pair")
@@ -180,15 +184,18 @@ namespace MR {
         DEBUG("searching for suitable phase encoding data...");
         using namespace App;
 
+        const auto opt_table = get_options("import_pe_table");
         const auto opt_topup = get_options("import_pe_topup");
         const auto opt_eddy = get_options("import_pe_eddy");
-        if (opt_topup.size() + opt_eddy.size() > 1)
+        if (opt_table.size() + opt_topup.size() + opt_eddy.size() > 1)
           throw Exception("Cannot specify more than one command-line option"
                           " for importing phase encoding information from external file(s)");
 
         scheme_type result;
         try {
-          if (!opt_topup.empty())
+          if (!opt_table.empty())
+            result = load_table(opt_table[0][0], header);
+          else if (!opt_topup.empty())
             result = load_topup(opt_topup[0][0], header);
           else if (!opt_eddy.empty())
             result = load_eddy(opt_eddy[0][0], opt_eddy[0][1], header);
@@ -220,8 +227,13 @@ namespace MR {
         try {
           pe_scheme = parse_scheme(keyval, H);
         } catch (Exception& e) {
-          WARN("Unable to conform phase encoding information to image realignment "
-               " for image \"" + H.name() + "\"; erasing");
+          if ((keyval.find("PhaseEncodingDirection") != keyval.end()
+              && keyval["PhaseEncodingDirection"] != "variable")
+              || (keyval.find("pe_scheme") != keyval.end()
+              && keyval["pe_scheme"] != "variable")) {
+            WARN("Unable to conform phase encoding information to image realignment"
+                 " for image \"" + H.name() + "\"; erasing");
+          }
           clear_scheme(keyval);
           return;
         }
@@ -344,7 +356,11 @@ namespace MR {
 
         auto scheme = parse_scheme(header.keyval(), header);
 
-        auto opt = get_options("export_pe_topup");
+        auto opt = get_options("export_pe_table");
+        if (!opt.empty())
+          save_table(check(scheme), header, opt[0][0]);
+
+        opt = get_options("export_pe_topup");
         if (!opt.empty())
           save_topup(check(scheme), header, opt[0][0]);
 
@@ -355,7 +371,32 @@ namespace MR {
 
 
 
+      scheme_type load_table(const std::string& path, const Header& header) {
+        if (Path::has_suffix(header.name(), {".nii", ".nii.gz", ".img", ".mgh", "mgz"})) {
+          WARN("Note use of -import_pe_table in conjunction with MGH / NIfTI image"
+               " interprets phase encoding directions as being strictly with respect to image axes,"
+               " not with respect to the FSL internal convention;"
+               " consider if -import_pe_topup is more appropriate for your use case"
+               " (see: mrtrix.readthedocs.org/en/" MRTRIX_BASE_VERSION "/concepts/pe_scheme.html#reference-axes-for-phase-encoding-directions)");
+        }
+        const scheme_type PE = load_matrix(path);
+        check(PE, header);
+        // As with JSON import, need to query the header to discover if the
+        //   strides / transform were modified on image load to make the image
+        //   data appear approximately axial, in which case we need to apply the
+        //   same transforms to the phase encoding data on load
+        return transform_for_image_load(PE, header);
+      }
+
+
+
       scheme_type load_topup(const std::string& path, const Header& header) {
+        if (!Path::has_suffix(header.name(), {".nii", ".nii.gz", ".img", ".mgh", "mgz"})) {
+          WARN("Loading FSL topup format phase encoding information"
+               " accompanying image \"" + header.name() + "\" that is not MGH / NIfTI format"
+               " may be erroneous due to possible flipping of first image axis"
+               " (see: mrtrix.readthedocs.org/en/" MRTRIX_BASE_VERSION "/concepts/pe_scheme.html#reference-axes-for-phase-encoding-directions)");
+        }
         scheme_type PE = load_matrix(path);
         check(PE, header);
         // Flip of first image axis based on determinant of image transform
@@ -369,6 +410,12 @@ namespace MR {
 
 
       scheme_type load_eddy(const std::string& config_path, const std::string& index_path, const Header& header) {
+        if (!Path::has_suffix(header.name(), {".nii", ".nii.gz", ".img", ".mgh", "mgz"})) {
+          WARN("Loading FSL eddy format phase encoding information"
+               " accompanying image \"" + header.name() + "\" that is not MGH / NIfTI format"
+               " may be erroneous due to possible flipping of first image axis"
+               " (see: mrtrix.readthedocs.org/en/" MRTRIX_BASE_VERSION "/concepts/pe_scheme.html#reference-axes-for-phase-encoding-directions)");
+        }
         const Eigen::MatrixXd config = load_matrix(config_path);
         const Eigen::Array<int, Eigen::Dynamic, 1> indices = load_vector<int>(index_path);
         scheme_type PE = eddy2topup(config, indices);

core/metadata/phase_encoding.h

@@ -27,6 +27,7 @@
 #include "header.h"
 #include "metadata/bids.h"
 #include "types.h"
+#include "version.h"
 
 namespace MR {
   namespace Metadata {
@@ -118,7 +119,8 @@ namespace MR {
 
         if (Path::has_suffix(header.name(), {".mgh", ".mgz", ".nii", ".nii.gz", ".img"})) {
           WARN("External phase encoding table \"" + path + "\" for image \"" + header.name() + "\""
-               " may not be suitable for FSL topup; consider use of -export_pe_topup instead");
+               " may not be suitable for FSL topup; consider use of -export_pe_topup instead"
+               " (see: mrtrix.readthedocs.org/en/" MRTRIX_BASE_VERSION "/concepts/pe_scheme.html#reference-axes-for-phase-encoding-directions)");
           save_table(transform_for_nifti_write(PE, header), path, true);
         } else {
           save_table(PE, path, true);
@@ -133,10 +135,11 @@ namespace MR {
           throw Exception(e, "Cannot export phase-encoding table to file \"" + path + "\"");
         }
 
-        // TODO Should this check be in place?
-        if (!Path::has_suffix(header.name(), {".mgh", ".mgz", ".nii", ".nii.gz", ".img"}))
-          throw Exception("Only export phase encoding table to FSL topup format"
-                          " in conjunction with MGH / NIfTI format images");
+        if (!Path::has_suffix(header.name(), {".mgh", ".mgz", ".nii", ".nii.gz", ".img"})) {
+          WARN("Beware use of -export_pe_topup in conjunction image format other than MGH / NIfTI;"
+               " -export_pe_table may be more suitable"
+               " (see: mrtrix.readthedocs.org/en/" MRTRIX_BASE_VERSION "/concepts/pe_scheme.html#reference-axes-for-phase-encoding-directions)");
+        }
 
         scheme_type table = transform_for_nifti_write(PE, header);
         // The flipping of first axis based on the determinant of the image header transform
@@ -154,9 +157,12 @@ namespace MR {
                      const HeaderType& header,
                      const std::string& config_path,
                      const std::string& index_path) {
-        if (!Path::has_suffix(header.name(), {".mgh", ".mgz", ".nii", ".nii.gz", ".img"}))
-          throw Exception("Only export phase encoding table to FSL eddy format"
-                          " in conjunction with MGH / NIfTI format images");
+        if (!Path::has_suffix(header.name(), {".mgh", ".mgz", ".nii", ".nii.gz", ".img"})) {
+          WARN("Exporting phase encoding table to FSL eddy format"
+               " in conjunction with format other than MGH / NIfTI"
+               " risks erroneous interpretation due to possible flipping of first image axis"
+               " (see: mrtrix.readthedocs.org/en/" MRTRIX_BASE_VERSION "/concepts/pe_scheme.html#reference-axes-for-phase-encoding-directions)");
+        }
         scheme_type table = transform_for_nifti_write(PE, header);
         Axes::permutations_type order;
         const auto adjusted_transform = File::NIfTI::adjust_transform(header, order);