FIx compilation on Windows with GCC16

cpp/core/dwi/tractography/algorithms/iFOD1.h

@@ -195,7 +195,7 @@ class iFOD1 : public MethodBase {
             max_truncation = val / max_val;
         }
 
-        if (uniform(rng) < val / max_val) {
+        if (uniform(rng()) < val / max_val) {
           dir = new_dir;
           dir.normalize();
           pos += S.step_size * dir;

cpp/core/dwi/tractography/algorithms/iFOD2.h

@@ -245,7 +245,7 @@ class iFOD2 : public MethodBase {
           max_truncation = val / max_val;
       }
 
-      if (uniform(rng) < val / max_val) {
+      if (uniform(rng()) < val / max_val) {
         mean_sample_num += n;
         half_log_prob0 = last_half_log_probN;
         pos = positions[0];

cpp/core/dwi/tractography/algorithms/nulldist.h

@@ -68,7 +68,7 @@ class NullDist1 : public MethodBase {
     return std::nullopt;
   }
 
-  float get_metric(const Eigen::Vector3f &, const Eigen::Vector3f &) override { return uniform(rng); }
+  float get_metric(const Eigen::Vector3f &, const Eigen::Vector3f &) override { return uniform(rng()); }
 
 protected:
   const Shared &S;
@@ -129,7 +129,7 @@ class NullDist2 : public iFOD2 {
     iFOD2::truncate_track(tck, length_to_revert_from, revert_step);
   }
 
-  float get_metric(const Eigen::Vector3f &, const Eigen::Vector3f &) override { return uniform(rng); }
+  float get_metric(const Eigen::Vector3f &, const Eigen::Vector3f &) override { return uniform(rng()); }
 
 protected:
   const Shared &S;

cpp/core/dwi/tractography/algorithms/tensor_prob.h

@@ -89,7 +89,7 @@ class Tensor_Prob : public Tensor_Det {
 
       for (ssize_t i = 0; i < residuals.size(); ++i) {
         residuals[i] = residuals[i] ? (data[i] - std::exp(-residuals[i])) : float(0.0);
-        data[i] += uniform_int(rng) ? residuals[i] : -residuals[i];
+        data[i] += uniform_int(rng()) ? residuals[i] : -residuals[i];
       }
     }
 

cpp/core/dwi/tractography/rng.cpp

@@ -18,6 +18,9 @@
 
 namespace MR::DWI::Tractography {
 
-thread_local Math::RNG rng;
-
+Math::RNG &rng() {
+  static thread_local Math::RNG instance;
+  return instance;
 }
+
+} // namespace MR::DWI::Tractography

cpp/core/dwi/tractography/rng.h

@@ -20,7 +20,12 @@
 
 namespace MR::DWI::Tractography {
 
-//! thread-local, but globally accessible RNG to vastly simplify multi-threading
-extern thread_local Math::RNG rng;
+//! Thread-local, but globally accessible RNG to vastly simplify multi-threading.
+/*! Implemented as a Meyers-style accessor returning a reference to a
+ *  function-local \c thread_local instance, lazily constructed on first use.
+ *  This avoids a namespace-scope \c thread_local definition,
+ *  which triggers compilation failures under GCC 16 on Windows
+ *  due to changes in thread-local storage (TLS) handling. */
+Math::RNG &rng();
 
 } // namespace MR::DWI::Tractography

cpp/core/dwi/tractography/seeding/basic.cpp

@@ -34,7 +34,7 @@ Sphere::Sphere(std::string_view in)                                      //
 bool Sphere::get_seed(Eigen::Vector3f &p) const {
   std::uniform_real_distribution<float> uniform;
   do {
-    p = {2.0F * uniform(rng) - 1.0F, 2.0F * uniform(rng) - 1.0F, 2.0F * uniform(rng) - 1.0F};
+    p = {2.0F * uniform(rng()) - 1.0F, 2.0F * uniform(rng()) - 1.0F, 2.0F * uniform(rng()) - 1.0F};
   } while (p.squaredNorm() > 1.0F);
   p = pos + rad * p;
   return true;
@@ -50,14 +50,14 @@ SeedMask::SeedMask(const std::filesystem::path &in)
 bool SeedMask::get_seed(Eigen::Vector3f &p) const {
   auto seed = mask;
   do {
-    seed.index(0) = std::uniform_int_distribution<int>(0, mask.size(0) - 1)(rng);
-    seed.index(1) = std::uniform_int_distribution<int>(0, mask.size(1) - 1)(rng);
-    seed.index(2) = std::uniform_int_distribution<int>(0, mask.size(2) - 1)(rng);
+    seed.index(0) = std::uniform_int_distribution<int>(0, mask.size(0) - 1)(rng());
+    seed.index(1) = std::uniform_int_distribution<int>(0, mask.size(1) - 1)(rng());
+    seed.index(2) = std::uniform_int_distribution<int>(0, mask.size(2) - 1)(rng());
   } while (!seed.value());
   std::uniform_real_distribution<float> uniform;
-  p = {static_cast<float>(seed.index(0)) + uniform(rng) - 0.5F,
-       static_cast<float>(seed.index(1)) + uniform(rng) - 0.5F,
-       static_cast<float>(seed.index(2)) + uniform(rng) - 0.5F};
+  p = {static_cast<float>(seed.index(0)) + uniform(rng()) - 0.5F,
+       static_cast<float>(seed.index(1)) + uniform(rng()) - 0.5F,
+       static_cast<float>(seed.index(2)) + uniform(rng()) - 0.5F};
   p = (*mask.voxel2scanner) * p;
   return true;
 }
@@ -101,9 +101,9 @@ bool Random_per_voxel::get_seed(Eigen::Vector3f &p) const {
   }
 
   std::uniform_real_distribution<float> uniform;
-  p = {static_cast<float>(mask.index(0)) + uniform(rng) - 0.5F,
-       static_cast<float>(mask.index(1)) + uniform(rng) - 0.5F,
-       static_cast<float>(mask.index(2)) + uniform(rng) - 0.5F};
+  p = {static_cast<float>(mask.index(0)) + uniform(rng()) - 0.5F,
+       static_cast<float>(mask.index(1)) + uniform(rng()) - 0.5F,
+       static_cast<float>(mask.index(2)) + uniform(rng()) - 0.5F};
   p = (*mask.voxel2scanner) * p;
   return true;
 }
@@ -228,22 +228,25 @@ bool Rejection_per_voxel::get_seed(Eigen::Vector3f &p) const {
   float selector;
   Eigen::Vector3f pos;
   do {
-    pos = {
-        uniform(rng) * (interp.size(0) - 1), uniform(rng) * (interp.size(1) - 1), uniform(rng) * (interp.size(2) - 1)};
+    pos = {uniform(rng()) * (interp.size(0) - 1),
+           uniform(rng()) * (interp.size(1) - 1),
+           uniform(rng()) * (interp.size(2) - 1)};
     seed.voxel(pos);
-    selector = rng->Uniform() * max;
+    selector = rng().Uniform() * max;
   } while (seed.value() < selector);
   p = interp.voxel2scanner * pos;
 #else
   auto seed = image;
   float selector;
   do {
-    seed.index(0) = std::uniform_int_distribution<int>(0, image.size(0) - 1)(rng);
-    seed.index(1) = std::uniform_int_distribution<int>(0, image.size(1) - 1)(rng);
-    seed.index(2) = std::uniform_int_distribution<int>(0, image.size(2) - 1)(rng);
-    selector = uniform(rng) * max;
+    seed.index(0) = std::uniform_int_distribution<int>(0, image.size(0) - 1)(rng());
+    seed.index(1) = std::uniform_int_distribution<int>(0, image.size(1) - 1)(rng());
+    seed.index(2) = std::uniform_int_distribution<int>(0, image.size(2) - 1)(rng());
+    selector = uniform(rng()) * max;
   } while (seed.value() < selector);
-  p = {seed.index(0) + uniform(rng) - 0.5f, seed.index(1) + uniform(rng) - 0.5f, seed.index(2) + uniform(rng) - 0.5f};
+  p = {seed.index(0) + uniform(rng()) - 0.5f,
+       seed.index(1) + uniform(rng()) - 0.5f,
+       seed.index(2) + uniform(rng()) - 0.5f};
   p = voxel2scanner * p;
 #endif
   return true;
@@ -314,7 +317,7 @@ Random_coordinates::Random_coordinates(const std::filesystem::path &path) //
 }
 
 bool Random_coordinates::get_seed(Eigen::Vector3f &p) const {
-  p = coords.row(std::uniform_int_distribution<>(0, static_cast<int>(num_coordinates() - 1))(rng));
+  p = coords.row(std::uniform_int_distribution<>(0, static_cast<int>(num_coordinates() - 1))(rng()));
   return true;
 }
 
@@ -334,8 +337,8 @@ bool Rejection_per_coord::get_seed(Eigen::Vector3f &p) const {
   std::uniform_real_distribution<float> selector;
   ssize_t coordinate_index(-1);
   do {
-    coordinate_index = index_selector(rng);
-  } while (weights(coordinate_index) < selector(rng));
+    coordinate_index = index_selector(rng());
+  } while (weights(coordinate_index) < selector(rng()));
   p = coords.row(coordinate_index);
   return true;
 }

cpp/core/dwi/tractography/seeding/dynamic.cpp

@@ -135,7 +135,7 @@ Dynamic::Dynamic(const std::filesystem::path &in,
   // Pick a good fixel to use for testing / debugging
   std::uniform_real_distribution<float> uniform;
   do {
-    test_fixel = std::uniform_int_distribution<int>(1, fixels.size() - 1)(rng);
+    test_fixel = std::uniform_int_distribution<int>(1, fixels.size() - 1)(rng());
   } while (fixels[test_fixel].get_weight() < 1.0 || fixels[test_fixel].get_FOD() < 0.5);
 
   output_fixel_images("begin");
@@ -186,7 +186,7 @@ bool Dynamic::get_seed(Eigen::Vector3f &p, Eigen::Vector3f &d) {
   while (1) {
 
     ++this_attempts;
-    const size_t fixel_index = 1 + uniform_int(rng);
+    const size_t fixel_index = 1 + uniform_int(rng());
     Fixel &fixel = fixels[fixel_index];
     float seed_prob;
     if (fixel.can_update()) {
@@ -223,11 +223,11 @@ bool Dynamic::get_seed(Eigen::Vector3f &p, Eigen::Vector3f &d) {
       seed_prob = fixel.get_old_prob();
     }
 
-    if (seed_prob > uniform_float(rng)) {
+    if (seed_prob > uniform_float(rng())) {
 
       const Eigen::Vector3i &v(fixel.get_voxel());
       const Eigen::Vector3f vp(
-          v[0] + uniform_float(rng) - 0.5, v[1] + uniform_float(rng) - 0.5, v[2] + uniform_float(rng) - 0.5);
+          v[0] + uniform_float(rng()) - 0.5, v[1] + uniform_float(rng()) - 0.5, v[2] + uniform_float(rng()) - 0.5);
       p = transform.voxel2scanner.cast<float>() * vp;
 
       bool good_seed = !act;

cpp/core/dwi/tractography/seeding/gmwmi.cpp

@@ -52,7 +52,8 @@ bool GMWMI::perturb(Eigen::Vector3f &p, Interp &interp) const {
     plane_one = (normal.cross(Eigen::Vector3f(0.0f, 1.0f, 0.0f))).normalized();
   const Eigen::Vector3f plane_two((normal.cross(plane_one)).normalized());
   std::uniform_real_distribution<float> uniform;
-  p += ((uniform(rng) - 0.5f) * perturb_max_step * plane_one) + ((uniform(rng) - 0.5f) * perturb_max_step * plane_two);
+  p += ((uniform(rng()) - 0.5f) * perturb_max_step * plane_one) +
+       ((uniform(rng()) - 0.5f) * perturb_max_step * plane_two);
   return find_interface(p, interp);
 }
 

cpp/core/dwi/tractography/seeding/list.cpp

@@ -59,7 +59,7 @@ bool List::get_seed(Eigen::Vector3f &p, Eigen::Vector3f &d) {
 
     do {
       float incrementer = 0.0;
-      const float sample = uniform(rng) * total_volume;
+      const float sample = uniform(rng()) * total_volume;
       for (auto &i : seeders) {
         incrementer += i->vol();
         if (incrementer > sample)

cpp/core/dwi/tractography/tracking/exec.h

@@ -454,7 +454,7 @@ template <class Method> class Exec {
     }
       return;
     case ACT::sgm_trunc_t::RANDOM:
-      tck.resize(tck.size() - static_cast<size_t>(std::round(method.act().sgm_depth * method.uniform(rng))));
+      tck.resize(tck.size() - static_cast<size_t>(std::round(method.act().sgm_depth * method.uniform(rng()))));
       return;
     case ACT::sgm_trunc_t::ROULETTE: {
       const size_t sgm_start = tck.size() - method.act().sgm_depth;
@@ -483,7 +483,7 @@ template <class Method> class Exec {
         }
       }
       if (invalid_vertices.empty()) {
-        const default_type sample = method.uniform(rng) * total_sum;
+        const default_type sample = method.uniform(rng()) * total_sum;
         default_type cumulative_sum = 0.0;
         while (valid_vertices.size() > 1) {
           cumulative_sum += valid_vertices.front().value();
@@ -493,7 +493,7 @@ template <class Method> class Exec {
         }
         tck.resize(valid_vertices.front().index() + 1);
       } else {
-        const default_type sample = method.uniform(rng) * static_cast<default_type>(invalid_vertices.size());
+        const default_type sample = method.uniform(rng()) * static_cast<default_type>(invalid_vertices.size());
         const size_t truncation_vertex = invalid_vertices[std::trunc(sample)];
         tck.resize(truncation_vertex + 1);
       }

cpp/core/dwi/tractography/tracking/method.cpp

@@ -54,20 +54,20 @@ bool MethodBase::check_seed() {
 Eigen::Vector3f MethodBase::random_direction() {
   Eigen::Vector3f d;
   do {
-    d[0] = 2.0 * uniform(rng) - 1.0;
-    d[1] = 2.0 * uniform(rng) - 1.0;
-    d[2] = 2.0 * uniform(rng) - 1.0;
+    d[0] = 2.0 * uniform(rng()) - 1.0;
+    d[1] = 2.0 * uniform(rng()) - 1.0;
+    d[2] = 2.0 * uniform(rng()) - 1.0;
   } while (d.squaredNorm() > 1.0);
   d.normalize();
   return d;
 }
 
 Eigen::Vector3f MethodBase::random_direction(const float max_angle, const float sin_max_angle) {
-  float phi = 2.0 * Math::pi * uniform(rng);
+  float phi = 2.0 * Math::pi * uniform(rng());
   float theta;
   do {
-    theta = max_angle * uniform(rng);
-  } while (sin_max_angle * uniform(rng) > sin(theta));
+    theta = max_angle * uniform(rng());
+  } while (sin_max_angle * uniform(rng()) > sin(theta));
   return Eigen::Vector3f(sin(theta) * cos(phi), sin(theta) * sin(phi), cos(theta));
 }