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 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..c31bd9ed59 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" @@ -122,15 +121,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 () { @@ -138,71 +128,39 @@ 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 + // Either one file: .lst or .tck + // OR multiple .tck files Tractography::Properties properties; - size_t count = 0; - vector 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 (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::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); } - DEBUG ("estimated number of input tracks: " + str(count)); + properties.load_ROIs(); - load_rois (properties); + size_t count = properties.value_or_default("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,7 +171,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, @@ -221,7 +178,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/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/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/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..23a8044f4b 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, 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/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); } diff --git a/src/dwi/tractography/file.h b/src/dwi/tractography/file.h index 43c2e5712a..41afd3c59e 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) + { + 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 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,206 @@ 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 + vector tck_files, weights_files; + + + + + //! initialise current reader instance + void construct( const vector& files, Properties& prop ) + { + clear(); + + 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; + + // 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; + + Properties dummy; + open( tck_files[file_index], "tracks", dummy ); - Eigen::Matrix get_next_point () + 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..1420c8ad73 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 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..59e5b71db7 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,64 @@ namespace MR namespace Tractography { + + // JH: information about .tck or .tsf file + struct TrackFileInfo + { NOMEMALIGN + + 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 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 +103,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..93f6edddaa --- /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 = prior_rois.equal_range ("include"); + for (iter it = range.first; it != range.second; ++it) include.add (it->second); + range = prior_rois.equal_range ("exclude"); + for (iter it = range.first; it != range.second; ++it) exclude.add (it->second); + range = prior_rois.equal_range ("mask"); + for (iter it = range.first; it != range.second; ++it) mask.add (it->second); + + prior_rois.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"