Skip to content

Commit

Permalink
Fixes to bvecs import
Browse files Browse the repository at this point in the history
In order for a bvecs file to be properly read in, the information about how the image was transformed upon image load needs to be available. However if an attempt is made to perform such a read based on a Header that it has itself been constructed from an arbitrary ImageType, this information is lost. Therefore, reading from a bvecs file needs to be done using the Header class in which the corresponding NIfTI was originally opened.
In the process, some functionalities of the tracking of that transformation on image load have been made more explicit.
Also fixes some syntax errors relating to the manual back-porting of #2917 in #3011.
  • Loading branch information
Lestropie committed Sep 22, 2024
1 parent 01a09b8 commit 2be257c
Show file tree
Hide file tree
Showing 23 changed files with 224 additions and 192 deletions.
28 changes: 14 additions & 14 deletions cmd/amp2sh.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -201,8 +201,9 @@ class Amp2SH { MEMALIGN(Amp2SH)

void run ()
{
auto amp = Image<value_type>::open (argument[0]).with_direct_io (3);
Header header (amp);
Header header_in (Header::open (argument[0]));
Header header_out (header_in);
header_out.datatype() = DataType::Float32;

vector<size_t> bzeros, dwis;
Eigen::MatrixXd dirs;
Expand All @@ -213,8 +214,8 @@ void run ()
dirs = Math::Sphere::cartesian2spherical (dirs);
}
else {
auto hit = header.keyval().find ("directions");
if (hit != header.keyval().end()) {
auto hit = header_in.keyval().find ("directions");
if (hit != header_in.keyval().end()) {
vector<default_type> dir_vector;
for (auto line : split_lines (hit->second)) {
auto v = parse_floats (line);
Expand All @@ -225,35 +226,34 @@ void run ()
dirs(i/2, 0) = dir_vector[i];
dirs(i/2, 1) = dir_vector[i+1];
}
header.keyval()["basis_directions"] = hit->second;
header.keyval().erase (hit);
header_out.keyval()["basis_directions"] = hit->second;
header_out.keyval().erase (hit);
}
else {
auto grad = DWI::get_DW_scheme (amp);
auto grad = DWI::get_DW_scheme (header_in);
DWI::Shells shells (grad);
shells.select_shells (true, false, false);
if (shells.smallest().is_bzero())
bzeros = shells.smallest().get_volumes();
dwis = shells.largest().get_volumes();
dirs = DWI::gen_direction_matrix (grad, dwis);
DWI::stash_DW_scheme (header, grad);
DWI::stash_DW_scheme (header_out, grad);
}
}
Metadata::PhaseEncoding::clear_scheme (header.keyval());
Metadata::PhaseEncoding::clear_scheme (header_out.keyval());

auto sh2amp = DWI::compute_SH2amp_mapping (dirs, true, 8);


bool normalise = get_options ("normalise").size();
if (normalise && !bzeros.size())
throw Exception ("the normalise option is only available if the input data contains b=0 images.");

header_out.size (3) = sh2amp.cols();
Stride::set_from_command_line (header_out);

header.size (3) = sh2amp.cols();
Stride::set_from_command_line (header);
auto SH = Image<value_type>::create (argument[1], header);

auto amp = header_in.get_image<value_type>().with_direct_io (3);
Amp2SHCommon common (sh2amp, bzeros, dwis, normalise);
auto SH = Image<value_type>::create (argument[1], header_out);

opt = get_options ("rician");
if (opt.size()) {
Expand Down
21 changes: 11 additions & 10 deletions cmd/dwi2adc.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -80,11 +80,11 @@ class DWI2ADC { MEMALIGN(DWI2ADC)


void run () {
auto dwi = Header::open (argument[0]).get_image<value_type>();
auto grad = DWI::get_DW_scheme (dwi);
auto header_in = Header::open (argument[0]);
auto grad = DWI::get_DW_scheme (header_in);

size_t dwi_axis = 3;
while (dwi.size (dwi_axis) < 2)
while (header_in.size (dwi_axis) < 2)
++dwi_axis;
INFO ("assuming DW images are stored along axis " + str (dwi_axis));

Expand All @@ -96,14 +96,15 @@ void run () {

auto binv = Math::pinv (b);

Header header (dwi);
header.datatype() = DataType::Float32;
DWI::stash_DW_scheme (header, grad);
Metadata::PhaseEncoding::clear_scheme (header.keyval());
header.ndim() = 4;
header.size(3) = 2;
Header header_out (header_in);
header_out.datatype() = DataType::Float32;
DWI::stash_DW_scheme (header_out, grad);
Metadata::PhaseEncoding::clear_scheme (header_out.keyval());
header_out.ndim() = 4;
header_out.size(3) = 2;

auto adc = Image<value_type>::create (argument[1], header);
auto dwi = header_in.get_image<value_type>();
auto adc = Image<value_type>::create (argument[1], header_out);

ThreadedLoop ("computing ADC values", dwi, 0, 3)
.run (DWI2ADC (binv, dwi_axis), dwi, adc);
Expand Down
10 changes: 5 additions & 5 deletions cmd/dwi2mask.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -66,13 +66,13 @@ OPTIONS

void run () {

auto input = Image<float>::open (argument[0]).with_direct_io (3);
auto grad = DWI::get_DW_scheme (input);

if (input.ndim() != 4)
auto header = Header::open(argument[0]);
if (header.ndim() != 4)
throw Exception ("input DWI image must be 4D");
auto grad = DWI::get_DW_scheme (header);
auto input = header.get_image<float>().with_direct_io (3);

Filter::DWIBrainMask dwi_brain_mask_filter (input, grad);
Filter::DWIBrainMask dwi_brain_mask_filter (header, grad);
dwi_brain_mask_filter.set_message ("computing dwi brain mask");
auto temp_mask = Image<bool>::scratch (dwi_brain_mask_filter, "brain mask");
dwi_brain_mask_filter (input, temp_mask);
Expand Down
35 changes: 18 additions & 17 deletions cmd/dwi2tensor.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -202,52 +202,53 @@ inline Processor<MASKType, B0Type, DKTType, PredictType> processor (const Eigen:

void run ()
{
auto dwi = Header::open (argument[0]).get_image<value_type>();
auto grad = DWI::get_DW_scheme (dwi);
auto header_in = Header::open (argument[0]);
auto grad = DWI::get_DW_scheme (header_in);

Image<bool> mask;
auto opt = get_options ("mask");
if (opt.size()) {
mask = Image<bool>::open (opt[0][0]);
check_dimensions (dwi, mask, 0, 3);
check_dimensions (header_in, mask, 0, 3);
}

bool ols = get_options ("ols").size();

// depending on whether first (initialisation) loop should be considered an iteration
auto iter = get_option_value ("iter", DEFAULT_NITER);

Header header (dwi);
header.datatype() = DataType::Float32;
header.ndim() = 4;
DWI::stash_DW_scheme (header, grad);
Metadata::PhaseEncoding::clear_scheme (header.keyval());
Header header_out (header_in);
header_out.datatype() = DataType::Float32;
header_out.ndim() = 4;
DWI::stash_DW_scheme (header_out, grad);
Metadata::PhaseEncoding::clear_scheme (header_out.keyval());

Image<value_type> predict;
opt = get_options ("predicted_signal");
if (opt.size())
predict = Image<value_type>::create (opt[0][0], header);
predict = Image<value_type>::create (opt[0][0], header_out);

header.size(3) = 6;
auto dt = Image<value_type>::create (argument[1], header);
header_out.size(3) = 6;
auto dt = Image<value_type>::create (argument[1], header_out);

Image<value_type> b0;
opt = get_options ("b0");
if (opt.size()) {
header.ndim() = 3;
b0 = Image<value_type>::create (opt[0][0], header);
header_out.ndim() = 3;
b0 = Image<value_type>::create (opt[0][0], header_out);
}

Image<value_type> dkt;
opt = get_options ("dkt");
if (opt.size()) {
header.ndim() = 4;
header.size(3) = 15;
dkt = Image<value_type>::create (opt[0][0], header);
header_out.ndim() = 4;
header_out.size(3) = 15;
dkt = Image<value_type>::create (opt[0][0], header_out);
}

Eigen::MatrixXd b = -DWI::grad2bmatrix<double> (grad, opt.size()>0);
Eigen::MatrixXd b = -DWI::grad2bmatrix<double> (grad, dkt.valid());

auto dwi = Header::open (argument[0]).get_image<value_type>();
ThreadedLoop("computing tensors", dwi, 0, 3).run (processor (b, ols, iter, mask, b0, dkt, predict), dwi, dt);
}

23 changes: 12 additions & 11 deletions cmd/dwiextract.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -66,10 +66,10 @@ void usage ()

void run()
{
auto input_image = Image<float>::open (argument[0]);
if (input_image.ndim() < 4)
auto header_in = Header::open (argument[0]);
if (header_in.ndim() < 4)
throw Exception ("Epected input image to contain more than three dimensions");
auto grad = DWI::get_DW_scheme (input_image);
auto grad = DWI::get_DW_scheme (header_in);

// Want to support non-shell-like data if it's just a straight extraction
// of all dwis or all bzeros i.e. don't initialise the Shells class
Expand Down Expand Up @@ -101,7 +101,7 @@ void run()
}

auto opt = get_options ("pe");
const auto pe_scheme = Metadata::PhaseEncoding::get_scheme (input_image);
const auto pe_scheme = Metadata::PhaseEncoding::get_scheme (header_in);
if (opt.size()) {
if (!pe_scheme.rows())
throw Exception ("Cannot filter volumes by phase-encoding: No such information present");
Expand Down Expand Up @@ -134,24 +134,25 @@ void run()

std::sort (volumes.begin(), volumes.end());

Header header (input_image);
Stride::set_from_command_line (header);
header.size (3) = volumes.size();
Header header_out (header_in);
Stride::set_from_command_line (header_out);
header_out.size (3) = volumes.size();

Eigen::MatrixXd new_grad (volumes.size(), grad.cols());
for (size_t i = 0; i < volumes.size(); i++)
new_grad.row (i) = grad.row (volumes[i]);
DWI::set_DW_scheme (header, new_grad);
DWI::set_DW_scheme (header_out, new_grad);

if (pe_scheme.rows()) {
Eigen::MatrixXd new_scheme (volumes.size(), pe_scheme.cols());
for (size_t i = 0; i != volumes.size(); ++i)
new_scheme.row(i) = pe_scheme.row (volumes[i]);
Metadata::PhaseEncoding::set_scheme (header.keyval(), new_scheme);
Metadata::PhaseEncoding::set_scheme (header_out.keyval(), new_scheme);
}

auto output_image = Image<float>::create (argument[1], header);
DWI::export_grad_commandline (header);
auto input_image = Image<float>::open (argument[0]);
auto output_image = Image<float>::create (argument[1], header_out);
DWI::export_grad_commandline (header_out);

auto input_volumes = Adapter::make<Adapter::Extract1D> (input_image, 3, volumes);
threaded_copy_with_progress_message ("extracting volumes", input_volumes, output_image);
Expand Down
24 changes: 12 additions & 12 deletions cmd/tckglobal.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -218,8 +218,7 @@ void run ()
using namespace DWI::Tractography::GT;

// Inputs -----------------------------------------------------------------------------

auto dwi = Image<float>::open(argument[0]).with_direct_io(3);
Header header_in = Header::open (argument[0]);

Properties properties;
properties.resp_WM = load_matrix<float> (argument[1]);
Expand All @@ -237,7 +236,7 @@ void run ()
opt = get_options("mask");
if (opt.size()) {
mask = Image<bool>::open(opt[0][0]);
check_dimensions(dwi, mask, 0, 3);
check_dimensions(header_in, mask, 0, 3);
}

// Parameters -------------------------------------------------------------------------
Expand Down Expand Up @@ -304,9 +303,10 @@ void run ()
if (opt.size())
stats.open_stream(opt[0][0]);

auto dwi = header_in.get_image<float>().with_direct_io(3);
ParticleGrid pgrid (dwi);

ExternalEnergyComputer* Eext = new ExternalEnergyComputer(stats, dwi, properties);
ExternalEnergyComputer* Eext = new ExternalEnergyComputer(stats, header_in, properties);
InternalEnergyComputer* Eint = new InternalEnergyComputer(stats, pgrid);
Eint->setConnPot(cpot);
EnergySumComputer* Esum = new EnergySumComputer(stats, Eint, properties.lam_int, Eext, properties.lam_ext / ( wmscale2 * properties.weight*properties.weight));
Expand Down Expand Up @@ -345,14 +345,14 @@ void run ()


// Save fiso, tod and eext
Header header (dwi);
header.datatype() = DataType::Float32;
Header header_out (dwi);
header_out.datatype() = DataType::Float32;

opt = get_options("fod");
if (opt.size()) {
INFO("Saving fODF image to file");
header.size(3) = Math::SH::NforL(properties.Lmax);
auto fODF = Image<float>::create (opt[0][0], header);
header_out.size(3) = Math::SH::NforL(properties.Lmax);
auto fODF = Image<float>::create (opt[0][0], header_out);
auto f = __copy_fod<float>(properties.Lmax, properties.weight, !get_options("noapo").size());
ThreadedLoop(Eext->getTOD(), 0, 3).run(f, Eext->getTOD(), fODF);
}
Expand All @@ -361,8 +361,8 @@ void run ()
if (opt.size()) {
if (properties.resp_ISO.size() > 0) {
INFO("Saving isotropic fractions to file");
header.size(3) = properties.resp_ISO.size();
auto Fiso = Image<float>::create (opt[0][0], header);
header_out.size(3) = properties.resp_ISO.size();
auto Fiso = Image<float>::create (opt[0][0], header_out);
threaded_copy(Eext->getFiso(), Fiso);
}
else {
Expand All @@ -373,8 +373,8 @@ void run ()
opt = get_options("eext");
if (opt.size()) {
INFO("Saving external energy to file");
header.ndim() = 3;
auto EextI = Image<float>::create (opt[0][0], header);
header_out.ndim() = 3;
auto EextI = Image<float>::create (opt[0][0], header_out);
threaded_copy(Eext->getEext(), EextI);
}

Expand Down
15 changes: 11 additions & 4 deletions core/axes.h
Original file line number Diff line number Diff line change
Expand Up @@ -19,6 +19,7 @@


#include <array>
#include <set>

#include "types.h"

Expand All @@ -31,15 +32,21 @@ namespace MR



using permutations_type = std::array<size_t, 3>;
using permutations_type = std::array<ssize_t, 3>;
using flips_type = std::array<bool, 3>;
class Shuffle {
class Shuffle { NOMEMALIGN
public:
Shuffle() : permutations ({0, 1, 2}), flips ({false, false, false}) {}
operator bool() const {
Shuffle() : permutations ({-1, -1, -1}), flips ({false, false, false}) {}
bool is_identity() const {
return (permutations[0] != 0 || permutations[1] != 1 || permutations[2] != 2 || //
flips[0] || flips[1] || flips[2]);
}
bool is_set() const {
return permutations != permutations_type{-1, -1, -1};
}
bool valid() const {
return std::set<ssize_t>(permutations.begin(), permutations.end()) == std::set<ssize_t>({0, 1, 2});
}
permutations_type permutations;
flips_type flips;
};
Expand Down
7 changes: 5 additions & 2 deletions core/dwi/gradient.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -115,6 +115,8 @@ namespace MR

Eigen::MatrixXd load_bvecs_bvals (const Header& header, const std::string& bvecs_path, const std::string& bvals_path)
{
assert (header.realignment().orig_transform().matrix().allFinite());

Eigen::MatrixXd bvals, bvecs;
try {
bvals = load_matrix<> (bvals_path);
Expand Down Expand Up @@ -167,8 +169,8 @@ namespace MR
// transform have negative determinant:
if (adjusted_transform.linear().determinant() > 0.0)
bvecs.row(0) = -bvecs.row(0);
save_matrix(bvecs, bvecs_path, KeyValues(), false);
save_matrix(grad.col(3), bvals_path, KeyValues(), false);
save_matrix (bvecs, bvecs_path, KeyValues(), false);
save_matrix (grad.col(3).transpose(), bvals_path, KeyValues(), false);
}


Expand Down Expand Up @@ -257,6 +259,7 @@ namespace MR
"(maximum scaling factor = " + str(max_scaling_factor) + ")");
}
}
assert (grad.allFinite());

// write the scheme as interpreted back into the header if:
// - vector normalisation effect is large, regardless of whether or not b-value scaling was applied
Expand Down
Loading

0 comments on commit 2be257c

Please sign in to comment.