Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

[FIX] Make changes for newer version of Eigen3 #40

Open
wants to merge 2 commits into
base: master
Choose a base branch
from
Open

[FIX] Make changes for newer version of Eigen3 #40

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
edf8c7f
update for eigen3 update
mattcieslak Apr 5, 2024
5f4c31b
Vector3
mattcieslak Apr 5, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
8 changes: 4 additions & 4 deletions cmd/amp2response.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -83,7 +83,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,
Expand All @@ -92,11 +92,11 @@ Eigen::Matrix<default_type, 3, 3> gen_rotation_matrix (const Eigen::Vector3& dir
// 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;
Expand Down Expand Up @@ -156,7 +156,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();
Expand Down
2 changes: 1 addition & 1 deletion cmd/fixel2tsf.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -95,7 +95,7 @@ void run ()
DWI::Tractography::Streamline<float> tck;

Transform transform (in_index_image);
Eigen::Vector3 voxel_pos;
Eigen::Vector3d voxel_pos;

while (reader (tck)) {
SetVoxelDir dixels;
Expand Down
12 changes: 6 additions & 6 deletions cmd/fixel2voxel.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -395,16 +395,16 @@ class DEC_unit : protected Base

void operator() (Image<uint32_t>& 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())
Expand All @@ -427,13 +427,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();
}
Expand All @@ -444,7 +444,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();
}
}
Expand Down
4 changes: 2 additions & 2 deletions cmd/fixelcfestats.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -202,7 +202,7 @@ void run() {
fixel2row[f] = f;
}

vector<Eigen::Vector3> positions (num_fixels);
vector<Eigen::Vector3d> positions (num_fixels);
vector<direction_type> directions (num_fixels);

const std::string output_fixel_directory = argument[5];
Expand All @@ -213,7 +213,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) {
Eigen::Vector3 vox ((default_type)index_image.index(0), (default_type)index_image.index(1), (default_type)index_image.index(2));
Eigen::Vector3d vox ((default_type)index_image.index(0), (default_type)index_image.index(1), (default_type)index_image.index(2));
vox = image_transform.voxel2scanner * vox;
index_image.index(3) = 1;
const index_type offset = index_image.value();
Expand Down
2 changes: 1 addition & 1 deletion cmd/labelconvert.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -149,7 +149,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;
Expand Down
4 changes: 2 additions & 2 deletions cmd/mrcolour.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -86,7 +86,7 @@ void run ()
{
Header H_in = Header::open (argument[0]);
const GUI::MRView::ColourMap::Entry colourmap = GUI::MRView::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");
Expand All @@ -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) {
Expand Down
4 changes: 2 additions & 2 deletions cmd/mrconvert.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -274,8 +274,8 @@ void permute_slice_direction (Header& H, const vector<int>& axes)
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);
}

Expand Down
6 changes: 3 additions & 3 deletions cmd/mrdegibbs.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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:
Expand Down Expand Up @@ -327,7 +327,7 @@ void run ()
auto slice_encoding_it = header.keyval().find ("SliceEncodingDirection");
if (slice_encoding_it != header.keyval().end()) {
try {
const Eigen::Vector3 slice_endocing_axis_onehot = Axes::id2dir (slice_encoding_it->second);
const Eigen::Vector3d slice_endocing_axis_onehot = Axes::id2dir (slice_encoding_it->second);
vector<size_t> auto_slice_axes = { 0, 0 };
if (slice_endocing_axis_onehot[0])
auto_slice_axes = { 1, 2 };
Expand Down
10 changes: 5 additions & 5 deletions cmd/mredit.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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]));
Expand Down Expand Up @@ -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
Expand All @@ -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)) {
Expand All @@ -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);
Expand Down
6 changes: 3 additions & 3 deletions cmd/mrtransform.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -453,7 +453,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);
Expand Down Expand Up @@ -488,13 +488,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;
Expand Down
16 changes: 8 additions & 8 deletions cmd/mtnormalise.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -99,7 +99,7 @@ struct PolyBasisFunction { MEMALIGN (PolyBasisFunction)

const int n_basis_vecs = 20;

FORCE_INLINE Eigen::MatrixXd operator () (const Eigen::Vector3& pos) {
FORCE_INLINE Eigen::MatrixXd operator () (const Eigen::Vector3d& pos) {
double x = pos[0];
double y = pos[1];
double z = pos[2];
Expand Down Expand Up @@ -134,7 +134,7 @@ struct PolyBasisFunction<0> { MEMALIGN (PolyBasisFunction<0>)
const int n_basis_vecs = 1;

PolyBasisFunction() {}
FORCE_INLINE Eigen::MatrixXd operator () (const Eigen::Vector3&) {
FORCE_INLINE Eigen::MatrixXd operator () (const Eigen::Vector3d&) {
Eigen::MatrixXd basis(n_basis_vecs, 1);
basis(0) = 1.0;
return basis;
Expand All @@ -146,7 +146,7 @@ template <>
struct PolyBasisFunction<1> { MEMALIGN (PolyBasisFunction<1>)
const int n_basis_vecs = 4;

FORCE_INLINE Eigen::MatrixXd operator () (const Eigen::Vector3& pos) {
FORCE_INLINE Eigen::MatrixXd operator () (const Eigen::Vector3d& pos) {
double x = pos[0];
double y = pos[1];
double z = pos[2];
Expand All @@ -164,7 +164,7 @@ template <>
struct PolyBasisFunction<2> { MEMALIGN (PolyBasisFunction<2>)
const int n_basis_vecs = 10;

FORCE_INLINE Eigen::MatrixXd operator () (const Eigen::Vector3& pos) {
FORCE_INLINE Eigen::MatrixXd operator () (const Eigen::Vector3d& pos) {
double x = pos[0];
double y = pos[1];
double z = pos[2];
Expand Down Expand Up @@ -466,8 +466,8 @@ void run_primitive () {
uint32_t index = 0;
for (auto i = Loop (0, 3) (mask, combined_tissue); i; ++i) {
if (mask.value()) {
Eigen::Vector3 vox (mask.index(0), mask.index(1), mask.index(2));
Eigen::Vector3 pos = transform.voxel2scanner * vox;
Eigen::Vector3d vox (mask.index(0), mask.index(1), mask.index(2));
Eigen::Vector3d pos = transform.voxel2scanner * vox;
norm_field_basis.row (index) = basis_function (pos).col(0);

double sum = 0.0;
Expand All @@ -483,8 +483,8 @@ void run_primitive () {

// Generate normalisation field in the log domain
for (auto i = Loop (0, 3) (norm_field_log); i; ++i) {
Eigen::Vector3 vox (norm_field_log.index(0), norm_field_log.index(1), norm_field_log.index(2));
Eigen::Vector3 pos = transform.voxel2scanner * vox;
Eigen::Vector3d vox (norm_field_log.index(0), norm_field_log.index(1), norm_field_log.index(2));
Eigen::Vector3d pos = transform.voxel2scanner * vox;
norm_field_log.value() = basis_function (pos).col(0).dot (norm_field_weights.col(0));
}

Expand Down
12 changes: 6 additions & 6 deletions cmd/tck2fixel.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -41,7 +41,7 @@ class TrackProcessor { MEMALIGN (TrackProcessor)
using SetVoxelDir = DWI::Tractography::Mapping::SetVoxelDir;

TrackProcessor (Image<uint32_t>& 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) ,
Expand All @@ -63,7 +63,7 @@ class TrackProcessor { MEMALIGN (TrackProcessor)
uint32_t last_index = first_index + num_fibres;
uint32_t 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 (uint32_t j = first_index; j < last_index; ++j) {
const float dp = abs (dir.dot (fixel_directions[j]));
if (dp > largest_dp) {
Expand All @@ -83,7 +83,7 @@ class TrackProcessor { MEMALIGN (TrackProcessor)

private:
Image<uint32_t> fixel_indexer;
const vector<Eigen::Vector3>& fixel_directions;
const vector<Eigen::Vector3d>& fixel_directions;
vector<uint16_t>& fixel_TDI;
const float angular_threshold_dp;
};
Expand Down Expand Up @@ -132,8 +132,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);
Expand All @@ -143,7 +143,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;
uint32_t offset = index_image.value();
size_t fixel_index = 0;
Expand Down
2 changes: 1 addition & 1 deletion cmd/tckconvert.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -331,7 +331,7 @@ class PLYWriter: public WriterInterface<float> { MEMALIGN(PLYWriter)
}

// 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);
Expand Down
14 changes: 7 additions & 7 deletions cmd/transformcompose.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -26,16 +26,16 @@ 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;
};


class Warp : public TransformBase { MEMALIGN(Warp)
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
Expand All @@ -52,8 +52,8 @@ class Linear : public TransformBase { MEMALIGN(Linear)
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;
}

Expand Down Expand Up @@ -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);
Expand Down
Loading