dwicat: Support different voxel grids

cmd/mraverageheader.cpp

@@ -26,10 +26,9 @@ using namespace MR;
 using namespace App;
 using namespace Registration;
 
-default_type PADDING_DEFAULT = 0.0;
-
-enum RESOLUTION { MAX, MEAN };
-const char *resolution_choices[] = {"max", "mean", nullptr};
+const default_type PADDING_DEFAULT = 0.0;
+const avgspace_voxspacing_t SPACING_DEFAULT_VALUE = avgspace_voxspacing_t::MEAN_NEAREST;
+const std::string SPACING_DEFAULT_STRING = "mean_nearest";
 
 // clang-format off
 void usage() {
@@ -38,22 +37,32 @@ void usage() {
 
   SYNOPSIS = "Calculate the average (unbiased) coordinate space of all input images";
 
+  DESCRIPTION
+  +"The voxel spacings of the calculated average header can be determined from the set of "
+   "input images in one of four different ways, "
+   "which may be more or less appropriate in different use cases. "
+   "These options vary in two key ways: "
+   "1. Projected voxel spacing of the input image in the direction of the average header axes "
+   "versus the voxel spacing of the input image axis that is closest to the average header axis; "
+   "2. Selecting the minimum of these spacings across input images to maintain maximal "
+   "spatial resolution versus the mean across images to produce an unbiased average.";
+
   ARGUMENTS
   + Argument ("input", "the input image(s).").type_image_in().allow_multiple()
   + Argument ("output", "the output image").type_image_out();
 
   OPTIONS
-  + Option ("padding", " boundary box padding in voxels."
-                       " Default: " + str(PADDING_DEFAULT))
+  + Option ("padding",
+            "boundary box padding in voxels."
+            " Default: " + str(PADDING_DEFAULT))
     + Argument ("value").type_float(0.0)
-
-  + Option ("resolution", "subsampling of template compared to smallest voxel size in any input image."
-                          "Valid options are: "
-                          "- 'mean': unbiased but loss of resolution for individual images possible; "
-                          "- 'max': smallest voxel size of any input image defines the resolution."
-                          " Default: mean")
-    + Argument ("type").type_choice (resolution_choices)
-
+  + Option ("spacing",
+            "Method for determination of voxel spacings based on"
+            " the set of input images and the average header axes"
+            " (see Description)."
+            " Valid options are: " + join(avgspace_voxspacing_choices, ",") + ";"
+            " default = " + SPACING_DEFAULT_STRING)
+    + Argument("type").type_choice(avgspace_voxspacing_choices)
   + Option ("fill", "set the intensity in the first volume of the average space to 1")
   + DataType::options();
 
@@ -64,14 +73,12 @@ void run() {
 
   const size_t num_inputs = argument.size() - 1;
 
-  auto opt = get_options("padding");
-  const default_type p = opt.size() ? default_type(opt[0][0]) : PADDING_DEFAULT;
+  const default_type p = get_option_value("padding", PADDING_DEFAULT);
   auto padding = Eigen::Matrix<default_type, 4, 1>(p, p, p, 1.0);
   INFO("padding in template voxels: " + str(padding.transpose().head<3>()));
 
-  opt = get_options("resolution");
-  const int resolution = opt.size() ? int(opt[0][0]) : 1;
-  INFO("template voxel subsampling: " + str(resolution));
+  auto opt = get_options("spacing");
+  const avgspace_voxspacing_t spacing = opt.empty() ? SPACING_DEFAULT_VALUE : avgspace_voxspacing_t(int(opt[0][0]));
 
   bool fill = get_options("fill").size();
 
@@ -92,7 +99,7 @@ void run() {
   }
 
   std::vector<Eigen::Transform<default_type, 3, Eigen::Projective>> transform_header_with;
-  auto H = compute_minimum_average_header(headers_in, transform_header_with, resolution, padding);
+  auto H = compute_minimum_average_header(headers_in, transform_header_with, spacing, padding);
   H.datatype() = DataType::Bit;
   if (fill) {
     H.ndim() = dim;

cmd/transformcalc.cpp

@@ -70,7 +70,7 @@ void usage() {
 
   + Example ("Calculate the transformation matrix from an original image"
              " and an image with modified header",
-             "transformcalc mov mapmovhdr header output",
+             "transformcalc orig_image modified_image header output",
              "")
 
   + Example ("Calculate the average affine matrix of a set of input matrices",

core/axes.cpp

@@ -69,15 +69,26 @@ Eigen::Vector3d id2dir(const std::string &id) {
     throw Exception("Malformed image axis identifier: \"" + id + "\"");
 }
 
-void get_permutation_to_make_axial(const transform_type &T, std::array<size_t, 3> &perm, std::array<bool, 3> &flip) {
+void get_shuffle_to_make_axial(const transform_type &T, std::array<size_t, 3> &perm, std::array<bool, 3> &flip) {
+  perm = closest(T.matrix().topLeftCorner<3, 3>());
+  // Figure out whether any of the rows of the transform point in the
+  //   opposite direction to the MRtrix convention
+  flip[perm[0]] = T(0, perm[0]) < 0.0;
+  flip[perm[1]] = T(1, perm[1]) < 0.0;
+  flip[perm[2]] = T(2, perm[2]) < 0.0;
+}
+
+std::array<size_t, 3> closest(const Eigen::Matrix3d &M) {
+  std::array<size_t, 3> result{3, 3, 3};
   // Find which row of the transform is closest to each scanner axis
-  decltype(T.matrix().topLeftCorner<3, 3>())::Index index;
-  T.matrix().topLeftCorner<3, 3>().row(0).cwiseAbs().maxCoeff(&index);
-  perm[0] = index;
-  T.matrix().topLeftCorner<3, 3>().row(1).cwiseAbs().maxCoeff(&index);
-  perm[1] = index;
-  T.matrix().topLeftCorner<3, 3>().row(2).cwiseAbs().maxCoeff(&index);
-  perm[2] = index;
+  Eigen::Matrix3d::Index index(0);
+  M.row(0).cwiseAbs().maxCoeff(&index);
+  result[0] = index;
+  M.row(1).cwiseAbs().maxCoeff(&index);
+  result[1] = index;
+  M.row(2).cwiseAbs().maxCoeff(&index);
+  result[2] = index;
+  assert(result[0] < 3 && result[1] < 3 && result[2] < 3);
 
   // Disambiguate permutations
   auto not_any_of = [](size_t a, size_t b) -> size_t {
@@ -89,19 +100,15 @@ void get_permutation_to_make_axial(const transform_type &T, std::array<size_t, 3
     assert(0);
     return std::numeric_limits<size_t>::max();
   };
-  if (perm[0] == perm[1])
-    perm[1] = not_any_of(perm[0], perm[2]);
-  if (perm[0] == perm[2])
-    perm[2] = not_any_of(perm[0], perm[1]);
-  if (perm[1] == perm[2])
-    perm[2] = not_any_of(perm[0], perm[1]);
-  assert(perm[0] != perm[1] && perm[1] != perm[2] && perm[2] != perm[0]);
+  if (result[0] == result[1])
+    result[1] = not_any_of(result[0], result[2]);
+  if (result[0] == result[2])
+    result[2] = not_any_of(result[0], result[1]);
+  if (result[1] == result[2])
+    result[2] = not_any_of(result[0], result[1]);
+  assert(result[0] != result[1] && result[1] != result[2] && result[2] != result[0]);
 
-  // Figure out whether any of the rows of the transform point in the
-  //   opposite direction to the MRtrix convention
-  flip[perm[0]] = T(0, perm[0]) < 0.0;
-  flip[perm[1]] = T(1, perm[1]) < 0.0;
-  flip[perm[2]] = T(2, perm[2]) < 0.0;
+  return result;
 }
 
 } // namespace Axes

core/axes.h

@@ -34,7 +34,10 @@ Eigen::Vector3d id2dir(const std::string &);
 
 //! determine the axis permutations and flips necessary to make an image
 //!   appear approximately axial
-void get_permutation_to_make_axial(const transform_type &T, std::array<size_t, 3> &perm, std::array<bool, 3> &flip);
+void get_shuffle_to_make_axial(const transform_type &T, std::array<size_t, 3> &perm, std::array<bool, 3> &flip);
+
+//! determine which vectors of a 3x3 transform are closest to the three axis indices
+std::array<size_t, 3> closest(const Eigen::Matrix3d &M);
 
 } // namespace Axes
 } // namespace MR

core/header.cpp

@@ -611,7 +611,7 @@ void Header::sanitise_transform() {
 
 void Header::realign_transform() {
   // find which row of the transform is closest to each scanner axis:
-  Axes::get_permutation_to_make_axial(transform(), realign_perm_, realign_flip_);
+  Axes::get_shuffle_to_make_axial(transform(), realign_perm_, realign_flip_);
 
   // check if image is already near-axial, return if true:
   if (realign_perm_[0] == 0 && realign_perm_[1] == 1 && realign_perm_[2] == 2 && !realign_flip_[0] &&

core/math/average_space.cpp

@@ -15,6 +15,7 @@
  */
 
 #include "math/average_space.h"
+#include "axes.h"
 
 namespace MR {
 namespace Math {
@@ -57,13 +58,24 @@ double matrix_average(std::vector<Eigen::MatrixXd> const &mat_in, Eigen::MatrixX
 
 namespace MR {
 
-FORCE_INLINE Eigen::Matrix<default_type, 8, 4> get_cuboid_corners(const Eigen::Matrix<default_type, 4, 1> &xyz1) {
+const char *const avgspace_voxspacing_choices[]{
+    "min_projection", "mean_projection", "min_nearest", "mean_nearest", nullptr};
+
+FORCE_INLINE Eigen::Matrix<default_type, 8, 4> get_cuboid_corners(const Eigen::Matrix<default_type, 4, 1> &xyz_sizes) {
   Eigen::Matrix<default_type, 8, 4> corners;
   // MatrixType faces = MatrixType(8,4);
   // faces << 1, 2, 3, 4,   2, 6, 7, 3,   4, 3, 7, 8,   1, 5, 8, 4,   1, 2, 6, 5,   5, 6, 7, 8;
-  corners << 0.0, 0.0, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0, 1.0, 1.0, 0.0, 1.0, 1.0, 0.0, 0.0, 1.0, 0.0, 0.0, 1.0, 1.0, 0.0,
-      1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 0.0, 1.0, 1.0;
-  corners *= xyz1.asDiagonal();
+  // clang-format off
+  corners << 0.0, 0.0, 0.0, 1.0,
+             0.0, 1.0, 0.0, 1.0,
+             1.0, 1.0, 0.0, 1.0,
+             1.0, 0.0, 0.0, 1.0,
+             0.0, 0.0, 1.0, 1.0,
+             0.0, 1.0, 1.0, 1.0,
+             1.0, 1.0, 1.0, 1.0,
+             1.0, 0.0, 1.0, 1.0;
+  // clang-format on
+  corners *= xyz_sizes.asDiagonal();
   return corners;
 }
 
@@ -163,7 +175,7 @@ void compute_average_voxel2scanner(
     const std::vector<Header> &input_headers,
     const Eigen::Matrix<default_type, 4, 1> &padding,
     const std::vector<Eigen::Transform<default_type, 3, Eigen::Projective>> &transform_header_with,
-    int voxel_subsampling) {
+    const avgspace_voxspacing_t voxel_spacing_calculation) {
   const size_t num_images = input_headers.size();
   DEBUG("compute_average_voxel2scanner num_images:" + str(num_images));
   std::vector<Eigen::Transform<default_type, 3, Eigen::Projective>> transformation_matrices;
@@ -200,38 +212,41 @@ void compute_average_voxel2scanner(
   for (size_t itrafo = 0; itrafo < num_images; itrafo++) {
     Eigen::MatrixXd Other = ScannerSpaceAxis3 * transformation_matrices[itrafo].rotation().transpose();
     Eigen::Quaterniond quat = rot_match_coordinates(ScannerSpaceAxis6, Other);
-    quaternions.col(itrafo) =
-        quat.coeffs(); // internally the coefficients are stored in the following order: [x, y, z, w]
+    // internally the coefficients are stored in the following order: [x, y, z, w]
+    quaternions.col(itrafo) = quat.coeffs();
   }
-  Eigen::SelfAdjointEigenSolver<Eigen::MatrixXd> es(quaternions * quaternions.transpose(), Eigen::ComputeEigenvectors);
+  const Eigen::SelfAdjointEigenSolver<Eigen::MatrixXd> es(quaternions * quaternions.transpose(),
+                                                          Eigen::ComputeEigenvectors);
   Eigen::Quaterniond average_quat;
   average_quat.coeffs() = es.eigenvectors().col(3).transpose();
   average_quat.normalize();
-  const Eigen::Matrix3d Raverage(average_quat.toRotationMatrix());
+  const Eigen::Matrix3d Raverage(average_quat.toRotationMatrix().transpose());
   Eigen::MatrixXd rot_vox_size = Eigen::MatrixXd(3, num_images);
   for (size_t itrafo = 0; itrafo < num_images; itrafo++) {
-    Eigen::MatrixXd M = (Raverage * transformation_matrices[itrafo].linear()).cwiseAbs();
-    if (voxel_subsampling == 0)
-      rot_vox_size.col(itrafo) = (Raverage * transformation_matrices[itrafo].linear()).cwiseAbs().diagonal();
-    else
-      rot_vox_size.col(itrafo) = (Raverage * transformation_matrices[itrafo].linear()).cwiseAbs().colwise().sum();
+    const Eigen::Matrix3d M = Raverage.transpose() * transformation_matrices[itrafo].linear();
+    switch (voxel_spacing_calculation) {
+    case avgspace_voxspacing_t::MIN_PROJECTION:
+      rot_vox_size.col(itrafo) = M.cwiseAbs().diagonal();
+      break;
+    case avgspace_voxspacing_t::MEAN_PROJECTION:
+      rot_vox_size.col(itrafo) = M.cwiseAbs().colwise().sum();
+      break;
+    case avgspace_voxspacing_t::MIN_NEAREST:
+    case avgspace_voxspacing_t::MEAN_NEAREST: {
+      std::array<size_t, 3> perm = Axes::closest(M);
+      for (size_t axis = 0; axis != 3; ++axis)
+        rot_vox_size(axis, itrafo) = input_headers[itrafo].spacing(perm[axis]);
+    } break;
+    }
   }
 
-  switch (voxel_subsampling) {
-  case 0: // maximum voxel size determined by minimum of all projected input image voxel sizes, sampling at or above
-          // Nyquist limit
-    projected_voxel_sizes = rot_vox_size.rowwise().minCoeff();
-    break;
-  case 1: // maximum voxel size determined by mean of all projected input image voxel sizes
-    projected_voxel_sizes = rot_vox_size.rowwise().mean();
-    break;
-  default:
-    assert(0 && "compute_average_voxel2scanner: invalid voxel_subsampling option");
-    break;
-  }
+  projected_voxel_sizes = (voxel_spacing_calculation == avgspace_voxspacing_t::MIN_PROJECTION ||
+                           voxel_spacing_calculation == avgspace_voxspacing_t::MIN_NEAREST)
+                              ? Eigen::Vector3d(rot_vox_size.rowwise().minCoeff())
+                              : Eigen::Vector3d(rot_vox_size.rowwise().mean());
 
   Eigen::MatrixXd mat_avg = Eigen::MatrixXd::Zero(4, 4);
-  mat_avg.block<3, 3>(0, 0) = Raverage.transpose();
+  mat_avg.block<3, 3>(0, 0) = Raverage;
   mat_avg.col(0) *= projected_voxel_sizes(0);
   mat_avg.col(1) *= projected_voxel_sizes(1);
   mat_avg.col(2) *= projected_voxel_sizes(2);
@@ -265,7 +280,7 @@ void compute_average_voxel2scanner(
 Header compute_minimum_average_header(
     const std::vector<Header> &input_headers,
     const std::vector<Eigen::Transform<default_type, 3, Eigen::Projective>> &transform_header_with,
-    int voxel_subsampling,
+    const avgspace_voxspacing_t voxel_spacing_calculation,
     Eigen::Matrix<default_type, 4, 1> padding) {
   Eigen::Transform<default_type, 3, Eigen::Projective> average_v2s_trafo;
   Eigen::Vector3d average_space_voxel_extent;
@@ -276,7 +291,7 @@ Header compute_minimum_average_header(
                                 input_headers,
                                 padding,
                                 transform_header_with,
-                                voxel_subsampling);
+                                voxel_spacing_calculation);
 
   // create average space header
   Header header_out;

core/math/average_space.h

@@ -32,14 +32,17 @@ double matrix_average(std::vector<Eigen::MatrixXd> const &mat_in, Eigen::MatrixX
 
 namespace MR {
 
-Eigen::Matrix<default_type, 8, 4> get_cuboid_corners(const Eigen::Matrix<default_type, 4, 1> &xzx1);
+extern const char *const avgspace_voxspacing_choices[];
+enum class avgspace_voxspacing_t { MIN_PROJECTION, MEAN_PROJECTION, MIN_NEAREST, MEAN_NEAREST };
+
+Eigen::Matrix<default_type, 8, 4> get_cuboid_corners(const Eigen::Matrix<default_type, 4, 1> &xyz_sizes);
 Eigen::Matrix<default_type, 8, 4>
 get_bounding_box(const Header &header, const Eigen::Transform<default_type, 3, Eigen::Projective> &voxel2scanner);
 
 Header compute_minimum_average_header(
     const std::vector<Header> &input_headers,
     const std::vector<Eigen::Transform<default_type, 3, Eigen::Projective>> &transform_header_with,
-    int voxel_subsampling = 1,
+    avgspace_voxspacing_t voxel_spacing_calculation = avgspace_voxspacing_t::MEAN_NEAREST,
     Eigen::Matrix<default_type, 4, 1> padding = Eigen::Matrix<default_type, 4, 1>(1.0, 1.0, 1.0, 1.0));
 
 template <class ImageType1, class ImageType2>
@@ -51,10 +54,10 @@ Header compute_minimum_average_header(
     Eigen::Transform<default_type, 3, Eigen::Projective> transform_2 =
         Eigen::Transform<default_type, 3, Eigen::Projective>::Identity(),
     Eigen::Matrix<default_type, 4, 1> padding = Eigen::Matrix<default_type, 4, 1>(1.0, 1.0, 1.0, 1.0),
-    int voxel_subsampling = 1) {
+    const avgspace_voxspacing_t voxel_spacing_calculation = avgspace_voxspacing_t::MEAN_NEAREST) {
   std::vector<Eigen::Transform<default_type, 3, Eigen::Projective>> init_transforms{transform_1, transform_2};
   std::vector<Header> headers{im1, im2};
-  return compute_minimum_average_header(headers, init_transforms, voxel_subsampling, padding);
+  return compute_minimum_average_header(headers, init_transforms, voxel_spacing_calculation, padding);
 }
 } // namespace MR
 #endif

docs/reference/commands/dwicat.rst

@@ -21,13 +21,19 @@ Usage
 Description
 -----------
 
-This script concatenates two or more 4D DWI series, accounting for the fact that there may be differences in intensity scaling between those series. This intensity scaling is corrected by determining scaling factors that will make the overall image intensities in the b=0 volumes of each series approximately equivalent.
+This script concatenates two or more 4D DWI series, accounting for various detrimental confounds that may affect such an operation. Each of those confounds are described in separate paragraphs below.
+
+There may be differences in intensity scaling between the input series. This intensity scaling is corrected by determining scaling factors that will make the overall image intensities in the b=0 volumes of each series approximately equivalent. This operation is only appropriate if the sequence parameters that influence the contrast of the b=0 image volumes are identical.
+
+Concatenation of DWI series defined on different voxel grids. If the voxel grids of the input DWI series are not precisely identical, then it may not be possible to do a simple concatenation operation. In this scenario the script will determine the appropriate way to combine the input series, ideally only manipulating header transformations and avoiding image interpolation if possible.
 
 Options
 -------
 
 - **-mask image** Provide a binary mask within which image intensities will be matched
 
+- **-nointensity** Do not perform intensity matching based on b=0 volumes
+
 Additional standard options for Python scripts
 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
 

docs/reference/commands/mraverageheader.rst

@@ -18,12 +18,17 @@ Usage
 -  *input*: the input image(s).
 -  *output*: the output image
 
+Description
+-----------
+
+The voxel spacings of the calculated average header can be determined from the set of input images in one of four different ways, which may be more or less appropriate in different use cases. These options vary in two key ways: 1. Projected voxel spacing of the input image in the direction of the average header axes versus the voxel spacing of the input image axis that is closest to the average header axis; 2. Selecting the minimum of these spacings across input images to maintain maximal spatial resolution versus the mean across images to produce an unbiased average.
+
 Options
 -------
 
--  **-padding value**  boundary box padding in voxels. Default: 0
+-  **-padding value** boundary box padding in voxels. Default: 0
 
--  **-resolution type** subsampling of template compared to smallest voxel size in any input image.Valid options are: - 'mean': unbiased but loss of resolution for individual images possible; - 'max': smallest voxel size of any input image defines the resolution. Default: mean
+-  **-spacing type** Method for determination of voxel spacings based on the set of input images and the average header axes (see Description).Valid options are: min_projection,mean_projection,min_nearest,mean_nearest; default = mean_nearest
 
 -  **-fill** set the intensity in the first volume of the average space to 1
 

docs/reference/commands/transformcalc.rst

@@ -36,7 +36,7 @@ Example usages
 
 -   *Calculate the transformation matrix from an original image and an image with modified header*::
 
-        $ transformcalc mov mapmovhdr header output
+        $ transformcalc orig_image modified_image header output
 
 -   *Calculate the average affine matrix of a set of input matrices*::