diff --git a/cmd/tckglobal.cpp b/cmd/tckglobal.cpp index 4cf47f7642..246fd4904e 100644 --- a/cmd/tckglobal.cpp +++ b/cmd/tckglobal.cpp @@ -303,17 +303,13 @@ void run () if (opt.size()) stats.open_stream(opt[0][0]); - auto dwi = header_in.get_image().with_direct_io(3); - ParticleGrid pgrid (dwi); - + ParticleGrid pgrid (header_in); ExternalEnergyComputer* Eext = new ExternalEnergyComputer(stats, header_in, 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. - - + MHSampler mhs (header_in, properties, stats, pgrid, Esum, mask); // All EnergyComputers are recursively destroyed upon destruction of mhs, except for the shared data. INFO("Start MH sampler"); Thread::run (Thread::multi(mhs), "MH sampler"); @@ -345,7 +341,7 @@ void run () // Save fiso, tod and eext - Header header_out (dwi); + Header header_out (header_in); header_out.datatype() = DataType::Float32; opt = get_options("fod"); diff --git a/core/spinlock.h b/core/spinlock.h new file mode 100644 index 0000000000..7973cfb5d6 --- /dev/null +++ b/core/spinlock.h @@ -0,0 +1,31 @@ +/* Copyright (c) 2008-2025 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/. + */ + +#include +#include + +namespace MR { + + class SpinLock + { NOMEMALIGN + public: + void lock() noexcept { while (m_.test_and_set(std::memory_order_acquire)); } + void unlock() noexcept { m_.clear(std::memory_order_release); } + protected: + std::atomic_flag m_ = ATOMIC_FLAG_INIT; + }; + +} diff --git a/src/dwi/tractography/GT/externalenergy.cpp b/src/dwi/tractography/GT/externalenergy.cpp index aa29a87ed1..797666cace 100644 --- a/src/dwi/tractography/GT/externalenergy.cpp +++ b/src/dwi/tractography/GT/externalenergy.cpp @@ -28,7 +28,7 @@ namespace MR { namespace Tractography { namespace GT { - ExternalEnergyComputer::ExternalEnergyComputer(Stats& stat, Header& dwiheader, const Properties& props) + ExternalEnergyComputer::ExternalEnergyComputer(Stats& stat, Header &dwiheader, const Properties& props) : EnergyComputer(stat), dwi(dwiheader.get_image().with_direct_io(3)), T(Transform(dwiheader).scanner2voxel), diff --git a/src/dwi/tractography/GT/gt.h b/src/dwi/tractography/GT/gt.h index 5f9af4ba79..8ee3803c47 100644 --- a/src/dwi/tractography/GT/gt.h +++ b/src/dwi/tractography/GT/gt.h @@ -108,6 +108,7 @@ namespace MR { } double getTint() const { + std::lock_guard lock (mutex); return Tint; } @@ -118,10 +119,12 @@ namespace MR { double getEextTotal() const { + std::lock_guard lock (mutex); return EextTot; } double getEintTotal() const { + std::lock_guard lock (mutex); return EintTot; } @@ -197,7 +200,7 @@ namespace MR { protected: - std::mutex mutex; + mutable std::mutex mutex; double Text, Tint; double EextTot, EintTot; double alpha; diff --git a/src/dwi/tractography/GT/internalenergy.cpp b/src/dwi/tractography/GT/internalenergy.cpp index f52288d2bc..eabf0d7c39 100644 --- a/src/dwi/tractography/GT/internalenergy.cpp +++ b/src/dwi/tractography/GT/internalenergy.cpp @@ -21,8 +21,8 @@ namespace MR { namespace DWI { namespace Tractography { namespace GT { - - + + double InternalEnergyComputer::stageConnect(const ParticleEnd& pe1, ParticleEnd &pe2) { // new @@ -37,31 +37,34 @@ namespace MR { } return dEint / stats.getTint(); } - - - + + + void InternalEnergyComputer::scanNeighbourhood(const Particle* p, const int alpha0, const double currTemp) { neighbourhood.resize(1); normalization = 1.0; - + Point_t ep = p->getEndPoint(alpha0); + if (pGrid.isoutofbounds(ep)) + return; size_t x, y, z; pGrid.pos2xyz(ep, x, y, z); - + float tolerance2 = Particle::L * Particle::L; // distance threshold (particle length), hard coded float costheta = Math::sqrt1_2; // angular threshold (45 degrees), hard coded ParticleEnd pe; float d1, d2, d, ct; - + for (int i = -1; i <= 1; i++) { for (int j = -1; j <= 1; j++) { for (int k = -1; k <= 1; k++) { - const ParticleGrid::ParticleVectorType* pvec = pGrid.at(x+i, y+j, z+k); - if (pvec == NULL) + const ParticleGrid::ParticleContainer* pvec = pGrid.at(x+i, y+j, z+k); + if (!pvec) continue; - - for (ParticleGrid::ParticleVectorType::const_iterator it = pvec->begin(); it != pvec->end(); ++it) + + std::lock_guard lock (pvec->mutex); + for (ParticleGrid::ParticleContainer::const_iterator it = pvec->begin(); it != pvec->end(); ++it) { pe.par = *it; if (pe.par == p) @@ -81,14 +84,14 @@ namespace MR { neighbourhood.push_back(pe); } } - + } } } - + } - - + + ParticleEnd InternalEnergyComputer::pickNeighbour() { double sum = 0.0; @@ -103,8 +106,8 @@ namespace MR { } return pe; } - - + + } } } diff --git a/src/dwi/tractography/GT/mhsampler.cpp b/src/dwi/tractography/GT/mhsampler.cpp index 9a342516de..bcd5fac2f0 100644 --- a/src/dwi/tractography/GT/mhsampler.cpp +++ b/src/dwi/tractography/GT/mhsampler.cpp @@ -60,7 +60,7 @@ namespace MR { void MHSampler::birth() { - //TRACE; + //std::cerr << 'b'; stats.incN('b'); Point_t pos; @@ -85,7 +85,7 @@ namespace MR { void MHSampler::death() { - //TRACE; + //std::cerr << 'd'; stats.incN('d'); Particle* par; @@ -111,7 +111,7 @@ namespace MR { void MHSampler::randshift() { - //TRACE; + //std::cerr << 'r'; stats.incN('r'); Particle* par; @@ -143,7 +143,7 @@ namespace MR { void MHSampler::optshift() { - //TRACE; + //std::cerr << 'o'; stats.incN('o'); Particle* par; @@ -176,7 +176,7 @@ namespace MR { void MHSampler::connect() // TODO Current implementation does not prevent loops. { - //TRACE; + //std::cerr << 'c'; stats.incN('c'); Particle* par; diff --git a/src/dwi/tractography/GT/mhsampler.h b/src/dwi/tractography/GT/mhsampler.h index f79de46ede..97ec07cadc 100644 --- a/src/dwi/tractography/GT/mhsampler.h +++ b/src/dwi/tractography/GT/mhsampler.h @@ -17,6 +17,7 @@ #ifndef __gt_mhsampler_h__ #define __gt_mhsampler_h__ +#include "header.h" #include "image.h" #include "transform.h" @@ -40,77 +41,77 @@ namespace MR { class MHSampler { MEMALIGN(MHSampler) public: - MHSampler(const Image& dwi, Properties &p, Stats &s, ParticleGrid &pgrid, + MHSampler(const Header& dwiheader, Properties &p, Stats &s, ParticleGrid &pgrid, EnergyComputer* e, Image& m) - : props(p), stats(s), pGrid(pgrid), E(e), T(dwi), - dims{size_t(dwi.size(0)), size_t(dwi.size(1)), size_t(dwi.size(2))}, - mask(m), lock(make_shared>(5*Particle::L)), + : props(p), stats(s), pGrid(pgrid), E(e), T(dwiheader), + dims{size_t(dwiheader.size(0)), size_t(dwiheader.size(1)), size_t(dwiheader.size(2))}, + mask(m), lock(make_shared>(std::max(5*Particle::L, float(2*pGrid.spacing())))), sigpos(Particle::L / 8.), sigdir(0.2) { DEBUG("Initialise Metropolis Hastings sampler."); } - + MHSampler(const MHSampler& other) - : props(other.props), stats(other.stats), pGrid(other.pGrid), E(other.E->clone()), + : props(other.props), stats(other.stats), pGrid(other.pGrid), E(other.E->clone()), T(other.T), dims(other.dims), mask(other.mask), lock(other.lock), rng_uniform(), rng_normal(), sigpos(other.sigpos), sigdir(other.sigdir) { DEBUG("Copy Metropolis Hastings sampler."); } - + ~MHSampler() { delete E; } - + void execute(); - + void next(); - + void birth(); void death(); void randshift(); void optshift(); void connect(); - - + + protected: - + Properties& props; Stats& stats; ParticleGrid& pGrid; EnergyComputer* E; // Polymorphic copy requires call to EnergyComputer::clone(), hence references or smart pointers won't do. - + Transform T; vector dims; Image mask; - + std::shared_ptr< SpatialLock > lock; Math::RNG::Uniform rng_uniform; Math::RNG::Normal rng_normal; float sigpos, sigdir; - - + + Point_t getRandPosInMask(); - + bool inMask(const Point_t p); - + Point_t getRandDir(); - + void moveRandom(const Particle* par, Point_t& pos, Point_t& dir); - + bool moveOptimal(const Particle* par, Point_t& pos, Point_t& dir) const; - + inline double calcShiftProb(const Particle* par, const Point_t& pos, const Point_t& dir) const { Point_t Dpos = par->getPosition() - pos; Point_t Ddir = par->getDirection() - dir; return gaussian_pdf(Dpos, sigpos) * gaussian_pdf(Ddir, sigdir); } - + inline double gaussian_pdf(const Point_t& x, double sigma) const { return std::exp( -x.squaredNorm() / (2*sigma) ) / std::sqrt( 2*Math::pi * sigma*sigma); } - - + + }; - + } } diff --git a/src/dwi/tractography/GT/particle.h b/src/dwi/tractography/GT/particle.h index eb6e733d76..80918ce1f7 100644 --- a/src/dwi/tractography/GT/particle.h +++ b/src/dwi/tractography/GT/particle.h @@ -18,28 +18,27 @@ #define __gt_particle_h__ #include "types.h" - - +#include "spinlock.h" namespace MR { namespace DWI { namespace Tractography { namespace GT { - + using Point_t = Eigen::Vector3f; - + /** * A particle is a segment of a track and consists of a position and a direction. */ class Particle { MEMALIGN(Particle) public: - + // Particle length static float L; - + // Constructors and destructor -------------------------------------------------- - + Particle() { predecessor = nullptr; @@ -47,177 +46,226 @@ namespace MR { visited = false; alive = false; } - + Particle(const Point_t& p, const Point_t& d) { init(p, d); } - + ~Particle() { finalize(); } - + inline void init(const Point_t& p, const Point_t& d) { - setPosition(p); - setDirection(d); + const std::lock_guard lock (spinlock); + pos = p; + dir = d; predecessor = nullptr; successor = nullptr; visited = false; alive = true; } - + inline void finalize() { + const std::lock_guard lock (spinlock); if (predecessor) - removePredecessor(); + removePredecessor_nolock(); if (successor) - removeSuccessor(); + removeSuccessor_nolock(); alive = false; } // disable copy and assignment Particle(const Particle&) = delete; Particle& operator=(const Particle&) = delete; - + // default move constructor and assignment - Particle(Particle&&) = default; - Particle& operator=(Particle&&) = default; + Particle(Particle&&) = delete; + Particle& operator=(Particle&&) = delete; + - // Getters and setters ---------------------------------------------------------- - + Point_t getPosition() const { + const std::lock_guard lock (spinlock); return pos; } - + void setPosition(const Point_t& p) { + const std::lock_guard lock (spinlock); pos = p; } - + Point_t getDirection() const { + const std::lock_guard lock (spinlock); return dir; } - + void setDirection(const Point_t& d) { + const std::lock_guard lock (spinlock); dir = d; dir.normalize(); } - + Point_t getEndPoint(const int a) const { + const std::lock_guard lock (spinlock); return (pos + a*L*dir); } - + bool hasPredecessor() const { + const std::lock_guard lock (spinlock); return (predecessor != nullptr); } - + Particle* getPredecessor() const { + const std::lock_guard lock (spinlock); return predecessor; } - + void connectPredecessor(Particle* p1, const int a1) { + const std::lock_guard lock (spinlock); assert(p1 != nullptr); - setPredecessor(p1); + setPredecessor_nolock(p1); if (a1 == 1) p1->setSuccessor(this); if (a1 == -1) p1->setPredecessor(this); } - + + void setPredecessor(Particle* p1) + { + const std::lock_guard lock (spinlock); + setPredecessor_nolock (p1); + } + + void setSuccessor(Particle* p1) + { + const std::lock_guard lock (spinlock); + setSuccessor_nolock (p1); + } + void removePredecessor() { - assert(predecessor != nullptr); - assert(predecessor->predecessor == this || predecessor->successor == this); - if (predecessor->predecessor == this) - predecessor->predecessor = nullptr; - if (predecessor->successor == this) - predecessor->successor = nullptr; - predecessor = nullptr; + const std::lock_guard lock (spinlock); + removePredecessor_nolock(); } - + + void removeSuccessor() + { + const std::lock_guard lock (spinlock); + removeSuccessor_nolock(); + } + bool hasSuccessor() const { + const std::lock_guard lock (spinlock); return (successor != nullptr); } - + Particle* getSuccessor() const { + const std::lock_guard lock (spinlock); return successor; } - + void connectSuccessor(Particle* p1, const int a1) { assert(p1 != nullptr); - setSuccessor(p1); + const std::lock_guard lock (spinlock); + setSuccessor_nolock(p1); if (a1 == 1) p1->setSuccessor(this); if (a1 == -1) p1->setPredecessor(this); } - - void removeSuccessor() - { - assert(successor != nullptr); - assert(successor->predecessor == this || successor->successor == this); - if (successor->predecessor == this) - successor->predecessor = nullptr; - if (successor->successor == this) - successor->successor = nullptr; - successor = nullptr; - } - + bool isVisited() const { + const std::lock_guard lock (spinlock); return visited; } - + void setVisited(const bool v) { + const std::lock_guard lock (spinlock); visited = v; } - + bool isAlive() const { + const std::lock_guard lock (spinlock); return alive; } - + protected: - + + mutable SpinLock spinlock; Point_t pos, dir; Particle* predecessor; Particle* successor; bool visited; bool alive; - - void setPredecessor(Particle* p1) + + void setPredecessor_nolock(Particle* p1) { + if (predecessor == p1) + return; if (predecessor) - removePredecessor(); + removePredecessor_nolock(); predecessor = p1; } - - void setSuccessor(Particle* p1) + + void setSuccessor_nolock(Particle* p1) { + if (successor == p1) + return; if (successor) - removeSuccessor(); + removeSuccessor_nolock(); successor = p1; } - + + void removePredecessor_nolock() + { + assert(predecessor != nullptr); + assert(predecessor != this); + assert(predecessor->predecessor == this || predecessor->successor == this); + const std::lock_guard lock (predecessor->spinlock); + if (predecessor->predecessor == this) + predecessor->predecessor = nullptr; + if (predecessor->successor == this) + predecessor->successor = nullptr; + predecessor = nullptr; + } + + void removeSuccessor_nolock() + { + assert(successor != nullptr); + assert(successor != this); + assert(successor->predecessor == this || successor->successor == this); + const std::lock_guard lock (successor->spinlock); + if (successor->predecessor == this) + successor->predecessor = nullptr; + if (successor->successor == this) + successor->successor = nullptr; + successor = nullptr; + } + }; - - + + /** * Small data structure that refers to one end of a particle. * It is used to represent candidate neighbours of a given particle, @@ -230,9 +278,9 @@ namespace MR { float e_conn; double p_suc; }; - - - + + + } } } diff --git a/src/dwi/tractography/GT/particlegrid.cpp b/src/dwi/tractography/GT/particlegrid.cpp index 638442cb02..8e15e4f848 100644 --- a/src/dwi/tractography/GT/particlegrid.cpp +++ b/src/dwi/tractography/GT/particlegrid.cpp @@ -25,18 +25,20 @@ namespace MR { ParticleGrid::ParticleGrid(const Header& H) { DEBUG("Initialise particle grid."); - dims[0] = Math::ceil( H.size(0) * H.spacing(0) / (2.0*Particle::L) ); - dims[1] = Math::ceil( H.size(1) * H.spacing(1) / (2.0*Particle::L) ); - dims[2] = Math::ceil( H.size(2) * H.spacing(2) / (2.0*Particle::L) ); + // define (isotropic) grid spacing + default_type vox = std::min({H.spacing(0), H.spacing(1), H.spacing(2)}); + grid_spacing = std::max(2.0 * Particle::L, vox); + + // set grid dimensions + dims[0] = Math::ceil( (H.size(0)-1) * H.spacing(0) / grid_spacing ) + 1; + dims[1] = Math::ceil( (H.size(1)-1) * H.spacing(1) / grid_spacing ) + 1; + dims[2] = Math::ceil( (H.size(2)-1) * H.spacing(2) / grid_spacing ) + 1; grid.resize(dims[0]*dims[1]*dims[2]); // Initialise scanner-to-grid transform - Eigen::DiagonalMatrix newspacing (2.0*Particle::L, 2.0*Particle::L, 2.0*Particle::L); - Eigen::Vector3d shift (H.spacing(0)/2.0 - Particle::L, - H.spacing(1)/2.0 - Particle::L, - H.spacing(2)/2.0 - Particle::L); - T_s2g = H.transform() * newspacing; - T_s2g = T_s2g.inverse().translate(shift); + Eigen::DiagonalMatrix newspacing (grid_spacing, grid_spacing, grid_spacing); + transform_type T_g2s = H.transform() * newspacing; + T_s2g = T_g2s.inverse(); } @@ -44,6 +46,7 @@ namespace MR { { Particle* p = pool.create(pos, dir); size_t gidx = pos2idx(pos); + std::lock_guard lock (mutex); grid[gidx].push_back(p); } @@ -51,7 +54,8 @@ namespace MR { { size_t gidx0 = pos2idx(p->getPosition()); size_t gidx1 = pos2idx(pos); - grid[gidx0].erase(std::remove(grid[gidx0].begin(), grid[gidx0].end(), p), grid[gidx0].end()); + std::lock_guard lock (mutex); + grid[gidx0].remove(p); p->setPosition(pos); p->setDirection(dir); grid[gidx1].push_back(p); @@ -60,7 +64,10 @@ namespace MR { void ParticleGrid::remove(Particle* p) { size_t gidx0 = pos2idx(p->getPosition()); - grid[gidx0].erase(std::remove(grid[gidx0].begin(), grid[gidx0].end(), p), grid[gidx0].end()); + { + std::lock_guard lock (mutex); + grid[gidx0].remove (p);// (std::remove(grid[gidx0].begin(), grid[gidx0].end(), p), grid[gidx0].end()); + } pool.destroy(p); } @@ -70,7 +77,7 @@ namespace MR { pool.clear(); } - const ParticleGrid::ParticleVectorType* ParticleGrid::at(const ssize_t x, const ssize_t y, const ssize_t z) const + const ParticleGrid::ParticleContainer* ParticleGrid::at(const ssize_t x, const ssize_t y, const ssize_t z) const { if ((x < 0) || (size_t(x) >= dims[0]) || (y < 0) || (size_t(y) >= dims[1]) || (z < 0) || (size_t(z) >= dims[2])) // out of bounds return nullptr; @@ -86,7 +93,7 @@ namespace MR { int alpha = 0; vector track; // Loop through all unvisited particles - for (ParticleVectorType& gridvox : grid) + for (ParticleContainer& gridvox : grid) { for (Particle* par0 : gridvox) { @@ -126,7 +133,7 @@ namespace MR { } } // Free all particle locks - for (ParticleVectorType& gridvox : grid) { + for (ParticleContainer& gridvox : grid) { for (Particle* par : gridvox) { par->setVisited(false); } diff --git a/src/dwi/tractography/GT/particlegrid.h b/src/dwi/tractography/GT/particlegrid.h index 17fe625d36..6343246ee5 100644 --- a/src/dwi/tractography/GT/particlegrid.h +++ b/src/dwi/tractography/GT/particlegrid.h @@ -18,6 +18,7 @@ #define __gt_particlegrid_h__ #include +#include #include "header.h" #include "transform.h" @@ -40,7 +41,34 @@ namespace MR { { MEMALIGN(ParticleGrid) public: - using ParticleVectorType = vector; + class ParticleContainer + { NOMEMALIGN + public: + ParticleContainer() = default; + ParticleContainer (const ParticleContainer& other) : particles (other.particles) { } + + void push_back (Particle* p) { + std::lock_guard lock (mutex); + particles.push_back (p); + } + void remove (Particle* p) { + std::lock_guard lock (mutex); + particles.erase (std::remove (particles.begin(), particles.end(), p), particles.end()); + } + + using iterator = std::deque::iterator; + using const_iterator = std::deque::const_iterator; + + iterator begin() { return particles.begin(); } + iterator end() { return particles.end(); } + const_iterator begin() const { return particles.begin(); } + const_iterator end() const { return particles.end(); } + + mutable std::mutex mutex; + + private: + std::deque particles; + }; ParticleGrid(const Header&); ParticleGrid(const ParticleGrid&) = delete; @@ -62,7 +90,7 @@ namespace MR { void clear(); - const ParticleVectorType* at(const ssize_t x, const ssize_t y, const ssize_t z) const; + const ParticleContainer* at(const ssize_t x, const ssize_t y, const ssize_t z) const; inline Particle* getRandom() { return pool.random(); @@ -74,10 +102,11 @@ namespace MR { protected: std::mutex mutex; ParticlePool pool; - vector grid; + vector grid; Math::RNG rng; transform_type T_s2g; size_t dims[3]; + default_type grid_spacing; inline size_t pos2idx(const Point_t& pos) const @@ -88,6 +117,13 @@ namespace MR { } public: + inline bool isoutofbounds(const Point_t& pos) const + { + Point_t gpos = T_s2g.cast() * pos; + return (gpos[0] <= -0.5) || (gpos[1] <= -0.5) || (gpos[2] <= -0.5) || + (gpos[0] >= dims[0]-0.5) || (gpos[1] >= dims[1]-0.5) || (gpos[2] >= dims[2]-0.5); + } + inline void pos2xyz(const Point_t& pos, size_t& x, size_t& y, size_t& z) const { Point_t gpos = T_s2g.cast() * pos; @@ -97,6 +133,10 @@ namespace MR { z = Math::round(gpos[2]); } + inline default_type spacing () { + return grid_spacing; + } + protected: inline size_t xyz2idx(const size_t x, const size_t y, const size_t z) const { diff --git a/src/dwi/tractography/GT/particlepool.h b/src/dwi/tractography/GT/particlepool.h index e58c6f3fdd..a3e6a4673c 100644 --- a/src/dwi/tractography/GT/particlepool.h +++ b/src/dwi/tractography/GT/particlepool.h @@ -73,6 +73,7 @@ namespace MR { * @brief Return number of Particles in the pool. */ inline size_t size() const { + std::lock_guard lock (mutex); return pool.size() - avail.size(); } @@ -104,7 +105,7 @@ namespace MR { } protected: - std::mutex mutex; + mutable std::mutex mutex; deque pool; std::stack > avail; Math::RNG rng; diff --git a/src/dwi/tractography/GT/spatiallock.h b/src/dwi/tractography/GT/spatiallock.h index 7046e8cb2b..bfd52b32c5 100644 --- a/src/dwi/tractography/GT/spatiallock.h +++ b/src/dwi/tractography/GT/spatiallock.h @@ -43,6 +43,7 @@ namespace MR { SpatialLock(const value_type tx, const value_type ty, const value_type tz) : _tx(tx), _ty(ty), _tz(tz) { } ~SpatialLock() { + std::lock_guard lock(mutex); lockcentres.clear(); } @@ -61,6 +62,10 @@ namespace MR { { NOMEMALIGN public: Guard(SpatialLock& l) : lock(l), idx(-1) { } + Guard(const Guard&) = delete; + Guard& operator= (const Guard&) = delete; + Guard(Guard&& other) noexcept = delete; + Guard& operator= (Guard&&) = delete; ~Guard() { if (idx >= 0) @@ -113,6 +118,7 @@ namespace MR { } void unlock(const size_t idx) { + std::lock_guard lock(mutex); lockcentres[idx].second = false; } diff --git a/testing/binaries/tests/tckglobal b/testing/binaries/tests/tckglobal new file mode 100644 index 0000000000..cf445e8770 --- /dev/null +++ b/testing/binaries/tests/tckglobal @@ -0,0 +1,3 @@ +dwiextract dwi.mif - | tckglobal - response.txt tmp.tck -fod tmp_fod.mif -eext tmp_eext.mif -etrend tmp_etrend.txt -niter 100000 -force +dwiextract dwi.mif - | tckglobal - response.txt tmp.tck -mask mask.mif -niter 100000 -force +tckglobal dwi2fod/msmt/dwi.mif dwi2fod/msmt/wm.txt tmp.tck -riso dwi2fod/msmt/gm.txt -riso dwi2fod/msmt/csf.txt -fiso tmp_fiso.mif -eext tmp_eext.mif -niter 100000 -force