From a3641ec4caed2b51258fe8028f7898445c2c347e Mon Sep 17 00:00:00 2001 From: Robert Smith Date: Tue, 22 Mar 2016 17:06:52 +1100 Subject: [PATCH] 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 | 34 ++++++--- .../tractography/connectome/connectome.cpp | 5 +- src/dwi/tractography/connectome/connectome.h | 4 - src/dwi/tractography/editing/editing.cpp | 3 + src/dwi/tractography/editing/worker.cpp | 15 +++- src/dwi/tractography/editing/worker.h | 5 +- src/dwi/tractography/mapping/mapper.h | 2 + src/dwi/tractography/resample.cpp | 74 ++++++++++++++++++- src/dwi/tractography/resample.h | 32 +++++--- 9 files changed, 141 insertions(+), 33 deletions(-) diff --git a/cmd/tckedit.cpp b/cmd/tckedit.cpp index 90417b461f..343cd115f6 100644 --- a/cmd/tckedit.cpp +++ b/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(properties["step_size"])); - else - step_size = to(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(properties["step_size"])); + else + step_size = to(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 ( diff --git a/src/dwi/tractography/connectome/connectome.cpp b/src/dwi/tractography/connectome/connectome.cpp index 4708a75cfe..d4858ad438 100644 --- a/src/dwi/tractography/connectome/connectome.cpp +++ b/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 }; diff --git a/src/dwi/tractography/connectome/connectome.h b/src/dwi/tractography/connectome/connectome.h index 174e356715..191d084d59 100644 --- a/src/dwi/tractography/connectome/connectome.h +++ b/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&); diff --git a/src/dwi/tractography/editing/editing.cpp b/src/dwi/tractography/editing/editing.cpp index 942bba93cc..3c5f8c2931 100644 --- a/src/dwi/tractography/editing/editing.cpp +++ b/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"); diff --git a/src/dwi/tractography/editing/worker.cpp b/src/dwi/tractography/editing/worker.cpp index 8358e88775..542fd325b2 100644 --- a/src/dwi/tractography/editing/worker.cpp +++ b/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 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; diff --git a/src/dwi/tractography/editing/worker.h b/src/dwi/tractography/editing/worker.h index 8f7b613f08..0eb5371b9c 100644 --- a/src/dwi/tractography/editing/worker.h +++ b/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 diff --git a/src/dwi/tractography/mapping/mapper.h b/src/dwi/tractography/mapping/mapper.h index ddbcfd0472..89f49cc8bc 100644 --- a/src/dwi/tractography/mapping/mapper.h +++ b/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 diff --git a/src/dwi/tractography/resample.cpp b/src/dwi/tractography/resample.cpp index 93c06dff69..b906b28c49 100644 --- a/src/dwi/tractography/resample.cpp +++ b/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& tck) const + { + Math::Hermite interp (0.1); + std::vector 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; + } + + + + + + } } } diff --git a/src/dwi/tractography/resample.h b/src/dwi/tractography/resample.h index 3a87100461..7cdf9dadcb 100644 --- a/src/dwi/tractography/resample.h +++ b/src/dwi/tractography/resample.h @@ -19,7 +19,6 @@ #include -#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&) const; + + private: + float step_size; + + }; + + + }