[med-svn] [Git][med-team/mrtrix3][master] Fix compatiblity with Eigen 3.4
Julien Lamy (@lamy-guest)
gitlab at salsa.debian.org
Wed Dec 22 13:49:22 GMT 2021
Julien Lamy pushed to branch master at Debian Med / mrtrix3
Commits:
7e63f69c by Julien Lamy at 2021-12-22T14:48:50+01:00
Fix compatiblity with Eigen 3.4
- - - - -
3 changed files:
- debian/changelog
- + debian/patches/eigen-3.4.patch
- debian/patches/series
Changes:
=====================================
debian/changelog
=====================================
@@ -1,7 +1,11 @@
mrtrix3 (3.0.3-2) UNRELEASED; urgency=medium
+ [ Andreas Tille ]
* Team upload.
* Fix watch file
+
+ [ Julien Lamy ]
+ * Fix compatibility with Eigen 3.4 Closes: #1002203
-- Andreas Tille <tille at debian.org> Wed, 22 Dec 2021 11:09:05 +0100
=====================================
debian/patches/eigen-3.4.patch
=====================================
@@ -0,0 +1,2752 @@
+Description: Fix compatibility with Eigen 3.4
+Origin: https://github.com/MRtrix3/mrtrix3
+Index: mrtrix3/cmd/amp2response.cpp
+===================================================================
+--- mrtrix3.orig/cmd/amp2response.cpp
++++ mrtrix3/cmd/amp2response.cpp
+@@ -82,7 +82,7 @@ void usage ()
+
+
+
+-Eigen::Matrix<default_type, 3, 3> gen_rotation_matrix (const Eigen::Vector3& dir)
++Eigen::Matrix<default_type, 3, 3> gen_rotation_matrix (const Eigen::Vector3d& dir)
+ {
+ static Math::RNG::Normal<default_type> rng;
+ // Generates a matrix that will rotate a unit vector into a new frame of reference,
+@@ -91,11 +91,11 @@ Eigen::Matrix<default_type, 3, 3> gen_ro
+ // Here the other two axes are determined at random (but both are orthogonal to the FOD peak direction)
+ Eigen::Matrix<default_type, 3, 3> R;
+ R (2, 0) = dir[0]; R (2, 1) = dir[1]; R (2, 2) = dir[2];
+- Eigen::Vector3 vec2 (rng(), rng(), rng());
++ Eigen::Vector3d vec2 (rng(), rng(), rng());
+ vec2 = dir.cross (vec2);
+ vec2.normalize();
+ R (0, 0) = vec2[0]; R (0, 1) = vec2[1]; R (0, 2) = vec2[2];
+- Eigen::Vector3 vec3 = dir.cross (vec2);
++ Eigen::Vector3d vec3 = dir.cross (vec2);
+ vec3.normalize();
+ R (1, 0) = vec3[0]; R (1, 1) = vec3[1]; R (1, 2) = vec3[2];
+ return R;
+@@ -155,7 +155,7 @@ class Accumulator { MEMALIGN(Accumulator
+ ++count;
+
+ // Grab the fibre direction
+- Eigen::Vector3 fibre_dir;
++ Eigen::Vector3d fibre_dir;
+ for (dir_image.index(3) = 0; dir_image.index(3) != 3; ++dir_image.index(3))
+ fibre_dir[dir_image.index(3)] = dir_image.value();
+ fibre_dir.normalize();
+Index: mrtrix3/cmd/fixel2tsf.cpp
+===================================================================
+--- mrtrix3.orig/cmd/fixel2tsf.cpp
++++ mrtrix3/cmd/fixel2tsf.cpp
+@@ -101,7 +101,7 @@ void run ()
+ DWI::Tractography::TrackScalar<float> scalars;
+
+ const Transform transform (in_index_image);
+- Eigen::Vector3 voxel_pos;
++ Eigen::Vector3d voxel_pos;
+
+ while (reader (tck)) {
+ SetVoxelDir dixels;
+Index: mrtrix3/cmd/fixel2voxel.cpp
+===================================================================
+--- mrtrix3.orig/cmd/fixel2voxel.cpp
++++ mrtrix3/cmd/fixel2voxel.cpp
+@@ -398,16 +398,16 @@ class DEC_unit : protected Base
+
+ void operator() (Image<index_type>& index, Image<float>& out)
+ {
+- Eigen::Vector3 sum_dec = {0.0, 0.0, 0.0};
++ Eigen::Vector3d sum_dec = {0.0, 0.0, 0.0};
+ if (vol.valid()) {
+ for (auto f = Base::Loop (index) (data, vol, dir); f; ++f) {
+ if (!f.padding())
+- sum_dec += Eigen::Vector3 (abs (dir.row(1)[0]), abs (dir.row(1)[1]), abs (dir.row(1)[2])) * data.value() * vol.value();
++ sum_dec += Eigen::Vector3d (abs (dir.row(1)[0]), abs (dir.row(1)[1]), abs (dir.row(1)[2])) * data.value() * vol.value();
+ }
+ } else {
+ for (auto f = Base::Loop (index) (data, dir); f; ++f) {
+ if (!f.padding())
+- sum_dec += Eigen::Vector3 (abs (dir.row(1)[0]), abs (dir.row(1)[1]), abs (dir.row(1)[2])) * data.value();
++ sum_dec += Eigen::Vector3d (abs (dir.row(1)[0]), abs (dir.row(1)[1]), abs (dir.row(1)[2])) * data.value();
+ }
+ }
+ if ((sum_dec.array() != 0.0).any())
+@@ -430,13 +430,13 @@ class DEC_scaled : protected Base
+
+ void operator() (FixelIndexType& index, Image<float>& out)
+ {
+- Eigen::Vector3 sum_dec = {0.0, 0.0, 0.0};
++ Eigen::Vector3d sum_dec = {0.0, 0.0, 0.0};
+ default_type sum_value = 0.0;
+ if (vol.valid()) {
+ default_type sum_volume = 0.0;
+ for (auto f = Base::Loop (index) (data, vol, dir); f; ++f) {
+ if (!f.padding()) {
+- sum_dec += Eigen::Vector3 (abs (dir.row(1)[0]), abs (dir.row(1)[1]), abs (dir.row(1)[2])) * data.value() * vol.value();
++ sum_dec += Eigen::Vector3d (abs (dir.row(1)[0]), abs (dir.row(1)[1]), abs (dir.row(1)[2])) * data.value() * vol.value();
+ sum_volume += vol.value();
+ sum_value += vol.value() * data.value();
+ }
+@@ -447,7 +447,7 @@ class DEC_scaled : protected Base
+ } else {
+ for (auto f = Base::Loop (index) (data, dir); f; ++f) {
+ if (!f.padding()) {
+- sum_dec += Eigen::Vector3 (abs (dir.row(1)[0]), abs (dir.row(1)[1]), abs (dir.row(1)[2])) * data.value();
++ sum_dec += Eigen::Vector3d (abs (dir.row(1)[0]), abs (dir.row(1)[1]), abs (dir.row(1)[2])) * data.value();
+ sum_value += data.value();
+ }
+ }
+Index: mrtrix3/cmd/labelconvert.cpp
+===================================================================
+--- mrtrix3.orig/cmd/labelconvert.cpp
++++ mrtrix3/cmd/labelconvert.cpp
+@@ -150,7 +150,7 @@ void run ()
+ Transform transform (out);
+ Interp::Nearest<decltype(in_spine)> nearest (in_spine);
+ for (auto l = Loop (out) (out); l; ++l) {
+- Eigen::Vector3 p (out.index (0), out.index (1), out.index (2));
++ Eigen::Vector3d p (out.index (0), out.index (1), out.index (2));
+ p = transform.voxel2scanner * p;
+ if (nearest.scanner (p) && nearest.value())
+ out.value() = spine_index;
+Index: mrtrix3/cmd/labelstats.cpp
+===================================================================
+--- mrtrix3.orig/cmd/labelstats.cpp
++++ mrtrix3/cmd/labelstats.cpp
+@@ -74,7 +74,7 @@ void run ()
+ coms.conservativeResizeLike (matrix_type::Zero (value, 3));
+ masses.conservativeResizeLike (vector_type::Zero (value));
+ }
+- coms.row(value-1) += Eigen::Vector3 (image.index(0), image.index(1), image.index(2));
++ coms.row(value-1) += Eigen::Vector3d (image.index(0), image.index(1), image.index(2));
+ masses[value-1]++;
+ }
+ }
+Index: mrtrix3/cmd/mrcentroid.cpp
+===================================================================
+--- mrtrix3.orig/cmd/mrcentroid.cpp
++++ mrtrix3/cmd/mrcentroid.cpp
+@@ -59,18 +59,18 @@ void run ()
+ check_dimensions (image, mask);
+ }
+
+- Eigen::Vector3 com (0.0, 0.0, 0.0);
++ Eigen::Vector3d com (0.0, 0.0, 0.0);
+ default_type mass = 0.0;
+ if (mask.valid()) {
+ for (auto l = Loop(image) (image, mask); l; ++l) {
+ if (mask.value()) {
+- com += Eigen::Vector3 (image.index(0), image.index(1), image.index(2)) * image.value();
++ com += Eigen::Vector3d (image.index(0), image.index(1), image.index(2)) * image.value();
+ mass += image.value();
+ }
+ }
+ } else {
+ for (auto l = Loop(image) (image); l; ++l) {
+- com += Eigen::Vector3 (image.index(0), image.index(1), image.index(2)) * image.value();
++ com += Eigen::Vector3d (image.index(0), image.index(1), image.index(2)) * image.value();
+ mass += image.value();
+ }
+ }
+Index: mrtrix3/cmd/mrcolour.cpp
+===================================================================
+--- mrtrix3.orig/cmd/mrcolour.cpp
++++ mrtrix3/cmd/mrcolour.cpp
+@@ -86,7 +86,7 @@ void run ()
+ {
+ Header H_in = Header::open (argument[0]);
+ const ColourMap::Entry colourmap = ColourMap::maps[argument[1]];
+- Eigen::Vector3 fixed_colour (NaN, NaN, NaN);
++ Eigen::Vector3d fixed_colour (NaN, NaN, NaN);
+ if (colourmap.is_colour) {
+ if (!(H_in.ndim() == 3 || (H_in.ndim() == 4 && H_in.size(3) == 1)))
+ throw Exception ("For applying a fixed colour, command expects a 3D image as input");
+@@ -96,7 +96,7 @@ void run ()
+ const auto values = parse_floats (opt[0][0]);
+ if (values.size() != 3)
+ throw Exception ("Target colour must be specified as a comma-separated list of three values");
+- fixed_colour = Eigen::Vector3 (values.data());
++ fixed_colour = Eigen::Vector3d (values.data());
+ if (fixed_colour.minCoeff() < 0.0)
+ throw Exception ("Values for fixed colour provided via -colour option cannot be negative");
+ } else if (colourmap.is_rgb) {
+Index: mrtrix3/cmd/mrconvert.cpp
+===================================================================
+--- mrtrix3.orig/cmd/mrconvert.cpp
++++ mrtrix3/cmd/mrconvert.cpp
+@@ -281,8 +281,8 @@ void permute_slice_direction (Header& H,
+ auto it = H.keyval().find ("SliceEncodingDirection");
+ if (it == H.keyval().end())
+ return;
+- const Eigen::Vector3 orig_dir = Axes::id2dir (it->second);
+- const Eigen::Vector3 new_dir (orig_dir[axes[0]], orig_dir[axes[1]], orig_dir[axes[2]]);
++ const Eigen::Vector3d orig_dir = Axes::id2dir (it->second);
++ const Eigen::Vector3d new_dir (orig_dir[axes[0]], orig_dir[axes[1]], orig_dir[axes[2]]);
+ it->second = Axes::dir2id (new_dir);
+ }
+
+Index: mrtrix3/cmd/mrdegibbs.cpp
+===================================================================
+--- mrtrix3.orig/cmd/mrdegibbs.cpp
++++ mrtrix3/cmd/mrdegibbs.cpp
+@@ -122,12 +122,12 @@ class ComputeSlice
+ assign_pos_of (pos, outer_axes).to (in, out);
+
+ for (auto l = Loop (slice_axes) (in); l; ++l)
+- im1 (in.index(X), in.index(Y)) = cdouble (in.value(), 0.0);
++ im1 (ssize_t(in.index(X)), ssize_t(in.index(Y))) = cdouble (in.value(), 0.0);
+
+ unring_2d ();
+
+ for (auto l = Loop (slice_axes) (out); l; ++l)
+- out.value() = im1 (out.index(X), out.index(Y)).real();
++ out.value() = im1 (ssize_t(out.index(X)), ssize_t(out.index(Y))).real();
+ }
+
+ private:
+@@ -331,7 +331,7 @@ void run ()
+ auto slice_encoding_it = header.keyval().find ("SliceEncodingDirection");
+ if (slice_encoding_it != header.keyval().end()) {
+ try {
+- const Eigen::Vector3 slice_encoding_axis_onehot = Axes::id2dir (slice_encoding_it->second);
++ const Eigen::Vector3d slice_encoding_axis_onehot = Axes::id2dir (slice_encoding_it->second);
+ vector<size_t> auto_slice_axes = { 0, 0 };
+ if (slice_encoding_axis_onehot[0])
+ auto_slice_axes = { 1, 2 };
+Index: mrtrix3/cmd/mredit.cpp
+===================================================================
+--- mrtrix3.orig/cmd/mredit.cpp
++++ mrtrix3/cmd/mredit.cpp
+@@ -78,7 +78,7 @@ class Vox : public Eigen::Array3i
+ { MEMALIGN (Vox)
+ public:
+ using Eigen::Array3i::Array3i;
+- Vox (const Eigen::Vector3& p) :
++ Vox (const Eigen::Vector3d& p) :
+ Eigen::Array3i { int(std::round (p[0])), int(std::round (p[1])), int(std::round (p[2])) } { }
+ bool operator< (const Vox& i) const {
+ return (i[0] == (*this)[0] ? (i[1] == (*this)[1] ? (i[2] < (*this)[2]) : (i[1] < (*this)[1])) : (i[0] < (*this)[0]));
+@@ -146,12 +146,12 @@ void run ()
+ operation_count += opt.size();
+ for (auto s : opt) {
+ const auto position = parse_floats (s[0]);
+- Eigen::Vector3 centre_scannerspace (position[0], position[1], position[2]);
++ Eigen::Vector3d centre_scannerspace (position[0], position[1], position[2]);
+ const default_type radius = s[1];
+ const float value = s[2];
+ if (position.size() != 3)
+ throw Exception ("Centre of sphere must be defined using 3 comma-separated values");
+- Eigen::Vector3 centre_voxelspace (centre_scannerspace);
++ Eigen::Vector3d centre_voxelspace (centre_scannerspace);
+ if (scanner)
+ centre_voxelspace = transform.scanner2voxel * centre_scannerspace;
+ else
+@@ -164,7 +164,7 @@ void run ()
+ while (to_expand.size()) {
+ const Vox v (to_expand.back());
+ to_expand.pop_back();
+- const Eigen::Vector3 v_scanner = transform.voxel2scanner * v.matrix().cast<default_type>();
++ const Eigen::Vector3d v_scanner = transform.voxel2scanner * v.matrix().cast<default_type>();
+ const default_type distance = (v_scanner - centre_scannerspace).norm();
+ if (distance < radius) {
+ if (!is_out_of_bounds (H, v)) {
+@@ -190,7 +190,7 @@ void run ()
+ if (position.size() != H.ndim())
+ throw Exception ("Image has " + str(H.ndim()) + " dimensions, but -voxel option position " + std::string(v[0]) + " provides only " + str(position.size()) + " coordinates");
+ if (scanner) {
+- Eigen::Vector3 p (position[0], position[1], position[2]);
++ Eigen::Vector3d p (position[0], position[1], position[2]);
+ p = transform.scanner2voxel * p;
+ const Vox voxel (p);
+ assign_pos_of (voxel).to (out);
+Index: mrtrix3/cmd/mrregister.cpp
+===================================================================
+--- mrtrix3.orig/cmd/mrregister.cpp
++++ mrtrix3/cmd/mrregister.cpp
+@@ -368,7 +368,7 @@ void run () {
+ bool init_rigid_matrix_set = false;
+ if (opt.size()) {
+ init_rigid_matrix_set = true;
+- Eigen::Vector3 centre;
++ Eigen::Vector3d centre;
+ transform_type rigid_transform = load_transform (opt[0][0], centre);
+ rigid.set_transform (rigid_transform);
+ if (!std::isfinite(centre(0))) {
+@@ -518,7 +518,7 @@ void run () {
+ throw Exception ("you cannot initialise with -affine_init_matrix since a rigid registration is being performed");
+
+ init_affine_matrix_set = true;
+- Eigen::Vector3 centre;
++ Eigen::Vector3d centre;
+ transform_type affine_transform = load_transform (opt[0][0], centre);
+ affine.set_transform (affine_transform);
+ if (!std::isfinite(centre(0))) {
+Index: mrtrix3/cmd/mrtransform.cpp
+===================================================================
+--- mrtrix3.orig/cmd/mrtransform.cpp
++++ mrtrix3/cmd/mrtransform.cpp
+@@ -471,7 +471,7 @@ void run ()
+ "therefore should not be used to reorient directions / diffusion gradients");
+ }
+ for (ssize_t n = 0; n < grad.rows(); ++n) {
+- Eigen::Vector3 grad_vector = grad.block<1,3>(n,0);
++ Eigen::Vector3d grad_vector = grad.block<1,3>(n,0);
+ grad.block<1,3>(n,0) = rotation * grad_vector;
+ }
+ DWI::set_DW_scheme (output_header, grad);
+@@ -507,13 +507,13 @@ void run ()
+ }
+ if (result.cols() == 2) {
+ Eigen::Matrix<default_type, 2, 1> azel (v.data());
+- Eigen::Vector3 dir;
++ Eigen::Vector3d dir;
+ Math::Sphere::spherical2cartesian (azel, dir);
+ dir = rotation * dir;
+ Math::Sphere::cartesian2spherical (dir, azel);
+ result.row (l) = azel;
+ } else {
+- const Eigen::Vector3 dir = rotation * Eigen::Vector3 (v.data());
++ const Eigen::Vector3d dir = rotation * Eigen::Vector3d (v.data());
+ result.row (l) = dir;
+ }
+ std::stringstream s;
+Index: mrtrix3/cmd/mtnormalise.cpp
+===================================================================
+--- mrtrix3.orig/cmd/mtnormalise.cpp
++++ mrtrix3/cmd/mtnormalise.cpp
+@@ -141,7 +141,7 @@ struct PolyBasisFunction { MEMALIGN (Pol
+
+ const int n_basis_vecs;
+
+- FORCE_INLINE Eigen::VectorXd operator () (const Eigen::Vector3& pos) const {
++ FORCE_INLINE Eigen::VectorXd operator () (const Eigen::Vector3d& pos) const {
+ double x = pos[0];
+ double y = pos[1];
+ double z = pos[2];
+@@ -231,8 +231,8 @@ Eigen::MatrixXd initialise_basis (IndexT
+ const uint32_t idx = index.value();
+ if (idx != std::numeric_limits<uint32_t>::max()) {
+ assert (idx < basis.rows());
+- Eigen::Vector3 vox (index.index(0), index.index(1), index.index(2));
+- Eigen::Vector3 pos = transform.voxel2scanner * vox;
++ Eigen::Vector3d vox (index.index(0), index.index(1), index.index(2));
++ Eigen::Vector3d pos = transform.voxel2scanner * vox;
+ basis.row(idx) = basis_function (pos);
+ }
+ }
+@@ -413,8 +413,8 @@ ImageType compute_full_field (int order,
+
+ struct FieldWriter { NOMEMALIGN
+ void operator() (ImageType& field) const {
+- Eigen::Vector3 vox (field.index(0), field.index(1), field.index(2));
+- Eigen::Vector3 pos = transform.voxel2scanner * vox;
++ Eigen::Vector3d vox (field.index(0), field.index(1), field.index(2));
++ Eigen::Vector3d pos = transform.voxel2scanner * vox;
+ field.value() = std::exp (basis_function (pos).dot (field_coeffs));
+ }
+
+Index: mrtrix3/cmd/peaks2fixel.cpp
+===================================================================
+--- mrtrix3.orig/cmd/peaks2fixel.cpp
++++ mrtrix3/cmd/peaks2fixel.cpp
+@@ -45,12 +45,12 @@ void usage ()
+
+
+
+-vector<Eigen::Vector3> get (Image<float>& data)
++vector<Eigen::Vector3d> get (Image<float>& data)
+ {
+ data.index(3) = 0;
+- vector<Eigen::Vector3> result;
++ vector<Eigen::Vector3d> result;
+ while (data.index(3) < data.size(3)) {
+- Eigen::Vector3 direction;
++ Eigen::Vector3d direction;
+ for (size_t axis = 0; axis != 3; ++axis) {
+ direction[axis] = data.value();
+ data.index(3)++;
+Index: mrtrix3/cmd/tck2fixel.cpp
+===================================================================
+--- mrtrix3.orig/cmd/tck2fixel.cpp
++++ mrtrix3/cmd/tck2fixel.cpp
+@@ -44,7 +44,7 @@ class TrackProcessor { MEMALIGN (TrackPr
+ using SetVoxelDir = DWI::Tractography::Mapping::SetVoxelDir;
+
+ TrackProcessor (Image<index_type>& fixel_indexer,
+- const vector<Eigen::Vector3>& fixel_directions,
++ const vector<Eigen::Vector3d>& fixel_directions,
+ vector<uint16_t>& fixel_TDI,
+ const float angular_threshold):
+ fixel_indexer (fixel_indexer) ,
+@@ -66,7 +66,7 @@ class TrackProcessor { MEMALIGN (TrackPr
+ index_type last_index = first_index + num_fibres;
+ index_type closest_fixel_index = 0;
+ float largest_dp = 0.0;
+- const Eigen::Vector3 dir (i->get_dir().normalized());
++ const Eigen::Vector3d dir (i->get_dir().normalized());
+ for (index_type j = first_index; j < last_index; ++j) {
+ const float dp = abs (dir.dot (fixel_directions[j]));
+ if (dp > largest_dp) {
+@@ -86,7 +86,7 @@ class TrackProcessor { MEMALIGN (TrackPr
+
+ private:
+ Image<index_type> fixel_indexer;
+- const vector<Eigen::Vector3>& fixel_directions;
++ const vector<Eigen::Vector3d>& fixel_directions;
+ vector<uint16_t>& fixel_TDI;
+ const float angular_threshold_dp;
+ };
+@@ -135,8 +135,8 @@ void run ()
+
+ const float angular_threshold = get_option_value ("angle", DEFAULT_ANGLE_THRESHOLD);
+
+- vector<Eigen::Vector3> positions (num_fixels);
+- vector<Eigen::Vector3> directions (num_fixels);
++ vector<Eigen::Vector3d> positions (num_fixels);
++ vector<Eigen::Vector3d> directions (num_fixels);
+
+ const std::string output_fixel_folder = argument[2];
+ Fixel::copy_index_and_directions_file (input_fixel_folder, output_fixel_folder);
+@@ -146,7 +146,7 @@ void run ()
+ // Load template fixel directions
+ Transform image_transform (index_image);
+ for (auto i = Loop ("loading template fixel directions and positions", index_image, 0, 3)(index_image); i; ++i) {
+- const Eigen::Vector3 vox ((default_type)index_image.index(0), (default_type)index_image.index(1), (default_type)index_image.index(2));
++ const Eigen::Vector3d vox ((default_type)index_image.index(0), (default_type)index_image.index(1), (default_type)index_image.index(2));
+ index_image.index(3) = 1;
+ index_type offset = index_image.value();
+ index_type fixel_index = 0;
+Index: mrtrix3/cmd/tckconvert.cpp
+===================================================================
+--- mrtrix3.orig/cmd/tckconvert.cpp
++++ mrtrix3/cmd/tckconvert.cpp
+@@ -331,7 +331,7 @@ class PLYWriter: public WriterInterface<
+ }
+
+ // calculate centroid
+- Eigen::Vector3 centroid(coord.row(0).mean(), coord.row(1).mean(), coord.row(2).mean());
++ Eigen::Vector3d centroid(coord.row(0).mean(), coord.row(1).mean(), coord.row(2).mean());
+
+ // subtract centroid
+ coord.row(0).array() -= centroid(0);
+Index: mrtrix3/cmd/transformcompose.cpp
+===================================================================
+--- mrtrix3.orig/cmd/transformcompose.cpp
++++ mrtrix3/cmd/transformcompose.cpp
+@@ -26,7 +26,7 @@ using namespace App;
+ class TransformBase { MEMALIGN(TransformBase)
+ public:
+ virtual ~TransformBase(){}
+- virtual Eigen::Vector3 transform_point (const Eigen::Vector3& input) = 0;
++ virtual Eigen::Vector3d transform_point (const Eigen::Vector3d& input) = 0;
+ };
+
+
+@@ -34,8 +34,8 @@ class Warp : public TransformBase { MEMA
+ public:
+ Warp (Image<default_type>& in) : interp (in) {}
+
+- Eigen::Vector3 transform_point (const Eigen::Vector3 &input) {
+- Eigen::Vector3 output;
++ Eigen::Vector3d transform_point (const Eigen::Vector3d &input) {
++ Eigen::Vector3d output;
+ if (interp.scanner (input))
+ output = interp.row(3);
+ else
+@@ -52,8 +52,8 @@ class Linear : public TransformBase { ME
+ public:
+ Linear (const transform_type& transform) : transform (transform) {}
+
+- Eigen::Vector3 transform_point (const Eigen::Vector3 &input) {
+- Eigen::Vector3 output = transform * input;
++ Eigen::Vector3d transform_point (const Eigen::Vector3d &input) {
++ Eigen::Vector3d output = transform * input;
+ return output;
+ }
+
+@@ -158,11 +158,11 @@ void run ()
+
+ Transform template_transform (output);
+ for (auto i = Loop ("composing transformations", output, 0, 3) (output); i ; ++i) {
+- Eigen::Vector3 voxel ((default_type) output.index(0),
++ Eigen::Vector3d voxel ((default_type) output.index(0),
+ (default_type) output.index(1),
+ (default_type) output.index(2));
+
+- Eigen::Vector3 position = template_transform.voxel2scanner * voxel;
++ Eigen::Vector3d position = template_transform.voxel2scanner * voxel;
+ ssize_t index = transform_list.size() - 1;
+ while (index >= 0) {
+ position = transform_list[index]->transform_point (position);
+Index: mrtrix3/cmd/warpinit.cpp
+===================================================================
+--- mrtrix3.orig/cmd/warpinit.cpp
++++ mrtrix3/cmd/warpinit.cpp
+@@ -65,7 +65,7 @@ void run ()
+ Transform transform (header);
+
+ auto func = [&transform](Image<float>& image) {
+- Eigen::Vector3 voxel_pos ((float)image.index(0), (float)image.index(1), (float)image.index(2));
++ Eigen::Vector3d voxel_pos ((float)image.index(0), (float)image.index(1), (float)image.index(2));
+ Eigen::Vector3f scanner_pos = (transform.voxel2scanner * voxel_pos).cast<float>();
+ for (auto l = Loop (3) (image); l; ++l)
+ image.value() = scanner_pos[image.index(3)];
+Index: mrtrix3/core/adapter/reslice.h
+===================================================================
+--- mrtrix3.orig/core/adapter/reslice.h
++++ mrtrix3/core/adapter/reslice.h
+@@ -149,12 +149,12 @@ namespace MR
+ OS[0] = OS[1] = OS[2] = 1;
+ }
+ else {
+- Vector3 y = direct_transform * Vector3 (0.0, 0.0, 0.0);
+- Vector3 x0 = direct_transform * Vector3 (1.0, 0.0, 0.0);
++ Vector3d y = direct_transform * Vector3d (0.0, 0.0, 0.0);
++ Vector3d x0 = direct_transform * Vector3d (1.0, 0.0, 0.0);
+ OS[0] = std::ceil ((1.0-std::numeric_limits<default_type>::epsilon()) * (y-x0).norm());
+- x0 = direct_transform * Vector3 (0.0, 1.0, 0.0);
++ x0 = direct_transform * Vector3d (0.0, 1.0, 0.0);
+ OS[1] = std::ceil ((1.0-std::numeric_limits<default_type>::epsilon()) * (y-x0).norm());
+- x0 = direct_transform * Vector3 (0.0, 0.0, 1.0);
++ x0 = direct_transform * Vector3d (0.0, 0.0, 1.0);
+ OS[2] = std::ceil ((1.0-std::numeric_limits<default_type>::epsilon()) * (y-x0).norm());
+ }
+ }
+@@ -194,9 +194,9 @@ namespace MR
+ value_type value () {
+ using namespace Eigen;
+ if (oversampling) {
+- Vector3 d (x[0]+from[0], x[1]+from[1], x[2]+from[2]);
++ Vector3d d (x[0]+from[0], x[1]+from[1], x[2]+from[2]);
+ default_type sum (0.0);
+- Vector3 s;
++ Vector3d s;
+ for (uint32_t z = 0; z < OS[2]; ++z) {
+ s[2] = d[2] + z*inc[2];
+ for (uint32_t y = 0; y < OS[1]; ++y) {
+@@ -210,7 +210,7 @@ namespace MR
+ }
+ return normalise<value_type> (sum, norm);
+ }
+- interp.voxel (direct_transform * Vector3 (x[0], x[1], x[2]));
++ interp.voxel (direct_transform * Vector3d (x[0], x[1], x[2]));
+ return interp.value();
+ }
+
+Index: mrtrix3/core/adapter/warp.h
+===================================================================
+--- mrtrix3.orig/core/adapter/warp.h
++++ mrtrix3/core/adapter/warp.h
+@@ -90,7 +90,7 @@ namespace MR
+
+
+ value_type value () {
+- Eigen::Vector3 pos = get_position();
++ Eigen::Vector3d pos = get_position();
+ if (std::isnan(pos[0]) || std::isnan(pos[1]) || std::isnan(pos[2]))
+ return value_when_out_of_bounds;
+ interp.scanner (pos);
+Index: mrtrix3/core/axes.cpp
+===================================================================
+--- mrtrix3.orig/core/axes.cpp
++++ mrtrix3/core/axes.cpp
+@@ -27,7 +27,7 @@ namespace MR
+
+
+
+- std::string dir2id (const Eigen::Vector3& axis)
++ std::string dir2id (const Eigen::Vector3d& axis)
+ {
+ if (axis[0] == -1) {
+ assert (!axis[1]); assert (!axis[2]); return "i-";
+@@ -48,7 +48,7 @@ namespace MR
+
+
+
+- Eigen::Vector3 id2dir (const std::string& id)
++ Eigen::Vector3d id2dir (const std::string& id)
+ {
+ if (id == "i-")
+ return { -1, 0, 0 };
+Index: mrtrix3/core/axes.h
+===================================================================
+--- mrtrix3.orig/core/axes.h
++++ mrtrix3/core/axes.h
+@@ -36,8 +36,8 @@ namespace MR
+ * phase-encoding direction between a 3-vector (e.g.
+ * [0 1 0] ) and a NIfTI axis identifier (e.g. 'i-')
+ */
+- std::string dir2id (const Eigen::Vector3&);
+- Eigen::Vector3 id2dir (const std::string&);
++ std::string dir2id (const Eigen::Vector3d&);
++ Eigen::Vector3d id2dir (const std::string&);
+
+
+
+Index: mrtrix3/core/file/dicom/image.h
+===================================================================
+--- mrtrix3.orig/core/file/dicom/image.h
++++ mrtrix3/core/file/dicom/image.h
+@@ -54,7 +54,7 @@ namespace MR {
+ }
+
+ size_t acq_dim[2], dim[2], series_num, instance, acq, sequence, echo_index;
+- Eigen::Vector3 position_vector, orientation_x, orientation_y, orientation_z, G;
++ Eigen::Vector3d position_vector, orientation_x, orientation_y, orientation_z, G;
+ default_type distance, pixel_size[2], slice_thickness, slice_spacing, scale_slope, scale_intercept, bvalue;
+ size_t data, bits_alloc, data_size, frame_offset;
+ std::string filename, image_type;
+@@ -99,7 +99,7 @@ namespace MR {
+ else {
+ if (!orientation_x.allFinite() || !orientation_y.allFinite())
+ throw Exception ("slice orientation information missing from DICOM header!");
+- Eigen::Vector3 normal = orientation_x.cross (orientation_y);
++ Eigen::Vector3d normal = orientation_x.cross (orientation_y);
+ if (normal.dot (orientation_z) < 0.0)
+ orientation_z = -normal;
+ else
+Index: mrtrix3/core/file/json_utils.cpp
+===================================================================
+--- mrtrix3.orig/core/file/json_utils.cpp
++++ mrtrix3/core/file/json_utils.cpp
+@@ -160,8 +160,8 @@ namespace MR
+ auto slice_encoding_it = header.keyval().find ("SliceEncodingDirection");
+ if (slice_encoding_it != header.keyval().end()) {
+ if (do_realign) {
+- const Eigen::Vector3 orig_dir (Axes::id2dir (slice_encoding_it->second));
+- Eigen::Vector3 new_dir;
++ const Eigen::Vector3d orig_dir (Axes::id2dir (slice_encoding_it->second));
++ Eigen::Vector3d new_dir;
+ for (size_t axis = 0; axis != 3; ++axis)
+ new_dir[axis] = flip[perm[axis]] ? -orig_dir[perm[axis]] : orig_dir[perm[axis]];
+ slice_encoding_it->second = Axes::dir2id (new_dir);
+@@ -293,8 +293,8 @@ namespace MR
+ }
+ auto slice_encoding_it = H_adj.keyval().find ("SliceEncodingDirection");
+ if (slice_encoding_it != H_adj.keyval().end()) {
+- const Eigen::Vector3 orig_dir (Axes::id2dir (slice_encoding_it->second));
+- Eigen::Vector3 new_dir;
++ const Eigen::Vector3d orig_dir (Axes::id2dir (slice_encoding_it->second));
++ Eigen::Vector3d new_dir;
+ for (size_t axis = 0; axis != 3; ++axis)
+ new_dir[axis] = flip[axis] ? orig_dir[order[axis]] : -orig_dir[order[axis]];
+ slice_encoding_it->second = Axes::dir2id (new_dir);
+Index: mrtrix3/core/filter/gradient.h
+===================================================================
+--- mrtrix3.orig/core/filter/gradient.h
++++ mrtrix3/core/filter/gradient.h
+@@ -144,7 +144,7 @@ namespace MR
+ if (wrt_scanner) {
+ Transform transform (in);
+ for (auto l = Loop(0,3) (out); l; ++l)
+- out.row(3) = transform.image2scanner.linear() * Eigen::Vector3 (out.row(3));
++ out.row(3) = transform.image2scanner.linear() * Eigen::Vector3d (out.row(3));
+ }
+ }
+ }
+Index: mrtrix3/core/filter/resize.h
+===================================================================
+--- mrtrix3.orig/core/filter/resize.h
++++ mrtrix3/core/filter/resize.h
+@@ -78,7 +78,7 @@ namespace MR
+ if (voxel_size.size() != 3)
+ throw Exception ("the voxel size must be defined using a value for all three dimensions.");
+
+- Eigen::Vector3 original_extent;
++ Eigen::Vector3d original_extent;
+ for (size_t j = 0; j < 3; ++j) {
+ if (voxel_size[j] <= 0.0)
+ throw Exception ("the voxel size must be larger than zero");
+Index: mrtrix3/core/header.cpp
+===================================================================
+--- mrtrix3.orig/core/header.cpp
++++ mrtrix3/core/header.cpp
+@@ -727,8 +727,8 @@ namespace MR
+ // header, that's also necessary to update here
+ auto slice_encoding_it = keyval().find ("SliceEncodingDirection");
+ if (slice_encoding_it != keyval().end()) {
+- const Eigen::Vector3 orig_dir (Axes::id2dir (slice_encoding_it->second));
+- Eigen::Vector3 new_dir;
++ const Eigen::Vector3d orig_dir (Axes::id2dir (slice_encoding_it->second));
++ Eigen::Vector3d new_dir;
+ for (size_t axis = 0; axis != 3; ++axis)
+ new_dir[axis] = orig_dir[realign_perm_[axis]] * (realign_flip_[realign_perm_[axis]] ? -1.0 : 1.0);
+ slice_encoding_it->second = Axes::dir2id (new_dir);
+Index: mrtrix3/core/image_helpers.h
+===================================================================
+--- mrtrix3.orig/core/image_helpers.h
++++ mrtrix3/core/image_helpers.h
+@@ -452,8 +452,8 @@ namespace MR
+ if (!dimensions_match(in1, in2, 0, 3))
+ return false;
+
+- const Eigen::Vector3 vs1 (in1.spacing(0), in1.spacing(1), in1.spacing(2));
+- const Eigen::Vector3 vs2 (in2.spacing(0), in2.spacing(1), in2.spacing(2));
++ const Eigen::Vector3d vs1 (in1.spacing(0), in1.spacing(1), in1.spacing(2));
++ const Eigen::Vector3d vs2 (in2.spacing(0), in2.spacing(1), in2.spacing(2));
+
+ Eigen::MatrixXd voxel_coord = Eigen::MatrixXd::Zero(4,4);
+ voxel_coord.row(3).fill(1.0);
+Index: mrtrix3/core/interp/base.h
+===================================================================
+--- mrtrix3.orig/core/interp/base.h
++++ mrtrix3/core/interp/base.h
+@@ -208,10 +208,10 @@ namespace MR
+ }
+
+ template <class VectorType>
+- Eigen::Vector3 intravoxel_offset (const VectorType& pos) {
++ Eigen::Vector3d intravoxel_offset (const VectorType& pos) {
+ if (set_out_of_bounds (pos))
+- return Eigen::Vector3 (NaN, NaN, NaN);
+- return Eigen::Vector3 (pos[0]-std::floor (pos[0]), pos[1]-std::floor (pos[1]), pos[2]-std::floor (pos[2]));
++ return Eigen::Vector3d (NaN, NaN, NaN);
++ return Eigen::Vector3d (pos[0]-std::floor (pos[0]), pos[1]-std::floor (pos[1]), pos[2]-std::floor (pos[2]));
+ }
+
+ };
+Index: mrtrix3/core/interp/cubic.h
+===================================================================
+--- mrtrix3.orig/core/interp/cubic.h
++++ mrtrix3/core/interp/cubic.h
+@@ -80,7 +80,7 @@ namespace MR
+
+ protected:
+ SplineType H[3];
+- Eigen::Vector3 P;
++ Eigen::Vector3d P;
+
+ ssize_t clamp (ssize_t x, ssize_t dim) const {
+ if (x < 0) return 0;
+@@ -120,7 +120,7 @@ namespace MR
+ /*! See file interp/base.h for details. */
+ template <class VectorType>
+ bool voxel (const VectorType& pos) {
+- Eigen::Vector3 f = Base<ImageType>::intravoxel_offset (pos);
++ Eigen::Vector3d f = Base<ImageType>::intravoxel_offset (pos);
+ if (Base<ImageType>::out_of_bounds)
+ return false;
+ P = pos;
+@@ -244,7 +244,7 @@ namespace MR
+ /*! See file interp/base.h for details. */
+ template <class VectorType>
+ bool voxel (const VectorType& pos) {
+- Eigen::Vector3 f = Base<ImageType>::intravoxel_offset (pos);
++ Eigen::Vector3d f = Base<ImageType>::intravoxel_offset (pos);
+ if (Base<ImageType>::out_of_bounds)
+ return false;
+ P = pos;
+@@ -405,7 +405,7 @@ namespace MR
+ /*! See file interp/base.h for details. */
+ template <class VectorType>
+ bool voxel (const VectorType& pos) {
+- Eigen::Vector3 f = Base<ImageType>::intravoxel_offset (pos);
++ Eigen::Vector3d f = Base<ImageType>::intravoxel_offset (pos);
+ if (Base<ImageType>::out_of_bounds)
+ return false;
+ P = pos;
+Index: mrtrix3/core/interp/linear.h
+===================================================================
+--- mrtrix3.orig/core/interp/linear.h
++++ mrtrix3/core/interp/linear.h
+@@ -107,7 +107,7 @@ namespace MR
+
+ protected:
+ const coef_type zero, eps;
+- Eigen::Vector3 P;
++ Eigen::Vector3d P;
+
+ ssize_t clamp (ssize_t x, ssize_t dim) const {
+ if (x < 0) return 0;
+@@ -148,7 +148,7 @@ namespace MR
+ /*! See file interp/base.h for details. */
+ template <class VectorType>
+ bool voxel (const VectorType& pos) {
+- Eigen::Vector3 f = Base<ImageType>::intravoxel_offset (pos);
++ Eigen::Vector3d f = Base<ImageType>::intravoxel_offset (pos);
+ if (Base<ImageType>::out_of_bounds)
+ return false;
+ P = pos;
+@@ -276,7 +276,7 @@ namespace MR
+ /*! See file interp/base.h for details. */
+ template <class VectorType>
+ bool voxel (const VectorType& pos) {
+- Eigen::Vector3 f = Base<ImageType>::intravoxel_offset (pos.template cast<default_type>());
++ Eigen::Vector3d f = Base<ImageType>::intravoxel_offset (pos.template cast<default_type>());
+ if (Base<ImageType>::out_of_bounds)
+ return false;
+ P = pos;
+@@ -436,7 +436,7 @@ namespace MR
+ /*! See file interp/base.h for details. */
+ template <class VectorType>
+ bool voxel (const VectorType& pos) {
+- Eigen::Vector3 f = Base<ImageType>::intravoxel_offset (pos.template cast<default_type>());
++ Eigen::Vector3d f = Base<ImageType>::intravoxel_offset (pos.template cast<default_type>());
+ if (Base<ImageType>::out_of_bounds)
+ return false;
+ P = pos;
+Index: mrtrix3/core/phase_encoding.h
+===================================================================
+--- mrtrix3.orig/core/phase_encoding.h
++++ mrtrix3/core/phase_encoding.h
+@@ -116,7 +116,7 @@ namespace MR
+ erase ("TotalReadoutTime");
+ } else {
+ erase ("pe_scheme");
+- const Eigen::Vector3 dir { PE(0, 0), PE(0, 1), PE(0, 2) };
++ const Eigen::Vector3d dir { PE(0, 0), PE(0, 1), PE(0, 2) };
+ header.keyval()["PhaseEncodingDirection"] = Axes::dir2id (dir);
+ if (PE.cols() >= 4)
+ header.keyval()["TotalReadoutTime"] = str(PE(0, 3), 3);
+Index: mrtrix3/core/types.h
+===================================================================
+--- mrtrix3.orig/core/types.h
++++ mrtrix3/core/types.h
+@@ -322,11 +322,6 @@ namespace std
+
+ }
+
+-namespace Eigen {
+- using Vector3 = Matrix<MR::default_type, 3, 1>;
+- using Vector4 = Matrix<MR::default_type, 4, 1>;
+-}
+-
+ #endif
+
+
+Index: mrtrix3/src/dwi/directions/set.cpp
+===================================================================
+--- mrtrix3.orig/src/dwi/directions/set.cpp
++++ mrtrix3/src/dwi/directions/set.cpp
+@@ -103,7 +103,7 @@ namespace MR {
+ Vertex (const Set& set, const index_type index, const bool inverse) :
+ dir (set[index] * (inverse ? -1.0 : 1.0)),
+ index (index) { }
+- const Eigen::Vector3 dir;
++ const Eigen::Vector3d dir;
+ const index_type index; // Indexes the underlying direction set
+ };
+
+@@ -115,7 +115,7 @@ namespace MR {
+ dist (std::max ( { vertices[one].dir.dot (normal), vertices[two].dir.dot (normal), vertices[three].dir.dot (normal) } ) ) { }
+ bool includes (const index_type i) const { return (indices[0] == i || indices[1] == i || indices[2] == i); }
+ const std::array<index_type,3> indices; // Indexes the vertices vector
+- const Eigen::Vector3 normal;
++ const Eigen::Vector3d normal;
+ const default_type dist;
+ };
+
+@@ -327,7 +327,7 @@ namespace MR {
+
+
+
+- index_type FastLookupSet::select_direction (const Eigen::Vector3& p) const
++ index_type FastLookupSet::select_direction (const Eigen::Vector3d& p) const
+ {
+
+ const size_t grid_index = dir2gridindex (p);
+@@ -349,7 +349,7 @@ namespace MR {
+
+
+
+- index_type FastLookupSet::select_direction_slow (const Eigen::Vector3& p) const
++ index_type FastLookupSet::select_direction_slow (const Eigen::Vector3d& p) const
+ {
+
+ index_type dir = 0;
+@@ -417,7 +417,7 @@ namespace MR {
+ case 3: el += el_grid_step; break;
+ }
+
+- const Eigen::Vector3 p (cos(az) * sin(el), sin(az) * sin(el), cos (el));
++ const Eigen::Vector3d p (cos(az) * sin(el), sin(az) * sin(el), cos (el));
+ const index_type nearest_dir = select_direction_slow (p);
+ bool dir_present = false;
+ for (vector<index_type>::const_iterator d = grid_lookup[i].begin(); !dir_present && d != grid_lookup[i].end(); ++d)
+@@ -453,7 +453,7 @@ namespace MR {
+
+
+
+- size_t FastLookupSet::dir2gridindex (const Eigen::Vector3& p) const
++ size_t FastLookupSet::dir2gridindex (const Eigen::Vector3d& p) const
+ {
+
+ const default_type azimuth = atan2(p[1], p[0]);
+@@ -477,7 +477,7 @@ namespace MR {
+ size_t error_count = 0;
+ const size_t checks = 1000000;
+ for (size_t i = 0; i != checks; ++i) {
+- Eigen::Vector3 p (normal(rng), normal(rng), normal(rng));
++ Eigen::Vector3d p (normal(rng), normal(rng), normal(rng));
+ p.normalize();
+ if (select_direction (p) != select_direction_slow (p))
+ ++error_count;
+Index: mrtrix3/src/dwi/directions/set.h
+===================================================================
+--- mrtrix3.orig/src/dwi/directions/set.h
++++ mrtrix3/src/dwi/directions/set.h
+@@ -84,7 +84,7 @@ namespace MR {
+ }
+
+ size_t size () const { return unit_vectors.size(); }
+- const Eigen::Vector3& get_dir (const size_t i) const { assert (i < size()); return unit_vectors[i]; }
++ const Eigen::Vector3d& get_dir (const size_t i) const { assert (i < size()); return unit_vectors[i]; }
+ const vector<index_type>& get_adj_dirs (const size_t i) const { assert (i < size()); return adj_dirs[i]; }
+ bool dirs_are_adjacent (const index_type one, const index_type two) const {
+ assert (one < size());
+@@ -98,13 +98,13 @@ namespace MR {
+
+ index_type get_min_linkage (const index_type one, const index_type two) const;
+
+- const vector<Eigen::Vector3>& get_dirs() const { return unit_vectors; }
+- const Eigen::Vector3& operator[] (const size_t i) const { assert (i < size()); return unit_vectors[i]; }
++ const vector<Eigen::Vector3d>& get_dirs() const { return unit_vectors; }
++ const Eigen::Vector3d& operator[] (const size_t i) const { assert (i < size()); return unit_vectors[i]; }
+
+
+ protected:
+
+- vector<Eigen::Vector3> unit_vectors;
++ vector<Eigen::Vector3d> unit_vectors;
+ vector< vector<index_type> > adj_dirs; // Note: not self-inclusive
+
+
+@@ -179,7 +179,7 @@ namespace MR {
+ az_begin (that.az_begin),
+ el_begin (that.el_begin) { }
+
+- index_type select_direction (const Eigen::Vector3&) const;
++ index_type select_direction (const Eigen::Vector3d&) const;
+
+
+
+@@ -192,11 +192,11 @@ namespace MR {
+
+ FastLookupSet ();
+
+- index_type select_direction_slow (const Eigen::Vector3&) const;
++ index_type select_direction_slow (const Eigen::Vector3d&) const;
+
+ void initialise();
+
+- size_t dir2gridindex (const Eigen::Vector3&) const;
++ size_t dir2gridindex (const Eigen::Vector3d&) const;
+
+ void test_lookup() const;
+
+Index: mrtrix3/src/dwi/fmls.cpp
+===================================================================
+--- mrtrix3.orig/src/dwi/fmls.cpp
++++ mrtrix3/src/dwi/fmls.cpp
+@@ -264,7 +264,7 @@ namespace MR {
+
+ // Revise multiple peaks if present
+ for (size_t peak_index = 0; peak_index != i->num_peaks(); ++peak_index) {
+- Eigen::Vector3 newton_peak_dir = i->get_peak_dir (peak_index); // to be updated by subsequent Math::SH::get_peak() call
++ Eigen::Vector3d newton_peak_dir = i->get_peak_dir (peak_index); // to be updated by subsequent Math::SH::get_peak() call
+ const default_type newton_peak_value = Math::SH::get_peak (in, lmax, newton_peak_dir, &(*precomputer));
+ if (std::isfinite (newton_peak_value) && newton_peak_dir.allFinite()) {
+
+Index: mrtrix3/src/dwi/fmls.h
+===================================================================
+--- mrtrix3.orig/src/dwi/fmls.h
++++ mrtrix3/src/dwi/fmls.h
+@@ -92,13 +92,13 @@ namespace MR
+ assert ((value <= 0.0 && neg) || (value >= 0.0 && !neg));
+ mask[bin] = true;
+ values[bin] = value;
+- const Eigen::Vector3& dir = mask.get_dirs()[bin];
++ const Eigen::Vector3d& dir = mask.get_dirs()[bin];
+ const default_type multiplier = (mean_dir.dot (dir)) > 0.0 ? 1.0 : -1.0;
+ mean_dir += dir * multiplier * abs(value) * weight;
+ integral += abs (value * weight);
+ }
+
+- void revise_peak (const size_t index, const Eigen::Vector3& revised_peak_dir, const default_type revised_peak_value)
++ void revise_peak (const size_t index, const Eigen::Vector3d& revised_peak_dir, const default_type revised_peak_value)
+ {
+ assert (!neg);
+ assert (index < num_peaks());
+@@ -142,8 +142,8 @@ namespace MR
+ const Eigen::Array<default_type, Eigen::Dynamic, 1>& get_values() const { return values; }
+ default_type get_max_peak_value() const { return max_peak_value; }
+ size_t num_peaks() const { return peak_dirs.size(); }
+- const Eigen::Vector3& get_peak_dir (const size_t i) const { assert (i < num_peaks()); return peak_dirs[i]; }
+- const Eigen::Vector3& get_mean_dir() const { return mean_dir; }
++ const Eigen::Vector3d& get_peak_dir (const size_t i) const { assert (i < num_peaks()); return peak_dirs[i]; }
++ const Eigen::Vector3d& get_mean_dir() const { return mean_dir; }
+ default_type get_integral() const { return integral; }
+ bool is_negative() const { return neg; }
+
+@@ -152,8 +152,8 @@ namespace MR
+ DWI::Directions::Mask mask;
+ Eigen::Array<default_type, Eigen::Dynamic, 1> values;
+ default_type max_peak_value;
+- vector<Eigen::Vector3> peak_dirs;
+- Eigen::Vector3 mean_dir;
++ vector<Eigen::Vector3d> peak_dirs;
++ Eigen::Vector3d mean_dir;
+ default_type integral;
+ bool neg;
+
+Index: mrtrix3/src/dwi/tractography/GT/particlegrid.h
+===================================================================
+--- mrtrix3.orig/src/dwi/tractography/GT/particlegrid.h
++++ mrtrix3/src/dwi/tractography/GT/particlegrid.h
+@@ -53,7 +53,7 @@ namespace MR {
+
+ // Initialise scanner-to-grid transform
+ Eigen::DiagonalMatrix<default_type, 3> newspacing (2.0*Particle::L, 2.0*Particle::L, 2.0*Particle::L);
+- Eigen::Vector3 shift (image.spacing(0)/2.0 - Particle::L,
++ Eigen::Vector3d shift (image.spacing(0)/2.0 - Particle::L,
+ image.spacing(1)/2.0 - Particle::L,
+ image.spacing(2)/2.0 - Particle::L);
+ T_s2g = image.transform() * newspacing;
+Index: mrtrix3/src/dwi/tractography/SIFT/model_base.h
+===================================================================
+--- mrtrix3.orig/src/dwi/tractography/SIFT/model_base.h
++++ mrtrix3/src/dwi/tractography/SIFT/model_base.h
+@@ -74,7 +74,7 @@ namespace MR
+ weight (1.0),
+ dir () { }
+
+- FixelBase (const default_type amp, const Eigen::Vector3& d) :
++ FixelBase (const default_type amp, const Eigen::Vector3d& d) :
+ FOD (amp),
+ TD (0.0),
+ weight (1.0),
+@@ -91,7 +91,7 @@ namespace MR
+ default_type get_FOD() const { return FOD; }
+ default_type get_TD() const { return TD; }
+ default_type get_weight() const { return weight; }
+- const Eigen::Vector3& get_dir() const { return dir; }
++ const Eigen::Vector3d& get_dir() const { return dir; }
+
+ void scale_FOD (const default_type factor) { FOD *= factor; }
+ void set_weight (const default_type w) { weight = w; }
+@@ -105,7 +105,7 @@ namespace MR
+
+ protected:
+ default_type FOD, TD, weight;
+- Eigen::Vector3 dir;
++ Eigen::Vector3d dir;
+
+ default_type get_cost_unweighted (const default_type mu) const { return Math::pow2 (get_diff (mu)); }
+
+Index: mrtrix3/src/dwi/tractography/connectome/tck2nodes.cpp
+===================================================================
+--- mrtrix3.orig/src/dwi/tractography/connectome/tck2nodes.cpp
++++ mrtrix3/src/dwi/tractography/connectome/tck2nodes.cpp
+@@ -31,8 +31,8 @@ namespace Connectome {
+
+ node_t Tck2nodes_end_voxels::select_node (const Tractography::Streamline<>& tck, Image<node_t>& v, const bool end) const
+ {
+- const Eigen::Vector3 p ((end ? tck.back() : tck.front()).cast<default_type>());
+- const Eigen::Vector3 v_float (transform->scanner2voxel * p);
++ const Eigen::Vector3d p ((end ? tck.back() : tck.front()).cast<default_type>());
++ const Eigen::Vector3d v_float (transform->scanner2voxel * p);
+ for (size_t axis = 0; axis != 3; ++axis)
+ v.index(axis) = std::round (v_float[axis]);
+ if (is_out_of_bounds (v))
+@@ -72,14 +72,14 @@ node_t Tck2nodes_radial::select_node (co
+ default_type min_dist = max_dist;
+ node_t node = 0;
+
+- const Eigen::Vector3 p = (end ? tck.back() : tck.front()).cast<default_type>();
+- const Eigen::Vector3 v_float = transform->scanner2voxel * p;
++ const Eigen::Vector3d p = (end ? tck.back() : tck.front()).cast<default_type>();
++ const Eigen::Vector3d v_float = transform->scanner2voxel * p;
+ const voxel_type centre { int(std::round (v_float[0])), int(std::round (v_float[1])), int(std::round (v_float[2])) };
+
+ for (vector<voxel_type>::const_iterator offset = radial_search.begin(); offset != radial_search.end(); ++offset) {
+
+ const voxel_type this_voxel (centre + *offset);
+- const Eigen::Vector3 p_voxel (transform->voxel2scanner * this_voxel.matrix().cast<default_type>());
++ const Eigen::Vector3d p_voxel (transform->voxel2scanner * this_voxel.matrix().cast<default_type>());
+ const default_type dist ((p - p_voxel).norm());
+
+ if (dist > min_dist + 2*max_add_dist)
+@@ -112,7 +112,7 @@ node_t Tck2nodes_revsearch::select_node
+ default_type dist = 0.0;
+
+ for (int index = start_index; index != midpoint_index; index += step) {
+- const Eigen::Vector3 v_float = transform->scanner2voxel * tck[index].cast<default_type>();
++ const Eigen::Vector3d v_float = transform->scanner2voxel * tck[index].cast<default_type>();
+ const voxel_type voxel { int(std::round (v_float[0])), int(std::round (v_float[1])), int(std::round (v_float[2])) };
+ assign_pos_of (voxel).to (v);
+ if (!is_out_of_bounds (v)) {
+@@ -137,15 +137,15 @@ node_t Tck2nodes_forwardsearch::select_n
+ // Start by defining the endpoint and the tangent at the endpoint
+ const int index = end ? (tck.size() - 1) : 0;
+ const int step = end ? -1 : 1;
+- const Eigen::Vector3 p = tck[index].cast<default_type>();
+- Eigen::Vector3 t;
++ const Eigen::Vector3d p = tck[index].cast<default_type>();
++ Eigen::Vector3d t;
+
+ // Heuristic for determining the tangent at the streamline endpoint
+ if (tck.size() > 2) {
+- const Eigen::Vector3 second_last_step = (tck[index+step] - tck[index+2*step]).cast<default_type>();
+- const Eigen::Vector3 last_step = (tck[index] - tck[index+step]).cast<default_type>();
++ const Eigen::Vector3d second_last_step = (tck[index+step] - tck[index+2*step]).cast<default_type>();
++ const Eigen::Vector3d last_step = (tck[index] - tck[index+step]).cast<default_type>();
+ const default_type length_ratio = last_step.norm() / second_last_step.norm();
+- const Eigen::Vector3 curvature (last_step - (second_last_step * length_ratio));
++ const Eigen::Vector3d curvature (last_step - (second_last_step * length_ratio));
+ t = last_step + curvature;
+ } else {
+ t = (tck[index] - tck[index+step]).cast<default_type>();
+@@ -159,7 +159,7 @@ node_t Tck2nodes_forwardsearch::select_n
+
+ // Voxel containing streamline endpoint not guaranteed to be appropriate
+ // Should it be tested anyway? Probably
+- const Eigen::Vector3 vp (transform->scanner2voxel * p);
++ const Eigen::Vector3d vp (transform->scanner2voxel * p);
+ const voxel_type voxel { int(std::round (vp[0])), int(std::round (vp[1])), int(std::round (vp[2])) };
+ if (is_out_of_bounds (v, voxel))
+ return 0;
+@@ -203,11 +203,11 @@ node_t Tck2nodes_forwardsearch::select_n
+
+
+
+-default_type Tck2nodes_forwardsearch::get_cf (const Eigen::Vector3& p, const Eigen::Vector3& t, const voxel_type& v) const
++default_type Tck2nodes_forwardsearch::get_cf (const Eigen::Vector3d& p, const Eigen::Vector3d& t, const voxel_type& v) const
+ {
+- const Eigen::Vector3 vfloat { default_type(v[0]), default_type(v[1]), default_type(v[2]) };
+- const Eigen::Vector3 vp (transform->voxel2scanner * vfloat);
+- Eigen::Vector3 offset (vp - p);
++ const Eigen::Vector3d vfloat { default_type(v[0]), default_type(v[1]), default_type(v[2]) };
++ const Eigen::Vector3d vp (transform->voxel2scanner * vfloat);
++ Eigen::Vector3d offset (vp - p);
+ const default_type dist = offset.norm();
+ offset.normalize();
+ const default_type angle = std::acos (t.dot (offset));
+@@ -234,7 +234,7 @@ void Tck2nodes_all_voxels::select_nodes
+ {
+ std::set<node_t> result;
+ for (Streamline<>::const_iterator p = tck.begin(); p != tck.end(); ++p) {
+- const Eigen::Vector3 v_float = transform->scanner2voxel * p->cast<default_type>();
++ const Eigen::Vector3d v_float = transform->scanner2voxel * p->cast<default_type>();
+ const voxel_type voxel (int(std::round (v_float[0])), int(std::round (v_float[1])), int(std::round (v_float[2])));
+ assign_pos_of (voxel).to (v);
+ if (!is_out_of_bounds (v)) {
+Index: mrtrix3/src/dwi/tractography/connectome/tck2nodes.h
+===================================================================
+--- mrtrix3.orig/src/dwi/tractography/connectome/tck2nodes.h
++++ mrtrix3/src/dwi/tractography/connectome/tck2nodes.h
+@@ -204,7 +204,7 @@ class Tck2nodes_forwardsearch : public T
+ const default_type max_dist;
+ const default_type angle_limit;
+
+- default_type get_cf (const Eigen::Vector3&, const Eigen::Vector3&, const voxel_type&) const;
++ default_type get_cf (const Eigen::Vector3d&, const Eigen::Vector3d&, const voxel_type&) const;
+
+ };
+
+Index: mrtrix3/src/dwi/tractography/mapping/gaussian/mapper.h
+===================================================================
+--- mrtrix3.orig/src/dwi/tractography/mapping/gaussian/mapper.h
++++ mrtrix3/src/dwi/tractography/mapping/gaussian/mapper.h
+@@ -103,10 +103,10 @@ namespace MR {
+ template <class Cont> void voxelise_precise (const Streamline<>&, Cont&) const;
+ template <class Cont> void voxelise_ends (const Streamline<>&, Cont&) const;
+
+- inline void add_to_set (SetVoxel& , const Eigen::Vector3i&, const Eigen::Vector3&, const default_type, const default_type) const;
+- inline void add_to_set (SetVoxelDEC&, const Eigen::Vector3i&, const Eigen::Vector3&, const default_type, const default_type) const;
+- inline void add_to_set (SetDixel& , const Eigen::Vector3i&, const Eigen::Vector3&, const default_type, const default_type) const;
+- inline void add_to_set (SetVoxelTOD&, const Eigen::Vector3i&, const Eigen::Vector3&, const default_type, const default_type) const;
++ inline void add_to_set (SetVoxel& , const Eigen::Vector3i&, const Eigen::Vector3d&, const default_type, const default_type) const;
++ inline void add_to_set (SetVoxelDEC&, const Eigen::Vector3i&, const Eigen::Vector3d&, const default_type, const default_type) const;
++ inline void add_to_set (SetDixel& , const Eigen::Vector3i&, const Eigen::Vector3d&, const default_type, const default_type) const;
++ inline void add_to_set (SetVoxelTOD&, const Eigen::Vector3i&, const Eigen::Vector3d&, const default_type, const default_type) const;
+
+ // Convenience function to convert from streamline position index to a linear-interpolated
+ // factor value (TrackMapperTWI member field factors[] only contains one entry per pre-upsampled point)
+@@ -128,7 +128,7 @@ namespace MR {
+ for (size_t i = 0; i != last; ++i) {
+ vox = round (scanner2voxel * tck[i]);
+ if (check (vox, info)) {
+- const Eigen::Vector3 dir ((tck[i+1] - tck[prev]).cast<default_type>().normalized());
++ const Eigen::Vector3d dir ((tck[i+1] - tck[prev]).cast<default_type>().normalized());
+ const default_type factor = tck_index_to_factor (i);
+ add_to_set (output, vox, dir, 1.0, factor);
+ }
+@@ -137,7 +137,7 @@ namespace MR {
+
+ vox = round (scanner2voxel * tck[last]);
+ if (check (vox, info)) {
+- const Eigen::Vector3 dir ((tck[last] - tck[prev]).cast<default_type>().normalized());
++ const Eigen::Vector3d dir ((tck[last] - tck[prev]).cast<default_type>().normalized());
+ const default_type factor = tck_index_to_factor (last);
+ add_to_set (output, vox, dir, 1.0f, factor);
+ }
+@@ -227,7 +227,7 @@ namespace MR {
+ }
+
+ length += (p_prev - p_voxel_exit).norm();
+- Eigen::Vector3 traversal_vector = (p_voxel_exit - p_voxel_entry).cast<default_type>().normalized();
++ Eigen::Vector3d traversal_vector = (p_voxel_exit - p_voxel_entry).cast<default_type>().normalized();
+ if (traversal_vector.allFinite() && check (this_voxel, info)) {
+ const default_type index_voxel_exit = default_type(p) + mu;
+ const size_t mean_tck_index = std::round (0.5 * (index_voxel_entry + index_voxel_exit));
+@@ -247,7 +247,7 @@ namespace MR {
+ for (size_t end = 0; end != 2; ++end) {
+ const Eigen::Vector3i vox = round (scanner2voxel * (end ? tck.back() : tck.front()));
+ if (check (vox, info)) {
+- const Eigen::Vector3 dir = (end ? (tck[tck.size()-1] - tck[tck.size()-2]) : (tck[0] - tck[1])).cast<default_type>().normalized();
++ const Eigen::Vector3d dir = (end ? (tck[tck.size()-1] - tck[tck.size()-2]) : (tck[0] - tck[1])).cast<default_type>().normalized();
+ const default_type factor = (end ? factors.back() : factors.front());
+ add_to_set (out, vox, dir, 1.0, factor);
+ }
+@@ -256,21 +256,21 @@ namespace MR {
+
+
+
+- inline void TrackMapper::add_to_set (SetVoxel& out, const Eigen::Vector3i& v, const Eigen::Vector3& d, const default_type l, const default_type f) const
++ inline void TrackMapper::add_to_set (SetVoxel& out, const Eigen::Vector3i& v, const Eigen::Vector3d& d, const default_type l, const default_type f) const
+ {
+ out.insert (v, l, f);
+ }
+- inline void TrackMapper::add_to_set (SetVoxelDEC& out, const Eigen::Vector3i& v, const Eigen::Vector3& d, const default_type l, const default_type f) const
++ inline void TrackMapper::add_to_set (SetVoxelDEC& out, const Eigen::Vector3i& v, const Eigen::Vector3d& d, const default_type l, const default_type f) const
+ {
+ out.insert (v, d, l, f);
+ }
+- inline void TrackMapper::add_to_set (SetDixel& out, const Eigen::Vector3i& v, const Eigen::Vector3& d, const default_type l, const default_type f) const
++ inline void TrackMapper::add_to_set (SetDixel& out, const Eigen::Vector3i& v, const Eigen::Vector3d& d, const default_type l, const default_type f) const
+ {
+ assert (dixel_plugin);
+ const size_t bin = (*dixel_plugin) (d);
+ out.insert (v, bin, l, f);
+ }
+- inline void TrackMapper::add_to_set (SetVoxelTOD& out, const Eigen::Vector3i& v, const Eigen::Vector3& d, const default_type l, const default_type f) const
++ inline void TrackMapper::add_to_set (SetVoxelTOD& out, const Eigen::Vector3i& v, const Eigen::Vector3d& d, const default_type l, const default_type f) const
+ {
+ assert (tod_plugin);
+ VoxelTOD::vector_type sh;
+Index: mrtrix3/src/dwi/tractography/mapping/gaussian/voxel.h
+===================================================================
+--- mrtrix3.orig/src/dwi/tractography/mapping/gaussian/voxel.h
++++ mrtrix3/src/dwi/tractography/mapping/gaussian/voxel.h
+@@ -80,18 +80,18 @@ namespace MR {
+ public:
+ VoxelDEC () : Base (), VoxelAddon () { }
+ VoxelDEC (const Eigen::Vector3i& V) : Base (V), VoxelAddon () { }
+- VoxelDEC (const Eigen::Vector3i& V, const Eigen::Vector3& d) : Base (V, d), VoxelAddon () { }
+- VoxelDEC (const Eigen::Vector3i& V, const Eigen::Vector3& d, const default_type l) : Base (V, d, l), VoxelAddon () { }
+- VoxelDEC (const Eigen::Vector3i& V, const Eigen::Vector3& d, const default_type l, const default_type f) : Base (V, d, l), VoxelAddon (f) { }
++ VoxelDEC (const Eigen::Vector3i& V, const Eigen::Vector3d& d) : Base (V, d), VoxelAddon () { }
++ VoxelDEC (const Eigen::Vector3i& V, const Eigen::Vector3d& d, const default_type l) : Base (V, d, l), VoxelAddon () { }
++ VoxelDEC (const Eigen::Vector3i& V, const Eigen::Vector3d& d, const default_type l, const default_type f) : Base (V, d, l), VoxelAddon (f) { }
+
+ VoxelDEC& operator= (const VoxelDEC& V) { Base::operator= (V); VoxelAddon::operator= (V); return (*this); }
+ void operator+= (const default_type) const { assert (0); }
+- void operator+= (const Eigen::Vector3&) const { assert (0); }
++ void operator+= (const Eigen::Vector3d&) const { assert (0); }
+ bool operator== (const VoxelDEC& V) const { return Base::operator== (V); }
+ bool operator< (const VoxelDEC& V) const { return Base::operator< (V); }
+
+- void add (const Eigen::Vector3&, const default_type) const { assert (0); }
+- void add (const Eigen::Vector3& i, const default_type l, const default_type f) const { Base::add (i, l); VoxelAddon::operator+= (f); }
++ void add (const Eigen::Vector3d&, const default_type) const { assert (0); }
++ void add (const Eigen::Vector3d& i, const default_type l, const default_type f) const { Base::add (i, l); VoxelAddon::operator+= (f); }
+ void normalize() const { VoxelAddon::normalize (get_length()); Base::normalize(); }
+
+ };
+@@ -191,7 +191,7 @@ namespace MR {
+
+ using VoxType = VoxelDEC;
+
+- inline void insert (const Eigen::Vector3i& v, const Eigen::Vector3& d, const default_type l, const default_type f)
++ inline void insert (const Eigen::Vector3i& v, const Eigen::Vector3d& d, const default_type l, const default_type f)
+ {
+ const VoxelDEC temp (v, d, l, f);
+ iterator existing = std::set<VoxelDEC>::find (temp);
+Index: mrtrix3/src/dwi/tractography/mapping/mapper.cpp
+===================================================================
+--- mrtrix3.orig/src/dwi/tractography/mapping/mapper.cpp
++++ mrtrix3/src/dwi/tractography/mapping/mapper.cpp
+@@ -254,7 +254,7 @@ void TrackMapperTWI::load_factors (const
+ if (contrast != CURVATURE)
+ throw Exception ("Unsupported contrast in function TrackMapperTWI::load_factors()");
+
+- vector<Eigen::Vector3> tangents;
++ vector<Eigen::Vector3d> tangents;
+ tangents.reserve (tck.size());
+
+ // Would like to be able to manipulate the length over which the tangent calculation is affected
+@@ -271,7 +271,7 @@ void TrackMapperTWI::load_factors (const
+ step_sizes.reserve (tck.size());
+
+ for (size_t i = 0; i != tck.size(); ++i) {
+- Eigen::Vector3 this_tangent;
++ Eigen::Vector3d this_tangent;
+ if (i == 0)
+ this_tangent = ((tck[1] - tck[0] ).cast<default_type>().normalized());
+ else if (i == tck.size() - 1)
+@@ -322,7 +322,7 @@ void TrackMapperTWI::load_factors (const
+ // Smooth both the tangent vectors and the principal normal vectors according to a Gaussuan kernel
+ // Remember: tangent vectors are unit length, but for principal normal vectors length must be preserved!
+
+- vector<Eigen::Vector3> smoothed_tangents;
++ vector<Eigen::Vector3d> smoothed_tangents;
+ smoothed_tangents.reserve (tangents.size());
+
+ static const default_type gaussian_theta = CURVATURE_TRACK_SMOOTHING_FWHM / (2.0 * sqrt (2.0 * log (2.0)));
+@@ -330,7 +330,7 @@ void TrackMapperTWI::load_factors (const
+
+ for (size_t i = 0; i != tck.size(); ++i) {
+
+- Eigen::Vector3 this_tangent (0.0, 0.0, 0.0);
++ Eigen::Vector3d this_tangent (0.0, 0.0, 0.0);
+
+ for (size_t j = 0; j != tck.size(); ++j) {
+ const default_type distance = spline_distances (i, j);
+Index: mrtrix3/src/dwi/tractography/mapping/mapper.h
+===================================================================
+--- mrtrix3.orig/src/dwi/tractography/mapping/mapper.h
++++ mrtrix3/src/dwi/tractography/mapping/mapper.h
+@@ -160,11 +160,11 @@ namespace MR {
+ virtual void postprocess (const Streamline<>& tck, SetVoxelExtras& out) const { }
+
+ // Used by voxelise() and voxelise_precise() to increment the relevant set
+- inline void add_to_set (SetVoxel& , const Eigen::Vector3i&, const Eigen::Vector3&, const default_type) const;
+- inline void add_to_set (SetVoxelDEC&, const Eigen::Vector3i&, const Eigen::Vector3&, const default_type) const;
+- inline void add_to_set (SetVoxelDir&, const Eigen::Vector3i&, const Eigen::Vector3&, const default_type) const;
+- inline void add_to_set (SetDixel& , const Eigen::Vector3i&, const Eigen::Vector3&, const default_type) const;
+- inline void add_to_set (SetVoxelTOD&, const Eigen::Vector3i&, const Eigen::Vector3&, const default_type) const;
++ inline void add_to_set (SetVoxel& , const Eigen::Vector3i&, const Eigen::Vector3d&, const default_type) const;
++ inline void add_to_set (SetVoxelDEC&, const Eigen::Vector3i&, const Eigen::Vector3d&, const default_type) const;
++ inline void add_to_set (SetVoxelDir&, const Eigen::Vector3i&, const Eigen::Vector3d&, const default_type) const;
++ inline void add_to_set (SetDixel& , const Eigen::Vector3i&, const Eigen::Vector3d&, const default_type) const;
++ inline void add_to_set (SetVoxelTOD&, const Eigen::Vector3i&, const Eigen::Vector3d&, const default_type) const;
+
+ DWI::Tractography::Resampling::Upsampler upsampler;
+
+@@ -183,7 +183,7 @@ namespace MR {
+ for (auto i = tck.cbegin(); i != last; ++i) {
+ vox = round (scanner2voxel * (*i));
+ if (check (vox, info)) {
+- const Eigen::Vector3 dir = (*(i+1) - *prev).cast<default_type>().normalized();
++ const Eigen::Vector3d dir = (*(i+1) - *prev).cast<default_type>().normalized();
+ if (dir.allFinite() && !dir.isZero())
+ add_to_set (output, vox, dir, 1.0);
+ }
+@@ -192,7 +192,7 @@ namespace MR {
+
+ vox = round (scanner2voxel * (*last));
+ if (check (vox, info)) {
+- const Eigen::Vector3 dir = (*last - *prev).cast<default_type>().normalized();
++ const Eigen::Vector3d dir = (*last - *prev).cast<default_type>().normalized();
+ if (dir.allFinite() && !dir.isZero())
+ add_to_set (output, vox, dir, 1.0);
+ }
+@@ -282,7 +282,7 @@ namespace MR {
+ }
+
+ length += (p_prev - p_voxel_exit).norm();
+- const Eigen::Vector3 traversal_vector = (p_voxel_exit - p_voxel_entry).cast<default_type>().normalized();
++ const Eigen::Vector3d traversal_vector = (p_voxel_exit - p_voxel_entry).cast<default_type>().normalized();
+ if (std::isfinite (traversal_vector[0]) && check (this_voxel, info))
+ add_to_set (out, this_voxel, traversal_vector, length);
+
+@@ -300,13 +300,13 @@ namespace MR {
+ if (tck.size() == 1) {
+ const auto vox = round (scanner2voxel * tck.front());
+ if (check (vox, info))
+- add_to_set (out, vox, Eigen::Vector3(NaN, NaN, NaN), 1.0);
++ add_to_set (out, vox, Eigen::Vector3d(NaN, NaN, NaN), 1.0);
+ return;
+ }
+ for (size_t end = 0; end != 2; ++end) {
+ const auto vox = round (scanner2voxel * (end ? tck.back() : tck.front()));
+ if (check (vox, info)) {
+- Eigen::Vector3 dir { NaN, NaN, NaN };
++ Eigen::Vector3d dir { NaN, NaN, NaN };
+ if (tck.size() > 1)
+ dir = (end ? (tck[tck.size()-1] - tck[tck.size()-2]) : (tck[0] - tck[1])).cast<default_type>().normalized();
+ add_to_set (out, vox, dir, 1.0);
+@@ -318,25 +318,25 @@ namespace MR {
+
+
+ // These are inlined to make as fast as possible
+- inline void TrackMapperBase::add_to_set (SetVoxel& out, const Eigen::Vector3i& v, const Eigen::Vector3& d, const default_type l) const
++ inline void TrackMapperBase::add_to_set (SetVoxel& out, const Eigen::Vector3i& v, const Eigen::Vector3d& d, const default_type l) const
+ {
+ out.insert (v, l);
+ }
+- inline void TrackMapperBase::add_to_set (SetVoxelDEC& out, const Eigen::Vector3i& v, const Eigen::Vector3& d, const default_type l) const
++ inline void TrackMapperBase::add_to_set (SetVoxelDEC& out, const Eigen::Vector3i& v, const Eigen::Vector3d& d, const default_type l) const
+ {
+ out.insert (v, d, l);
+ }
+- inline void TrackMapperBase::add_to_set (SetVoxelDir& out, const Eigen::Vector3i& v, const Eigen::Vector3& d, const default_type l) const
++ inline void TrackMapperBase::add_to_set (SetVoxelDir& out, const Eigen::Vector3i& v, const Eigen::Vector3d& d, const default_type l) const
+ {
+ out.insert (v, d, l);
+ }
+- inline void TrackMapperBase::add_to_set (SetDixel& out, const Eigen::Vector3i& v, const Eigen::Vector3& d, const default_type l) const
++ inline void TrackMapperBase::add_to_set (SetDixel& out, const Eigen::Vector3i& v, const Eigen::Vector3d& d, const default_type l) const
+ {
+ assert (dixel_plugin);
+ const DWI::Directions::index_type bin = (*dixel_plugin) (d);
+ out.insert (v, bin, l);
+ }
+- inline void TrackMapperBase::add_to_set (SetVoxelTOD& out, const Eigen::Vector3i& v, const Eigen::Vector3& d, const default_type l) const
++ inline void TrackMapperBase::add_to_set (SetVoxelTOD& out, const Eigen::Vector3i& v, const Eigen::Vector3d& d, const default_type l) const
+ {
+ assert (tod_plugin);
+ Eigen::Matrix<default_type, Eigen::Dynamic, 1> sh;
+Index: mrtrix3/src/dwi/tractography/mapping/mapper_plugins.cpp
+===================================================================
+--- mrtrix3.orig/src/dwi/tractography/mapping/mapper_plugins.cpp
++++ mrtrix3/src/dwi/tractography/mapping/mapper_plugins.cpp
+@@ -65,7 +65,7 @@ namespace MR {
+ if (statistic == ENDS_CORR) {
+
+ for (; index >= 0 && index < ssize_t(tck.size()); index += step) {
+- const Eigen::Vector3 p = interp.scanner2voxel * tck[index].cast<default_type>();
++ const Eigen::Vector3d p = interp.scanner2voxel * tck[index].cast<default_type>();
+ const Eigen::Array3i v ( { int(std::round (p[0])), int(std::round (p[1])), int(std::round (p[2])) } );
+ if (!is_out_of_bounds (backtrack_mask, v)) {
+ assign_pos_of (v, 0, 3).to (backtrack_mask);
+@@ -152,7 +152,7 @@ namespace MR {
+ if (interp.scanner (tck[index])) {
+ for (interp.index(3) = 0; interp.index(3) != interp.size(3); ++interp.index(3))
+ sh_coeffs[interp.index(3)] = interp.value();
+- const Eigen::Vector3 dir = (tck[(index == ssize_t(tck.size()-1)) ? index : (index+1)] - tck[index ? (index-1) : 0]).cast<default_type>().normalized();
++ const Eigen::Vector3d dir = (tck[(index == ssize_t(tck.size()-1)) ? index : (index+1)] - tck[index ? (index-1) : 0]).cast<default_type>().normalized();
+ factors.push_back (precomputer->value (sh_coeffs, dir));
+ } else {
+ factors.push_back (NaN);
+@@ -170,7 +170,7 @@ namespace MR {
+ for (interp.index(3) = 0; interp.index(3) != interp.size(3); ++interp.index(3))
+ sh_coeffs[interp.index(3)] = interp.value();
+ // Get the FOD amplitude along the streamline tangent
+- const Eigen::Vector3 dir = (tck[(i == tck.size()-1) ? i : (i+1)] - tck[i ? (i-1) : 0]).cast<default_type>().normalized();
++ const Eigen::Vector3d dir = (tck[(i == tck.size()-1) ? i : (i+1)] - tck[i ? (i-1) : 0]).cast<default_type>().normalized();
+ factors.push_back (precomputer->value (sh_coeffs, dir));
+ } else {
+ factors.push_back (NaN);
+Index: mrtrix3/src/dwi/tractography/mapping/mapper_plugins.h
+===================================================================
+--- mrtrix3.orig/src/dwi/tractography/mapping/mapper_plugins.h
++++ mrtrix3/src/dwi/tractography/mapping/mapper_plugins.h
+@@ -45,7 +45,7 @@ namespace MR {
+ dirs (directions) { }
+ DixelMappingPlugin (const DixelMappingPlugin& that) :
+ dirs (that.dirs) { }
+- DWI::Directions::index_type operator() (const Eigen::Vector3& d) const { return dirs.select_direction (d); }
++ DWI::Directions::index_type operator() (const Eigen::Vector3d& d) const { return dirs.select_direction (d); }
+ private:
+ const DWI::Directions::FastLookupSet& dirs;
+ };
+Index: mrtrix3/src/dwi/tractography/mapping/voxel.h
+===================================================================
+--- mrtrix3.orig/src/dwi/tractography/mapping/voxel.h
++++ mrtrix3/src/dwi/tractography/mapping/voxel.h
+@@ -47,7 +47,7 @@ namespace MR {
+ return (V[0] >= 0 && V[0] < H.size(0) && V[1] >= 0 && V[1] < H.size(1) && V[2] >= 0 && V[2] < H.size(2));
+ }
+
+- inline Eigen::Vector3 vec2DEC (const Eigen::Vector3& d)
++ inline Eigen::Vector3d vec2DEC (const Eigen::Vector3d& d)
+ {
+ return { abs(d[0]), abs(d[1]), abs(d[2]) };
+ }
+@@ -85,11 +85,11 @@ namespace MR {
+ Voxel (V),
+ colour (0.0, 0.0, 0.0) { }
+
+- VoxelDEC (const Eigen::Vector3i& V, const Eigen::Vector3& d) :
++ VoxelDEC (const Eigen::Vector3i& V, const Eigen::Vector3d& d) :
+ Voxel (V),
+ colour (vec2DEC (d)) { }
+
+- VoxelDEC (const Eigen::Vector3i& V, const Eigen::Vector3& d, const float l) :
++ VoxelDEC (const Eigen::Vector3i& V, const Eigen::Vector3d& d, const float l) :
+ Voxel (V, l),
+ colour (vec2DEC (d)) { }
+
+@@ -101,13 +101,13 @@ namespace MR {
+ bool operator< (const VoxelDEC& V) const { return Voxel::operator< (V); }
+
+ void normalize() const { colour.normalize(); Voxel::normalize(); }
+- void set_dir (const Eigen::Vector3& i) { colour = vec2DEC (i); }
+- void add (const Eigen::Vector3& i, const default_type l) const { Voxel::operator+= (l); colour += vec2DEC (i); }
+- void operator+= (const Eigen::Vector3& i) const { Voxel::operator+= (1.0); colour += vec2DEC (i); }
+- const Eigen::Vector3& get_colour() const { return colour; }
++ void set_dir (const Eigen::Vector3d& i) { colour = vec2DEC (i); }
++ void add (const Eigen::Vector3d& i, const default_type l) const { Voxel::operator+= (l); colour += vec2DEC (i); }
++ void operator+= (const Eigen::Vector3d& i) const { Voxel::operator+= (1.0); colour += vec2DEC (i); }
++ const Eigen::Vector3d& get_colour() const { return colour; }
+
+ private:
+- mutable Eigen::Vector3 colour;
++ mutable Eigen::Vector3d colour;
+
+ };
+
+@@ -127,11 +127,11 @@ namespace MR {
+ Voxel (V),
+ dir (0.0, 0.0, 0.0) { }
+
+- VoxelDir (const Eigen::Vector3i& V, const Eigen::Vector3& d) :
++ VoxelDir (const Eigen::Vector3i& V, const Eigen::Vector3d& d) :
+ Voxel (V),
+ dir (d) { }
+
+- VoxelDir (const Eigen::Vector3i& V, const Eigen::Vector3& d, const default_type l) :
++ VoxelDir (const Eigen::Vector3i& V, const Eigen::Vector3d& d, const default_type l) :
+ Voxel (V, l),
+ dir (d) { }
+
+@@ -142,13 +142,13 @@ namespace MR {
+ bool operator< (const VoxelDir& V) const { return Voxel::operator< (V); }
+
+ void normalize() const { dir.normalize(); Voxel::normalize(); }
+- void set_dir (const Eigen::Vector3& i) { dir = i; }
+- void add (const Eigen::Vector3& i, const default_type l) const { Voxel::operator+= (l); dir += i * (dir.dot(i) < 0.0 ? -1.0 : 1.0); }
+- void operator+= (const Eigen::Vector3& i) const { Voxel::operator+= (1.0); dir += i * (dir.dot(i) < 0.0 ? -1.0 : 1.0); }
+- const Eigen::Vector3& get_dir() const { return dir; }
++ void set_dir (const Eigen::Vector3d& i) { dir = i; }
++ void add (const Eigen::Vector3d& i, const default_type l) const { Voxel::operator+= (l); dir += i * (dir.dot(i) < 0.0 ? -1.0 : 1.0); }
++ void operator+= (const Eigen::Vector3d& i) const { Voxel::operator+= (1.0); dir += i * (dir.dot(i) < 0.0 ? -1.0 : 1.0); }
++ const Eigen::Vector3d& get_dir() const { return dir; }
+
+ private:
+- mutable Eigen::Vector3 dir;
++ mutable Eigen::Vector3d dir;
+
+ };
+
+@@ -317,12 +317,12 @@ namespace MR {
+ else
+ existing->add (v.get_colour(), v.get_length());
+ }
+- inline void insert (const Eigen::Vector3i& v, const Eigen::Vector3& d)
++ inline void insert (const Eigen::Vector3i& v, const Eigen::Vector3d& d)
+ {
+ const VoxelDEC temp (v, d);
+ insert (temp);
+ }
+- inline void insert (const Eigen::Vector3i& v, const Eigen::Vector3& d, const default_type l)
++ inline void insert (const Eigen::Vector3i& v, const Eigen::Vector3d& d, const default_type l)
+ {
+ const VoxelDEC temp (v, d, l);
+ insert (temp);
+@@ -344,12 +344,12 @@ namespace MR {
+ else
+ existing->add (v.get_dir(), v.get_length());
+ }
+- inline void insert (const Eigen::Vector3i& v, const Eigen::Vector3& d)
++ inline void insert (const Eigen::Vector3i& v, const Eigen::Vector3d& d)
+ {
+ const VoxelDir temp (v, d);
+ insert (temp);
+ }
+- inline void insert (const Eigen::Vector3i& v, const Eigen::Vector3& d, const default_type l)
++ inline void insert (const Eigen::Vector3i& v, const Eigen::Vector3d& d, const default_type l)
+ {
+ const VoxelDir temp (v, d, l);
+ insert (temp);
+Index: mrtrix3/src/dwi/tractography/mapping/writer.h
+===================================================================
+--- mrtrix3.orig/src/dwi/tractography/mapping/writer.h
++++ mrtrix3/src/dwi/tractography/mapping/writer.h
+@@ -277,8 +277,8 @@ namespace MR {
+
+
+ // Convenience functions for Directionally-Encoded Colour processing
+- Eigen::Vector3 get_dec ();
+- void set_dec (const Eigen::Vector3&);
++ Eigen::Vector3d get_dec ();
++ void set_dec (const Eigen::Vector3d&);
+
+ // Convenience functions for Track Orientation Distribution processing
+ void get_tod ( VoxelTOD::vector_type&);
+@@ -458,10 +458,10 @@ namespace MR {
+
+
+ template <typename value_type>
+- Eigen::Vector3 MapWriter<value_type>::get_dec ()
++ Eigen::Vector3d MapWriter<value_type>::get_dec ()
+ {
+ assert (type == DEC);
+- Eigen::Vector3 value;
++ Eigen::Vector3d value;
+ buffer.index(3) = 0; value[0] = buffer.value();
+ buffer.index(3)++; value[1] = buffer.value();
+ buffer.index(3)++; value[2] = buffer.value();
+@@ -469,7 +469,7 @@ namespace MR {
+ }
+
+ template <typename value_type>
+- void MapWriter<value_type>::set_dec (const Eigen::Vector3& value)
++ void MapWriter<value_type>::set_dec (const Eigen::Vector3d& value)
+ {
+ assert (type == DEC);
+ buffer.index(3) = 0; buffer.value() = value[0];
+Index: mrtrix3/src/fixel/filter/smooth.cpp
+===================================================================
+--- mrtrix3.orig/src/fixel/filter/smooth.cpp
++++ mrtrix3/src/fixel/filter/smooth.cpp
+@@ -46,7 +46,7 @@ namespace MR
+ fixel_positions.resize (matrix.size());
+ const Transform transform (index_image);
+ for (auto i = Loop (index_image, 0, 3) (index_image); i; ++i) {
+- const Eigen::Vector3 vox ((default_type)index_image.index(0), (default_type)index_image.index(1), (default_type)index_image.index(2));
++ const Eigen::Vector3d vox ((default_type)index_image.index(0), (default_type)index_image.index(1), (default_type)index_image.index(2));
+ const Eigen::Vector3f scanner = (transform.voxel2scanner * vox).cast<float>();
+ index_image.index(3) = 0;
+ const index_type count = index_image.value();
+Index: mrtrix3/src/fixel/matrix.cpp
+===================================================================
+--- mrtrix3.orig/src/fixel/matrix.cpp
++++ mrtrix3/src/fixel/matrix.cpp
+@@ -147,7 +147,7 @@ namespace MR
+ bool operator() (const DWI::Tractography::Streamline<>& tck,
+ vector<index_type>& out) const
+ {
+- using direction_type = Eigen::Vector3;
++ using direction_type = Eigen::Vector3d;
+ using SetVoxelDir = DWI::Tractography::Mapping::SetVoxelDir;
+
+ SetVoxelDir in;
+Index: mrtrix3/src/gui/dwi/renderer.cpp
+===================================================================
+--- mrtrix3.orig/src/gui/dwi/renderer.cpp
++++ mrtrix3/src/gui/dwi/renderer.cpp
+@@ -623,14 +623,14 @@ namespace MR
+ if (I == i) {
+
+ // Invert a direction if required
+- std::array<Eigen::Vector3, 3> d {{ dirs[i], dirs[j], dirs[k] }};
+- const Eigen::Vector3 mean_dir ((d[0]+d[1]+d[2]).normalized());
++ std::array<Eigen::Vector3d, 3> d {{ dirs[i], dirs[j], dirs[k] }};
++ const Eigen::Vector3d mean_dir ((d[0]+d[1]+d[2]).normalized());
+ for (size_t v = 0; v != 3; ++v) {
+ if (d[v].dot (mean_dir) < 0.0)
+ d[v] = -d[v];
+ }
+ // Conform to right hand rule
+- const Eigen::Vector3 normal (((d[1]-d[0]).cross (d[2]-d[1])).normalized());
++ const Eigen::Vector3d normal (((d[1]-d[0]).cross (d[2]-d[1])).normalized());
+ if (normal.dot (mean_dir) < 0.0)
+ indices_data.push_back ( {{GLint(i), GLint(k), GLint(j)}} );
+ else
+Index: mrtrix3/src/registration/linear.h
+===================================================================
+--- mrtrix3.orig/src/registration/linear.h
++++ mrtrix3/src/registration/linear.h
+@@ -442,7 +442,7 @@ namespace MR
+
+ // set control point coordinates inside +-1/3 of the midway_image size
+ {
+- Eigen::Vector3 ext (midway_image_header.spacing(0) / 6.0,
++ Eigen::Vector3d ext (midway_image_header.spacing(0) / 6.0,
+ midway_image_header.spacing(1) / 6.0,
+ midway_image_header.spacing(2) / 6.0);
+ for (size_t i = 0; i<3; ++i)
+Index: mrtrix3/src/registration/metric/cc_helper.h
+===================================================================
+--- mrtrix3.orig/src/registration/metric/cc_helper.h
++++ mrtrix3/src/registration/metric/cc_helper.h
+@@ -45,7 +45,7 @@ namespace MR
+ int nmax = extent[0] * extent[1] * extent[2];
+ Eigen::VectorXd n1 = Eigen::VectorXd(nmax);
+ Eigen::VectorXd n2 = Eigen::VectorXd(nmax);
+- Eigen::Vector3 pos;
++ Eigen::Vector3d pos;
+ for (auto l = Loop("precomputing cross correlation values") (im1_image); l; ++l) {
+ pos[0] = im1_image.index(0);
+ pos[1] = im1_image.index(1);
+Index: mrtrix3/src/registration/metric/cross_correlation.h
+===================================================================
+--- mrtrix3.orig/src/registration/metric/cross_correlation.h
++++ mrtrix3/src/registration/metric/cross_correlation.h
+@@ -39,7 +39,7 @@ namespace MR
+ const Iterator& iter,
+ Eigen::Matrix<default_type, Eigen::Dynamic, 1>& gradient) {
+
+- // const Eigen::Vector3 pos = Eigen::Vector3 (iter.index(0), iter.index(1), iter.index(2));
++ // const Eigen::Vector3d pos = Eigen::Vector3d (iter.index(0), iter.index(1), iter.index(2));
+
+ assert (params.processed_mask.valid());
+ assert (params.processed_image.valid());
+Index: mrtrix3/src/registration/metric/demons.h
+===================================================================
+--- mrtrix3.orig/src/registration/metric/demons.h
++++ mrtrix3/src/registration/metric/demons.h
+@@ -118,8 +118,8 @@ namespace MR
+ im1_update.row(3) = 0.0;
+ im2_update.row(3) = 0.0;
+ } else {
+- im1_update.row(3) = Eigen::Vector3 (speed * grad.array() / denominator);
+- im2_update.row(3) = -Eigen::Vector3 (im1_update.row(3));
++ im1_update.row(3) = Eigen::Vector3d (speed * grad.array() / denominator);
++ im2_update.row(3) = -Eigen::Vector3d (im1_update.row(3));
+ }
+ }
+
+Index: mrtrix3/src/registration/metric/difference_robust.h
+===================================================================
+--- mrtrix3.orig/src/registration/metric/difference_robust.h
++++ mrtrix3/src/registration/metric/difference_robust.h
+@@ -35,9 +35,9 @@ namespace MR
+
+ template <class Params>
+ default_type operator() (Params& params,
+- const Eigen::Vector3 im1_point,
+- const Eigen::Vector3 im2_point,
+- const Eigen::Vector3 midway_point,
++ const Eigen::Vector3d im1_point,
++ const Eigen::Vector3d im2_point,
++ const Eigen::Vector3d midway_point,
+ Eigen::Matrix<default_type, Eigen::Dynamic, 1>& gradient) {
+
+ assert (!this->weighted && "FIXME: set_weights not implemented for 3D metric");
+@@ -103,9 +103,9 @@ namespace MR
+
+ template <class Params>
+ default_type operator() (Params& params,
+- const Eigen::Vector3& im1_point,
+- const Eigen::Vector3& im2_point,
+- const Eigen::Vector3& midway_point,
++ const Eigen::Vector3d& im1_point,
++ const Eigen::Vector3d& im2_point,
++ const Eigen::Vector3d& midway_point,
+ Eigen::Matrix<default_type, Eigen::Dynamic, 1>& gradient) {
+
+ params.im1_image_interp->value_and_gradient_row_wrt_scanner (im1_values, im1_grad);
+Index: mrtrix3/src/registration/metric/linear_base.h
+===================================================================
+--- mrtrix3.orig/src/registration/metric/linear_base.h
++++ mrtrix3/src/registration/metric/linear_base.h
+@@ -56,9 +56,9 @@ namespace MR
+
+ template <class Params>
+ default_type operator() (Params& params,
+- const Eigen::Vector3& im1_point,
+- const Eigen::Vector3& im2_point,
+- const Eigen::Vector3& midway_point,
++ const Eigen::Vector3d& im1_point,
++ const Eigen::Vector3d& im2_point,
++ const Eigen::Vector3d& midway_point,
+ Eigen::Matrix<default_type, Eigen::Dynamic, 1>& gradient);
+
+ protected:
+Index: mrtrix3/src/registration/metric/local_cross_correlation.h
+===================================================================
+--- mrtrix3.orig/src/registration/metric/local_cross_correlation.h
++++ mrtrix3/src/registration/metric/local_cross_correlation.h
+@@ -36,7 +36,7 @@ namespace MR
+ void operator() (MaskType& mask, ImageType3& out) {
+ if (!mask.value())
+ return;
+- Eigen::Vector3 pos (mask.index(0), mask.index(1), mask.index(2));
++ Eigen::Vector3d pos (mask.index(0), mask.index(1), mask.index(2));
+ out.index(0) = pos[0];
+ out.index(1) = pos[1];
+ out.index(2) = pos[2];
+@@ -237,7 +237,7 @@ namespace MR
+ return 0.0;
+ }
+
+- const Eigen::Vector3 pos = Eigen::Vector3(default_type(iter.index(0)), default_type(iter.index(0)), default_type(iter.index(0)));
++ const Eigen::Vector3d pos = Eigen::Vector3d(default_type(iter.index(0)), default_type(iter.index(0)), default_type(iter.index(0)));
+ params.processed_image_interp->voxel(pos);
+ typename Params::Im1ValueType val1;
+ typename Params::Im2ValueType val2;
+@@ -266,9 +266,9 @@ namespace MR
+
+ // ITK:
+ // derivWRTImage[dim] = 2.0 * sFixedMoving / (sFixedFixed_sMovingMoving) * (fixedI - sFixedMoving / sMovingMoving * movingI) * movingImageGradient[dim];
+- Eigen::Vector3 derivWRTImage = - A_BC * ((val2 - A/B * val1) * grad1 - 0.0 * (val1 - A/C * val2) * grad2);
++ Eigen::Vector3d derivWRTImage = - A_BC * ((val2 - A/B * val1) * grad1 - 0.0 * (val1 - A/C * val2) * grad2);
+
+- const Eigen::Vector3 midway_point = midway_v2s * pos;
++ const Eigen::Vector3d midway_point = midway_v2s * pos;
+ const auto jacobian_vec = params.transformation.get_jacobian_vector_wrt_params (midway_point);
+ gradient.segment<4>(0) += derivWRTImage(0) * jacobian_vec;
+ gradient.segment<4>(4) += derivWRTImage(1) * jacobian_vec;
+Index: mrtrix3/src/registration/metric/mean_squared.h
+===================================================================
+--- mrtrix3.orig/src/registration/metric/mean_squared.h
++++ mrtrix3/src/registration/metric/mean_squared.h
+@@ -32,9 +32,9 @@ namespace MR
+ public:
+ template <class Params>
+ default_type operator() (Params& params,
+- const Eigen::Vector3& im1_point,
+- const Eigen::Vector3& im2_point,
+- const Eigen::Vector3& midway_point,
++ const Eigen::Vector3d& im1_point,
++ const Eigen::Vector3d& im2_point,
++ const Eigen::Vector3d& midway_point,
+ Eigen::Matrix<default_type, Eigen::Dynamic, 1>& gradient) {
+
+ assert (!this->weighted && "FIXME: set_weights not implemented for 3D metric");
+@@ -67,9 +67,9 @@ namespace MR
+ public:
+ template <class Params>
+ default_type operator() (Params& params,
+- const Eigen::Vector3& im1_point,
+- const Eigen::Vector3& im2_point,
+- const Eigen::Vector3& midway_point,
++ const Eigen::Vector3d& im1_point,
++ const Eigen::Vector3d& im2_point,
++ const Eigen::Vector3d& midway_point,
+ Eigen::Matrix<default_type, Eigen::Dynamic, 1>& gradient) {
+
+ assert (!this->weighted && "FIXME: set_weights not implemented for 3D metric");
+@@ -99,9 +99,9 @@ namespace MR
+
+ template <class Params>
+ default_type operator() (Params& params,
+- const Eigen::Vector3& im1_point,
+- const Eigen::Vector3& im2_point,
+- const Eigen::Vector3& midway_point,
++ const Eigen::Vector3d& im1_point,
++ const Eigen::Vector3d& im2_point,
++ const Eigen::Vector3d& midway_point,
+ Eigen::Matrix<default_type, Eigen::Dynamic, 1>& gradient) {
+
+ assert (!this->weighted && "FIXME: set_weights not implemented for 3D metric");
+@@ -139,9 +139,9 @@ namespace MR
+ public:
+ template <class Params>
+ default_type operator() (Params& params,
+- const Eigen::Vector3& im1_point,
+- const Eigen::Vector3& im2_point,
+- const Eigen::Vector3& midway_point,
++ const Eigen::Vector3d& im1_point,
++ const Eigen::Vector3d& im2_point,
++ const Eigen::Vector3d& midway_point,
+ Eigen::Matrix<default_type, Eigen::Dynamic, 1>& gradient) {
+
+ const ssize_t volumes = params.im1_image_interp->size(3);
+@@ -192,9 +192,9 @@ namespace MR
+
+ template <class Params>
+ default_type operator() (Params& params,
+- const Eigen::Vector3& im1_point,
+- const Eigen::Vector3& im2_point,
+- const Eigen::Vector3& midway_point,
++ const Eigen::Vector3d& im1_point,
++ const Eigen::Vector3d& im2_point,
++ const Eigen::Vector3d& midway_point,
+ Eigen::Matrix<default_type, Eigen::Dynamic, 1>& gradient) {
+ const ssize_t volumes = params.im1_image_interp->size(3);
+
+@@ -229,9 +229,9 @@ namespace MR
+ public:
+ template <class Params>
+ default_type operator() (Params& params,
+- const Eigen::Vector3& im1_point,
+- const Eigen::Vector3& im2_point,
+- const Eigen::Vector3& midway_point,
++ const Eigen::Vector3d& im1_point,
++ const Eigen::Vector3d& im2_point,
++ const Eigen::Vector3d& midway_point,
+ Eigen::Matrix<default_type, Eigen::Dynamic, 1>& gradient) {
+
+ const ssize_t volumes = params.im1_image_interp->size(3);
+@@ -302,9 +302,9 @@ namespace MR
+
+ template <class Params>
+ Eigen::Matrix<default_type, Eigen::Dynamic, 1> operator() (Params& params,
+- const Eigen::Vector3& im1_point,
+- const Eigen::Vector3& im2_point,
+- const Eigen::Vector3& midway_point,
++ const Eigen::Vector3d& im1_point,
++ const Eigen::Vector3d& im2_point,
++ const Eigen::Vector3d& midway_point,
+ Eigen::Matrix<default_type, Eigen::Dynamic, 1>& gradient) {
+
+ assert (volumes == params.im1_image_interp->size(3));
+Index: mrtrix3/src/registration/metric/params.h
+===================================================================
+--- mrtrix3.orig/src/registration/metric/params.h
++++ mrtrix3/src/registration/metric/params.h
+@@ -131,7 +131,7 @@ namespace MR
+ }
+
+ void update_control_points () {
+- const Eigen::Vector3 centre = transformation.get_centre();
++ const Eigen::Vector3d centre = transformation.get_centre();
+ control_points.resize(4, 4);
+ // tetrahedron centred at centre of midspace scaled by control_point_exent
+ control_points << 1.0, -1.0, -1.0, 1.0,
+@@ -181,10 +181,10 @@ namespace MR
+ im2_image, midway_image, trafo2, no_oversampling, NAN);
+
+ auto T = MR::Transform(midway_image).voxel2scanner;
+- Eigen::Vector3 midway_point, voxel_pos, im1_point, im2_point;
++ Eigen::Vector3d midway_point, voxel_pos, im1_point, im2_point;
+
+ for (auto i = Loop (midway_image) (check, im1_reslicer, im2_reslicer); i; ++i) {
+- voxel_pos = Eigen::Vector3 ((default_type)check.index(0), (default_type)check.index(1), (default_type)check.index(2));
++ voxel_pos = Eigen::Vector3d ((default_type)check.index(0), (default_type)check.index(1), (default_type)check.index(2));
+ midway_point = T * voxel_pos;
+
+ check.index(3) = 0;
+@@ -237,7 +237,7 @@ namespace MR
+ MR::copy_ptr<Im1MaskInterpolatorType> im1_mask_interp;
+ MR::copy_ptr<Im2MaskInterpolatorType> im2_mask_interp;
+ default_type loop_density;
+- Eigen::Vector3 control_point_exent;
++ Eigen::Vector3d control_point_exent;
+
+ bool robust_estimate_subset;
+ bool robust_estimate_use_score;
+Index: mrtrix3/src/registration/metric/thread_kernel.h
+===================================================================
+--- mrtrix3.orig/src/registration/metric/thread_kernel.h
++++ mrtrix3/src/registration/metric/thread_kernel.h
+@@ -125,12 +125,12 @@ namespace MR
+ typename cost_is_vector<U>::no = 0,
+ typename is_asymmetric<U>::no = 0) {
+
+- Eigen::Vector3 voxel_pos ((default_type)iter.index(0), (default_type)iter.index(1), (default_type)iter.index(2));
+- Eigen::Vector3 midway_point = voxel2scanner * voxel_pos;
++ Eigen::Vector3d voxel_pos ((default_type)iter.index(0), (default_type)iter.index(1), (default_type)iter.index(2));
++ Eigen::Vector3d midway_point = voxel2scanner * voxel_pos;
+
+
+
+- Eigen::Vector3 im2_point;
++ Eigen::Vector3d im2_point;
+ params.transformation.transform_half_inverse (im2_point, midway_point);
+ if (params.im2_mask_interp) {
+ params.im2_mask_interp->scanner (im2_point);
+@@ -143,7 +143,7 @@ namespace MR
+ return;
+ }
+
+- Eigen::Vector3 im1_point;
++ Eigen::Vector3d im1_point;
+ params.transformation.transform_half (im1_point, midway_point);
+ if (params.im1_mask_interp) {
+ params.im1_mask_interp->scanner (im1_point);
+@@ -176,8 +176,8 @@ namespace MR
+ typename cost_is_vector<U>::no = 0,
+ typename is_asymmetric<U>::yes = 0) {
+
+- Eigen::Vector3 voxel_pos ((default_type)iter.index(0), (default_type)iter.index(1), (default_type)iter.index(2));
+- Eigen::Vector3 im2_point = voxel2scanner * voxel_pos; // image 2 == midway_point == fixed image
++ Eigen::Vector3d voxel_pos ((default_type)iter.index(0), (default_type)iter.index(1), (default_type)iter.index(2));
++ Eigen::Vector3d im2_point = voxel2scanner * voxel_pos; // image 2 == midway_point == fixed image
+
+ // shift voxel position as evaluate iterates over a subset of the image
+ if (params.robust_estimate_subset) {
+@@ -212,7 +212,7 @@ namespace MR
+ return;
+ }
+
+- Eigen::Vector3 im1_point; // moving
++ Eigen::Vector3d im1_point; // moving
+ params.transformation.transform_half (im1_point, im2_point); // transform_half is full transformation, transform_half_inverse is identity
+ if (params.im1_mask_interp) {
+ params.im1_mask_interp->scanner (im1_point);
+@@ -240,11 +240,11 @@ namespace MR
+ typename cost_is_vector<U>::yes = 0,
+ typename is_asymmetric<U>::no = 0) {
+
+- Eigen::Vector3 voxel_pos ((default_type)iter.index(0), (default_type)iter.index(1), (default_type)iter.index(2));
++ Eigen::Vector3d voxel_pos ((default_type)iter.index(0), (default_type)iter.index(1), (default_type)iter.index(2));
+
+- Eigen::Vector3 midway_point = voxel2scanner * voxel_pos;
++ Eigen::Vector3d midway_point = voxel2scanner * voxel_pos;
+
+- Eigen::Vector3 im2_point;
++ Eigen::Vector3d im2_point;
+ params.transformation.transform_half_inverse (im2_point, midway_point);
+ if (params.im2_mask_interp) {
+ params.im2_mask_interp->scanner (im2_point);
+@@ -252,7 +252,7 @@ namespace MR
+ return;
+ }
+
+- Eigen::Vector3 im1_point;
++ Eigen::Vector3d im1_point;
+ params.transformation.transform_half (im1_point, midway_point);
+ if (params.im1_mask_interp) {
+ params.im1_mask_interp->scanner (im1_point);
+@@ -298,7 +298,7 @@ namespace MR
+ typename is_asymmetric<U>::no = 0) {
+ assert(params.processed_image.valid());
+
+- Eigen::Vector3 voxel_pos ((default_type)iter.index(0), (default_type)iter.index(1), (default_type)iter.index(2));
++ Eigen::Vector3d voxel_pos ((default_type)iter.index(0), (default_type)iter.index(1), (default_type)iter.index(2));
+
+ if (params.processed_mask.valid()){
+ assign_pos_of (iter, 0, 3).to (params.processed_mask);
+Index: mrtrix3/src/registration/transform/affine.cpp
+===================================================================
+--- mrtrix3.orig/src/registration/transform/affine.cpp
++++ mrtrix3/src/registration/transform/affine.cpp
+@@ -248,14 +248,14 @@ namespace MR
+ */
+
+
+- Eigen::Matrix<default_type, 4, 1> Affine::get_jacobian_vector_wrt_params (const Eigen::Vector3& p) const {
++ Eigen::Matrix<default_type, 4, 1> Affine::get_jacobian_vector_wrt_params (const Eigen::Vector3d& p) const {
+ Eigen::Matrix<default_type, 4, 1> jac;
+ jac.head(3) = p - centre;
+ jac(3) = 1.0;
+ return jac;
+ }
+
+- Eigen::MatrixXd Affine::get_jacobian_wrt_params (const Eigen::Vector3& p) const {
++ Eigen::MatrixXd Affine::get_jacobian_wrt_params (const Eigen::Vector3d& p) const {
+ Eigen::MatrixXd jacobian (3, 12);
+ jacobian.setZero();
+ const auto v = get_jacobian_vector_wrt_params(p);
+Index: mrtrix3/src/registration/transform/affine.h
+===================================================================
+--- mrtrix3.orig/src/registration/transform/affine.h
++++ mrtrix3/src/registration/transform/affine.h
+@@ -110,9 +110,9 @@ namespace MR
+ this->optimiser_weights << weights, weights, weights;
+ }
+
+- Eigen::Matrix<default_type, 4, 1> get_jacobian_vector_wrt_params (const Eigen::Vector3& p) const ;
++ Eigen::Matrix<default_type, 4, 1> get_jacobian_vector_wrt_params (const Eigen::Vector3d& p) const ;
+
+- Eigen::MatrixXd get_jacobian_wrt_params (const Eigen::Vector3& p) const ;
++ Eigen::MatrixXd get_jacobian_wrt_params (const Eigen::Vector3d& p) const ;
+
+ void set_parameter_vector (const Eigen::Matrix<ParameterType, Eigen::Dynamic, 1>& param_vector);
+
+Index: mrtrix3/src/registration/transform/base.h
+===================================================================
+--- mrtrix3.orig/src/registration/transform/base.h
++++ mrtrix3/src/registration/transform/base.h
+@@ -118,7 +118,7 @@ namespace MR
+ return number_of_parameters;
+ }
+
+- Eigen::Matrix<default_type, 4, 1> get_jacobian_vector_wrt_params (const Eigen::Vector3& p) const {
++ Eigen::Matrix<default_type, 4, 1> get_jacobian_vector_wrt_params (const Eigen::Vector3d& p) const {
+ throw Exception ("FIXME: get_jacobian_vector_wrt_params not implemented for this metric");
+ Eigen::Matrix<default_type, 4, 1> jac;
+ return jac;
+@@ -173,23 +173,23 @@ namespace MR
+ compute_halfspace_transformations();
+ }
+
+- const Eigen::Vector3 get_translation() const {
++ const Eigen::Vector3d get_translation() const {
+ return trafo.translation();
+ }
+
+- void set_centre_without_transform_update (const Eigen::Vector3& centre_in) {
++ void set_centre_without_transform_update (const Eigen::Vector3d& centre_in) {
+ centre = centre_in;
+ DEBUG ("centre: " + str(centre.transpose()));
+ }
+
+- void set_centre (const Eigen::Vector3& centre_in) {
++ void set_centre (const Eigen::Vector3d& centre_in) {
+ centre = centre_in;
+ DEBUG ("centre: " + str(centre.transpose()));
+ compute_offset();
+ compute_halfspace_transformations();
+ }
+
+- const Eigen::Vector3 get_centre() const {
++ const Eigen::Vector3d get_centre() const {
+ return centre;
+ }
+
+@@ -212,7 +212,7 @@ namespace MR
+ nonsymmetric = asym;
+ }
+
+- void set_offset (const Eigen::Vector3& offset_in) {
++ void set_offset (const Eigen::Vector3d& offset_in) {
+ trafo.translation() = offset_in;
+ compute_halfspace_transformations();
+ }
+@@ -280,7 +280,7 @@ namespace MR
+ Eigen::Transform<ParameterType, 3, Eigen::AffineCompact> trafo;
+ Eigen::Transform<ParameterType, 3, Eigen::AffineCompact> trafo_half;
+ Eigen::Transform<ParameterType, 3, Eigen::AffineCompact> trafo_half_inverse;
+- Eigen::Vector3 centre;
++ Eigen::Vector3d centre;
+ Eigen::VectorXd optimiser_weights;
+ bool nonsymmetric;
+
+Index: mrtrix3/src/registration/transform/initialiser.cpp
+===================================================================
+--- mrtrix3.orig/src/registration/transform/initialiser.cpp
++++ mrtrix3/src/registration/transform/initialiser.cpp
+@@ -36,8 +36,8 @@ namespace MR
+ const vector<MultiContrastSetting>& contrast_settings) {
+
+ CONSOLE ("initialising centre of rotation using centre of mass");
+- Eigen::Vector3 im1_centre_mass, im2_centre_mass;
+- Eigen::Vector3 im1_centre_mass_transformed, im2_centre_mass_transformed;
++ Eigen::Vector3d im1_centre_mass, im2_centre_mass;
++ Eigen::Vector3d im1_centre_mass_transformed, im2_centre_mass_transformed;
+
+ Image<default_type> bogus_mask;
+
+@@ -47,7 +47,7 @@ namespace MR
+ transform.transform_half_inverse (im1_centre_mass_transformed, im1_centre_mass);
+ transform.transform_half (im2_centre_mass_transformed, im2_centre_mass);
+
+- Eigen::Vector3 centre = (im1_centre_mass + im2_centre_mass) * 0.5;
++ Eigen::Vector3d centre = (im1_centre_mass + im2_centre_mass) * 0.5;
+ DEBUG("centre: " + str(centre.transpose()));
+ transform.set_centre_without_transform_update (centre);
+ }
+@@ -61,13 +61,13 @@ namespace MR
+ Registration::Transform::Init::LinearInitialisationParams& init) {
+
+ CONSOLE ("initialising centre of rotation using geometric centre");
+- Eigen::Vector3 im1_centre_scanner;
++ Eigen::Vector3d im1_centre_scanner;
+ get_geometric_centre (im1, im1_centre_scanner);
+
+- Eigen::Vector3 im2_centre_scanner;
++ Eigen::Vector3d im2_centre_scanner;
+ get_geometric_centre (im2, im2_centre_scanner);
+
+- Eigen::Vector3 centre = (im1_centre_scanner + im2_centre_scanner) / 2.0;
++ Eigen::Vector3d centre = (im1_centre_scanner + im2_centre_scanner) / 2.0;
+ DEBUG("centre: " + str(centre.transpose()));
+ transform.set_centre_without_transform_update (centre);
+ }
+@@ -82,14 +82,14 @@ namespace MR
+ Registration::Transform::Init::LinearInitialisationParams& init) {
+
+ CONSOLE ("initialising centre of rotation and translation using geometric centre");
+- Eigen::Vector3 im1_centre_scanner;
++ Eigen::Vector3d im1_centre_scanner;
+ get_geometric_centre (im1, im1_centre_scanner);
+
+- Eigen::Vector3 im2_centre_scanner;
++ Eigen::Vector3d im2_centre_scanner;
+ get_geometric_centre (im2, im2_centre_scanner);
+
+- Eigen::Vector3 translation = im1_centre_scanner - im2_centre_scanner;
+- Eigen::Vector3 centre = (im1_centre_scanner + im2_centre_scanner) / 2.0;
++ Eigen::Vector3d translation = im1_centre_scanner - im2_centre_scanner;
++ Eigen::Vector3d centre = (im1_centre_scanner + im2_centre_scanner) / 2.0;
+ transform.set_centre (centre);
+ transform.set_translation (translation);
+ }
+Index: mrtrix3/src/registration/transform/initialiser_helpers.cpp
+===================================================================
+--- mrtrix3.orig/src/registration/transform/initialiser_helpers.cpp
++++ mrtrix3/src/registration/transform/initialiser_helpers.cpp
+@@ -98,8 +98,8 @@ namespace MR
+ if (!mask.value())
+ return;
+ }
+- Eigen::Vector3 voxel_pos ((default_type)image.index(0), (default_type)image.index(1), (default_type)image.index(2));
+- Eigen::Vector3 scanner_pos = transform.voxel2scanner * voxel_pos;
++ Eigen::Vector3d voxel_pos ((default_type)image.index(0), (default_type)image.index(1), (default_type)image.index(2));
++ Eigen::Vector3d scanner_pos = transform.voxel2scanner * voxel_pos;
+
+ default_type xc = scanner_pos[0] - centre[0];
+ default_type yc = scanner_pos[1] - centre[1];
+@@ -165,8 +165,8 @@ namespace MR
+
+ if (!mask.valid()) {
+ for (auto i = Loop (0, 3)(image); i; ++i) {
+- Eigen::Vector3 voxel_pos ((default_type)image.index(0), (default_type)image.index(1), (default_type)image.index(2));
+- Eigen::Vector3 scanner_pos = transform.voxel2scanner * voxel_pos;
++ Eigen::Vector3d voxel_pos ((default_type)image.index(0), (default_type)image.index(1), (default_type)image.index(2));
++ Eigen::Vector3d scanner_pos = transform.voxel2scanner * voxel_pos;
+
+ default_type xc = scanner_pos[0] - centre[0];
+ default_type yc = scanner_pos[1] - centre[1];
+@@ -191,8 +191,8 @@ namespace MR
+ for (auto i = Loop (0, 3)(image, mask); i; ++i) {
+ if (mask.value() <= 0.0)
+ continue;
+- Eigen::Vector3 voxel_pos ((default_type)image.index(0), (default_type)image.index(1), (default_type)image.index(2));
+- Eigen::Vector3 scanner_pos = transform.voxel2scanner * voxel_pos;
++ Eigen::Vector3d voxel_pos ((default_type)image.index(0), (default_type)image.index(1), (default_type)image.index(2));
++ Eigen::Vector3d scanner_pos = transform.voxel2scanner * voxel_pos;
+ default_type val = image.value();
+
+ default_type xc = scanner_pos[0] - centre[0];
+@@ -272,7 +272,7 @@ namespace MR
+
+ void get_centre_of_mass (Image<default_type>& im,
+ Image<default_type>& mask,
+- Eigen::Vector3& centre_of_mass,
++ Eigen::Vector3d& centre_of_mass,
+ const vector<MultiContrastSetting>& contrast_settings) {
+ centre_of_mass.setZero();
+ default_type mass (0.0);
+@@ -326,8 +326,8 @@ namespace MR
+ get_centre_of_mass (im1, init.init_translation.unmasked1 ? bogus_mask : mask1, im1_centre_of_mass, contrast_settings);
+ get_centre_of_mass (im2, init.init_translation.unmasked2 ? bogus_mask : mask2, im2_centre_of_mass, contrast_settings);
+
+- Eigen::Vector3 centre = (im1_centre_of_mass + im2_centre_of_mass) / 2.0;
+- Eigen::Vector3 translation = im1_centre_of_mass - im2_centre_of_mass;
++ Eigen::Vector3d centre = (im1_centre_of_mass + im2_centre_of_mass) / 2.0;
++ Eigen::Vector3d translation = im1_centre_of_mass - im2_centre_of_mass;
+ transform.set_centre_without_transform_update (centre);
+ transform.set_translation (translation);
+ #ifdef DEBUG_INIT
+@@ -375,8 +375,8 @@ namespace MR
+ // Transform tra1;
+ MR::Transform T1 (im1);
+ MR::Transform T2 (im2);
+- Eigen::Vector3 c1 = T1.scanner2voxel * im1_centre_of_mass;
+- Eigen::Vector3 c2 = T2.scanner2voxel * im2_centre_of_mass;
++ Eigen::Vector3d c1 = T1.scanner2voxel * im1_centre_of_mass;
++ Eigen::Vector3d c2 = T2.scanner2voxel * im2_centre_of_mass;
+ VEC(c1)
+ VEC(c2)
+ im1_moments.index(0) = std::round(c1[0]);
+@@ -405,9 +405,9 @@ namespace MR
+ void MomentsInitialiser::run () {
+ if (!calculate_eigenvectors(im1, im2, mask1, mask2)) {
+ WARN ("Image moments not successful. Using centre of mass instead.");
+- Eigen::Vector3 centre = (im1_centre_of_mass + im2_centre_of_mass) / 2.0;
++ Eigen::Vector3d centre = (im1_centre_of_mass + im2_centre_of_mass) / 2.0;
+ transform.set_centre (centre);
+- Eigen::Vector3 translation = im1_centre_of_mass - im2_centre_of_mass;
++ Eigen::Vector3d translation = im1_centre_of_mass - im2_centre_of_mass;
+ transform.set_translation (translation);
+ return;
+ }
+@@ -466,8 +466,8 @@ namespace MR
+ assert(abs(A.determinant() - 1.0) < 0.0001);
+ A = A.transpose().eval(); // A * im2_evec = im1_evec
+
+- Eigen::Vector3 centre = (im1_centre_of_mass + im2_centre_of_mass) / 2.0;
+- Eigen::Vector3 offset = im1_centre_of_mass - im2_centre_of_mass;
++ Eigen::Vector3d centre = (im1_centre_of_mass + im2_centre_of_mass) / 2.0;
++ Eigen::Vector3d offset = im1_centre_of_mass - im2_centre_of_mass;
+ transform.set_centre_without_transform_update (centre);
+
+ Eigen::Translation<default_type, 3> T_offset (offset), T_c2 (im2_centre_of_mass);
+@@ -530,8 +530,8 @@ namespace MR
+ Eigen::Matrix<default_type, 3, 1>& centre_of_mass) {
+ // centre of mass is calculated using only the zeroth order
+ // TODO multithread
+- Eigen::Vector3 voxel_pos = Eigen::Vector3::Zero();
+- Eigen::Vector3 scanner = Eigen::Vector3::Zero();
++ Eigen::Vector3d voxel_pos = Eigen::Vector3d::Zero();
++ Eigen::Vector3d scanner = Eigen::Vector3d::Zero();
+ MR::Transform im_transform (im);
+ default_type im_mass (0.0);
+
+Index: mrtrix3/src/registration/transform/initialiser_helpers.h
+===================================================================
+--- mrtrix3.orig/src/registration/transform/initialiser_helpers.h
++++ mrtrix3/src/registration/transform/initialiser_helpers.h
+@@ -37,7 +37,7 @@ namespace MR
+ template <class ImageType, class ValueType>
+ void get_geometric_centre (const ImageType& image, Eigen::Matrix<ValueType, 3, 1>& centre)
+ {
+- Eigen::Vector3 centre_voxel;
++ Eigen::Vector3d centre_voxel;
+ centre_voxel[0] = (static_cast<default_type>(image.size(0)) / 2.0) - 1.0;
+ centre_voxel[1] = (static_cast<default_type>(image.size(1)) / 2.0) - 1.0;
+ centre_voxel[2] = (static_cast<default_type>(image.size(2)) / 2.0) - 1.0;
+@@ -47,7 +47,7 @@ namespace MR
+
+ void get_centre_of_mass (Image<default_type>& im,
+ Image<default_type>& mask,
+- Eigen::Vector3& centre_of_mass,
++ Eigen::Vector3d& centre_of_mass,
+ const vector<MultiContrastSetting>& contrast_settings);
+
+ bool get_sorted_eigen_vecs_vals (const Eigen::Matrix<default_type, 3, 3>& mat,
+@@ -132,7 +132,7 @@ namespace MR
+ Image<default_type>& mask1;
+ Image<default_type>& mask2;
+ const vector<MultiContrastSetting>& contrast_settings;
+- Eigen::Vector3 im1_centre, im2_centre;
++ Eigen::Vector3d im1_centre, im2_centre;
+ Eigen::Matrix<default_type, 3, 1> im1_centre_of_mass, im2_centre_of_mass;
+ Eigen::Matrix<default_type, 3, 3> im1_covariance_matrix, im2_covariance_matrix;
+ Eigen::Matrix<default_type, Eigen::Dynamic, Eigen::Dynamic> im1_evec, im2_evec;
+Index: mrtrix3/src/registration/transform/rigid.cpp
+===================================================================
+--- mrtrix3.orig/src/registration/transform/rigid.cpp
++++ mrtrix3/src/registration/transform/rigid.cpp
+@@ -234,14 +234,14 @@ namespace MR
+ */
+
+
+- Eigen::Matrix<default_type, 4, 1> Rigid::get_jacobian_vector_wrt_params (const Eigen::Vector3& p) const {
++ Eigen::Matrix<default_type, 4, 1> Rigid::get_jacobian_vector_wrt_params (const Eigen::Vector3d& p) const {
+ Eigen::Matrix<default_type, 4, 1> jac;
+ jac.head(3) = p - centre;
+ jac(3) = 1.0;
+ return jac;
+ }
+
+- Eigen::MatrixXd Rigid::get_jacobian_wrt_params (const Eigen::Vector3& p) const {
++ Eigen::MatrixXd Rigid::get_jacobian_wrt_params (const Eigen::Vector3d& p) const {
+ Eigen::MatrixXd jacobian (3, 12);
+ jacobian.setZero();
+ const auto v = get_jacobian_vector_wrt_params(p);
+Index: mrtrix3/src/registration/transform/rigid.h
+===================================================================
+--- mrtrix3.orig/src/registration/transform/rigid.h
++++ mrtrix3/src/registration/transform/rigid.h
+@@ -97,9 +97,9 @@ namespace MR
+ this->optimiser_weights << weights, weights, weights;
+ }
+
+- Eigen::Matrix<default_type, 4, 1> get_jacobian_vector_wrt_params (const Eigen::Vector3& p) const ;
++ Eigen::Matrix<default_type, 4, 1> get_jacobian_vector_wrt_params (const Eigen::Vector3d& p) const ;
+
+- Eigen::MatrixXd get_jacobian_wrt_params (const Eigen::Vector3& p) const ;
++ Eigen::MatrixXd get_jacobian_wrt_params (const Eigen::Vector3d& p) const ;
+
+ void set_parameter_vector (const Eigen::Matrix<ParameterType, Eigen::Dynamic, 1>& param_vector);
+
+Index: mrtrix3/src/registration/warp/compose.h
+===================================================================
+--- mrtrix3.orig/src/registration/warp/compose.h
++++ mrtrix3/src/registration/warp/compose.h
+@@ -40,7 +40,7 @@ namespace MR
+
+ template <class InputDeformationFieldType, class OutputDeformationFieldType>
+ void operator() (InputDeformationFieldType& deform_input, OutputDeformationFieldType& deform_output) {
+- deform_output.row(3) = transform * Eigen::Vector3 (deform_input.row(3));
++ deform_output.row(3) = transform * Eigen::Vector3d (deform_input.row(3));
+ }
+
+ protected:
+@@ -58,8 +58,8 @@ namespace MR
+
+ template <class DisplacementFieldType, class DeformationFieldType>
+ void operator() (DisplacementFieldType& disp_input, DeformationFieldType& deform_output) {
+- Eigen::Vector3 voxel (disp_input.index(0), disp_input.index(1), disp_input.index(2));
+- deform_output.row(3) = linear_transform * (image_transform.voxel2scanner * voxel + Eigen::Vector3 (disp_input.row(3)));
++ Eigen::Vector3d voxel (disp_input.index(0), disp_input.index(1), disp_input.index(2));
++ deform_output.row(3) = linear_transform * (image_transform.voxel2scanner * voxel + Eigen::Vector3d (disp_input.row(3)));
+ }
+
+ protected:
+@@ -74,15 +74,15 @@ namespace MR
+
+
+ void operator() (Image<default_type>& disp_input1, Image<default_type>& disp_output) {
+- Eigen::Vector3 voxel ((default_type)disp_input1.index(0), (default_type)disp_input1.index(1), (default_type)disp_input1.index(2));
+- Eigen::Vector3 voxel_position = disp1_transform.voxel2scanner * voxel;
+- Eigen::Vector3 original_position = voxel_position + Eigen::Vector3(disp_input1.row(3));
++ Eigen::Vector3d voxel ((default_type)disp_input1.index(0), (default_type)disp_input1.index(1), (default_type)disp_input1.index(2));
++ Eigen::Vector3d voxel_position = disp1_transform.voxel2scanner * voxel;
++ Eigen::Vector3d original_position = voxel_position + Eigen::Vector3d(disp_input1.row(3));
+ disp2_interp.scanner (original_position);
+ if (!disp2_interp) {
+ disp_output.row(3) = disp_input1.row(3);
+ } else {
+- Eigen::Vector3 displacement (Eigen::Vector3(disp2_interp.row(3)).array() * step);
+- Eigen::Vector3 new_position = displacement + original_position;
++ Eigen::Vector3d displacement (Eigen::Vector3d(disp2_interp.row(3)).array() * step);
++ Eigen::Vector3d new_position = displacement + original_position;
+ disp_output.row(3) = new_position - voxel_position;
+ }
+ }
+@@ -106,18 +106,18 @@ namespace MR
+
+
+ void operator() (Image<default_type>& deform) {
+- Eigen::Vector3 voxel ((default_type)deform.index(0), (default_type)deform.index(1), (default_type)deform.index(2));
+- Eigen::Vector3 position = linear1 * voxel;
++ Eigen::Vector3d voxel ((default_type)deform.index(0), (default_type)deform.index(1), (default_type)deform.index(2));
++ Eigen::Vector3d position = linear1 * voxel;
+ deform1_interp.scanner (position);
+ if (!deform1_interp) {
+ deform.row(3) = out_of_bounds;
+ } else {
+- Eigen::Vector3 position2 = deform1_interp.row(3);
++ Eigen::Vector3d position2 = deform1_interp.row(3);
+ deform2_interp.scanner (position2);
+ if (!deform2_interp) {
+ deform.row(3) = out_of_bounds;
+ } else {
+- Eigen::Vector3 position3 = deform2_interp.row(3);
++ Eigen::Vector3d position3 = deform2_interp.row(3);
+ deform.row(3) = linear2 * position3;
+ }
+ }
+@@ -128,7 +128,7 @@ namespace MR
+ Interp::Linear<DeformationField2Type> deform1_interp;
+ Interp::Linear<DeformationField2Type> deform2_interp;
+ const transform_type linear2;
+- Eigen::Vector3 out_of_bounds;
++ Eigen::Vector3d out_of_bounds;
+ };
+
+
+@@ -162,7 +162,7 @@ namespace MR
+
+ default_type max_norm = 0.0;
+ auto max_norm_func = [&max_norm](Image<default_type>& update) {
+- default_type norm = Eigen::Vector3 (update.row(3)).norm();
++ default_type norm = Eigen::Vector3d (update.row(3)).norm();
+ if (norm > max_norm)
+ max_norm = norm;
+ };
+@@ -184,7 +184,7 @@ namespace MR
+ default_type scaled_step = step / scale_factor; // apply the step size and scale factor at once
+ ThreadedLoop (update).run (
+ [&scaled_step](Image<default_type>& update, Image<default_type>& scaled_update) {
+- scaled_update.row(3) = Eigen::Vector3 (update.row(3)) * scaled_step;
++ scaled_update.row(3) = Eigen::Vector3d (update.row(3)) * scaled_step;
+ }, update, *scaled_update);
+
+ // CONSOLE ("composing " + str(std::log2 (scale_factor)) + "times");
+Index: mrtrix3/src/registration/warp/convert.h
+===================================================================
+--- mrtrix3.orig/src/registration/warp/convert.h
++++ mrtrix3/src/registration/warp/convert.h
+@@ -34,8 +34,8 @@ namespace MR
+ void displacement2deformation (ImageType& input, ImageType& output) {
+ MR::Transform transform (input);
+ auto kernel = [&] (ImageType& input, ImageType& output) {
+- Eigen::Vector3 voxel ((default_type)input.index(0), (default_type)input.index(1), (default_type)input.index(2));
+- output.row(3) = (transform.voxel2scanner * voxel).template cast<typename ImageType::value_type> () + Eigen::Vector3 (input.row(3));
++ Eigen::Vector3d voxel ((default_type)input.index(0), (default_type)input.index(1), (default_type)input.index(2));
++ output.row(3) = (transform.voxel2scanner * voxel).template cast<typename ImageType::value_type> () + Eigen::Vector3d (input.row(3));
+ };
+ ThreadedLoop (input, 0, 3).run (kernel, input, output);
+ }
+@@ -44,8 +44,8 @@ namespace MR
+ void deformation2displacement (ImageType& input, ImageType& output) {
+ MR::Transform transform (input);
+ auto kernel = [&] (ImageType& input, ImageType& output) {
+- Eigen::Vector3 voxel ((default_type)input.index(0), (default_type)input.index(1), (default_type)input.index(2));
+- output.row(3) = Eigen::Vector3(input.row(3)) - transform.voxel2scanner * voxel;
++ Eigen::Vector3d voxel ((default_type)input.index(0), (default_type)input.index(1), (default_type)input.index(2));
++ output.row(3) = Eigen::Vector3d(input.row(3)) - transform.voxel2scanner * voxel;
+ };
+ ThreadedLoop (input, 0, 3).run (kernel, input, output);
+ }
+Index: mrtrix3/src/registration/warp/invert.h
+===================================================================
+--- mrtrix3.orig/src/registration/warp/invert.h
++++ mrtrix3/src/registration/warp/invert.h
+@@ -47,9 +47,9 @@ namespace MR
+
+ void operator() (Image<default_type>& displacement_inverse)
+ {
+- Eigen::Vector3 voxel ((default_type)displacement_inverse.index(0), (default_type)displacement_inverse.index(1), (default_type)displacement_inverse.index(2));
+- Eigen::Vector3 truth = transform.voxel2scanner * voxel;
+- Eigen::Vector3 current = truth + Eigen::Vector3(displacement_inverse.row(3));
++ Eigen::Vector3d voxel ((default_type)displacement_inverse.index(0), (default_type)displacement_inverse.index(1), (default_type)displacement_inverse.index(2));
++ Eigen::Vector3d truth = transform.voxel2scanner * voxel;
++ Eigen::Vector3d current = truth + Eigen::Vector3d(displacement_inverse.row(3));
+
+ size_t iter = 0;
+ default_type error = std::numeric_limits<default_type>::max();
+@@ -62,10 +62,10 @@ namespace MR
+
+ private:
+
+- default_type update (Eigen::Vector3& current, const Eigen::Vector3& truth)
++ default_type update (Eigen::Vector3d& current, const Eigen::Vector3d& truth)
+ {
+ displacement.scanner (current);
+- Eigen::Vector3 discrepancy = truth - (current + Eigen::Vector3 (displacement.row(3)));
++ Eigen::Vector3d discrepancy = truth - (current + Eigen::Vector3d (displacement.row(3)));
+ current += discrepancy;
+ return discrepancy.dot (discrepancy);
+ }
+@@ -91,9 +91,9 @@ namespace MR
+
+ void operator() (Image<default_type>& inv_deform)
+ {
+- Eigen::Vector3 voxel ((default_type)inv_deform.index(0), (default_type)inv_deform.index(1), (default_type)inv_deform.index(2));
+- Eigen::Vector3 truth = transform.voxel2scanner * voxel;
+- Eigen::Vector3 current = inv_deform.row(3);
++ Eigen::Vector3d voxel ((default_type)inv_deform.index(0), (default_type)inv_deform.index(1), (default_type)inv_deform.index(2));
++ Eigen::Vector3d truth = transform.voxel2scanner * voxel;
++ Eigen::Vector3d current = inv_deform.row(3);
+
+ size_t iter = 0;
+ default_type error = std::numeric_limits<default_type>::max();
+@@ -106,10 +106,10 @@ namespace MR
+
+ private:
+
+- default_type update (Eigen::Vector3& current, const Eigen::Vector3& truth)
++ default_type update (Eigen::Vector3d& current, const Eigen::Vector3d& truth)
+ {
+ deform.scanner (current);
+- Eigen::Vector3 discrepancy = truth - Eigen::Vector3 (deform.row(3));
++ Eigen::Vector3d discrepancy = truth - Eigen::Vector3d (deform.row(3));
+ current += discrepancy;
+ return discrepancy.dot (discrepancy);
+ }
+Index: mrtrix3/src/surface/algo/image2mesh.h
+===================================================================
+--- mrtrix3.orig/src/surface/algo/image2mesh.h
++++ mrtrix3/src/surface/algo/image2mesh.h
+@@ -122,7 +122,7 @@ namespace MR
+ const auto existing = vox2vertindex.find (voxels[in_vertex]);
+ if (existing == vox2vertindex.end()) {
+ triangle_vertices[out_vertex] = vertices.size();
+- Eigen::Vector3 pos_voxelspace (default_type(voxels[in_vertex][0]) - 0.5, default_type(voxels[in_vertex][1]) - 0.5, default_type(voxels[in_vertex][2]) - 0.5);
++ Eigen::Vector3d pos_voxelspace (default_type(voxels[in_vertex][0]) - 0.5, default_type(voxels[in_vertex][1]) - 0.5, default_type(voxels[in_vertex][2]) - 0.5);
+ vertices.push_back (transform.voxel2scanner * pos_voxelspace);
+ } else {
+ triangle_vertices[out_vertex] = existing->second;
+Index: mrtrix3/src/surface/algo/mesh2image.cpp
+===================================================================
+--- mrtrix3.orig/src/surface/algo/mesh2image.cpp
++++ mrtrix3/src/surface/algo/mesh2image.cpp
+@@ -47,7 +47,7 @@ namespace MR
+
+ // For speed, want the vertex data to be in voxel positions
+ Mesh mesh;
+- vector<Eigen::Vector3> polygon_normals;
++ vector<Eigen::Vector3d> polygon_normals;
+
+ // For every edge voxel, stores those polygons that may intersect the voxel
+ using Vox2Poly = std::map< Vox, vector<size_t> >;
+@@ -126,13 +126,13 @@ namespace MR
+ mesh.load_quad_vertices (vertices, poly_index - mesh.num_triangles());
+
+ // Test whether or not the two objects can be separated via projection onto an axis
+- auto separating_axis = [&] (const Eigen::Vector3& axis) -> bool {
++ auto separating_axis = [&] (const Eigen::Vector3d& axis) -> bool {
+ default_type voxel_low = std::numeric_limits<default_type>::infinity();
+ default_type voxel_high = -std::numeric_limits<default_type>::infinity();
+ default_type poly_low = std::numeric_limits<default_type>::infinity();
+ default_type poly_high = -std::numeric_limits<default_type>::infinity();
+
+- static const Eigen::Vector3 voxel_offsets[8] = { { -0.5, -0.5, -0.5 },
++ static const Eigen::Vector3d voxel_offsets[8] = { { -0.5, -0.5, -0.5 },
+ { -0.5, -0.5, 0.5 },
+ { -0.5, 0.5, -0.5 },
+ { -0.5, 0.5, 0.5 },
+@@ -142,7 +142,7 @@ namespace MR
+ { 0.5, 0.5, 0.5 } };
+
+ for (size_t i = 0; i != 8; ++i) {
+- const Eigen::Vector3 v (vox.matrix().cast<default_type>() + voxel_offsets[i]);
++ const Eigen::Vector3d v (vox.matrix().cast<default_type>() + voxel_offsets[i]);
+ const default_type projection = axis.dot (v);
+ voxel_low = std::min (voxel_low, projection);
+ voxel_high = std::max (voxel_high, projection);
+@@ -163,7 +163,7 @@ namespace MR
+ // All cross-products between voxel and polygon edges
+ // Polygon normal
+ for (size_t i = 0; i != 3; ++i) {
+- Eigen::Vector3 axis (0.0, 0.0, 0.0);
++ Eigen::Vector3d axis (0.0, 0.0, 0.0);
+ axis[i] = 1.0;
+ if (separating_axis (axis))
+ return false;
+@@ -228,7 +228,7 @@ namespace MR
+ for (adj_voxel[1] = centre_voxel[1]-1; adj_voxel[1] <= centre_voxel[1]+1; ++adj_voxel[1]) {
+ for (adj_voxel[0] = centre_voxel[0]-1; adj_voxel[0] <= centre_voxel[0]+1; ++adj_voxel[0]) {
+ if (!is_out_of_bounds (H, adj_voxel) && (adj_voxel - centre_voxel).any()) {
+- const Eigen::Vector3 offset (adj_voxel.cast<default_type>().matrix() - mesh.vert(i));
++ const Eigen::Vector3d offset (adj_voxel.cast<default_type>().matrix() - mesh.vert(i));
+ const default_type dp_normal = offset.dot (mesh.norm(i));
+ const default_type offset_on_plane = (offset - (mesh.norm(i) * dp_normal)).norm();
+ assign_pos_of (adj_voxel).to (sum_distances);
+@@ -376,13 +376,13 @@ namespace MR
+ class Pipe
+ { NOMEMALIGN
+ public:
+- Pipe (const Mesh& mesh, const vector<Eigen::Vector3>& polygon_normals) :
++ Pipe (const Mesh& mesh, const vector<Eigen::Vector3d>& polygon_normals) :
+ mesh (mesh),
+ polygon_normals (polygon_normals)
+
+ {
+ // Generate a set of points within this voxel that need to be tested individually
+- offsets_to_test.reset(new vector<Eigen::Vector3>());
++ offsets_to_test.reset(new vector<Eigen::Vector3d>());
+ offsets_to_test->reserve (pve_nsamples);
+ for (size_t x_idx = 0; x_idx != pve_os_ratio; ++x_idx) {
+ const default_type x = -0.5 + ((default_type(x_idx) + 0.5) / default_type(pve_os_ratio));
+@@ -404,7 +404,7 @@ namespace MR
+ size_t inside_mesh_count = 0;
+ for (vector<Vertex>::const_iterator i_p = offsets_to_test->begin(); i_p != offsets_to_test->end(); ++i_p) {
+ Vertex p (*i_p);
+- p += Eigen::Vector3 (voxel[0], voxel[1], voxel[2]);
++ p += Eigen::Vector3d (voxel[0], voxel[1], voxel[2]);
+
+ default_type best_min_edge_distance_on_plane = -std::numeric_limits<default_type>::infinity();
+ bool best_result_inside = false;
+@@ -412,7 +412,7 @@ namespace MR
+
+ // Only test against those polygons that are near this voxel
+ for (vector<size_t>::const_iterator polygon_index = in.second.begin(); polygon_index != in.second.end(); ++polygon_index) {
+- const Eigen::Vector3& n (polygon_normals[*polygon_index]);
++ const Eigen::Vector3d& n (polygon_normals[*polygon_index]);
+
+ const size_t polygon_num_vertices = (*polygon_index < mesh.num_triangles()) ? 3 : 4;
+ VertexList v;
+@@ -516,9 +516,9 @@ namespace MR
+
+ private:
+ const Mesh& mesh;
+- const vector<Eigen::Vector3>& polygon_normals;
++ const vector<Eigen::Vector3d>& polygon_normals;
+
+- std::shared_ptr<vector<Eigen::Vector3>> offsets_to_test;
++ std::shared_ptr<vector<Eigen::Vector3d>> offsets_to_test;
+
+ };
+
+Index: mrtrix3/src/surface/filter/vertex_transform.cpp
+===================================================================
+--- mrtrix3.orig/src/surface/filter/vertex_transform.cpp
++++ mrtrix3/src/surface/filter/vertex_transform.cpp
+@@ -93,7 +93,7 @@ namespace MR
+ case transform_t::FS2REAL:
+ vector<size_t> axes( 3 );
+ auto M = File::NIfTI::adjust_transform( header, axes );
+- Eigen::Vector3 cras( 3, 1 );
++ Eigen::Vector3d cras( 3, 1 );
+ for ( size_t i = 0; i < 3; i++ )
+ {
+ cras[ i ] = M( i, 3 );
+Index: mrtrix3/src/surface/mesh.cpp
+===================================================================
+--- mrtrix3.orig/src/surface/mesh.cpp
++++ mrtrix3/src/surface/mesh.cpp
+@@ -293,7 +293,7 @@ namespace MR
+ warn_attribute = true;
+
+ triangles.push_back ( vector<uint32_t> { uint32_t(vertices.size()-3), uint32_t(vertices.size()-2), uint32_t(vertices.size()-1) } );
+- const Eigen::Vector3 computed_normal = Surface::normal (*this, triangles.back());
++ const Eigen::Vector3d computed_normal = Surface::normal (*this, triangles.back());
+ if (computed_normal.dot (normal.cast<default_type>()) < 0.0)
+ warn_right_hand_rule = true;
+ if (abs (computed_normal.dot (normal.cast<default_type>())) < 0.99)
+@@ -357,7 +357,7 @@ namespace MR
+ throw Exception ("Error parsing STL file " + Path::basename (path) + ": facet ended with " + str(vertex_index) + " vertices");
+ triangles.push_back ( vector<uint32_t> { uint32_t(vertices.size()-3), uint32_t(vertices.size()-2), uint32_t(vertices.size()-1) } );
+ vertex_index = 0;
+- const Eigen::Vector3 computed_normal = Surface::normal (*this, triangles.back());
++ const Eigen::Vector3d computed_normal = Surface::normal (*this, triangles.back());
+ if (computed_normal.dot (normal) < 0.0)
+ warn_right_hand_rule = true;
+ if (abs (computed_normal.dot (normal)) < 0.99)
+@@ -720,7 +720,7 @@ namespace MR
+ out.write (reinterpret_cast<const char*>(&count), sizeof(uint32_t));
+ const uint16_t attribute_byte_count = 0;
+ for (TriangleList::const_iterator i = triangles.begin(); i != triangles.end(); ++i) {
+- const Eigen::Vector3 n (normal (*this, *i));
++ const Eigen::Vector3d n (normal (*this, *i));
+ const float n_temp[3] { float(n[0]), float(n[1]), float(n[2]) };
+ out.write (reinterpret_cast<const char*>(&n_temp[0]), 3 * sizeof(float));
+ for (size_t v = 0; v != 3; ++v) {
+@@ -737,7 +737,7 @@ namespace MR
+ File::OFStream out (path);
+ out << "solid \n";
+ for (TriangleList::const_iterator i = triangles.begin(); i != triangles.end(); ++i) {
+- const Eigen::Vector3 n (normal (*this, *i));
++ const Eigen::Vector3d n (normal (*this, *i));
+ out << "facet normal " << str (n[0]) << " " << str (n[1]) << " " << str (n[2]) << "\n";
+ out << " outer loop\n";
+ for (size_t v = 0; v != 3; ++v) {
+Index: mrtrix3/src/surface/types.h
+===================================================================
+--- mrtrix3.orig/src/surface/types.h
++++ mrtrix3/src/surface/types.h
+@@ -29,7 +29,7 @@ namespace MR
+ {
+
+
+- using Vertex = Eigen::Vector3;
++ using Vertex = Eigen::Vector3d;
+ using VertexList = vector<Vertex>;
+ using Triangle = Polygon<3>;
+ using TriangleList = vector<Triangle>;
+@@ -41,7 +41,7 @@ namespace MR
+ public:
+ using Eigen::Array3i::Array3i;
+ Vox () : Eigen::Array3i (-1, -1, -1) { }
+- Vox (const Eigen::Vector3& p) : Eigen::Array3i (int(std::round (p[0])), int(std::round (p[1])), int(std::round (p[2]))) { }
++ Vox (const Eigen::Vector3d& p) : Eigen::Array3i (int(std::round (p[0])), int(std::round (p[1])), int(std::round (p[2]))) { }
+ bool operator< (const Vox& i) const
+ {
+ return ((*this)[2] == i[2] ? (((*this)[1] == i[1]) ? ((*this)[0] < i[0]) : ((*this)[1] < i[1])) : ((*this)[2] < i[2]));
+Index: mrtrix3/testing/cmd/testing_diff_mesh.cpp
+===================================================================
+--- mrtrix3.orig/testing/cmd/testing_diff_mesh.cpp
++++ mrtrix3/testing/cmd/testing_diff_mesh.cpp
+@@ -84,12 +84,12 @@ void run ()
+
+ for (size_t i = 0; i != in1.num_triangles(); ++i) {
+ // Explicitly load the vertex data
+- std::array<Eigen::Vector3, 3> v1;
++ std::array<Eigen::Vector3d, 3> v1;
+ for (size_t vertex = 0; vertex != 3; ++vertex)
+ v1[vertex] = in1.vert(in1.tri(i)[vertex]);
+ bool match_found = false;
+ for (size_t j = 0; j != in2.num_triangles() && !match_found; ++j) {
+- std::array<Eigen::Vector3, 3> v2;
++ std::array<Eigen::Vector3d, 3> v2;
+ for (size_t vertex = 0; vertex != 3; ++vertex)
+ v2[vertex] = in2.vert (in2.tri(j)[vertex]);
+ bool all_vertices_matched = true;
+@@ -111,12 +111,12 @@ void run ()
+ }
+
+ for (size_t i = 0; i != in1.num_quads(); ++i) {
+- std::array<Eigen::Vector3, 4> v1;
++ std::array<Eigen::Vector3d, 4> v1;
+ for (size_t vertex = 0; vertex != 4; ++vertex)
+ v1[vertex] = in1.vert (in1.quad(i)[vertex]);
+ bool match_found = false;
+ for (size_t j = 0; j != in2.num_quads() && !match_found; ++j) {
+- std::array<Eigen::Vector3, 4> v2;
++ std::array<Eigen::Vector3d, 4> v2;
+ for (size_t vertex = 0; vertex != 4; ++vertex)
+ v2[vertex] = in2.vert (in2.quad(j)[vertex]);
+ bool all_vertices_matched = true;
+Index: mrtrix3/testing/cmd/testing_diff_peaks.cpp
+===================================================================
+--- mrtrix3.orig/testing/cmd/testing_diff_peaks.cpp
++++ mrtrix3/testing/cmd/testing_diff_peaks.cpp
+@@ -68,7 +68,7 @@ void run ()
+ .run ([&tol] (decltype(in1)& a, decltype(in2)& b)
+ {
+ for (size_t i = 0; i != size_t(a.size(3)); i += 3) {
+- Eigen::Vector3 veca, vecb;
++ Eigen::Vector3d veca, vecb;
+ for (size_t axis = 0; axis != 3; ++axis) {
+ a.index(3) = b.index(3) = i + axis;
+ veca[axis] = a.value();
+Index: mrtrix3/core/eigen_plugins/dense_base.h
+===================================================================
+--- mrtrix3.orig/core/eigen_plugins/dense_base.h
++++ mrtrix3/core/eigen_plugins/dense_base.h
+@@ -19,7 +19,7 @@ template <class ImageType> \
+ inline Derived& operator ARG (const MR::Helper::ConstRow<ImageType>& row) { \
+ this->resize (row.image.size(row.axis),1); \
+ for (row.image.index(row.axis) = 0; row.image.index (row.axis) < row.image.size (row.axis); ++row.image.index (row.axis)) \
+- this->operator() (row.image.index (row.axis), 0) ARG row.image.value(); \
++ this->operator() (ssize_t (row.image.index (row.axis)), 0) ARG row.image.value(); \
+ return derived(); \
+ }
+
+Index: mrtrix3/core/eigen_plugins/matrix.h
+===================================================================
+--- mrtrix3.orig/core/eigen_plugins/matrix.h
++++ mrtrix3/core/eigen_plugins/matrix.h
+@@ -21,7 +21,7 @@ template <class ImageType> Matrix (const
+ template <class ImageType> inline Matrix& operator ARG (const MR::Helper::ConstRow<ImageType>& row) { \
+ this->resize (row.image.size(row.axis),1); \
+ for (row.image.index(row.axis) = 0; row.image.index (row.axis) < row.image.size (row.axis); ++row.image.index (row.axis)) \
+- this->operator() (row.image.index (row.axis), 0) ARG row.image.value(); \
++ this->operator() (ssize_t (row.image.index (row.axis)), 0) ARG row.image.value(); \
+ return *this; \
+ }
+
=====================================
debian/patches/series
=====================================
@@ -3,3 +3,4 @@
debian_build_cfg
python3.patch
mrtrix_version.patch
+eigen-3.4.patch
View it on GitLab: https://salsa.debian.org/med-team/mrtrix3/-/commit/7e63f69ccc938aa3482a50006f9529c7f86d3aea
--
View it on GitLab: https://salsa.debian.org/med-team/mrtrix3/-/commit/7e63f69ccc938aa3482a50006f9529c7f86d3aea
You're receiving this email because of your account on salsa.debian.org.
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://alioth-lists.debian.net/pipermail/debian-med-commit/attachments/20211222/d9eace09/attachment-0001.htm>
More information about the debian-med-commit
mailing list