tckedit: New option -resample

This option will re-sample a streamline to a requested step size, but will always retain the two endpoints.

cmd/tckedit.cpp

@@ -69,6 +69,7 @@ void usage ()
   + Option ("test_ends_only", "only test the ends of each streamline against the provided include/exclude ROIs")
 
   // TODO Input weights with multiple input files currently not supported
+  + OptionGroup ("Options for handling streamline weights")
   + Tractography::TrackWeightsInOption
   + Tractography::TrackWeightsOutOption;
 
@@ -86,16 +87,23 @@ void usage ()
 
 
 
-void update_output_step_size (Tractography::Properties& properties, const int upsample_ratio, const int downsample_ratio)
+void update_output_step_size (Tractography::Properties& properties, const int upsample_ratio, const int downsample_ratio, const float forced_step_size)
 {
-  if (upsample_ratio == 1 && downsample_ratio == 1)
-    return;
-  float step_size = 0.0;
-  if (properties.find ("output_step_size") == properties.end())
-    step_size = (properties.find ("step_size") == properties.end() ? 0.0 : to<float>(properties["step_size"]));
-  else
-    step_size = to<float>(properties["output_step_size"]);
-  properties["output_step_size"] = str(step_size * float(downsample_ratio) / float(upsample_ratio));
+  if (std::isfinite (forced_step_size)) {
+    properties["output_step_size"] = str(forced_step_size);
+  } else {
+    if (upsample_ratio == 1 && downsample_ratio == 1)
+      return;
+    float step_size = 0.0;
+    if (properties.find ("output_step_size") == properties.end())
+      step_size = (properties.find ("step_size") == properties.end() ? 0.0 : to<float>(properties["step_size"]));
+    else
+      step_size = to<float>(properties["output_step_size"]);
+    properties["output_step_size"] = str(step_size * float(downsample_ratio) / float(upsample_ratio));
+  }
+  auto downsample = properties.find ("downsample_factor");
+  if (downsample != properties.end())
+    properties.erase (downsample);
 }
 
 
@@ -167,6 +175,10 @@ void run ()
   const int upsample   = opt.size() ? int(opt[0][0]) : 1;
   opt = get_options ("downsample");
   const int downsample = opt.size() ? int(opt[0][0]) : 1;
+  opt = get_options ("resample");
+  if (opt.size() && (upsample != 1 || downsample != 1))
+    throw Exception ("Cannot combine -resample with -upsample or -downsample (order is ambiguous)");
+  const float step_size = opt.size() ? float(opt[0][0]) : NaN;
   const bool inverse = get_options ("inverse").size();
   const bool test_ends_only = get_options ("test_ends_only").size();
   const bool out_ends_only = get_options ("out_ends_only").size();
@@ -178,11 +190,11 @@ void run ()
   const size_t skip   = opt.size() ? size_t(opt[0][0]) : 0;
 
   Loader loader (input_file_list);
-  Worker worker (properties, upsample, downsample, inverse, test_ends_only);
+  Worker worker (properties, upsample, downsample, step_size, inverse, test_ends_only);
   // This needs to be run AFTER creation of the Worker class
   // (worker needs to be able to set max & min number of points based on step size in input file,
   //  receiver needs "output_step_size" field to have been updated before file creation)
-  update_output_step_size (properties, upsample, downsample);
+  update_output_step_size (properties, upsample, downsample, step_size);
   Receiver receiver (output_path, properties, number, skip, out_ends_only);
 
   Thread::run_queue (

src/dwi/tractography/connectome/connectome.cpp

@@ -30,12 +30,11 @@ namespace Connectome {
 
 
 
-const char* modes[] = { "assignment_end_voxels", "assignment_radial_search", "assignment_reverse_search", "assignment_forward_search", "assignment_all_voxels", NULL };
+using namespace App;
 
-const char* metrics[] = { "count", "meanlength", "invlength", "invnodevolume", "invlength_invnodevolume", "mean_scalar", NULL };
 
 
-using namespace App;
+const char* modes[] = { "assignment_end_voxels", "assignment_radial_search", "assignment_reverse_search", "assignment_forward_search", "assignment_all_voxels", NULL };
 
 
 

src/dwi/tractography/connectome/connectome.h

@@ -46,11 +46,7 @@ class Metric;
 
 
 
-extern const char* metrics[];
 extern const char* modes[];
-
-
-
 extern const App::OptionGroup AssignmentOptions;
 Tck2nodes_base* load_assignment_mode (Image<node_t>&);
 

src/dwi/tractography/editing/editing.cpp

@@ -53,6 +53,9 @@ const OptionGroup ResampleOption = OptionGroup ("Streamline resampling options")
                           "(decreases required storage space)")
     + Argument ("ratio").type_integer (1, 1, 1e6)
 
+  + Option ("resample", "re-sample the streamline to a desired step size (in mm)")
+    + Argument ("step_size").type_float (0.0, 0.0, 1e6)
+
   + Option ("out_ends_only", "only output the two endpoints of each streamline");
 
 

src/dwi/tractography/editing/worker.cpp

@@ -33,12 +33,13 @@ namespace MR {
           out.weight = in.weight;
 
           if (!thresholds (in)) {
-            // Want to test thresholds before wasting time on upsampling; but if -inverse is set,
-            //   still need to apply both the upsampler and downsampler before writing to output
+            // Want to test thresholds before wasting time on resampling; but if -inverse is set,
+            //   still need to apply resampling before writing to output
             if (inverse) {
               std::vector<Eigen::Vector3f> tck (in);
               upsampler (tck);
               downsampler (tck);
+              resampler (tck);
               tck.swap (out);
             }
             return true;
@@ -60,6 +61,7 @@ namespace MR {
                 if (properties.exclude.contains (p)) {
                   if (inverse) {
                     downsampler (tck);
+                    resampler (tck);
                     tck.swap (out);
                   }
                   return true;
@@ -71,6 +73,7 @@ namespace MR {
                 if (properties.exclude.contains (p)) {
                   if (inverse) {
                     downsampler (tck);
+                    resampler (tck);
                     tck.swap (out);
                   }
                   return true;
@@ -83,6 +86,7 @@ namespace MR {
               if (!i) {
                 if (inverse) {
                   downsampler (tck);
+                  resampler (tck);
                   tck.swap (out);
                 }
                 return true;
@@ -113,9 +117,11 @@ namespace MR {
             if (cropped_tracks.empty())
               return true;
 
-            // Apply downsampler independently to each
-            for (auto& i : cropped_tracks)
+            // Apply downsampler / resampler independently to each
+            for (auto& i : cropped_tracks) {
               downsampler (i);
+              resampler (i);
+            }
 
             if (cropped_tracks.size() == 1) {
               cropped_tracks[0].swap (out);
@@ -136,6 +142,7 @@ namespace MR {
 
             if (!inverse) {
               downsampler (tck);
+              resampler (tck);
               tck.swap (out);
             }
             return true;

src/dwi/tractography/editing/worker.h

@@ -38,10 +38,11 @@ namespace MR {
         {
 
           public:
-            Worker (Tractography::Properties& p, const size_t upsample_ratio, const size_t downsample_ratio, const bool inv, const bool end) :
+            Worker (Tractography::Properties& p, const size_t upsample_ratio, const size_t downsample_ratio, const float step_size, const bool inv, const bool end) :
               properties (p),
               upsampler (upsample_ratio),
               downsampler (downsample_ratio),
+              resampler (step_size),
               inverse (inv),
               ends_only (end),
               thresholds (p),
@@ -51,6 +52,7 @@ namespace MR {
               properties (that.properties),
               upsampler (that.upsampler),
               downsampler (that.downsampler),
+              resampler (that.resampler),
               inverse (that.inverse),
               ends_only (that.ends_only),
               thresholds (that.thresholds),
@@ -64,6 +66,7 @@ namespace MR {
             const Tractography::Properties& properties;
             Upsampler upsampler;
             Downsampler downsampler;
+            Resampler resampler;
             const bool inverse, ends_only;
 
             class Thresholds

src/dwi/tractography/mapping/mapper.h

@@ -34,6 +34,8 @@
 #include "dwi/tractography/mapping/twi_stats.h"
 #include "dwi/tractography/mapping/voxel.h"
 
+#include "math/hermite.h"
+
 
 // Didn't bother making this a command-line option, since curvature contrast results were very poor regardless of smoothing
 #define CURVATURE_TRACK_SMOOTHING_FWHM 10.0 // In mm

src/dwi/tractography/resample.cpp

@@ -15,6 +15,8 @@
 
 #include "dwi/tractography/resample.h"
 
+#include "math/hermite.h"
+
 
 namespace MR {
   namespace DWI {
@@ -70,9 +72,9 @@ namespace MR {
       {
         if (!M.rows() || in.size() < 2)
           return false;
-        const size_t s = in.size();
         // Abandoned curvature-based extrapolation - badly posed when step size is not guaranteed to be consistent,
         //   and probably makes little difference anyways
+        const size_t s = in.size();
         in.insert    (in.begin(), in[ 0 ] + (in[1] - in[0]));
         in.push_back (            in[ s ] + (in[s] - in[s-1]));
         for (size_t i = 0; i != 3; ++i) {
@@ -98,6 +100,9 @@ namespace MR {
 
 
 
+
+
+
       bool Downsampler::operator() (Tracking::GeneratedTrack& tck) const
       {
         if (ratio <= 1 || tck.empty())
@@ -141,6 +146,73 @@ namespace MR {
 
 
 
+      bool Resampler::operator() (std::vector<Eigen::Vector3f>& tck) const
+      {
+        Math::Hermite<float> interp (0.1);
+        std::vector<Eigen::Vector3f> output;
+        // Extensions required to enable Hermite interpolation in last streamline segment at either end
+        const size_t s = tck.size();
+        tck.insert    (tck.begin(), tck[ 0 ] + (tck[1] - tck[0]));
+        tck.push_back (             tck[ s ] + (tck[s] - tck[s-1]));
+        const ssize_t midpoint = tck.size()/2;
+        output.push_back (tck[midpoint]);
+        // Generate from the midpoint to the start, reverse, then generate from midpoint to the end
+        for (ssize_t step = -1; step <= 1; step += 2) {
+
+          ssize_t index = midpoint;
+          float mu_lower = 0.0f;
+
+          // Loop to generate points
+          do {
+
+            // If we don't have to step along the input track, can keep the mu from the previous
+            //   interpolation point as the lower bound
+            while (index > 1 && index < ssize_t(tck.size()-2) && (output.back() - tck[index+step]).norm() < step_size) {
+              index += step;
+              mu_lower = 0.0f;
+            }
+            // Always preserve the termination points, regardless of resampling
+            if (index == 1) {
+              output.push_back (tck[1]);
+              std::reverse (output.begin(), output.end());
+            } else if (index == ssize_t(tck.size()-2)) {
+              output.push_back (tck[s]);
+            } else {
+
+              // Perform binary search
+              Eigen::Vector3f p_lower = tck[index], p, p_upper = tck[index+step];
+              float mu_upper = 1.0f;
+              float mu = 0.5 * (mu_lower + mu_upper);
+              do {
+                mu = 0.5 * (mu_lower + mu_upper);
+                interp.set (mu);
+                p = interp.value (tck[index-step], tck[index], tck[index+step], tck[index+2*step]);
+                if ((p - output.back()).norm() < step_size) {
+                  mu_lower = mu;
+                  p_lower = p;
+                } else {
+                  mu_upper = mu;
+                  p_upper = p;
+                }
+              } while ((p_upper - p_lower).norm() > 0.001 * step_size);
+              output.push_back (p);
+
+            }
+
+            // Loop until an endpoint has been added
+          } while (index > 1 && index < ssize_t(tck.size()-2));
+
+        }
+
+        std::swap (tck, output);
+        return true;
+      }
+
+
+
+
+
+
     }
   }
 }

src/dwi/tractography/resample.h

@@ -19,7 +19,6 @@
 
 #include <vector>
 
-#include "math/hermite.h"
 #include "dwi/tractography/tracking/generated_track.h"
 
 
@@ -66,14 +65,6 @@ namespace MR {
 
 
 
-
-
-
-
-
-
-
-
       class Downsampler
       {
 
@@ -95,6 +86,29 @@ namespace MR {
 
 
 
+      class Resampler
+      {
+
+        public:
+          Resampler () :
+              step_size (0.0) { }
+
+          Resampler (const float ss) :
+              step_size (ss) { }
+
+          void set_step_size (const float ss) { step_size = ss; }
+          float get_step_size() const { return step_size; }
+          bool valid() const { return step_size; }
+
+          bool operator() (std::vector<Eigen::Vector3f>&) const;
+
+        private:
+          float step_size;
+
+      };
+
+
+
 
 
     }