Supporting multiple .tck files

.gitignore

@@ -17,6 +17,7 @@
 *.hdr
 *.qrc
 *.rste
+*.tar
 config
 build.*
 /configure.log
@@ -51,3 +52,4 @@ testing/bin/
 testing/src/project_version.cpp
 testing/src/project_version.h
 .__tmp.log
+.vscode

cmd/afdconnectivity.cpp

@@ -156,7 +156,7 @@ value_type AFDConnectivity::get (const std::string& path)
 
   Tractography::Properties properties;
   Tractography::Reader<value_type> reader (path, properties);
-  const size_t track_count = (properties.find ("count") == properties.end() ? 0 : to<size_t>(properties["count"]));
+  const size_t track_count = properties.value_or_default<size_t>("count",0);
   DWI::Tractography::Mapping::TrackLoader loader (reader, track_count, "summing apparent fibre density within track");
 
   // If WBFT is provided, this is the sum of (volume/length) across streamlines

cmd/fixel2tsf.cpp

@@ -86,7 +86,7 @@ void run ()
   float angular_threshold = get_option_value ("angle", DEFAULT_ANGULAR_THRESHOLD);
   const float angular_threshold_dp = cos (angular_threshold * (Math::pi / 180.0));
 
-  const size_t num_tracks = properties["count"].empty() ? 0 : to<int> (properties["count"]);
+  const size_t num_tracks = properties.value_or_default<size_t>("count",0);
 
   DWI::Tractography::Mapping::TrackMapperBase mapper (in_index_image);
   mapper.set_use_precise_mapping (true);

cmd/fixelcfestats.cpp

@@ -285,7 +285,7 @@ void run() {
   DWI::Tractography::Properties properties;
   DWI::Tractography::Reader<float> track_file (track_filename, properties);
   // Read in tracts, and compute whole-brain fixel-fixel connectivity
-  const size_t num_tracks = properties["count"].empty() ? 0 : to<size_t> (properties["count"]);
+  const size_t num_tracks = properties.value_or_default<size_t>("count",0);
   if (!num_tracks)
     throw Exception ("no tracks found in input file");
   if (num_tracks < 1000000)

cmd/tck2connectome.cpp

@@ -148,7 +148,7 @@ void execute (Image<node_t>& node_image, const node_t max_node_index, const std:
   Tractography::Reader<float> reader (argument[0], properties);
 
   // Initialise classes in preparation for multi-threading
-  Mapping::TrackLoader loader (reader, properties["count"].empty() ? 0 : to<size_t>(properties["count"]), "Constructing connectome");
+  Mapping::TrackLoader loader (reader, properties.value_or_default<size_t>("count",0), "Constructing connectome");
   Tractography::Connectome::Mapper mapper (*tck2nodes, metric);
   Tractography::Connectome::Matrix<T> connectome (max_node_index, statistic, vector_output, track_assignments);
 

cmd/tck2fixel.cpp

@@ -159,7 +159,7 @@ void run ()
   DWI::Tractography::Properties properties;
   DWI::Tractography::Reader<float> track_file (track_filename, properties);
   // Read in tracts, and compute whole-brain fixel-fixel connectivity
-  const size_t num_tracks = properties["count"].empty() ? 0 : to<int> (properties["count"]);
+  const size_t num_tracks = properties.value_or_default<size_t>("count",0);
   if (!num_tracks)
     throw Exception ("no tracks found in input file");
 

cmd/tckdfc.cpp

@@ -306,7 +306,7 @@ void run ()
     // TODO Constructor for properties using the file path?
     Tractography::Reader<float> tck_file (tck_path, properties);
   }
-  const size_t num_tracks = properties["count"].empty() ? 0 : to<size_t> (properties["count"]);
+  const size_t num_tracks = properties.value_or_default<size_t>("count",0);
 
   Image<float> fmri_image (Image<float>::open (argument[1]).with_direct_io(3));
 

cmd/tckedit.cpp

@@ -28,7 +28,6 @@
 #include "dwi/tractography/weights.h"
 
 #include "dwi/tractography/editing/editing.h"
-#include "dwi/tractography/editing/loader.h"
 #include "dwi/tractography/editing/receiver.h"
 #include "dwi/tractography/editing/worker.h"
 
@@ -122,87 +121,46 @@ void usage ()
 }
 
 
-void erase_if_present (Tractography::Properties& p, const std::string s)
-{
-  auto i = p.find (s);
-  if (i != p.end())
-    p.erase (i);
-}
-
-
-
 void run ()
 {
 
   const size_t num_inputs = argument.size() - 1;
   const std::string output_path = argument[num_inputs];
 
   // Make sure configuration is sensible
+  if (num_inputs == 0)
+    throw Exception("Insufficient number of inputs");
   if (get_options("tck_weights_in").size() && num_inputs > 1)
     throw Exception ("Cannot use per-streamline weighting with multiple input files");
 
-  // Get the consensus streamline properties from among the multiple input files
+  // Either one file: .lst or .tck
+  // OR multiple .tck files 
   Tractography::Properties properties;
-  size_t count = 0;
-  vector<std::string> input_file_list;
-
-  for (size_t file_index = 0; file_index != num_inputs; ++file_index) {
-
-    input_file_list.push_back (argument[file_index]);
-
-    Properties p;
-    Reader<float> (argument[file_index], p);
-
-    for (const auto& i : p.comments) {
-      bool present = false;
-      for (const auto& j: properties.comments)
-        if ( (present = (i == j)) )
-          break;
-      if (!present)
-        properties.comments.push_back (i);
-    }
-
-    for (const auto& i : p.prior_rois) {
-      const auto potential_matches = properties.prior_rois.equal_range (i.first);
-      bool present = false;
-      for (auto j = potential_matches.first; !present && j != potential_matches.second; ++j)
-        present = (i.second == j->second);
-      if (!present)
-        properties.prior_rois.insert (i);
-    }
-
-    size_t this_count = 0, this_total_count = 0;
-
-    for (const auto& i : p) {
-      if (i.first == "count") {
-        this_count = to<float> (i.second);
-      } else if (i.first == "total_count") {
-        this_total_count += to<float> (i.second);
-      } else {
-        auto existing = properties.find (i.first);
-        if (existing == properties.end())
-          properties.insert (i);
-        else if (i.second != existing->second)
-          existing->second = "variable";
-      }
-    }
-
-    count += this_count;
+  Tractography::Reader<float> *reader = nullptr;
+
+  if ( num_inputs == 1 ) {
+    reader = new Tractography::Reader<float>(argument[0], properties);
+  } else {
+    vector<std::string> input_file_list;
+    for (size_t file_index = 0; file_index < num_inputs; ++file_index)
+      input_file_list.push_back (argument[file_index]);
 
+    reader = new Tractography::Reader<float>(input_file_list, properties);
   }
 
-  DEBUG ("estimated number of input tracks: " + str(count));
+  properties.load_ROIs();
 
-  load_rois (properties);
+  size_t count = properties.value_or_default<size_t>("count",0);
+  DEBUG ("estimated number of input tracks: " + str(count));
 
   // Some properties from tracking may be overwritten by this editing process
   // Due to the potential use of masking, we have no choice but to clear the
   //   properties class of any fields that would otherwise propagate through
   //   and be applied as part of this editing
-  erase_if_present (properties, "min_dist");
-  erase_if_present (properties, "max_dist");
-  erase_if_present (properties, "min_weight");
-  erase_if_present (properties, "max_weight");
+  properties.erase_if_present("min_dist");
+  properties.erase_if_present("max_dist");
+  properties.erase_if_present("min_weight");
+  properties.erase_if_present("max_weight");
   Editing::load_properties (properties);
 
   // Parameters that the worker threads need to be aware of, but do not appear in Properties
@@ -213,15 +171,14 @@ void run ()
   const size_t number = get_option_value ("number", size_t(0));
   const size_t skip   = get_option_value ("skip",   size_t(0));
 
-  Loader loader (input_file_list);
   Worker worker (properties, inverse, 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)
   Receiver receiver (output_path, properties, number, skip);
 
   Thread::run_queue (
-      loader,
+      *reader,
       Thread::batch (Streamline<>()),
       Thread::multi (worker),
       Thread::batch (Streamline<>()),

cmd/tckinfo.cpp

@@ -52,20 +52,19 @@ void run ()
     std::cout << "***********************************\n";
     std::cout << "  Tracks file: \"" << argument[i] << "\"\n";
 
-    for (Tractography::Properties::iterator i = properties.begin(); i != properties.end(); ++i) {
-      std::string S (i->first + ':');
+    for (auto& p: properties) {
+      std::string S (p.first + ':');
       S.resize (22, ' ');
-      std::cout << "    " << S << i->second << "\n";
+      std::cout << "    " << S << p.second << "\n";
     }
 
-    if (properties.comments.size()) {
+    if (! properties.comments.empty()) {
       std::cout << "    Comments:             ";
-      for (vector<std::string>::iterator i = properties.comments.begin(); i != properties.comments.end(); ++i)
-        std::cout << (i == properties.comments.begin() ? "" : "                       ") << *i << "\n";
+      std::cout << join( properties.comments, "\n                          " ) << "\n";
     }
 
-    for (std::multimap<std::string,std::string>::const_iterator i = properties.prior_rois.begin(); i != properties.prior_rois.end(); ++i)
-      std::cout << "    ROI:                  " << i->first << " " << i->second << "\n";
+    for (auto& r: properties.prior_rois)
+      std::cout << "    ROI:                  " << r.first << " " << r.second << "\n";
 
 
 

cmd/tckmap.cpp

@@ -267,7 +267,7 @@ void run () {
   Tractography::Properties properties;
   Tractography::Reader<float> file (argument[0], properties);
 
-  const size_t num_tracks = properties["count"].empty() ? 0 : to<size_t> (properties["count"]);
+  const size_t num_tracks = properties.value_or_default<size_t>("count",0);
 
   vector<default_type> voxel_size = get_option_value ("vox", vector<default_type>());
 

cmd/tckresample.cpp

@@ -116,7 +116,7 @@ void run ()
 
   const std::unique_ptr<Resampling::Base> resampler (Resampling::get_resampler());
 
-  const float old_step_size = get_step_size (properties);
+  const float old_step_size = properties.get_step_size();
   if (!std::isfinite (old_step_size)) {
     INFO ("Do not have step size information from input track file");
   }
@@ -132,9 +132,7 @@ void run ()
   }
   properties["output_step_size"] = std::isfinite (new_step_size) ? str(new_step_size) : "variable";
 
-  auto downsample = properties.find ("downsample_factor");
-  if (downsample != properties.end())
-    properties.erase (downsample);
+  properties.erase_if_present("downsample_factor");
 
   Worker worker (resampler);
   Receiver receiver (argument[1], properties);

cmd/tcksample.cpp

@@ -465,9 +465,7 @@ void run ()
   if (nointerp && precise)
     throw Exception ("Option -nointerp and -precise are mutually exclusive");
   const interp_type interp = nointerp ? interp_type::NEAREST : (precise ? interp_type::PRECISE : interp_type::LINEAR);
-  const size_t num_tracks = properties.find("count") == properties.end() ?
-                            0 :
-                            to<size_t>(properties["count"]);
+  const size_t num_tracks = properties.value_or_default<size_t>("count",0);
 
   if (statistic == stat_tck::NONE && interp == interp_type::PRECISE)
     throw Exception ("Precise streamline mapping may only be used with per-streamline statistics");

cmd/tckstats.cpp

@@ -120,11 +120,10 @@ void run ()
     Tractography::Properties properties;
     Tractography::Reader<float> reader (argument[0], properties);
 
-    if (properties.find ("count") != properties.end())
-      header_count = to<size_t> (properties["count"]);
+    header_count = properties.value_or_default<size_t>("count",0);
 
     if (!get_options ("explicit").size()) {
-      step_size = get_step_size (properties);
+      step_size = properties.get_step_size();
       if (!std::isfinite (step_size) || !step_size) {
         INFO ("Streamline step size undefined in header; lengths will be calculated manually");
         if (get_options ("histogram").size()) {

cmd/tcktransform.cpp

@@ -44,23 +44,6 @@ void usage ()
 using value_type = float;
 using TrackType = Tractography::Streamline<value_type>;
 
-
-
-class Loader { MEMALIGN(Loader)
-  public:
-    Loader (const std::string& file) : reader (file, properties) {}
-
-    bool operator() (TrackType& item) {
-      return reader (item);
-    }
-
-    Tractography::Properties properties;
-  protected:
-    Tractography::Reader<value_type> reader;
-};
-
-
-
 class Warper { MEMALIGN(Warper)
   public:
     Warper (const Image<value_type>& warp) :
@@ -120,15 +103,16 @@ class Writer { MEMALIGN(Writer)
 
 void run ()
 {
-  Loader loader (argument[0]);
+  Tractography::Properties properties;
+  Tractography::Reader<value_type> reader ( argument[0], properties );
 
   auto data = Image<value_type>::open (argument[1]).with_direct_io (3);
   Warper warper (data);
 
-  Writer writer (argument[2], loader.properties);
+  Writer writer (argument[2], properties);
 
   Thread::run_queue (
-      loader,
+      reader,
       Thread::batch (TrackType(), 1024),
       Thread::multi (warper),
       Thread::batch (TrackType(), 1024),

core/app.cpp

@@ -1156,7 +1156,7 @@ namespace MR
         }
         if (i.arg->type == ArgDirectoryOut)
           check_overwrite (text);
-        if (i.arg->type == TracksIn && !Path::has_suffix (text, ".tck"))
+        if (i.arg->type == TracksIn && !Path::has_suffix (text, ".tck") && !Path::has_suffix (text, ".lst"))
           throw Exception ("input file \"" + text + "\" is not a valid track file");
         if (i.arg->type == TracksOut && !Path::has_suffix (text, ".tck"))
           throw Exception ("output track file \"" + text + "\" must use the .tck suffix");
@@ -1184,7 +1184,7 @@ namespace MR
           }
           if (arg.type == ArgDirectoryOut)
             check_overwrite (text);
-          if (arg.type == TracksIn && !Path::has_suffix (text, ".tck"))
+          if (arg.type == TracksIn && !Path::has_suffix (text, ".tck") && !Path::has_suffix (text, ".lst"))
             throw Exception ("input file \"" + text + "\" for option \"-" + std::string(i.opt->id) + "\" is not a valid track file");
           if (arg.type == TracksOut && !Path::has_suffix (text, ".tck"))
             throw Exception ("output track file \"" + text + "\" for option \"-" + std::string(i.opt->id) + "\" must use the .tck suffix");

core/command.h

@@ -18,7 +18,10 @@
 #define __command_h__
 
 
-#include <xmmintrin.h>
+#ifdef FLUSH_TO_ZERO
+# include <xmmintrin.h>
+#endif
+
 #include "app.h"
 #include "exec_version.h"
 #ifdef MRTRIX_PROJECT

core/file/utils.h

@@ -205,6 +205,25 @@ namespace MR
 
 
 
+    inline size_t readlines(const std::string& name, vector<std::string>& lines)
+    {
+      std::string line;
+      std::ifstream in(name.c_str());
+      size_t count = 0;
+
+      if (!in)
+        throw Exception("cannot open file \"" + name + "\": " + strerror(errno));
+
+      while (getline(in,line)) {
+        lines.push_back(line);
+        ++count;
+      }
+
+      in.close();
+      return count;
+    }
+
+
 
     inline std::string create_tempfile (int64_t size = 0, const char* suffix = NULL)
     {