tckglobal: Fix operation

Changes in 2be257c made as part of #3011 resulted in there being two separate attempts to construct an Image class from the same input header.
Fixes #3171.

Author: Lestropie

cmd/tckglobal.cpp

@@ -303,15 +303,13 @@ void run ()
   if (opt.size())
     stats.open_stream(opt[0][0]);
 
-  auto dwi = header_in.get_image<float>().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");
@@ -345,7 +343,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");

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<float>().with_direct_io(3)),
             T(Transform(dwiheader).scanner2voxel),

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<float>& dwi, Properties &p, Stats &s, ParticleGrid &pgrid, 
+          MHSampler(const Header& dwiheader, Properties &p, Stats &s, ParticleGrid &pgrid,
                     EnergyComputer* e, Image<bool>& 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<SpatialLock<float>>(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<SpatialLock<float>>(5*Particle::L)),
               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<size_t> dims;
           Image<bool> mask;
-          
+
           std::shared_ptr< SpatialLock<float> > lock;
           Math::RNG::Uniform<float> rng_uniform;
           Math::RNG::Normal<float> 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);
           }
-          
-          
+
+
         };
-        
+
 
       }
     }