dwi2adc: Change IVIM interface

cpp/cmd/dwi2adc.cpp

@@ -25,6 +25,10 @@
 using namespace MR;
 using namespace App;
 
+using value_type = float;
+
+constexpr value_type ivim_cutoff_default = value_type(120);
+
 // clang-format off
 void usage ()
 {
@@ -33,137 +37,175 @@ void usage ()
   SYNOPSIS = "Calculate ADC and/or IVIM parameters.";
 
   DESCRIPTION
-  + "By default, the command will estimate the Apparent Diffusion Coefficient (ADC) "
-    "using the isotropic mono-exponential model: S(b) = S(0) * exp(-D * b). "
-    "The output consists of 2 volumes, respectively S(0) and D."
+  + "The command estimates the Apparent Diffusion Coefficient (ADC) "
+    "using the isotropic mono-exponential model: S(b) = S(0) * exp(-ADC * b). "
+    "The value of S(0) can be optionally exported using command-line option -szero."
 
   + "When using the -ivim option, the command will additionally estimate the "
     "Intra-Voxel Incoherent Motion (IVIM) parameters f and D', i.e., the perfusion "
     "fraction and the pseudo-diffusion coefficient. IVIM assumes a bi-exponential "
     "model: S(b) = S(0) * ((1-f) * exp(-D * b) + f * exp(-D' * b)). This command "
     "adopts a 2-stage fitting strategy, in which the ADC is first estimated based on "
-    "the DWI data with b > cutoff, and the other parameters are estimated subsequently. "
-    "The output consists of 4 volumes, respectively S(0), D, f, and D'."
-
-  + "Note that this command ignores the gradient orientation entirely. This approach is "
-    "therefore only suited for mean DWI (trace-weighted) images or for DWI data collected "
-    "with isotropically-distributed gradient directions.";
+    "the DWI data with b > cutoff, and the other parameters are estimated subsequently."
+
+  + "Note that this command ignores the gradient orientation entirely. "
+    "If a conventional DWI series is provided as input, "
+    "all volumes will contribute equally toward the model fit "
+    "irrespective of direction of diffusion sensitisation; "
+    "DWI data should therefore ideally consist of isotropically-distributed gradient directions.";
+    "The approach can alternative be applied to mean DWI (trace-weighted) images.";
 
   ARGUMENTS
-    + Argument ("input", "the input image.").type_image_in ()
-    + Argument ("output", "the output image.").type_image_out ();
+    + Argument ("input", "the input image").type_image_in()
+    + Argument ("output", "the output ADC image").type_image_out();
 
   OPTIONS
-    + Option ("ivim", "also estimate IVIM parameters in 2-stage fit.")
+    + Option("szero",
+             "export image of S(0); "
+             "that is, the model-estimated signal intensity in the absence of diffusion weighting")
+      + Argument("image").type_image_out()
+
+    + Option ("ivim", "also estimate IVIM parameters in 2-stage fit, "
+                      "yielding two images encoding signal fraction and diffusivity respectively of perfusion1 component")
+      + Argument("fraction").type_image_out()
+      + Argument("diffusivity").type_image_out()
 
-    + Option ("cutoff", "minimum b-value for ADC estimation in IVIM fit (default = 120 s/mm^2).")
-    +   Argument ("bval").type_integer (0, 1000)
+    + Option ("cutoff", "minimum b-value for ADC estimation in IVIM fit "
+                        "(default = " + str(ivim_cutoff_default) + " s/mm^2).")
+    +   Argument ("bval").type_float (0.0, 1000.0)
 
     + DWI::GradImportOptions();
-  
+
   REFERENCES
   + "Le Bihan, D.; Breton, E.; Lallemand, D.; Aubin, M.L.; Vignaud, J.; Laval-Jeantet, M. "
     "Separation of diffusion and perfusion in intravoxel incoherent motion MR imaging. "
     "Radiology, 1988, 168, 497–505."
-    
+
   + "Jalnefjord, O.; Andersson, M.; Montelius; M.; Starck, G.; Elf, A.; Johanson, V.; Svensson, J.; Ljungberg, M. "
     "Comparison of methods for estimation of the intravoxel incoherent motion (IVIM) "
     "diffusion coefficient (D) and perfusion fraction (f). "
     "Magn Reson Mater Phy, 2018, 31, 715–723.";
 }
 // clang-format on
 
-using value_type = float;
-
 class DWI2ADC {
 public:
-  DWI2ADC(const Eigen::VectorXd &bvals, size_t dwi_axis, bool ivim, int cutoff)
-      : bvals(bvals), dwi(bvals.size()), adc(2), dwi_axis(dwi_axis), ivim(ivim), cutoff(cutoff) {
-    Eigen::MatrixXd b(bvals.size(), 2);
+  DWI2ADC(const Header &header, const Eigen::VectorXd &bvals, const size_t dwi_axis)
+      : H(header),
+        bvals(bvals),
+        dwi(bvals.size()),
+        b(bvals.size(), 2),
+        logszero_and_adc(2),
+        dwi_axis(dwi_axis),
+        ivim_cutoff(std::numeric_limits<value_type>::signaling_NaN()) {
     for (ssize_t i = 0; i < b.rows(); ++i) {
       b(i, 0) = 1.0;
       b(i, 1) = -bvals(i);
     }
     binv = Math::pinv(b);
-    if (ivim) {
-      // select volumes with b-value > cutoff
-      for (ssize_t j = 0; j < bvals.size(); j++) {
-        if (bvals[j] > cutoff)
-          idx.push_back(j);
-      }
-      const Eigen::MatrixXd bsub = b(idx, Eigen::all);
-      bsubinv = Math::pinv(bsub);
+  }
+
+  void set_bzero_path(const std::string &path) { szero_image = Image<value_type>::create(path, H); }
+
+  void initialise_ivim(const std::string &ivimfrac_path, const std::string &ivimdiff_path, const value_type cutoff) {
+    ivimfrac_image = Image<value_type>::create(ivimfrac_path, H);
+    ivimdiff_image = Image<value_type>::create(ivimdiff_path, H);
+    ivim_cutoff = cutoff;
+    // select volumes with b-value > cutoff
+    for (ssize_t j = 0; j < bvals.size(); j++) {
+      if (bvals[j] > cutoff)
+        ivim_suprathresh_idx.push_back(j);
     }
+    const Eigen::MatrixXd bsub = b(ivim_suprathresh_idx, Eigen::all);
+    bsubinv = Math::pinv(bsub);
   }
 
-  template <class DWIType, class ADCType> void operator()(DWIType &dwi_image, ADCType &adc_image) {
+  void operator()(Image<value_type> &dwi_image, Image<value_type> &adc_image) {
     for (auto l = Loop(dwi_axis)(dwi_image); l; ++l) {
-      value_type val = dwi_image.value();
+      const value_type val = dwi_image.value();
       dwi[dwi_image.index(dwi_axis)] = val > 1.0e-12 ? std::log(val) : 1.0e-12;
     }
 
-    if (ivim) {
-      dwisub = dwi(idx);
-      adc = bsubinv * dwisub;
+    if (std::isnan(ivim_cutoff)) {
+      logszero_and_adc = binv * dwi;
     } else {
-      adc = binv * dwi;
+      dwisub = dwi(ivim_suprathresh_idx);
+      logszero_and_adc = bsubinv * dwisub;
+    }
+
+    adc_image.value() = logszero_and_adc[1];
+
+    if (std::isnan(ivim_cutoff)) {
+      if (szero_image.valid()) {
+        assign_pos_of(adc_image).to(szero_image);
+        szero_image.value() = std::exp(logszero_and_adc[0]);
+      }
+      return;
     }
 
-    adc_image.index(3) = 0;
-    adc_image.value() = std::exp(adc[0]);
-    adc_image.index(3) = 1;
-    adc_image.value() = adc[1];
-
-    if (ivim) {
-      const double A = std::exp(adc[0]);
-      const double D = adc[1];
-      Eigen::VectorXd logS = adc[0] - D * bvals.array();
-      Eigen::VectorXd logdiff = (dwi.array() > logS.array()).select(dwi, logS);
-      logdiff.array() += Eigen::log(1 - Eigen::exp(-(dwi - logS).array().abs()));
-      adc = binv * logdiff;
-      const double C = std::exp(adc[0]);
-      const double Dstar = adc[1];
-      const double S0 = A + C;
-      const double f = C / S0;
-      adc_image.index(3) = 0;
-      adc_image.value() = S0;
-      adc_image.index(3) = 2;
-      adc_image.value() = f;
-      adc_image.index(3) = 3;
-      adc_image.value() = Dstar;
+    const double A = std::exp(logszero_and_adc[0]);
+    const double D = logszero_and_adc[1];
+    const Eigen::VectorXd logS = logszero_and_adc[0] - D * bvals.array();
+    Eigen::VectorXd logdiff = (dwi.array() > logS.array()).select(dwi, logS);
+    logdiff.array() += Eigen::log(1 - Eigen::exp(-(dwi - logS).array().abs()));
+    logszero_and_adc = binv * logdiff;
+    const double C = std::exp(logszero_and_adc[0]);
+    const double Dstar = logszero_and_adc[1];
+    const double S0 = A + C;
+    const double f = C / S0;
+    if (szero_image.valid()) {
+      assign_pos_of(adc_image).to(szero_image);
+      szero_image.value() = S0;
     }
+    assign_pos_of(adc_image).to(ivimfrac_image, ivimdiff_image);
+    ivimfrac_image.value() = f;
+    ivimdiff_image.value() = Dstar;
   }
 
 private:
-  Eigen::VectorXd bvals, dwi, dwisub, adc;
-  Eigen::MatrixXd binv, bsubinv;
-  std::vector<size_t> idx;
+  const Header H;
+  Eigen::VectorXd bvals, dwi, dwisub, logszero_and_adc;
+  Eigen::MatrixXd b, binv;
   const size_t dwi_axis;
-  const bool ivim;
-  const int cutoff;
+
+  // Optional export
+  Image<value_type> szero_image;
+
+  // Members relating to IVIM operation
+  std::vector<size_t> ivim_suprathresh_idx;
+  Eigen::MatrixXd bsubinv;
+  Image<value_type> ivimfrac_image;
+  Image<value_type> ivimdiff_image;
+  value_type ivim_cutoff;
 };
 
 void run() {
-  auto dwi = Header::open(argument[0]).get_image<value_type>();
-  auto grad = DWI::get_DW_scheme(dwi);
+  auto H_in = Header::open(argument[0]);
+  auto grad = DWI::get_DW_scheme(H_in);
 
   size_t dwi_axis = 3;
-  while (dwi.size(dwi_axis) < 2)
+  while (H_in.size(dwi_axis) < 2)
     ++dwi_axis;
   INFO("assuming DW images are stored along axis " + str(dwi_axis));
 
-  auto opt = get_options("ivim");
-  bool ivim = opt.size();
-  int bmin = get_option_value("cutoff", 120);
+  Header H_out(H_in);
+  H_out.datatype() = DataType::Float32;
+  H_out.datatype().set_byte_order_native();
+  H_out.ndim() = 3;
+  DWI::stash_DW_scheme(H_out, grad);
+  PhaseEncoding::clear_scheme(H_out);
+
+  DWI2ADC functor(H_out, grad.col(3), dwi_axis);
 
-  Header header(dwi);
-  header.datatype() = DataType::Float32;
-  DWI::stash_DW_scheme(header, grad);
-  PhaseEncoding::clear_scheme(header);
-  header.ndim() = 4;
-  header.size(3) = ivim ? 4 : 2;
+  auto opt = get_options("szero");
+  if (opt.size())
+    functor.set_bzero_path(opt[0][0]);
+  opt = get_options("ivim");
+  if (opt.size())
+    functor.initialise_ivim(opt[0][0], opt[0][1], get_option_value("cutoff", ivim_cutoff_default));
 
-  auto adc = Image<value_type>::create(argument[1], header);
+  auto dwi = H_in.get_image<value_type>();
+  auto adc = Image<value_type>::create(argument[1], H_out);
 
-  ThreadedLoop("computing ADC values", dwi, 0, 3).run(DWI2ADC(grad.col(3), dwi_axis, ivim, bmin), dwi, adc);
+  ThreadedLoop("computing ADC values", H_out).run(functor, dwi, adc);
 }

docs/reference/commands/dwi2adc.rst

@@ -15,22 +15,24 @@ Usage
 
     dwi2adc [ options ]  input output
 
--  *input*: the input image.
--  *output*: the output image.
+-  *input*: the input image
+-  *output*: the output ADC image
 
 Description
 -----------
 
-By default, the command will estimate the Apparent Diffusion Coefficient (ADC) using the isotropic mono-exponential model: S(b) = S(0) * exp(-D * b). The output consists of 2 volumes, respectively S(0) and D.
+The command estimates the Apparent Diffusion Coefficient (ADC) using the isotropic mono-exponential model: S(b) = S(0) * exp(-ADC * b). The value of S(0) can be optionally exported using command-line option -szero.
 
-When using the -ivim option, the command will additionally estimate the Intra-Voxel Incoherent Motion (IVIM) parameters f and D', i.e., the perfusion fraction and the pseudo-diffusion coefficient. IVIM assumes a bi-exponential model: S(b) = S(0) * ((1-f) * exp(-D * b) + f * exp(-D' * b)). This command adopts a 2-stage fitting strategy, in which the ADC is first estimated based on the DWI data with b > cutoff, and the other parameters are estimated subsequently. The output consists of 4 volumes, respectively S(0), D, f, and D'.
+When using the -ivim option, the command will additionally estimate the Intra-Voxel Incoherent Motion (IVIM) parameters f and D', i.e., the perfusion fraction and the pseudo-diffusion coefficient. IVIM assumes a bi-exponential model: S(b) = S(0) * ((1-f) * exp(-D * b) + f * exp(-D' * b)). This command adopts a 2-stage fitting strategy, in which the ADC is first estimated based on the DWI data with b > cutoff, and the other parameters are estimated subsequently.
 
-Note that this command ignores the gradient orientation entirely. This approach is therefore only suited for mean DWI (trace-weighted) images or for DWI data collected with isotropically-distributed gradient directions.
+Note that this command ignores the gradient orientation entirely. If a conventional DWI series is provided as input, all volumes will contribute equally toward the model fit irrespective of direction of diffusion sensitisation; DWI data should therefore ideally consist of isotropically-distributed gradient directions.
 
 Options
 -------
 
--  **-ivim** also estimate IVIM parameters in 2-stage fit.
+-  **-szero image** export image of S(0); that is, the model-estimated signal intensity in the absence of diffusion weighting
+
+-  **-ivim fraction diffusivity** also estimate IVIM parameters in 2-stage fit, yielding two images encoding signal fraction and diffusivity respectively of perfusion1 component
 
 -  **-cutoff bval** minimum b-value for ADC estimation in IVIM fit (default = 120 s/mm^2).
 

testing/binaries/tests/dwi2adc/default

@@ -1,5 +1,18 @@
 #!/bin/bash
 # Basic test of default operation of dwi2adc command
-# COmpares output to that generated from a prior software version
-dwi2adc dwi2adc/in.mif - | \
+# Compares output to that generated from a prior software version
+#
+# Note that in the pre-generated reference output test data,
+#   dwi2adc created a single output image,
+#   where the first volume was S(0) and the second volume was ADC.
+# The updated interface involves writing a compulsory ADC image,
+#   and an optional S(0) volume written to a different image
+#   using explicit command-line option -szero.
+# In order to test the output of this command against the pre-existing reference
+#   without necessitating the re-generation of reference test data,
+#   these are explicitly concatenated before comparing to the reference data
+
+dwi2adc dwi2adc/in.mif tmp_adc.mif -szero tmp_szero.mif -force
+
+mrcat tmp_szero.mif tmp_adc.mif -axis 3 - | \
 testing_diff_image - dwi2adc/out.mif -frac 1e-5