Merge pull request MRtrix3#1507 from MRtrix3/streamline_length_precise

Precise quantification of streamlines lengths

cmd/tckgen.cpp

@@ -55,7 +55,7 @@ void usage ()
 
   SYNOPSIS = "Perform streamlines tractography";
 
-  DESCRIPTION 
+  DESCRIPTION
     + "By default, tckgen produces a fixed number of streamlines, by attempting "
       "to seed from new random locations until the target number of "
       "streamlines have been selected (in other words, after all inclusion & "
@@ -129,9 +129,15 @@ void usage ()
       "voxel for each streamline. These data are then sampled via trilinear "
       "interpolation at each streamline step, the diffusion tensor model is fitted, "
       "and the streamline follows the orientation of the principal eigenvector of "
-      "that tensor.";
+      "that tensor."
 
-  REFERENCES 
+    + "Note that the behaviour of the -angle option varies slightly depending on the "
+      "order of integration: for any first-order method, this angle corresponds to the "
+      "deviation in streamline trajectory per step; for higher-order methods, this "
+      "corresponds to the change in underlying fibre orientation between the start and "
+      "end points of each step.";
+
+  REFERENCES
    + "References based on streamlines algorithm used:"
 
    + "* FACT:\n"

cmd/tckresample.cpp

@@ -45,9 +45,18 @@ void usage ()
   SYNOPSIS = "Resample each streamline in a track file to a new set of vertices";
 
   DESCRIPTION
-  + "This may be either increasing or decreasing the number of samples along "
-    "each streamline, or changing the positions of the samples according to "
-    "some specified trajectory."
+  + "It is necessary to specify precisely ONE of the command-line options for "
+    "controlling how this resampling takes place; this may be either increasing "
+    "or decreasing the number of samples along each streamline, or may involve "
+    "changing the positions of the samples according to some specified trajectory."
+
+  + "Note that because the length of a streamline is calculated based on the sums of "
+    "distances between adjacent vertices, resampling a streamline to a new set of "
+    "vertices will typically change the quantified length of that streamline; the "
+    "magnitude of the difference will typically depend on the discrepancy in the "
+    "number of vertices, with less vertices leading to a shorter length (due to "
+    "taking chordal lengths of curved trajectories)."
+
 
   + DWI::Tractography::preserve_track_order_desc;
 

cmd/tckstats.cpp

@@ -48,7 +48,7 @@ void usage ()
 
   AUTHOR = "Robert E. Smith ([email protected])";
 
-  SYNOPSIS = "Calculate statistics on streamlines length";
+  SYNOPSIS = "Calculate statistics on streamlines lengths";
 
   ARGUMENTS
   + Argument ("tracks_in", "the input track file").type_tracks_in();
@@ -66,9 +66,6 @@ void usage ()
   + Option ("dump", "dump the streamlines lengths to a text file")
     + Argument ("path").type_file_out()
 
-  + Option ("explicit", "explicitly calculate the length of each streamline, "
-                        "ignoring any step size information present in the header")
-
   + Option ("ignorezero", "do not generate a warning if the track file contains streamlines with zero length")
 
   + Tractography::TrackWeightsInOption;
@@ -111,8 +108,9 @@ void run ()
   size_t count = 0, header_count = 0;
   float min_length = std::numeric_limits<float>::infinity();
   float max_length = -std::numeric_limits<float>::infinity();
-  double sum_lengths = 0.0, sum_weights = 0.0;
-  vector<double> histogram;
+  size_t empty_streamlines = 0, zero_length_streamlines = 0;
+  default_type sum_lengths = 0.0, sum_weights = 0.0;
+  vector<default_type> histogram;
   vector<LW> all_lengths;
   all_lengths.reserve (header_count);
 
@@ -123,14 +121,9 @@ void run ()
     if (properties.find ("count") != properties.end())
       header_count = to<size_t> (properties["count"]);
 
-    if (!get_options ("explicit").size()) {
-      step_size = get_step_size (properties);
-      if (!std::isfinite (step_size) || !step_size) {
-        INFO ("Streamline step size undefined in header; lengths will be calculated manually");
-        if (get_options ("histogram").size()) {
-          WARN ("Do not have streamline step size with which to construct histogram; histogram will be generated using 1mm bin widths");
-        }
-      }
+    step_size = get_step_size (properties);
+    if ((!std::isfinite (step_size) || !step_size) && get_options ("histogram").size()) {
+      WARN ("Do not have streamline step size with which to bin histogram; histogram will be generated using 1mm bin widths");
     }
 
     std::unique_ptr<File::OFStream> dump;
@@ -142,7 +135,7 @@ void run ()
     Streamline<> tck;
     while (reader (tck)) {
       ++count;
-      const float length = std::isfinite (step_size) ? tck.calc_length (step_size) : tck.calc_length();
+      const float length = Tractography::length (tck);
       if (std::isfinite (length)) {
         min_length = std::min (min_length, length);
         max_length = std::max (max_length, length);
@@ -153,15 +146,28 @@ void run ()
         while (histogram.size() <= index)
           histogram.push_back (0.0);
         histogram[index] += tck.weight;
+        if (!length)
+          ++zero_length_streamlines;
+      } else {
+        ++empty_streamlines;
       }
       if (dump)
         (*dump) << length << "\n";
       ++progress;
     }
   }
 
-  if (histogram.size() && histogram.front() && !get_options ("ignorezero").size())
-    WARN ("read " + str(histogram.front()) + " zero-length tracks");
+  if (!get_options ("ignorezero").size() && (empty_streamlines || zero_length_streamlines)) {
+    std::string s ("read");
+    if (empty_streamlines) {
+      s += " " + str(empty_streamlines) + " empty streamlines";
+      if (zero_length_streamlines)
+        s += " and";
+    }
+    if (zero_length_streamlines)
+      s += " " + str(zero_length_streamlines) + " streamlines with zero length (one vertex only)";
+    WARN (s);
+  }
   if (count != header_count)
     WARN ("expected " + str(header_count) + " tracks according to header; read " + str(count));
   if (!std::isfinite (min_length))
@@ -177,7 +183,7 @@ void run ()
       // Perform a weighted median calculation
       std::sort (all_lengths.begin(), all_lengths.end());
       size_t median_index = 0;
-      double sum = sum_weights - all_lengths[0].get_weight();
+      default_type sum = sum_weights - all_lengths[0].get_weight();
       while (sum > 0.5 * sum_weights) { sum -= all_lengths[++median_index].get_weight(); }
       median_length = all_lengths[median_index].get_length();
     } else {
@@ -187,10 +193,10 @@ void run ()
     median_length = NaN;
   }
 
-  double stdev = 0.0;
+  default_type ssd = 0.0;
   for (vector<LW>::const_iterator i = all_lengths.begin(); i != all_lengths.end(); ++i)
-    stdev += i->get_weight() * Math::pow2 (i->get_length() - mean_length);
-  stdev = sum_weights ? (std::sqrt (stdev / (((count - 1) / float(count)) * sum_weights))) : NaN;
+    ssd += i->get_weight() * Math::pow2 (i->get_length() - mean_length);
+  const float stdev = sum_weights ? (std::sqrt (ssd / (((count - 1) / default_type(count)) * sum_weights))) : NaN;
 
   vector<std::string> fields;
   auto opt = get_options ("output");

docs/reference/commands/tckgen.rst

@@ -41,6 +41,8 @@ Below is a list of available tracking algorithms, the input image data that they
 
 - Tensor_Prob: A probabilistic algorithm that takes as input a 4D diffusion-weighted image (DWI) series. Within each image voxel, a residual bootstrap is performed to obtain a unique realisation of the DWI data in that voxel for each streamline. These data are then sampled via trilinear interpolation at each streamline step, the diffusion tensor model is fitted, and the streamline follows the orientation of the principal eigenvector of that tensor.
 
+Note that the behaviour of the -angle option varies slightly depending on the order of integration: for any first-order method, this angle corresponds to the deviation in streamline trajectory per step; for higher-order methods, this corresponds to the change in underlying fibre orientation between the start and end points of each step.
+
 Options
 -------
 
@@ -53,7 +55,7 @@ Streamlines tractography options
 
 -  **-step size** set the step size of the algorithm in mm (default is 0.1 x voxelsize; for iFOD2: 0.5 x voxelsize).
 
--  **-angle theta** set the maximum angle between successive steps (default is 90deg x stepsize / voxelsize).
+-  **-angle theta** set the maximum angle between successive steps: default is chosen such that minimum radius of curvature is one voxel; interpretation depends on order of integration (see Description).
 
 -  **-minlength value** set the minimum length of any track in mm (default is 5 x voxelsize without ACT, 2 x voxelsize with ACT).
 

docs/reference/commands/tckresample.rst

@@ -21,7 +21,9 @@ Usage
 Description
 -----------
 
-This may be either increasing or decreasing the number of samples along each streamline, or changing the positions of the samples according to some specified trajectory.
+It is necessary to specify precisely ONE of the command-line options for controlling how this resampling takes place; this may be either increasing or decreasing the number of samples along each streamline, or may involve changing the positions of the samples according to some specified trajectory.
+
+Note that because the length of a streamline is calculated based on the sums of distances between adjacent vertices, resampling a streamline to a new set of vertices will typically change the quantified length of that streamline; the magnitude of the difference will typically depend on the discrepancy in the number of vertices, with less vertices leading to a shorter length (due to taking chordal lengths of curved trajectories).
 
 Note that if multi-threading is used in this command, the ordering of tracks in the output file is unlikely to match the order of the incoming data. If your application explicitly requires that the order of tracks not change, you should run this command with the option -nthreads 0.
 

docs/reference/commands/tckstats.rst

@@ -6,7 +6,7 @@ tckstats
 Synopsis
 --------
 
-Calculate statistics on streamlines length
+Calculate statistics on streamlines lengths
 
 Usage
 --------
@@ -26,8 +26,6 @@ Options
 
 -  **-dump path** dump the streamlines lengths to a text file
 
--  **-explicit** explicitly calculate the length of each streamline, ignoring any step size information present in the header
-
 -  **-ignorezero** do not generate a warning if the track file contains streamlines with zero length
 
 -  **-tck_weights_in path** specify a text scalar file containing the streamline weights

docs/reference/commands_list.rst

@@ -204,7 +204,7 @@ List of MRtrix3 commands
     :ref:`tcksample`, "Sample values of an associated image along tracks"
     :ref:`tcksift2`, "Successor to the SIFT method; instead of removing streamlines, use an EM framework to find an appropriate cross-section multiplier for each streamline"
     :ref:`tcksift`, "Filter a whole-brain fibre-tracking data set such that the streamline densities match the FOD lobe integrals"
-    :ref:`tckstats`, "Calculate statistics on streamlines length"
+    :ref:`tckstats`, "Calculate statistics on streamlines lengths"
     :ref:`tcktransform`, "Apply a spatial transformation to a tracks file"
     :ref:`tensor2metric`, "Generate maps of tensor-derived parameters"
     :ref:`transformcalc`, "Perform calculations on linear transformation matrices"

src/dwi/tractography/algorithms/calibrator.h

@@ -42,7 +42,7 @@ namespace MR {
           for (ssize_t j = -extent; j <= extent; ++j) {
             float x = i + 0.5*j, y = SQRT_3_OVER_2*j;
             float n = Math::pow2(x) + Math::pow2(y);
-            if (n > maxR) 
+            if (n > maxR)
               continue;
             n = spacing * std::sqrt (n);
             float z = std::cos (n);
@@ -62,24 +62,25 @@ namespace MR {
         };
       }
 
-      template <class Method> 
-        void calibrate (Method& method) 
+      template <class Method>
+        void calibrate (Method& method)
         {
           typename Method::Calibrate calibrate_func (method);
           const float sqrt3 = std::sqrt (3.0);
+          const float max_angle = std::isfinite (method.S.max_angle_ho) ? method.S.max_angle_ho : method.S.max_angle_1o;
 
           vector<Pair> amps;
-          for (float el = 0.0; el < method.S.max_angle; el += 0.001) {
+          for (float el = 0.0; el < max_angle; el += 0.001) {
             amps.push_back (Pair (el, calibrate_func (el)));
             if (!std::isfinite (amps.back().amp) || amps.back().amp <= 0.0) break;
           }
           float zero = amps.back().el;
 
           float N_min = Inf, theta_min = NaN, ratio = NaN;
           for (size_t i = 1; i < amps.size(); ++i) {
-            float N = Math::pow2 (method.S.max_angle);
+            float N = Math::pow2 (max_angle);
             float Ns =  N * (1.0+amps[0].amp/amps[i].amp)/(2.0*Math::pow2(zero));
-            auto dirs = direction_grid (method.S.max_angle + amps[i].el, sqrt3*amps[i].el);
+            auto dirs = direction_grid (max_angle + amps[i].el, sqrt3*amps[i].el);
             N = Ns+dirs.size();
             //std::cout << amps[i].el << " " << amps[i].amp << " " << Ns << " " << dirs.size() << " " << Ns+dirs.size() << "\n";
             if (N > 0.0 && N < N_min) {
@@ -89,10 +90,10 @@ namespace MR {
             }
           }
 
-          method.calibrate_list = direction_grid (method.S.max_angle+theta_min, sqrt3*theta_min);
+          method.calibrate_list = direction_grid (max_angle+theta_min, sqrt3*theta_min);
           method.calibrate_ratio = ratio;
 
-          INFO ("rejection sampling will use " + str (method.calibrate_list.size()) 
+          INFO ("rejection sampling will use " + str (method.calibrate_list.size())
               + " directions with a ratio of " + str (method.calibrate_ratio) + " (predicted number of samples per step = " + str (N_min) + ")");
         }
 

src/dwi/tractography/algorithms/fact.h

@@ -52,17 +52,19 @@ namespace MR
           if (rk4)
             throw Exception ("4th-order Runge-Kutta integration not valid for FACT algorithm");
 
-          set_step_size (0.1f);
+          set_step_size (rk4 ? 0.5f : 0.1f, false);
+          set_num_points();
           set_cutoff (TCKGEN_DEFAULT_CUTOFF_FIXEL);
+
           // If user specifies the angle threshold manually, want to enforce this as-is at each step
           // If it's calculated automatically, it needs to be corrected for the fact that the permissible
           //   angle per step has been calculated within set_step_size(), but FACT will not curve at each
           //   step; only at the voxel transitions.
-          if (!App::get_options("angle").size()) {
-            max_angle *= vox() / step_size;
-            cos_max_angle = std::cos (max_angle);
+          if (!App::get_options ("angle").size()) {
+            max_angle_1o = std::min (max_angle_1o * vox() / step_size, float(Math::pi_2));
+            cos_max_angle_1o = std::cos (max_angle_1o);
           }
-          dot_threshold = std::cos (max_angle);
+          dot_threshold = std::cos (max_angle_1o);
 
           properties["method"] = "FACT";
         }
@@ -113,7 +115,7 @@ namespace MR
         const float max_norm = do_next (dir);
 
         if (max_norm < S.threshold)
-          return BAD_SIGNAL;
+          return MODEL;
 
         pos += S.step_size * dir;
         return CONTINUE;

src/dwi/tractography/algorithms/iFOD1.h

@@ -44,7 +44,7 @@ namespace MR
           SharedBase (diff_path, property_set),
           lmax (Math::SH::LforN (source.size(3))),
           max_trials (TCKGEN_DEFAULT_MAX_TRIALS_PER_STEP),
-          sin_max_angle (std::sin (max_angle)),
+          sin_max_angle_1o (std::sin (max_angle_1o)),
           mean_samples (0.0),
           mean_truncations (0.0),
           max_max_truncation (0.0),
@@ -57,15 +57,14 @@ namespace MR
             throw Exception ("Algorithm iFOD1 expects as input a spherical harmonic (SH) image");
           }
 
-          set_step_size (0.1f);
+          set_step_size (rk4 ? 0.5f : 0.1f, rk4);
+          // max_angle needs to be set because it influences the cone in which FOD amplitudes are sampled
           if (rk4) {
-            max_angle = 0.5 * max_angle_rk4;
-            INFO ("minimum radius of curvature = " + str(step_size / (max_angle_rk4 / (0.5 * Math::pi))) + " mm");
-          } else {
-            INFO ("minimum radius of curvature = " + str(step_size / ( 2.0 * sin (max_angle / 2.0))) + " mm");
+            max_angle_1o = 0.5 * max_angle_ho;
+            cos_max_angle_1o = std::cos (max_angle_1o);
           }
-          sin_max_angle = std::sin (max_angle);
-
+          sin_max_angle_1o = std::sin (max_angle_1o);
+          set_num_points();
           set_cutoff (TCKGEN_DEFAULT_CUTOFF_FOD);
 
           properties["method"] = "iFOD1";
@@ -101,7 +100,7 @@ namespace MR
         }
 
         size_t lmax, max_trials;
-        float sin_max_angle;
+        float sin_max_angle_1o;
         Math::SH::PrecomputedAL<float> precomputer;
 
         private:
@@ -212,7 +211,7 @@ namespace MR
           }
         }
 
-        return BAD_SIGNAL;
+        return MODEL;
       }
 
 
@@ -238,7 +237,7 @@ namespace MR
         );
       }
 
-      Eigen::Vector3f rand_dir (const Eigen::Vector3f& d) { return (random_direction (d, S.max_angle, S.sin_max_angle)); }
+      Eigen::Vector3f rand_dir (const Eigen::Vector3f& d) { return (random_direction (d, S.max_angle_1o, S.sin_max_angle_1o)); }