From 8a54b7f1c7db97f31f8a516d4852e05ceb7c2a94 Mon Sep 17 00:00:00 2001 From: Jonathan H Date: Thu, 14 Feb 2019 14:20:22 +0000 Subject: [PATCH 01/11] Minor housekeeping: - add dependency explicitly (previously included implicitly via properties->roi->image->header), - and remove duplicate include. --- src/dwi/tractography/file_base.h | 1 + src/dwi/tractography/roi.h | 1 - 2 files changed, 1 insertion(+), 1 deletion(-) diff --git a/src/dwi/tractography/file_base.h b/src/dwi/tractography/file_base.h index 904e6f5027..582b626281 100644 --- a/src/dwi/tractography/file_base.h +++ b/src/dwi/tractography/file_base.h @@ -20,6 +20,7 @@ #include #include "types.h" +#include "datatype.h" #include "file/key_value.h" #include "file/ofstream.h" #include "file/path.h" diff --git a/src/dwi/tractography/roi.h b/src/dwi/tractography/roi.h index fd89e759c4..ce7a9b8112 100644 --- a/src/dwi/tractography/roi.h +++ b/src/dwi/tractography/roi.h @@ -18,7 +18,6 @@ #include "app.h" #include "image.h" -#include "image.h" #include "interp/linear.h" #include "math/rng.h" From c8cc8e753213d30b8d9c7a387dd5b7859a3b9117 Mon Sep 17 00:00:00 2001 From: Jonathan H Date: Thu, 14 Feb 2019 23:57:30 +0000 Subject: [PATCH 02/11] Add new readlines utility, parsing lines of file to vector of strings. --- core/file/utils.h | 19 +++++++++++++++++++ 1 file changed, 19 insertions(+) diff --git a/core/file/utils.h b/core/file/utils.h index 6cc94615db..93ed886f0f 100644 --- a/core/file/utils.h +++ b/core/file/utils.h @@ -172,6 +172,25 @@ namespace MR + inline size_t readlines(const std::string& name, std::vector& 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) { From 97a24c20eb6574e17761f78e8a47139b21319fa5 Mon Sep 17 00:00:00 2001 From: Jonathan H Date: Fri, 15 Feb 2019 10:18:45 +0000 Subject: [PATCH 03/11] Enhance Tractography::Properties, move implementations to .cpp file, and adapt commands accordingly. --- cmd/afdconnectivity.cpp | 2 +- cmd/fixel2tsf.cpp | 2 +- cmd/fixelcfestats.cpp | 2 +- cmd/tck2connectome.cpp | 2 +- cmd/tck2fixel.cpp | 2 +- cmd/tckdfc.cpp | 2 +- cmd/tckedit.cpp | 17 ++---- cmd/tckinfo.cpp | 15 +++-- cmd/tckmap.cpp | 2 +- cmd/tckresample.cpp | 6 +- cmd/tcksample.cpp | 4 +- cmd/tckstats.cpp | 5 +- src/dwi/tractography/properties.cpp | 63 +++++++++++++++++++ src/dwi/tractography/properties.h | 93 ++++++++++++----------------- 14 files changed, 125 insertions(+), 92 deletions(-) diff --git a/cmd/afdconnectivity.cpp b/cmd/afdconnectivity.cpp index 47c5ba479d..d42d1a62c5 100644 --- a/cmd/afdconnectivity.cpp +++ b/cmd/afdconnectivity.cpp @@ -155,7 +155,7 @@ value_type AFDConnectivity::get (const std::string& path) Tractography::Properties properties; Tractography::Reader reader (path, properties); - const size_t track_count = (properties.find ("count") == properties.end() ? 0 : to(properties["count"])); + const size_t track_count = properties.value_or_default("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 diff --git a/cmd/fixel2tsf.cpp b/cmd/fixel2tsf.cpp index dc949e2b38..cfd2606205 100644 --- a/cmd/fixel2tsf.cpp +++ b/cmd/fixel2tsf.cpp @@ -85,7 +85,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 (properties["count"]); + const size_t num_tracks = properties.value_or_default("count",0); DWI::Tractography::Mapping::TrackMapperBase mapper (in_index_image); mapper.set_use_precise_mapping (true); diff --git a/cmd/fixelcfestats.cpp b/cmd/fixelcfestats.cpp index 3ef898a091..cf298e6f21 100644 --- a/cmd/fixelcfestats.cpp +++ b/cmd/fixelcfestats.cpp @@ -284,7 +284,7 @@ void run() { DWI::Tractography::Properties properties; DWI::Tractography::Reader 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 (properties["count"]); + const size_t num_tracks = properties.value_or_default("count",0); if (!num_tracks) throw Exception ("no tracks found in input file"); if (num_tracks < 1000000) diff --git a/cmd/tck2connectome.cpp b/cmd/tck2connectome.cpp index d90d596561..5f7e9f8756 100644 --- a/cmd/tck2connectome.cpp +++ b/cmd/tck2connectome.cpp @@ -110,7 +110,7 @@ void execute (Image& node_image, const node_t max_node_index, const std: Tractography::Reader reader (argument[0], properties); // Initialise classes in preparation for multi-threading - Mapping::TrackLoader loader (reader, properties["count"].empty() ? 0 : to(properties["count"]), "Constructing connectome"); + Mapping::TrackLoader loader (reader, properties.value_or_default("count",0), "Constructing connectome"); Tractography::Connectome::Mapper mapper (*tck2nodes, metric); Tractography::Connectome::Matrix connectome (max_node_index, statistic, vector_output, track_assignments); diff --git a/cmd/tck2fixel.cpp b/cmd/tck2fixel.cpp index 1b18af82f3..9c4e1bbf9c 100644 --- a/cmd/tck2fixel.cpp +++ b/cmd/tck2fixel.cpp @@ -158,7 +158,7 @@ void run () DWI::Tractography::Properties properties; DWI::Tractography::Reader 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 (properties["count"]); + const size_t num_tracks = properties.value_or_default("count",0); if (!num_tracks) throw Exception ("no tracks found in input file"); diff --git a/cmd/tckdfc.cpp b/cmd/tckdfc.cpp index cf08daf1ea..84f9c3d515 100644 --- a/cmd/tckdfc.cpp +++ b/cmd/tckdfc.cpp @@ -305,7 +305,7 @@ void run () // TODO Constructor for properties using the file path? Tractography::Reader tck_file (tck_path, properties); } - const size_t num_tracks = properties["count"].empty() ? 0 : to (properties["count"]); + const size_t num_tracks = properties.value_or_default("count",0); Image fmri_image (Image::open (argument[1]).with_direct_io(3)); diff --git a/cmd/tckedit.cpp b/cmd/tckedit.cpp index 37ada0e6f8..7e321cadf8 100644 --- a/cmd/tckedit.cpp +++ b/cmd/tckedit.cpp @@ -86,15 +86,6 @@ 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 () { @@ -156,10 +147,10 @@ void run () // 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 diff --git a/cmd/tckinfo.cpp b/cmd/tckinfo.cpp index 1c14e0ef17..f3db8bd12a 100644 --- a/cmd/tckinfo.cpp +++ b/cmd/tckinfo.cpp @@ -51,20 +51,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::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::const_iterator i = properties.roi.begin(); i != properties.roi.end(); ++i) - std::cout << " ROI: " << i->first << " " << i->second << "\n"; + for (auto& r: properties.roi) + std::cout << " ROI: " << r.first << " " << r.second << "\n"; diff --git a/cmd/tckmap.cpp b/cmd/tckmap.cpp index 2dcc510d51..d89f0e340b 100644 --- a/cmd/tckmap.cpp +++ b/cmd/tckmap.cpp @@ -266,7 +266,7 @@ void run () { Tractography::Properties properties; Tractography::Reader file (argument[0], properties); - const size_t num_tracks = properties["count"].empty() ? 0 : to (properties["count"]); + const size_t num_tracks = properties.value_or_default("count",0); vector voxel_size = get_option_value ("vox", vector()); diff --git a/cmd/tckresample.cpp b/cmd/tckresample.cpp index c9731e9adb..7afce17b17 100644 --- a/cmd/tckresample.cpp +++ b/cmd/tckresample.cpp @@ -115,7 +115,7 @@ void run () const std::unique_ptr 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"); } @@ -131,9 +131,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); diff --git a/cmd/tcksample.cpp b/cmd/tcksample.cpp index 87dbaabed7..8a41f79e9f 100644 --- a/cmd/tcksample.cpp +++ b/cmd/tcksample.cpp @@ -464,9 +464,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(properties["count"]); + const size_t num_tracks = properties.value_or_default("count",0); if (statistic == stat_tck::NONE && interp == interp_type::PRECISE) throw Exception ("Precise streamline mapping may only be used with per-streamline statistics"); diff --git a/cmd/tckstats.cpp b/cmd/tckstats.cpp index 8ec420d69d..6f7d0bfbd5 100644 --- a/cmd/tckstats.cpp +++ b/cmd/tckstats.cpp @@ -119,11 +119,10 @@ void run () Tractography::Properties properties; Tractography::Reader reader (argument[0], properties); - if (properties.find ("count") != properties.end()) - header_count = to (properties["count"]); + header_count = properties.value_or_default("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()) { diff --git a/src/dwi/tractography/properties.cpp b/src/dwi/tractography/properties.cpp index 10c26a30ea..d2721bd3ae 100644 --- a/src/dwi/tractography/properties.cpp +++ b/src/dwi/tractography/properties.cpp @@ -34,6 +34,69 @@ namespace MR { roi.clear(); } + + + + + void check_compatible (const Properties& a, const Properties& b, const std::string& type) + { + if ( a.get_step_size() != b.get_step_size() ) + WARN ("mismatch between " + type + " properties (step-size)"); + + for (auto& p: {"source","act","method","threshold"}) + if (a.has(p) && b.has(p) && a.get(p) != b.get(p)) + WARN ("mismatch between " + type + " properties (" + p + ")"); + } + + + + + void check_timestamps (const Properties& a, const Properties& b, const std::string& type) + { + auto stamp_a = a.find ("timestamp"); + auto stamp_b = b.find ("timestamp"); + if (stamp_a == a.end() || stamp_b == b.end()) + throw Exception ("unable to verify " + type + " pair: missing timestamp"); + if (stamp_a->second != stamp_b->second) + throw Exception ("invalid " + type + " combination - timestamps do not match"); + } + + + + + void check_counts (const Properties& a, const Properties& b, const std::string& type, bool abort_on_fail) + { + auto count_a = a.find ("count"); + auto count_b = b.find ("count"); + if ((count_a == a.end()) || (count_b == b.end())) { + std::string mesg = "unable to validate " + type + " pair: missing count field"; + if (abort_on_fail) throw Exception (mesg); + else WARN (mesg); + } + if (to(count_a->second) != to(count_b->second)) { + std::string mesg = type + " files do not contain same number of elements"; + if (abort_on_fail) throw Exception (mesg); + else WARN (mesg); + } + } + + + + + std::ostream& operator<< (std::ostream& stream, const Properties& P) + { + stream << "seeds: " << P.seeds; + stream << "include: " << P.include << ", exclude: " << P.exclude << ", mask: " << P.mask << ", dict: "; + for (auto& Pk: P) + stream << "[ " << Pk.first << ": " << Pk.second << " ], "; + + stream << "comments: "; + for (auto& c: P.comments) + stream << "\"" << c << "\", "; + + return stream; + } + } } } diff --git a/src/dwi/tractography/properties.h b/src/dwi/tractography/properties.h index 70c4736ee0..eb8369dd13 100644 --- a/src/dwi/tractography/properties.h +++ b/src/dwi/tractography/properties.h @@ -17,6 +17,7 @@ #define __dwi_tractography_properties_h__ #include +#include "mrtrix.h" #include "timer.h" #include "dwi/tractography/roi.h" #include "dwi/tractography/seeding/list.h" @@ -66,79 +67,63 @@ namespace MR roi.clear(); } - template void set (T& variable, const std::string& name) { - if ((*this)[name].empty()) (*this)[name] = str (variable); - else variable = to ((*this)[name]); + bool has(const std::string& name) const { + return this->find(name) != this->end(); } - void load_ROIs (); - }; - - - inline void check_timestamps (const Properties& a, const Properties& b, const std::string& type) - { - Properties::const_iterator stamp_a = a.find ("timestamp"); - Properties::const_iterator stamp_b = b.find ("timestamp"); - if (stamp_a == a.end() || stamp_b == b.end()) - throw Exception ("unable to verify " + type + " pair: missing timestamp"); - if (stamp_a->second != stamp_b->second) - throw Exception ("invalid " + type + " combination - timestamps do not match"); - } + std::string get(const std::string& name) const { + auto p = this->find(name); + if (p == this->end()) + throw Exception("Undefined property: \"" + name + "\""); + else + return p->second; + } + template + T value_or_default(const std::string& name, const T& defval) const { + auto p = this->find(name); + return p == this->end() ? defval : to(p->second); + } + template + void set (T& variable, const std::string& name) { + if ((*this)[name].empty()) (*this)[name] = str (variable); + else variable = to ((*this)[name]); + } - inline void check_counts (const Properties& a, const Properties& b, const std::string& type, bool abort_on_fail) - { - Properties::const_iterator count_a = a.find ("count"); - Properties::const_iterator count_b = b.find ("count"); - if ((count_a == a.end()) || (count_b == b.end())) { - std::string mesg = "unable to validate " + type + " pair: missing count field"; - if (abort_on_fail) throw Exception (mesg); - else WARN (mesg); - } - if (to(count_a->second) != to(count_b->second)) { - std::string mesg = type + " files do not contain same number of elements"; - if (abort_on_fail) throw Exception (mesg); - else WARN (mesg); - } - } + void erase_if_present(const std::string& name) { + auto p = this->find(name); + if (p != this->end()) + this->erase(p); + } + inline float get_step_size() const { + return value_or_default( "output_step_size", value_or_default( "step_size", NaN ) ); + } + void load_ROIs (); + }; - inline float get_step_size (const Properties& p) - { - for (size_t index = 0; index != 2; ++index) { - const std::string key (index ? "step_size" : "output_step_size"); - auto it = p.find (key); - if (it != p.end()) { - try { - return to (it->second); - } catch (...) { } - } - } - return NaN; + // JH: kept for now, but should be removed in favour of: p.get_step_size() + inline float get_step_size (const Properties& p) { + DEBUG( "Function get_step_size is deprecated. Call properties.get_step_size() instead." ); + return p.get_step_size(); } + // JH: unused for now + void check_compatible(const Properties& a, const Properties& b, const std::string& type); + void check_timestamps (const Properties& a, const Properties& b, const std::string& type); - inline std::ostream& operator<< (std::ostream& stream, const Properties& P) - { - stream << "seeds: " << P.seeds; - stream << "include: " << P.include << ", exclude: " << P.exclude << ", mask: " << P.mask << ", dict: "; - for (std::map::const_iterator i = P.begin(); i != P.end(); ++i) - stream << "[ " << i->first << ": " << i->second << " ], "; - stream << "comments: "; - for (vector::const_iterator i = P.comments.begin(); i != P.comments.end(); ++i) - stream << "\"" << *i << "\", "; - return (stream); - } + void check_counts (const Properties& a, const Properties& b, const std::string& type, bool abort_on_fail); + std::ostream& operator<< (std::ostream& stream, const Properties& P); } } From ac61b7cac81f59399742682fe9c7ed53ae5a2f94 Mon Sep 17 00:00:00 2001 From: Jonathan H Date: Thu, 14 Feb 2019 23:57:50 +0000 Subject: [PATCH 04/11] Introduce TrackFileInfo object to decouple parsing of header information from streamline-reading. Introduce properties_consensus method, based on tckedit's current implementation, to find a consensus between multiple header information. Extend Tractography::Reader to support multiple tracks files by default (and associated weights files). --- core/app.cpp | 4 +- src/dwi/tractography/file.h | 309 ++++++++++++++++++++--------- src/dwi/tractography/file_base.cpp | 142 ++++++++++--- src/dwi/tractography/file_base.h | 60 +++++- 4 files changed, 384 insertions(+), 131 deletions(-) diff --git a/core/app.cpp b/core/app.cpp index 51b77b9d85..6e8562d11a 100644 --- a/core/app.cpp +++ b/core/app.cpp @@ -1085,7 +1085,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"); @@ -1113,7 +1113,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"); diff --git a/src/dwi/tractography/file.h b/src/dwi/tractography/file.h index aebe2ad377..c0dcbeb9f8 100644 --- a/src/dwi/tractography/file.h +++ b/src/dwi/tractography/file.h @@ -24,6 +24,7 @@ #include "file/config.h" #include "file/key_value.h" #include "file/ofstream.h" +#include "file/utils.h" #include "dwi/tractography/file_base.h" #include "dwi/tractography/properties.h" #include "dwi/tractography/streamline.h" @@ -60,7 +61,7 @@ namespace MR virtual bool operator() (const Streamline&) = 0; virtual ~WriterInterface() { } }; - + //! A class to read streamlines data @@ -70,64 +71,50 @@ namespace MR public: //! open the \c file for reading and load header into \c properties - Reader (const std::string& file, Properties& properties) : - current_index (0) { - open (file, "tracks", properties); - auto opt = App::get_options ("tck_weights_in"); - if (opt.size()) { - weights_file.reset (new std::ifstream (str(opt[0][0]).c_str(), std::ios_base::in)); - if (!weights_file->good()) - throw Exception ("Unable to open streamlines weights file " + str(opt[0][0])); - } + Reader (const std::string& file, Properties& prop) + { + std::vector files; + if (Path::has_suffix(file, ".lst")){ + File::readlines( file, files ); + INFO("opening list of " + str(files.size()) + " .tck file(s)"); } + else + files.push_back(file); + + construct( files, prop ); + } + Reader(const std::vector& files, Properties& prop) + { construct(files, prop); } - //! fetch next track from file - bool operator() (Streamline& tck) { - tck.clear(); + ~Reader() + { close_current_files(false); clear(); } - if (!in.is_open()) - return false; - do { - auto p = get_next_point(); - if (std::isinf (p[0])) { - in.close(); - check_excess_weights(); - return false; - } - if (in.eof()) { - in.close(); - check_excess_weights(); - return false; - } - if (std::isnan (p[0])) { - tck.index = current_index++; - if (weights_file) { + //! reset member properties + void clear() + { + std::ifstream a, b; // dummy streams + + file_index = 0; + track_index = sub_count = 0; + dtype = DataType::Undefined; + in.swap(a); + inw.swap(b); + tck_files.clear(); + weights_files.clear(); + } - (*weights_file) >> tck.weight; - if (weights_file->fail()) { - WARN ("Streamline weights file contains less entries than .tck file; only read " + str(current_index-1) + " streamlines"); - in.close(); - tck.clear(); - return false; - } - } else { - tck.weight = 1.0; - } - return true; - } - tck.push_back (p); - } while (in.good()); + //! fetch next track from file + inline bool operator() (Streamline& tck) { + return get_next_track(tck); + } - in.close(); - return false; - } @@ -135,57 +122,201 @@ namespace MR using __ReaderBase__::in; using __ReaderBase__::dtype; - uint64_t current_index; - std::unique_ptr weights_file; + size_t file_index; + uint64_t track_index, sub_count; + std::ifstream inw; - //! takes care of byte ordering issues + std::vector tck_files, weights_files; + + + + + //! initialise current reader instance + void construct( const std::vector& files, Properties& prop ) + { + clear(); + + if (files.empty()) + return; + + // assign properties (also checks that input files are valid .tck files) + properties_consensus(files, "tracks", prop); + tck_files = files; + + // try to find associated weights + auto opt = App::get_options ("tck_weights_in"); + std::string weights; + + if (opt.size()) { + + weights = str(opt[0][0]); + + if ( files.size()>1 && !Path::has_suffix(weights,".lst") ) + throw Exception("When reading several .tck files, the weights file should be a list with extension .lst, containing one filename per line."); + + if (Path::has_suffix(weights,".lst")) + File::readlines(weights, weights_files); + else + weights_files.push_back(weights); + + // check that number of files matches + if (weights_files.size() != tck_files.size()) + throw Exception("Mismatch between number of .tck files and number of streamline weights files."); + } + } + + + + + //! close currently opened files + void close_current_files(bool check_excess=true) + { + // close tck file + close(); + + // close weights file + float temp; + if (inw.is_open()) { + if (check_excess) { + inw >> temp; + if (!inw.fail()) + WARN ("Streamline weights file contains more entries than .tck file"); + } + + inw.close(); + } + + } + + + + + //! open next file from file list + bool open_next_files() + { + close_current_files(false); + + if (file_index >= tck_files.size()) + return false; - Eigen::Matrix get_next_point () + Properties dummy; + open( tck_files[file_index], "tracks", dummy ); + + if (!weights_files.empty()) { - using namespace ByteOrder; - switch (dtype()) { - case DataType::Float32LE: - { - float p[3]; - in.read ((char*) p, sizeof (p)); - return { ValueType(LE(p[0])), ValueType(LE(p[1])), ValueType(LE(p[2])) }; - } - case DataType::Float32BE: - { - float p[3]; - in.read ((char*) p, sizeof (p)); - return { ValueType(BE(p[0])), ValueType(BE(p[1])), ValueType(BE(p[2])) }; - } - case DataType::Float64LE: - { - double p[3]; - in.read ((char*) p, sizeof (p)); - return { ValueType(LE(p[0])), ValueType(LE(p[1])), ValueType(LE(p[2])) }; - } - case DataType::Float64BE: - { - double p[3]; - in.read ((char*) p, sizeof (p)); - return { ValueType(BE(p[0])), ValueType(BE(p[1])), ValueType(BE(p[2])) }; + inw.open( weights_files[file_index].c_str(), std::ios_base::in ); + if (!inw) + throw Exception ("Unable to open streamlines weights file " + weights_files[file_index]); + } + + ++file_index; + sub_count = 0; + return true; + } + + + + + //! get next track in currently opened file + bool get_next_track( Streamline& tck ) + { + tck.clear(); + + if ( !in.is_open() && !open_next_files() ) + return false; + + // keep adding points to current streamline + do { + + auto p = get_next_point(); + + // if we hit eof or Inf value, call again to open next file + if ( in.eof() || std::isinf(p[0]) ) { + close_current_files(); + return get_next_track(tck); + } + + // end of streamline: set tck index and weight, and return + if (std::isnan(p[0])) { + + tck.index = track_index; + ++track_index; + ++sub_count; + + if (inw.is_open()) { + + inw >> tck.weight; + if (inw.fail()) { + --track_index; + --sub_count; + close_current_files(false); + + WARN ("Streamline weights file contains fewer entries than .tck file; only read " + str(sub_count) + " streamlines."); + return get_next_track(tck); } - default: - assert (0); - break; + + } else { + tck.weight = 1.0; + } + + return true; } - return { NaN, NaN, NaN }; - } - //! Check that the weights file does not contain excess entries - void check_excess_weights() + // otherwise add point and iterate + tck.push_back (p); + + } while ( in.good() ); + + WARN ("Streamline-parsing from tracks file ended prematurely."); + close_current_files(); + return get_next_track(tck); + } + + + + + + //! takes care of byte ordering issues + Eigen::Matrix get_next_point () { - if (!weights_file) - return; - float temp; - (*weights_file) >> temp; - if (!weights_file->fail()) - WARN ("Streamline weights file contains more entries than .tck file"); + using namespace ByteOrder; + switch (dtype()) { + case DataType::Float32LE: + { + float p[3]; + in.read ((char*) p, sizeof (p)); + return { ValueType(LE(p[0])), ValueType(LE(p[1])), ValueType(LE(p[2])) }; + } + case DataType::Float32BE: + { + float p[3]; + in.read ((char*) p, sizeof (p)); + return { ValueType(BE(p[0])), ValueType(BE(p[1])), ValueType(BE(p[2])) }; + } + case DataType::Float64LE: + { + double p[3]; + in.read ((char*) p, sizeof (p)); + return { ValueType(LE(p[0])), ValueType(LE(p[1])), ValueType(LE(p[2])) }; + } + case DataType::Float64BE: + { + double p[3]; + in.read ((char*) p, sizeof (p)); + return { ValueType(BE(p[0])), ValueType(BE(p[1])), ValueType(BE(p[2])) }; + } + default: + assert (0); + break; + } + return { NaN, NaN, NaN }; } + + + + + //! copy construction explicitly disabled Reader (const Reader&) = delete; }; diff --git a/src/dwi/tractography/file_base.cpp b/src/dwi/tractography/file_base.cpp index f0fb3505e1..89c27de51f 100644 --- a/src/dwi/tractography/file_base.cpp +++ b/src/dwi/tractography/file_base.cpp @@ -21,63 +21,139 @@ namespace MR { namespace Tractography { - void __ReaderBase__::open (const std::string& file, const std::string& type, Properties& properties) + void __ReaderBase__::open (const std::string& file, const std::string& type, Properties& prop) { - properties.clear(); - dtype = DataType::Undefined; + TrackFileInfo info( file, type, prop ); + dtype = info.dtype; + + in.open (info.path.c_str(), std::ios::in | std::ios::binary); + if (!in) + throw Exception ("error opening " + type + " data file \"" + info.path + "\": " + strerror(errno)); + in.seekg (info.offset); + } + + - const std::string firstline ("mrtrix " + type); - File::KeyValue kv (file, firstline.c_str()); - std::string data_file; + void TrackFileInfo::check_dtype() + { + if (dtype == DataType::Undefined) + throw Exception ("no datatype specified for tracks file \"" + path + "\""); + if (dtype != DataType::Float32LE && dtype != DataType::Float32BE && + dtype != DataType::Float64LE && dtype != DataType::Float64BE) + throw Exception ("only supported datatype for tracks file are " + "Float32LE, Float32BE, Float64LE & Float64BE (in " + type + " file \"" + path + "\")"); + } + + + + + void TrackFileInfo::parse( const std::string& file, const std::string& file_type, Properties& prop ) + { + clear(); + prop.clear(); + type = file_type; + + // initialise KeyValue instance to parse metadata + const std::string magic ("mrtrix " + type); + File::KeyValue kv (file, magic.c_str()); + + // parse key/value pairs while (kv.next()) { const std::string key = lowercase (kv.key()); if (key == "roi") { try { vector V (split (kv.value(), " \t", true, 2)); - properties.roi.insert (std::pair (V[0], V[1])); + prop.roi.insert (std::pair (V[0], V[1])); } catch (...) { WARN ("invalid ROI specification in " + type + " file \"" + file + "\" - ignored"); } } - else if (key == "comment") properties.comments.push_back (kv.value()); - else if (key == "file") data_file = kv.value(); - else if (key == "datatype") dtype = DataType::parse (kv.value()); - else properties[kv.key()] = kv.value(); + else if (key == "comment") prop.comments.push_back (kv.value()); + else if (key == "file") path = kv.value(); + else if (key == "datatype") dtype = DataType::parse (kv.value()); + else prop[kv.key()] = kv.value(); } - if (dtype == DataType::Undefined) - throw Exception ("no datatype specified for tracks file \"" + file + "\""); - if (dtype != DataType::Float32LE && dtype != DataType::Float32BE && - dtype != DataType::Float64LE && dtype != DataType::Float64BE) - throw Exception ("only supported datatype for tracks file are " - "Float32LE, Float32BE, Float64LE & Float64BE (in " + type + " file \"" + file + "\")"); + // check field "datatype" + check_dtype(); - if (data_file.empty()) - throw Exception ("missing \"files\" specification for " + type + " file \"" + file + "\""); + // field "file" should be defined + if (path.empty()) + throw Exception ("missing \"file\" specification for " + type + " file \"" + file + "\""); - std::istringstream files_stream (data_file); - std::string fname; - files_stream >> fname; - int64_t offset = 0; - if (files_stream.good()) { - try { files_stream >> offset; } + // its value should be: "filename offset" (no spaces in filename!) + std::istringstream file_val (path); + std::string path_val; + file_val >> path_val; + + if (file_val.good()) { + try { file_val >> offset; } catch (...) { throw Exception ("invalid offset specified for file \"" - + fname + "\" in " + type + " file \"" + file + "\""); + + path_val + "\" in " + type + " file \"" + file + "\""); } } - if (fname != ".") - fname = Path::join (Path::dirname (file), fname); + // value "." interpreted as "this file" + if (path_val != ".") + path = Path::join (Path::dirname (file), path_val); else - fname = file; + path = file; + } - in.open (fname.c_str(), std::ios::in | std::ios::binary); - if (!in) - throw Exception ("error opening " + type + " data file \"" + fname + "\": " + strerror(errno)); - in.seekg (offset); + + + void properties_consensus( const std::vector& files, const std::string& type, Properties& prop ) + { + // ensure there is at least one file + size_t n_files = files.size(); + if (n_files == 0) + throw Exception("empty list of files for consensus!"); + + // read first file + TrackFileInfo info_ref( files[0], type, prop ); + + // create comments buffer and summary counts + std::unordered_set comments_buffer( prop.comments.begin(), prop.comments.end() ); + uint64_t count = prop.value_or_default("count",0); + uint64_t total_count = prop.value_or_default("total_count",0); + + // iterate over remaining files + TrackFileInfo info_tmp; + Properties prop_tmp; + for ( size_t k=1; k < n_files; ++k ) { + + info_tmp.parse( files[k], type, prop_tmp ); + if ( info_tmp.dtype != info_ref.dtype ) + throw Exception("datatype mismatch between " + type + " files \"" + files[0] + "\" and \"" + files[k] + "\""); + + for (auto& c: prop_tmp.comments) + if ( comments_buffer.find(c) == comments_buffer.end() ) { + prop.comments.push_back(c); + comments_buffer.insert(c); + } + + for (auto& p: prop_tmp) + if (p.first == "count") + count += to(p.second); + else if (p.first == "total_count") + total_count += to(p.second); + else { + auto existing = prop.find (p.first); + if (existing == prop.end()) + prop.insert(p); + else if (p.second != existing->second) { + INFO("variable property \"" + p.first + "\" across " + type + " files"); + existing->second = "variable"; + } + } + } + + // final counts + prop["count"] = str(count); + prop["total_count"] = str(total_count); } } diff --git a/src/dwi/tractography/file_base.h b/src/dwi/tractography/file_base.h index 582b626281..ebdda04107 100644 --- a/src/dwi/tractography/file_base.h +++ b/src/dwi/tractography/file_base.h @@ -17,7 +17,7 @@ #define __dwi_tractography_file_base_h__ #include -#include +#include #include "types.h" #include "datatype.h" @@ -36,18 +36,63 @@ namespace MR namespace Tractography { + + // JH: information about .tck or .tsf file + struct TrackFileInfo + { + std::string path, type; + uint64_t offset; + DataType dtype; + + TrackFileInfo() + { clear(); } + + TrackFileInfo( const std::string& file, const std::string& type, Properties& prop ) + { clear(); parse(file,type,prop); } + + + void clear() { + dtype = DataType::Undefined; + offset = 0; + path = ""; + type = ""; + } + + + inline bool empty() + { return type.empty() || path.empty(); } + + + void check_dtype(); + + + void parse( const std::string& file, const std::string& type, Properties& prop ); + + + }; + + + + // combine properties across several .tck files + void properties_consensus( const std::vector& files, const std::string& type, Properties& prop ); + + + //! \cond skip class __ReaderBase__ { NOMEMALIGN public: - ~__ReaderBase__ () { - if (in.is_open()) - in.close(); - } - void open (const std::string& file, const std::string& firstline, Properties& properties); + ~__ReaderBase__ () + { close(); } - void close () { in.close(); } + + void close () + { if (in.is_open()) in.close(); } + + + void open (const std::string& file, const std::string& type, Properties& prop); + protected: @@ -56,6 +101,7 @@ namespace MR }; + template class __WriterBase__ { NOMEMALIGN From 88fd979add14b128cfc0d3e0bf4d4f9114480a29 Mon Sep 17 00:00:00 2001 From: Jonathan H Date: Sun, 17 Feb 2019 11:25:21 +0000 Subject: [PATCH 05/11] Fix spacing --- cmd/tckinfo.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/cmd/tckinfo.cpp b/cmd/tckinfo.cpp index f3db8bd12a..4da91300e3 100644 --- a/cmd/tckinfo.cpp +++ b/cmd/tckinfo.cpp @@ -59,7 +59,7 @@ void run () if (! properties.comments.empty()) { std::cout << " Comments: "; - std::cout << join( properties.comments, "\n " ) << "\n"; + std::cout << join( properties.comments, "\n " ) << "\n"; } for (auto& r: properties.roi) From 7bf3343a3dae664d26348469a4d405d90195de40 Mon Sep 17 00:00:00 2001 From: Jonathan H Date: Wed, 20 Feb 2019 10:40:00 +0000 Subject: [PATCH 06/11] Merge changes proposed within dev branch --- cmd/afdconnectivity.cpp | 2 +- cmd/fixel2tsf.cpp | 2 +- cmd/fixelcfestats.cpp | 2 +- cmd/tck2connectome.cpp | 2 +- cmd/tck2fixel.cpp | 2 +- cmd/tckdfc.cpp | 2 +- cmd/tckedit.cpp | 17 +- cmd/tckinfo.cpp | 15 +- cmd/tckmap.cpp | 2 +- cmd/tckresample.cpp | 6 +- cmd/tcksample.cpp | 4 +- cmd/tckstats.cpp | 5 +- core/app.cpp | 4 +- core/command.h | 5 +- core/file/utils.h | 19 ++ src/dwi/tractography/file.h | 309 ++++++++++++++++++++-------- src/dwi/tractography/file_base.cpp | 142 ++++++++++--- src/dwi/tractography/file_base.h | 61 +++++- src/dwi/tractography/properties.cpp | 103 ++++++++++ src/dwi/tractography/properties.h | 92 ++++----- src/dwi/tractography/roi.h | 1 - 21 files changed, 573 insertions(+), 224 deletions(-) create mode 100644 src/dwi/tractography/properties.cpp diff --git a/cmd/afdconnectivity.cpp b/cmd/afdconnectivity.cpp index a3d31bf275..e210bb0167 100644 --- a/cmd/afdconnectivity.cpp +++ b/cmd/afdconnectivity.cpp @@ -156,7 +156,7 @@ value_type AFDConnectivity::get (const std::string& path) Tractography::Properties properties; Tractography::Reader reader (path, properties); - const size_t track_count = (properties.find ("count") == properties.end() ? 0 : to(properties["count"])); + const size_t track_count = properties.value_or_default("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 diff --git a/cmd/fixel2tsf.cpp b/cmd/fixel2tsf.cpp index e236a37d7f..d5ab20a1ee 100644 --- a/cmd/fixel2tsf.cpp +++ b/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 (properties["count"]); + const size_t num_tracks = properties.value_or_default("count",0); DWI::Tractography::Mapping::TrackMapperBase mapper (in_index_image); mapper.set_use_precise_mapping (true); diff --git a/cmd/fixelcfestats.cpp b/cmd/fixelcfestats.cpp index 9227c7a728..f2297b7be4 100644 --- a/cmd/fixelcfestats.cpp +++ b/cmd/fixelcfestats.cpp @@ -285,7 +285,7 @@ void run() { DWI::Tractography::Properties properties; DWI::Tractography::Reader 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 (properties["count"]); + const size_t num_tracks = properties.value_or_default("count",0); if (!num_tracks) throw Exception ("no tracks found in input file"); if (num_tracks < 1000000) diff --git a/cmd/tck2connectome.cpp b/cmd/tck2connectome.cpp index 39781dfdf7..c77d2bd18b 100644 --- a/cmd/tck2connectome.cpp +++ b/cmd/tck2connectome.cpp @@ -148,7 +148,7 @@ void execute (Image& node_image, const node_t max_node_index, const std: Tractography::Reader reader (argument[0], properties); // Initialise classes in preparation for multi-threading - Mapping::TrackLoader loader (reader, properties["count"].empty() ? 0 : to(properties["count"]), "Constructing connectome"); + Mapping::TrackLoader loader (reader, properties.value_or_default("count",0), "Constructing connectome"); Tractography::Connectome::Mapper mapper (*tck2nodes, metric); Tractography::Connectome::Matrix connectome (max_node_index, statistic, vector_output, track_assignments); diff --git a/cmd/tck2fixel.cpp b/cmd/tck2fixel.cpp index 570cdf392c..5e48ddb8dc 100644 --- a/cmd/tck2fixel.cpp +++ b/cmd/tck2fixel.cpp @@ -159,7 +159,7 @@ void run () DWI::Tractography::Properties properties; DWI::Tractography::Reader 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 (properties["count"]); + const size_t num_tracks = properties.value_or_default("count",0); if (!num_tracks) throw Exception ("no tracks found in input file"); diff --git a/cmd/tckdfc.cpp b/cmd/tckdfc.cpp index 2ba276f731..94bf2b810f 100644 --- a/cmd/tckdfc.cpp +++ b/cmd/tckdfc.cpp @@ -306,7 +306,7 @@ void run () // TODO Constructor for properties using the file path? Tractography::Reader tck_file (tck_path, properties); } - const size_t num_tracks = properties["count"].empty() ? 0 : to (properties["count"]); + const size_t num_tracks = properties.value_or_default("count",0); Image fmri_image (Image::open (argument[1]).with_direct_io(3)); diff --git a/cmd/tckedit.cpp b/cmd/tckedit.cpp index c2aeb404a9..56d3967753 100644 --- a/cmd/tckedit.cpp +++ b/cmd/tckedit.cpp @@ -122,15 +122,6 @@ 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 () { @@ -199,10 +190,10 @@ void run () // 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 diff --git a/cmd/tckinfo.cpp b/cmd/tckinfo.cpp index 8c44b99c7d..4cd6281058 100644 --- a/cmd/tckinfo.cpp +++ b/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::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::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"; diff --git a/cmd/tckmap.cpp b/cmd/tckmap.cpp index e092e72e22..e5de529cd0 100644 --- a/cmd/tckmap.cpp +++ b/cmd/tckmap.cpp @@ -267,7 +267,7 @@ void run () { Tractography::Properties properties; Tractography::Reader file (argument[0], properties); - const size_t num_tracks = properties["count"].empty() ? 0 : to (properties["count"]); + const size_t num_tracks = properties.value_or_default("count",0); vector voxel_size = get_option_value ("vox", vector()); diff --git a/cmd/tckresample.cpp b/cmd/tckresample.cpp index 81f3e77a27..cf084ad7cd 100644 --- a/cmd/tckresample.cpp +++ b/cmd/tckresample.cpp @@ -116,7 +116,7 @@ void run () const std::unique_ptr 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); diff --git a/cmd/tcksample.cpp b/cmd/tcksample.cpp index 5dbfa21d19..6f430e3606 100644 --- a/cmd/tcksample.cpp +++ b/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(properties["count"]); + const size_t num_tracks = properties.value_or_default("count",0); if (statistic == stat_tck::NONE && interp == interp_type::PRECISE) throw Exception ("Precise streamline mapping may only be used with per-streamline statistics"); diff --git a/cmd/tckstats.cpp b/cmd/tckstats.cpp index 311a45775b..2c061da448 100644 --- a/cmd/tckstats.cpp +++ b/cmd/tckstats.cpp @@ -120,11 +120,10 @@ void run () Tractography::Properties properties; Tractography::Reader reader (argument[0], properties); - if (properties.find ("count") != properties.end()) - header_count = to (properties["count"]); + header_count = properties.value_or_default("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()) { diff --git a/core/app.cpp b/core/app.cpp index 479e5183cf..b6ae21bb3c 100644 --- a/core/app.cpp +++ b/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"); diff --git a/core/command.h b/core/command.h index 7f47780480..0228883cc3 100644 --- a/core/command.h +++ b/core/command.h @@ -18,7 +18,10 @@ #define __command_h__ -#include +#ifdef FLUSH_TO_ZERO +# include +#endif + #include "app.h" #include "exec_version.h" #ifdef MRTRIX_PROJECT diff --git a/core/file/utils.h b/core/file/utils.h index ad392269f9..308c6f966d 100644 --- a/core/file/utils.h +++ b/core/file/utils.h @@ -205,6 +205,25 @@ namespace MR + inline size_t readlines(const std::string& name, std::vector& 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) { diff --git a/src/dwi/tractography/file.h b/src/dwi/tractography/file.h index 43c2e5712a..c5ec874124 100644 --- a/src/dwi/tractography/file.h +++ b/src/dwi/tractography/file.h @@ -25,6 +25,7 @@ #include "file/config.h" #include "file/key_value.h" #include "file/ofstream.h" +#include "file/utils.h" #include "dwi/tractography/file_base.h" #include "dwi/tractography/properties.h" #include "dwi/tractography/streamline.h" @@ -61,7 +62,7 @@ namespace MR virtual bool operator() (const Streamline&) = 0; virtual ~WriterInterface() { } }; - + //! A class to read streamlines data @@ -71,64 +72,50 @@ namespace MR public: //! open the \c file for reading and load header into \c properties - Reader (const std::string& file, Properties& properties) : - current_index (0) { - open (file, "tracks", properties); - auto opt = App::get_options ("tck_weights_in"); - if (opt.size()) { - weights_file.reset (new std::ifstream (str(opt[0][0]).c_str(), std::ios_base::in)); - if (!weights_file->good()) - throw Exception ("Unable to open streamlines weights file " + str(opt[0][0])); - } + Reader (const std::string& file, Properties& prop) + { + std::vector files; + if (Path::has_suffix(file, ".lst")){ + File::readlines( file, files ); + INFO("opening list of " + str(files.size()) + " .tck file(s)"); } + else + files.push_back(file); + + construct( files, prop ); + } + Reader(const std::vector& files, Properties& prop) + { construct(files, prop); } - //! fetch next track from file - bool operator() (Streamline& tck) { - tck.clear(); + ~Reader() + { close_current_files(false); clear(); } - if (!in.is_open()) - return false; - do { - auto p = get_next_point(); - if (std::isinf (p[0])) { - in.close(); - check_excess_weights(); - return false; - } - if (in.eof()) { - in.close(); - check_excess_weights(); - return false; - } - if (std::isnan (p[0])) { - tck.index = current_index++; - if (weights_file) { + //! reset member properties + void clear() + { + std::ifstream a, b; // dummy streams + + file_index = 0; + track_index = sub_count = 0; + dtype = DataType::Undefined; + in.swap(a); + inw.swap(b); + tck_files.clear(); + weights_files.clear(); + } - (*weights_file) >> tck.weight; - if (weights_file->fail()) { - WARN ("Streamline weights file contains less entries than .tck file; only read " + str(current_index-1) + " streamlines"); - in.close(); - tck.clear(); - return false; - } - } else { - tck.weight = 1.0; - } - return true; - } - tck.push_back (p); - } while (in.good()); + //! fetch next track from file + inline bool operator() (Streamline& tck) { + return get_next_track(tck); + } - in.close(); - return false; - } @@ -136,57 +123,201 @@ namespace MR using __ReaderBase__::in; using __ReaderBase__::dtype; - uint64_t current_index; - std::unique_ptr weights_file; + size_t file_index; + uint64_t track_index, sub_count; + std::ifstream inw; - //! takes care of byte ordering issues + std::vector tck_files, weights_files; + + + + + //! initialise current reader instance + void construct( const std::vector& files, Properties& prop ) + { + clear(); + + if (files.empty()) + return; + + // assign properties (also checks that input files are valid .tck files) + properties_consensus(files, "tracks", prop); + tck_files = files; + + // try to find associated weights + auto opt = App::get_options ("tck_weights_in"); + std::string weights; + + if (opt.size()) { + + weights = str(opt[0][0]); + + if ( files.size()>1 && !Path::has_suffix(weights,".lst") ) + throw Exception("When reading several .tck files, the weights file should be a list with extension .lst, containing one filename per line."); + + if (Path::has_suffix(weights,".lst")) + File::readlines(weights, weights_files); + else + weights_files.push_back(weights); + + // check that number of files matches + if (weights_files.size() != tck_files.size()) + throw Exception("Mismatch between number of .tck files and number of streamline weights files."); + } + } + + + + + //! close currently opened files + void close_current_files(bool check_excess=true) + { + // close tck file + close(); + + // close weights file + float temp; + if (inw.is_open()) { + if (check_excess) { + inw >> temp; + if (!inw.fail()) + WARN ("Streamline weights file contains more entries than .tck file"); + } + + inw.close(); + } + + } + + + + + //! open next file from file list + bool open_next_files() + { + close_current_files(false); + + if (file_index >= tck_files.size()) + return false; - Eigen::Matrix get_next_point () + Properties dummy; + open( tck_files[file_index], "tracks", dummy ); + + if (!weights_files.empty()) { - using namespace ByteOrder; - switch (dtype()) { - case DataType::Float32LE: - { - float p[3]; - in.read ((char*) p, sizeof (p)); - return { ValueType(LE(p[0])), ValueType(LE(p[1])), ValueType(LE(p[2])) }; - } - case DataType::Float32BE: - { - float p[3]; - in.read ((char*) p, sizeof (p)); - return { ValueType(BE(p[0])), ValueType(BE(p[1])), ValueType(BE(p[2])) }; - } - case DataType::Float64LE: - { - double p[3]; - in.read ((char*) p, sizeof (p)); - return { ValueType(LE(p[0])), ValueType(LE(p[1])), ValueType(LE(p[2])) }; - } - case DataType::Float64BE: - { - double p[3]; - in.read ((char*) p, sizeof (p)); - return { ValueType(BE(p[0])), ValueType(BE(p[1])), ValueType(BE(p[2])) }; + inw.open( weights_files[file_index].c_str(), std::ios_base::in ); + if (!inw) + throw Exception ("Unable to open streamlines weights file " + weights_files[file_index]); + } + + ++file_index; + sub_count = 0; + return true; + } + + + + + //! get next track in currently opened file + bool get_next_track( Streamline& tck ) + { + tck.clear(); + + if ( !in.is_open() && !open_next_files() ) + return false; + + // keep adding points to current streamline + do { + + auto p = get_next_point(); + + // if we hit eof or Inf value, call again to open next file + if ( in.eof() || std::isinf(p[0]) ) { + close_current_files(); + return get_next_track(tck); + } + + // end of streamline: set tck index and weight, and return + if (std::isnan(p[0])) { + + tck.index = track_index; + ++track_index; + ++sub_count; + + if (inw.is_open()) { + + inw >> tck.weight; + if (inw.fail()) { + --track_index; + --sub_count; + close_current_files(false); + + WARN ("Streamline weights file contains fewer entries than .tck file; only read " + str(sub_count) + " streamlines."); + return get_next_track(tck); } - default: - assert (0); - break; + + } else { + tck.weight = 1.0; + } + + return true; } - return { NaN, NaN, NaN }; - } - //! Check that the weights file does not contain excess entries - void check_excess_weights() + // otherwise add point and iterate + tck.push_back (p); + + } while ( in.good() ); + + WARN ("Streamline-parsing from tracks file ended prematurely."); + close_current_files(); + return get_next_track(tck); + } + + + + + + //! takes care of byte ordering issues + Eigen::Matrix get_next_point () { - if (!weights_file) - return; - float temp; - (*weights_file) >> temp; - if (!weights_file->fail()) - WARN ("Streamline weights file contains more entries than .tck file"); + using namespace ByteOrder; + switch (dtype()) { + case DataType::Float32LE: + { + float p[3]; + in.read ((char*) p, sizeof (p)); + return { ValueType(LE(p[0])), ValueType(LE(p[1])), ValueType(LE(p[2])) }; + } + case DataType::Float32BE: + { + float p[3]; + in.read ((char*) p, sizeof (p)); + return { ValueType(BE(p[0])), ValueType(BE(p[1])), ValueType(BE(p[2])) }; + } + case DataType::Float64LE: + { + double p[3]; + in.read ((char*) p, sizeof (p)); + return { ValueType(LE(p[0])), ValueType(LE(p[1])), ValueType(LE(p[2])) }; + } + case DataType::Float64BE: + { + double p[3]; + in.read ((char*) p, sizeof (p)); + return { ValueType(BE(p[0])), ValueType(BE(p[1])), ValueType(BE(p[2])) }; + } + default: + assert (0); + break; + } + return { NaN, NaN, NaN }; } + + + + + //! copy construction explicitly disabled Reader (const Reader&) = delete; }; diff --git a/src/dwi/tractography/file_base.cpp b/src/dwi/tractography/file_base.cpp index 85205646f0..699a056889 100644 --- a/src/dwi/tractography/file_base.cpp +++ b/src/dwi/tractography/file_base.cpp @@ -22,63 +22,139 @@ namespace MR { namespace Tractography { - void __ReaderBase__::open (const std::string& file, const std::string& type, Properties& properties) + void __ReaderBase__::open (const std::string& file, const std::string& type, Properties& prop) { - properties.clear(); - dtype = DataType::Undefined; + TrackFileInfo info( file, type, prop ); + dtype = info.dtype; + + in.open (info.path.c_str(), std::ios::in | std::ios::binary); + if (!in) + throw Exception ("error opening " + type + " data file \"" + info.path + "\": " + strerror(errno)); + in.seekg (info.offset); + } + + - const std::string firstline ("mrtrix " + type); - File::KeyValue kv (file, firstline.c_str()); - std::string data_file; + void TrackFileInfo::check_dtype() + { + if (dtype == DataType::Undefined) + throw Exception ("no datatype specified for tracks file \"" + path + "\""); + if (dtype != DataType::Float32LE && dtype != DataType::Float32BE && + dtype != DataType::Float64LE && dtype != DataType::Float64BE) + throw Exception ("only supported datatype for tracks file are " + "Float32LE, Float32BE, Float64LE & Float64BE (in " + type + " file \"" + path + "\")"); + } + + + + + void TrackFileInfo::parse( const std::string& file, const std::string& file_type, Properties& prop ) + { + clear(); + prop.clear(); + type = file_type; + + // initialise KeyValue instance to parse metadata + const std::string magic ("mrtrix " + type); + File::KeyValue kv (file, magic.c_str()); + + // parse key/value pairs while (kv.next()) { const std::string key = lowercase (kv.key()); if (key == "roi" || key == "prior_roi") { try { vector V (split (kv.value(), " \t", true, 2)); - properties.prior_rois.insert (std::pair (V[0], V[1])); + prop.prior_rois.insert (std::pair (V[0], V[1])); } catch (...) { WARN ("invalid ROI specification in " + type + " file \"" + file + "\" - ignored"); } } - else if (key == "comment") properties.comments.push_back (kv.value()); - else if (key == "file") data_file = kv.value(); - else if (key == "datatype") dtype = DataType::parse (kv.value()); - else properties[kv.key()] = kv.value(); + else if (key == "comment") prop.comments.push_back (kv.value()); + else if (key == "file") path = kv.value(); + else if (key == "datatype") dtype = DataType::parse (kv.value()); + else prop[kv.key()] = kv.value(); } - if (dtype == DataType::Undefined) - throw Exception ("no datatype specified for tracks file \"" + file + "\""); - if (dtype != DataType::Float32LE && dtype != DataType::Float32BE && - dtype != DataType::Float64LE && dtype != DataType::Float64BE) - throw Exception ("only supported datatype for tracks file are " - "Float32LE, Float32BE, Float64LE & Float64BE (in " + type + " file \"" + file + "\")"); + // check field "datatype" + check_dtype(); - if (data_file.empty()) - throw Exception ("missing \"files\" specification for " + type + " file \"" + file + "\""); + // field "file" should be defined + if (path.empty()) + throw Exception ("missing \"file\" specification for " + type + " file \"" + file + "\""); - std::istringstream files_stream (data_file); - std::string fname; - files_stream >> fname; - int64_t offset = 0; - if (files_stream.good()) { - try { files_stream >> offset; } + // its value should be: "filename offset" (no spaces in filename!) + std::istringstream file_val (path); + std::string path_val; + file_val >> path_val; + + if (file_val.good()) { + try { file_val >> offset; } catch (...) { throw Exception ("invalid offset specified for file \"" - + fname + "\" in " + type + " file \"" + file + "\""); + + path_val + "\" in " + type + " file \"" + file + "\""); } } - if (fname != ".") - fname = Path::join (Path::dirname (file), fname); + // value "." interpreted as "this file" + if (path_val != ".") + path = Path::join (Path::dirname (file), path_val); else - fname = file; + path = file; + } - in.open (fname.c_str(), std::ios::in | std::ios::binary); - if (!in) - throw Exception ("error opening " + type + " data file \"" + fname + "\": " + strerror(errno)); - in.seekg (offset); + + + void properties_consensus( const std::vector& files, const std::string& type, Properties& prop ) + { + // ensure there is at least one file + size_t n_files = files.size(); + if (n_files == 0) + throw Exception("empty list of files for consensus!"); + + // read first file + TrackFileInfo info_ref( files[0], type, prop ); + + // create comments buffer and summary counts + std::unordered_set comments_buffer( prop.comments.begin(), prop.comments.end() ); + uint64_t count = prop.value_or_default("count",0); + uint64_t total_count = prop.value_or_default("total_count",0); + + // iterate over remaining files + TrackFileInfo info_tmp; + Properties prop_tmp; + for ( size_t k=1; k < n_files; ++k ) { + + info_tmp.parse( files[k], type, prop_tmp ); + if ( info_tmp.dtype != info_ref.dtype ) + throw Exception("datatype mismatch between " + type + " files \"" + files[0] + "\" and \"" + files[k] + "\""); + + for (auto& c: prop_tmp.comments) + if ( comments_buffer.find(c) == comments_buffer.end() ) { + prop.comments.push_back(c); + comments_buffer.insert(c); + } + + for (auto& p: prop_tmp) + if (p.first == "count") + count += to(p.second); + else if (p.first == "total_count") + total_count += to(p.second); + else { + auto existing = prop.find (p.first); + if (existing == prop.end()) + prop.insert(p); + else if (p.second != existing->second) { + INFO("variable property \"" + p.first + "\" across " + type + " files"); + existing->second = "variable"; + } + } + } + + // final counts + prop["count"] = str(count); + prop["total_count"] = str(total_count); } } diff --git a/src/dwi/tractography/file_base.h b/src/dwi/tractography/file_base.h index 82dd240a31..9c41f68756 100644 --- a/src/dwi/tractography/file_base.h +++ b/src/dwi/tractography/file_base.h @@ -18,9 +18,10 @@ #define __dwi_tractography_file_base_h__ #include -#include +#include #include "types.h" +#include "datatype.h" #include "file/key_value.h" #include "file/ofstream.h" #include "file/path.h" @@ -36,18 +37,63 @@ namespace MR namespace Tractography { + + // JH: information about .tck or .tsf file + struct TrackFileInfo + { + std::string path, type; + uint64_t offset; + DataType dtype; + + TrackFileInfo() + { clear(); } + + TrackFileInfo( const std::string& file, const std::string& type, Properties& prop ) + { clear(); parse(file,type,prop); } + + + void clear() { + dtype = DataType::Undefined; + offset = 0; + path = ""; + type = ""; + } + + + inline bool empty() + { return type.empty() || path.empty(); } + + + void check_dtype(); + + + void parse( const std::string& file, const std::string& type, Properties& prop ); + + + }; + + + + // combine properties across several .tck files + void properties_consensus( const std::vector& files, const std::string& type, Properties& prop ); + + + //! \cond skip class __ReaderBase__ { NOMEMALIGN public: - ~__ReaderBase__ () { - if (in.is_open()) - in.close(); - } - void open (const std::string& file, const std::string& firstline, Properties& properties); + ~__ReaderBase__ () + { close(); } - void close () { in.close(); } + + void close () + { if (in.is_open()) in.close(); } + + + void open (const std::string& file, const std::string& type, Properties& prop); + protected: @@ -56,6 +102,7 @@ namespace MR }; + template class __WriterBase__ { NOMEMALIGN diff --git a/src/dwi/tractography/properties.cpp b/src/dwi/tractography/properties.cpp new file mode 100644 index 0000000000..d2721bd3ae --- /dev/null +++ b/src/dwi/tractography/properties.cpp @@ -0,0 +1,103 @@ +/* + * Copyright (c) 2008-2018 the MRtrix3 contributors. + * + * This Source Code Form is subject to the terms of the Mozilla Public + * License, v. 2.0. If a copy of the MPL was not distributed with this + * file, you can obtain one at http://mozilla.org/MPL/2.0/ + * + * MRtrix3 is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty + * of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. + * + * For more details, see http://www.mrtrix.org/ + */ + + +#include "dwi/tractography/properties.h" + + +namespace MR { + namespace DWI { + namespace Tractography { + + void Properties::load_ROIs () + { + using iter = std::multimap::const_iterator; + + std::pair range = roi.equal_range ("include"); + for (iter it = range.first; it != range.second; ++it) include.add (it->second); + range = roi.equal_range ("exclude"); + for (iter it = range.first; it != range.second; ++it) exclude.add (it->second); + range = roi.equal_range ("mask"); + for (iter it = range.first; it != range.second; ++it) mask.add (it->second); + + roi.clear(); + } + + + + + + void check_compatible (const Properties& a, const Properties& b, const std::string& type) + { + if ( a.get_step_size() != b.get_step_size() ) + WARN ("mismatch between " + type + " properties (step-size)"); + + for (auto& p: {"source","act","method","threshold"}) + if (a.has(p) && b.has(p) && a.get(p) != b.get(p)) + WARN ("mismatch between " + type + " properties (" + p + ")"); + } + + + + + void check_timestamps (const Properties& a, const Properties& b, const std::string& type) + { + auto stamp_a = a.find ("timestamp"); + auto stamp_b = b.find ("timestamp"); + if (stamp_a == a.end() || stamp_b == b.end()) + throw Exception ("unable to verify " + type + " pair: missing timestamp"); + if (stamp_a->second != stamp_b->second) + throw Exception ("invalid " + type + " combination - timestamps do not match"); + } + + + + + void check_counts (const Properties& a, const Properties& b, const std::string& type, bool abort_on_fail) + { + auto count_a = a.find ("count"); + auto count_b = b.find ("count"); + if ((count_a == a.end()) || (count_b == b.end())) { + std::string mesg = "unable to validate " + type + " pair: missing count field"; + if (abort_on_fail) throw Exception (mesg); + else WARN (mesg); + } + if (to(count_a->second) != to(count_b->second)) { + std::string mesg = type + " files do not contain same number of elements"; + if (abort_on_fail) throw Exception (mesg); + else WARN (mesg); + } + } + + + + + std::ostream& operator<< (std::ostream& stream, const Properties& P) + { + stream << "seeds: " << P.seeds; + stream << "include: " << P.include << ", exclude: " << P.exclude << ", mask: " << P.mask << ", dict: "; + for (auto& Pk: P) + stream << "[ " << Pk.first << ": " << Pk.second << " ], "; + + stream << "comments: "; + for (auto& c: P.comments) + stream << "\"" << c << "\", "; + + return stream; + } + + } + } +} + diff --git a/src/dwi/tractography/properties.h b/src/dwi/tractography/properties.h index b35a61b394..0ba02dc4b2 100644 --- a/src/dwi/tractography/properties.h +++ b/src/dwi/tractography/properties.h @@ -18,6 +18,7 @@ #define __dwi_tractography_properties_h__ #include +#include "mrtrix.h" #include "timer.h" #include "dwi/tractography/roi.h" #include "dwi/tractography/seeding/list.h" @@ -71,78 +72,63 @@ namespace MR comments.clear(); } - template void set (T& variable, const std::string& name) { - if ((*this)[name].empty()) (*this)[name] = str (variable); - else variable = to ((*this)[name]); + bool has(const std::string& name) const { + return this->find(name) != this->end(); } - }; - - - inline void check_timestamps (const Properties& a, const Properties& b, const std::string& type) - { - Properties::const_iterator stamp_a = a.find ("timestamp"); - Properties::const_iterator stamp_b = b.find ("timestamp"); - if (stamp_a == a.end() || stamp_b == b.end()) - throw Exception ("unable to verify " + type + " pair: missing timestamp"); - if (stamp_a->second != stamp_b->second) - throw Exception ("invalid " + type + " combination - timestamps do not match"); - } + std::string get(const std::string& name) const { + auto p = this->find(name); + if (p == this->end()) + throw Exception("Undefined property: \"" + name + "\""); + else + return p->second; + } + template + T value_or_default(const std::string& name, const T& defval) const { + auto p = this->find(name); + return p == this->end() ? defval : to(p->second); + } + template + void set (T& variable, const std::string& name) { + if ((*this)[name].empty()) (*this)[name] = str (variable); + else variable = to ((*this)[name]); + } - inline void check_counts (const Properties& a, const Properties& b, const std::string& type, bool abort_on_fail) - { - Properties::const_iterator count_a = a.find ("count"); - Properties::const_iterator count_b = b.find ("count"); - if ((count_a == a.end()) || (count_b == b.end())) { - std::string mesg = "unable to validate " + type + " pair: missing count field"; - if (abort_on_fail) throw Exception (mesg); - else WARN (mesg); - } - if (to(count_a->second) != to(count_b->second)) { - std::string mesg = type + " files do not contain same number of elements"; - if (abort_on_fail) throw Exception (mesg); - else WARN (mesg); - } - } + void erase_if_present(const std::string& name) { + auto p = this->find(name); + if (p != this->end()) + this->erase(p); + } + inline float get_step_size() const { + return value_or_default( "output_step_size", value_or_default( "step_size", NaN ) ); + } + void load_ROIs (); + }; - inline float get_step_size (const Properties& p) - { - for (size_t index = 0; index != 2; ++index) { - const std::string key (index ? "step_size" : "output_step_size"); - auto it = p.find (key); - if (it != p.end()) { - try { - return to (it->second); - } catch (...) { } - } - } - return NaN; + // JH: kept for now, but should be removed in favour of: p.get_step_size() + inline float get_step_size (const Properties& p) { + DEBUG( "Function get_step_size is deprecated. Call properties.get_step_size() instead." ); + return p.get_step_size(); } + // JH: unused for now + void check_compatible(const Properties& a, const Properties& b, const std::string& type); + void check_timestamps (const Properties& a, const Properties& b, const std::string& type); - inline std::ostream& operator<< (std::ostream& stream, const Properties& P) - { - stream << "seeds: " << P.seeds; - stream << "include: " << P.include << ", exclude: " << P.exclude << ", mask: " << P.mask << ", dict: "; - for (std::map::const_iterator i = P.begin(); i != P.end(); ++i) - stream << "[ " << i->first << ": " << i->second << " ], "; - stream << "comments: "; - for (vector::const_iterator i = P.comments.begin(); i != P.comments.end(); ++i) - stream << "\"" << *i << "\", "; - return (stream); - } + void check_counts (const Properties& a, const Properties& b, const std::string& type, bool abort_on_fail); + std::ostream& operator<< (std::ostream& stream, const Properties& P); } } diff --git a/src/dwi/tractography/roi.h b/src/dwi/tractography/roi.h index 10e4aca859..22fd8128cb 100644 --- a/src/dwi/tractography/roi.h +++ b/src/dwi/tractography/roi.h @@ -19,7 +19,6 @@ #include "app.h" #include "image.h" -#include "image.h" #include "interp/linear.h" #include "math/rng.h" From 1d8fe28ce32bd184f86f9d6ed276275cd57879fd Mon Sep 17 00:00:00 2001 From: Jonathan H Date: Wed, 20 Feb 2019 10:43:07 +0000 Subject: [PATCH 07/11] Fix syntax issue found during #1569 --- core/file/utils.h | 2 +- src/dwi/tractography/file.h | 8 ++++---- src/dwi/tractography/file_base.cpp | 2 +- src/dwi/tractography/file_base.h | 5 +++-- 4 files changed, 9 insertions(+), 8 deletions(-) diff --git a/core/file/utils.h b/core/file/utils.h index 93ed886f0f..2e0b05f7ba 100644 --- a/core/file/utils.h +++ b/core/file/utils.h @@ -172,7 +172,7 @@ namespace MR - inline size_t readlines(const std::string& name, std::vector& lines) + inline size_t readlines(const std::string& name, vector& lines) { std::string line; std::ifstream in(name.c_str()); diff --git a/src/dwi/tractography/file.h b/src/dwi/tractography/file.h index c0dcbeb9f8..d7a74f447e 100644 --- a/src/dwi/tractography/file.h +++ b/src/dwi/tractography/file.h @@ -73,7 +73,7 @@ namespace MR //! open the \c file for reading and load header into \c properties Reader (const std::string& file, Properties& prop) { - std::vector files; + vector files; if (Path::has_suffix(file, ".lst")){ File::readlines( file, files ); INFO("opening list of " + str(files.size()) + " .tck file(s)"); @@ -84,7 +84,7 @@ namespace MR construct( files, prop ); } - Reader(const std::vector& files, Properties& prop) + Reader(const vector& files, Properties& prop) { construct(files, prop); } ~Reader() @@ -126,13 +126,13 @@ namespace MR uint64_t track_index, sub_count; std::ifstream inw; - std::vector tck_files, weights_files; + vector tck_files, weights_files; //! initialise current reader instance - void construct( const std::vector& files, Properties& prop ) + void construct( const vector& files, Properties& prop ) { clear(); diff --git a/src/dwi/tractography/file_base.cpp b/src/dwi/tractography/file_base.cpp index 89c27de51f..711b7809c5 100644 --- a/src/dwi/tractography/file_base.cpp +++ b/src/dwi/tractography/file_base.cpp @@ -105,7 +105,7 @@ namespace MR { - void properties_consensus( const std::vector& files, const std::string& type, Properties& prop ) + void properties_consensus( const vector& files, const std::string& type, Properties& prop ) { // ensure there is at least one file size_t n_files = files.size(); diff --git a/src/dwi/tractography/file_base.h b/src/dwi/tractography/file_base.h index ebdda04107..bb762f1d24 100644 --- a/src/dwi/tractography/file_base.h +++ b/src/dwi/tractography/file_base.h @@ -39,7 +39,8 @@ namespace MR // JH: information about .tck or .tsf file struct TrackFileInfo - { + { NOMEMALIGN + std::string path, type; uint64_t offset; DataType dtype; @@ -74,7 +75,7 @@ namespace MR // combine properties across several .tck files - void properties_consensus( const std::vector& files, const std::string& type, Properties& prop ); + void properties_consensus( const vector& files, const std::string& type, Properties& prop ); From f4164ab6fde5123e6d667fca5bd9be93eb2f5313 Mon Sep 17 00:00:00 2001 From: Jonathan H Date: Wed, 20 Feb 2019 11:25:04 +0000 Subject: [PATCH 08/11] Change roi property name --- src/dwi/tractography/properties.cpp | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/dwi/tractography/properties.cpp b/src/dwi/tractography/properties.cpp index d2721bd3ae..93f6edddaa 100644 --- a/src/dwi/tractography/properties.cpp +++ b/src/dwi/tractography/properties.cpp @@ -24,14 +24,14 @@ namespace MR { { using iter = std::multimap::const_iterator; - std::pair range = roi.equal_range ("include"); + std::pair range = prior_rois.equal_range ("include"); for (iter it = range.first; it != range.second; ++it) include.add (it->second); - range = roi.equal_range ("exclude"); + range = prior_rois.equal_range ("exclude"); for (iter it = range.first; it != range.second; ++it) exclude.add (it->second); - range = roi.equal_range ("mask"); + range = prior_rois.equal_range ("mask"); for (iter it = range.first; it != range.second; ++it) mask.add (it->second); - roi.clear(); + prior_rois.clear(); } From a89e7680bd7fe671f2a82674fdd0bb5e1312bd88 Mon Sep 17 00:00:00 2001 From: Jonathan H Date: Fri, 22 Feb 2019 23:30:02 +0000 Subject: [PATCH 09/11] Remove unnecessary loader, track down potential bug with multi-threaded Editing::Worker/Receiver implementations. Use at instead of [] to throw exception if property does not exist. --- cmd/tckedit.cpp | 56 ++-------------- cmd/tcktransform.cpp | 24 ++----- src/dwi/tractography/editing/loader.h | 88 ------------------------- src/dwi/tractography/editing/worker.cpp | 47 ++++++------- 4 files changed, 31 insertions(+), 184 deletions(-) delete mode 100644 src/dwi/tractography/editing/loader.h diff --git a/cmd/tckedit.cpp b/cmd/tckedit.cpp index 56d3967753..5c7004b2b6 100644 --- a/cmd/tckedit.cpp +++ b/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" @@ -133,59 +132,17 @@ void run () throw Exception ("Cannot use per-streamline weighting with multiple input files"); // Get the consensus streamline properties from among the multiple input files - Tractography::Properties properties; - size_t count = 0; vector input_file_list; - - for (size_t file_index = 0; file_index != num_inputs; ++file_index) { - + for (size_t file_index = 0; file_index != num_inputs; ++file_index) input_file_list.push_back (argument[file_index]); - Properties p; - Reader (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 (i.second); - } else if (i.first == "total_count") { - this_total_count += to (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::Properties properties; + Tractography::Reader reader( input_file_list, properties ); + properties.load_ROIs(); + size_t count = properties.value_or_default("count",0); DEBUG ("estimated number of input tracks: " + str(count)); - load_rois (properties); - // 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 @@ -204,7 +161,6 @@ 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, @@ -212,7 +168,7 @@ void run () Receiver receiver (output_path, properties, number, skip); Thread::run_queue ( - loader, + reader, Thread::batch (Streamline<>()), Thread::multi (worker), Thread::batch (Streamline<>()), diff --git a/cmd/tcktransform.cpp b/cmd/tcktransform.cpp index 575991f4cd..5ba30d74d1 100644 --- a/cmd/tcktransform.cpp +++ b/cmd/tcktransform.cpp @@ -44,23 +44,6 @@ void usage () using value_type = float; using TrackType = Tractography::Streamline; - - -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 reader; -}; - - - class Warper { MEMALIGN(Warper) public: Warper (const Image& warp) : @@ -120,15 +103,16 @@ class Writer { MEMALIGN(Writer) void run () { - Loader loader (argument[0]); + Tractography::Properties properties; + Tractography::Reader reader ( argument[0], properties ); auto data = Image::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), diff --git a/src/dwi/tractography/editing/loader.h b/src/dwi/tractography/editing/loader.h deleted file mode 100644 index 43175f5889..0000000000 --- a/src/dwi/tractography/editing/loader.h +++ /dev/null @@ -1,88 +0,0 @@ -/* Copyright (c) 2008-2019 the MRtrix3 contributors. - * - * This Source Code Form is subject to the terms of the Mozilla Public - * License, v. 2.0. If a copy of the MPL was not distributed with this - * file, You can obtain one at http://mozilla.org/MPL/2.0/. - * - * Covered Software is provided under this License on an "as is" - * basis, without warranty of any kind, either expressed, implied, or - * statutory, including, without limitation, warranties that the - * Covered Software is free of defects, merchantable, fit for a - * particular purpose or non-infringing. - * See the Mozilla Public License v. 2.0 for more details. - * - * For more details, see http://www.mrtrix.org/. - */ - -#ifndef __dwi_tractography_editing_loader_h__ -#define __dwi_tractography_editing_loader_h__ - - -#include - -#include "memory.h" -#include "types.h" - -#include "dwi/tractography/file.h" -#include "dwi/tractography/properties.h" -#include "dwi/tractography/streamline.h" - - -namespace MR { - namespace DWI { - namespace Tractography { - namespace Editing { - - - - - class Loader - { MEMALIGN(Loader) - - public: - Loader (const vector& files) : - file_list (files), - dummy_properties (), - reader (new Reader<> (file_list[0], dummy_properties)), - file_index (0) { } - - bool operator() (Streamline<>&); - - - private: - const vector& file_list; - Properties dummy_properties; - std::unique_ptr > reader; - size_t file_index; - - }; - - - - bool Loader::operator() (Streamline<>& out) - { - out.clear(); - - if ((*reader) (out)) - return true; - - while (++file_index != file_list.size()) { - dummy_properties.clear(); - reader.reset (new Reader<> (file_list[file_index], dummy_properties)); - if ((*reader) (out)) - return true; - } - - return false; - - } - - - - - } - } - } -} - -#endif diff --git a/src/dwi/tractography/editing/worker.cpp b/src/dwi/tractography/editing/worker.cpp index 5340502ca3..f3427cd35e 100644 --- a/src/dwi/tractography/editing/worker.cpp +++ b/src/dwi/tractography/editing/worker.cpp @@ -41,7 +41,6 @@ namespace MR { // Assign to ROIs if (properties.include.size() || properties.exclude.size()) { - include_visited.assign (properties.include.size(), false); if (ends_only) { @@ -130,40 +129,36 @@ namespace MR { - Worker::Thresholds::Thresholds (Tractography::Properties& properties) : - max_length (std::numeric_limits::infinity()), - min_length (0.0f), - max_weight (std::numeric_limits::infinity()), - min_weight (0.0f), - step_size (get_step_size (properties)) + Worker::Thresholds::Thresholds (Tractography::Properties& properties) { - if (properties.find ("max_dist") != properties.end()) { - try { - max_length = to(properties["max_dist"]); - } catch (...) { } + static const float Inf = std::numeric_limits::infinity(); + + step_size = properties.get_step_size(); + + try { + max_length = to(properties.at("max_dist")); + } catch (...) { + max_length = Inf; } - if (properties.find ("min_dist") != properties.end()) { - try { - min_length = to(properties["min_dist"]); - } catch (...) { } + + try { + min_length = to(properties.at("min_dist")); + } catch (...) { + min_length = 0.0f; } if (std::isfinite (step_size)) { // User may set these values to a precise value, which may then fail due to floating-point // calculation of streamline length // Therefore throw a bit of error margin in here - float error_margin = 0.1; - if (properties.find ("downsample_factor") != properties.end()) - error_margin = 0.5 / to(properties["downsample_factor"]); + float error_margin = 0.5f / properties.value_or_default( "downsample_factor", 5.0f ); + max_length += error_margin * step_size; min_length -= error_margin * step_size; } - if (properties.find ("max_weight") != properties.end()) - max_weight = to(properties["max_weight"]); - - if (properties.find ("min_weight") != properties.end()) - min_weight = to(properties["min_weight"]); + max_weight = properties.value_or_default( "max_weight", Inf ); + min_weight = properties.value_or_default( "min_weight", 0.0f ); } @@ -181,11 +176,11 @@ namespace MR { bool Worker::Thresholds::operator() (const Streamline<>& in) const { - const float length = (std::isfinite (step_size) ? in.calc_length (step_size) : in.calc_length()); - return ((length <= max_length) && + const float length = std::isfinite (step_size) ? in.calc_length (step_size) : in.calc_length(); + return (length <= max_length) && (length >= min_length) && (in.weight <= max_weight) && - (in.weight >= min_weight)); + (in.weight >= min_weight); } From a1c62f6041319f3a5219a1c76c643ac28ebb1d36 Mon Sep 17 00:00:00 2001 From: Jonathan H Date: Sat, 23 Feb 2019 00:39:19 +0000 Subject: [PATCH 10/11] Prevent nested .lst files, or mixed lists of .lst and .tck files. --- cmd/tckedit.cpp | 24 +++++++++++++++++------- src/dwi/tractography/file.h | 5 +++++ 2 files changed, 22 insertions(+), 7 deletions(-) diff --git a/cmd/tckedit.cpp b/cmd/tckedit.cpp index 5c7004b2b6..c31bd9ed59 100644 --- a/cmd/tckedit.cpp +++ b/cmd/tckedit.cpp @@ -128,16 +128,26 @@ void run () 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 - vector input_file_list; - for (size_t file_index = 0; file_index != num_inputs; ++file_index) - input_file_list.push_back (argument[file_index]); - + // Either one file: .lst or .tck + // OR multiple .tck files Tractography::Properties properties; - Tractography::Reader reader( input_file_list, properties ); + Tractography::Reader *reader = nullptr; + + if ( num_inputs == 1 ) { + reader = new Tractography::Reader(argument[0], properties); + } else { + vector 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(input_file_list, properties); + } + properties.load_ROIs(); size_t count = properties.value_or_default("count",0); @@ -168,7 +178,7 @@ void run () Receiver receiver (output_path, properties, number, skip); Thread::run_queue ( - reader, + *reader, Thread::batch (Streamline<>()), Thread::multi (worker), Thread::batch (Streamline<>()), diff --git a/src/dwi/tractography/file.h b/src/dwi/tractography/file.h index d73c382088..41afd3c59e 100644 --- a/src/dwi/tractography/file.h +++ b/src/dwi/tractography/file.h @@ -140,6 +140,11 @@ namespace MR if (files.empty()) return; + // check that there are no .lst files + for (auto& file: files) + if (Path::has_suffix(file,".lst")) + throw Exception("Nested .lst files, or lists of files mixing .tck and .lst extensions, are not supported."); + // assign properties (also checks that input files are valid .tck files) properties_consensus(files, "tracks", prop); tck_files = files; From a5e1283293a5a38f37ffcdf46e1fa9596e8ecf0e Mon Sep 17 00:00:00 2001 From: Jonathan H Date: Tue, 26 Feb 2019 10:31:29 +0000 Subject: [PATCH 11/11] Ignore VSCode environment, and tar archives --- .gitignore | 2 ++ 1 file changed, 2 insertions(+) diff --git a/.gitignore b/.gitignore index 090f1149f8..5be2fab59d 100644 --- a/.gitignore +++ b/.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