Merge branch dev into 'file_mmap_destructor_no_longer_throws'

cmd/dirflip.cpp

@@ -38,7 +38,7 @@ void usage ()
   + "The orientations themselves are not affected, only their "
     "polarity; this is necessary to ensure near-optimal distribution of DW "
     "directions for eddy-current correction.";
-     
+
   ARGUMENTS
     + Argument ("in", "the input files for the directions.").type_file_in()
     + Argument ("out", "the output files for the directions.").type_file_out();
@@ -63,7 +63,7 @@ class Shared { MEMALIGN(Shared)
       progress ("optimising directions for eddy-currents", target_num_permutations),
       best_signs (directions.rows(), 1), best_eddy (std::numeric_limits<value_type>::max()) { }
 
-    bool update (value_type eddy, const vector<int>& signs) 
+    bool update (value_type eddy, const vector<int>& signs)
     {
       std::lock_guard<std::mutex> lock (mutex);
       if (eddy < best_eddy) {
@@ -99,7 +99,7 @@ class Shared { MEMALIGN(Shared)
     vector<int> best_signs;
     value_type best_eddy;
     std::mutex mutex;
-  
+
 };
 
 
@@ -114,7 +114,7 @@ class Processor { MEMALIGN(Processor)
       uniform (0, signs.size()-1) { }
 
     void execute () {
-      while (eval()); 
+      while (eval());
     }
 
 
@@ -128,8 +128,8 @@ class Processor { MEMALIGN(Processor)
       next_permutation();
 
       value_type eddy = 0.0;
-      for (size_t i = 0; i < signs.size(); ++i) 
-        for (size_t j = i+1; j < signs.size(); ++j) 
+      for (size_t i = 0; i < signs.size(); ++i)
+        for (size_t j = i+1; j < signs.size(); ++j)
           eddy += shared.eddy (i, j, signs);
 
       return shared.update (eddy, signs);
@@ -149,18 +149,20 @@ class Processor { MEMALIGN(Processor)
 
 
 
-void run () 
+void run ()
 {
   auto directions = DWI::Directions::load_cartesian (argument[0]);
 
   size_t num_permutations = get_option_value ("permutations", DEFAULT_PERMUTATIONS);
 
-  Shared eddy_shared (directions, num_permutations);
-  Thread::run (Thread::multi (Processor (eddy_shared)), "eval thread");
-
-  auto& signs = eddy_shared.get_best_signs();
+  vector<int> signs;
+  {
+    Shared eddy_shared (directions, num_permutations);
+    Thread::run (Thread::multi (Processor (eddy_shared)), "eval thread");
+    signs = eddy_shared.get_best_signs();
+  }
 
-  for (ssize_t n = 0; n < directions.rows(); ++n) 
+  for (ssize_t n = 0; n < directions.rows(); ++n)
     if (signs[n] < 0)
       directions.row(n) *= -1.0;
 

cmd/tckglobal.cpp

@@ -59,38 +59,39 @@ void usage ()
 {
 
   AUTHOR = "Daan Christiaens ([email protected])";
-  
+
   COPYRIGHT = "Copyright (C) 2015 KU Leuven, Dept. Electrical Engineering, ESAT/PSI,\n"
               "Herestraat 49 box 7003, 3000 Leuven, Belgium \n\n"
 
               "This is free software; see the source for copying conditions.\n"
               "There is NO warranty; not even for MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.";
 
   SYNOPSIS = "Multi-Shell Multi-Tissue Global Tractography";
-  
+
   DESCRIPTION
   + "This command will reconstruct the global white matter fibre tractogram that best "
     "explains the input DWI data, using a multi-tissue spherical convolution model."
-  
+
   + "Example use: "
-  
+
   + " $ tckglobal dwi.mif wmr.txt -riso csfr.txt -riso gmr.txt -mask mask.mif \n"
-    "   -niter 1e9 -fod fod.mif -fiso fiso.mif tracks.tck "
-
+    "   -niter 1e8 -fod fod.mif -fiso fiso.mif tracks.tck "
+
+
   + "in which dwi.mif is the input image, wmr.txt is an anisotropic, multi-shell response function for WM, "
     "and csfr.txt and gmr.txt are isotropic response functions for CSF and GM. The output tractogram is "
     "saved to tracks.tck. Optional output images fod.mif and fiso.mif contain the predicted WM fODF and "
     "isotropic tissue fractions of CSF and GM respectively, estimated as part of the global optimization "
     "and thus affected by spatial regularization.";
-  
+
   REFERENCES
   + "Christiaens, D.; Reisert, M.; Dhollander, T.; Sunaert, S.; Suetens, P. & Maes, F. " // Internal
     "Global tractography of multi-shell diffusion-weighted imaging data using a multi-tissue model. "
     "NeuroImage, 2015, 123, 89-101";
 
   ARGUMENTS
   + Argument ("source", "the image containing the raw DWI data.").type_image_in()
-  
+
   + Argument ("response", "the response of a track segment on the DWI signal.").type_file_in()
 
   + Argument ("tracks", "the output file containing the tracks generated.").type_tracks_out();
@@ -99,7 +100,7 @@ void usage ()
 
   OPTIONS
   + OptionGroup("Input options")
-  
+
   + Option ("grad", "specify the diffusion encoding scheme (required if not supplied in the header).")
     + Argument ("scheme").type_file_in()
 
@@ -111,13 +112,13 @@ void usage ()
 
 
   + OptionGroup("Parameters")
-  
+
   + Option ("lmax", "set the maximum harmonic order for the output series. (default = " + str(DEFAULT_LMAX) + ")")
     + Argument ("order").type_integer (2, 30)
 
   + Option ("length", "set the length of the particles (fibre segments). (default = " + str(DEFAULT_LENGTH, 2) + "mm)")
     + Argument ("size").type_float (1e-6)
-      
+
   + Option ("weight", "set the weight by which particles contribute to the model. (default = " + str(DEFAULT_WEIGHT, 2) + ")")
     + Argument ("w").type_float(1e-6, 1.0)
 
@@ -162,7 +163,7 @@ void usage ()
 
 
   + OptionGroup("Advanced parameters, if you really know what you're doing")
-  
+
   + Option ("balance", "balance internal and external energy. (default = " + str(DEFAULT_BALANCE, 2) + ")\n"
             "Negative values give more weight to the internal energy, positive to the external energy.")
     + Argument ("b").type_float(-100.0, 100.0)
@@ -191,7 +192,7 @@ void usage ()
 template<typename T>
 class __copy_fod { MEMALIGN(__copy_fod<T>)
   public:
-    __copy_fod (const int lmax, const double weight, const bool apodise) 
+    __copy_fod (const int lmax, const double weight, const bool apodise)
       : w(weight), a(apodise), apo (lmax), SH_in (Math::SH::NforL(lmax)), SH_out (SH_in.size()) { }
 
     void operator() (Image<T>& in, Image<T>& out) {
@@ -214,23 +215,23 @@ void run ()
 {
 
   using namespace DWI::Tractography::GT;
-  
+
   // Inputs -----------------------------------------------------------------------------
-  
+
   auto dwi = Image<float>::open(argument[0]).with_direct_io(3);
-  
+
   Properties properties;
   properties.resp_WM = load_matrix<float> (argument[1]);
   double wmscale2 = (properties.resp_WM(0,0)*properties.resp_WM(0,0))/M_4PI;
-  
+
   Eigen::VectorXf riso;
   auto opt = get_options("riso");
   for (auto popt : opt)
   {
     riso = load_vector<float>(popt[0]);
     properties.resp_ISO.push_back(riso);
   }
-  
+
   auto mask = Image<bool>();
   opt = get_options("mask");
   if (opt.size()) {
@@ -254,7 +255,7 @@ void run ()
   properties.lam_ext = 1.0;
   properties.lam_int = 1.0;
   properties.beta = get_option_value("beta", DEFAULT_BETA);
-  
+
   opt = get_options("balance");
   if (opt.size())
   {
@@ -263,7 +264,7 @@ void run ()
     properties.lam_ext = 2*b;
     properties.lam_int = 2*(1-b);
   }
-  
+
   opt = get_options("prob");
   if (opt.size())
   {
@@ -278,14 +279,14 @@ void run ()
       throw Exception("Specified list of proposal probabilities is invalid.");
     }
   }
-  
+
   uint64_t niter = get_option_value("niter", DEFAULT_NITER);
   double t0 = get_option_value("t0", DEFAULT_T0);
   double t1 = get_option_value("t1", DEFAULT_T1);
-  
+
   double mu = get_option_value("ppot", DEFAULT_PPOT);
   properties.ppot = mu * wmscale2 * properties.weight;
-  
+
   opt = get_options("lambda");
   if (opt.size()) {
     properties.lam_ext = 1.0;
@@ -294,36 +295,35 @@ void run ()
 
 
   // Prepare data structures ------------------------------------------------------------
-  
+
   INFO("Initialise data structures for global tractography.");
-  
+
   Stats stats (t0, t1, niter);
   opt = get_options("etrend");
   if (opt.size())
     stats.open_stream(opt[0][0]);
-  
+
   ParticleGrid pgrid (dwi);
-  
+
   ExternalEnergyComputer* Eext = new ExternalEnergyComputer(stats, dwi, properties);
   InternalEnergyComputer* Eint = new InternalEnergyComputer(stats, pgrid);
   Eint->setConnPot(cpot);
   EnergySumComputer* Esum = new EnergySumComputer(stats, Eint, properties.lam_int, Eext, properties.lam_ext / ( wmscale2 * properties.weight*properties.weight));
-  
+
   MHSampler mhs (dwi, properties, stats, pgrid, Esum, mask);   // All EnergyComputers are recursively destroyed upon destruction of mhs, except for the shared data.
-  
-  
+
+
   INFO("Start MH sampler");
-
-  auto t = Thread::run (Thread::multi(mhs), "MH sampler");
-  t.wait();
-
+
+  Thread::run (Thread::multi(mhs), "MH sampler");
+
   INFO("Final no. particles: " + std::to_string(pgrid.getTotalCount()));
   INFO("Final external energy: " + std::to_string(stats.getEextTotal()));
   INFO("Final internal energy: " + std::to_string(stats.getEintTotal()));
-  
-  
+
+
   // Copy results to output buffers -----------------------------------------------------
-  
+
   // Save the tracks (output)
   INFO("Saving tracks to file");
   MR::DWI::Tractography::Properties ftfileprops;
@@ -338,15 +338,15 @@ void run ()
   ftfileprops.comments.push_back("no. iterations = " + std::to_string((long long int) niter));
   ftfileprops.comments.push_back("T0 = " + std::to_string((long double) t0));
   ftfileprops.comments.push_back("T1 = " + std::to_string((long double) t1));
-  
+
   MR::DWI::Tractography::Writer<float> writer (argument[2], ftfileprops);
   pgrid.exportTracks(writer);
-  
-  
+
+
   // Save fiso, tod and eext
   Header header (dwi);
   header.datatype() = DataType::Float32;
-  
+
   opt = get_options("fod");
   if (opt.size()) {
     INFO("Saving fODF image to file");
@@ -355,7 +355,7 @@ void run ()
     auto f = __copy_fod<float>(properties.Lmax, properties.weight, !get_options("noapo").size());
     ThreadedLoop(Eext->getTOD(), 0, 3).run(f, Eext->getTOD(), fODF);
   }
-  
+
   opt = get_options("fiso");
   if (opt.size()) {
     if (properties.resp_ISO.size() > 0) {
@@ -368,16 +368,16 @@ void run ()
       WARN("Ignore saving file " + opt[0][0] + ", because no isotropic response functions were provided.");
     }
   }
-  
+
   opt = get_options("eext");
   if (opt.size()) {
     INFO("Saving external energy to file");
     header.ndim() = 3;
     auto EextI = Image<float>::create (opt[0][0], header);
     threaded_copy(Eext->getEext(), EextI);
   }
-  
-  
+
+
 
 
 }

core/algo/random_threaded_loop.h

@@ -17,6 +17,7 @@
 #define __algo_random_threaded_loop_h__
 
 #include "debug.h"
+#include "exception.h"
 #include "algo/loop.h"
 #include "algo/iterator.h"
 #include "thread.h"
@@ -178,8 +179,7 @@ namespace MR
               }
             } loop_thread = { shared, functor };
 
-            auto t = Thread::run (Thread::multi (loop_thread), "loop threads");
-            t.wait();
+            Thread::run (Thread::multi (loop_thread), "loop threads");
           }
 
 
@@ -194,6 +194,7 @@ namespace MR
               typename std::remove_reference<ImageType>::type...
                 > loop_thread (outer_loop.axes, inner_axes, functor, voxel_density, dimensions, vox...);
             run_outer (loop_thread, voxel_density, dimensions);
+            check_app_exit_code();
           }
 
       };