From 178f1ccbf83ccbcab2f12ef1167c0a6f3b57973f Mon Sep 17 00:00:00 2001 From: Robert Smith Date: Sun, 22 Sep 2024 14:39:46 +1000 Subject: [PATCH] Fixes to bvecs import 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. --- cmd/amp2sh.cpp | 28 +++++----- cmd/dwi2adc.cpp | 21 ++++---- cmd/dwi2mask.cpp | 10 ++-- cmd/dwi2tensor.cpp | 35 ++++++------ cmd/dwiextract.cpp | 23 ++++---- cmd/tckglobal.cpp | 24 ++++----- core/axes.h | 15 ++++-- core/dwi/gradient.cpp | 7 ++- core/file/nifti_utils.cpp | 4 +- core/filter/dwi_brain_mask.h | 3 +- core/header.cpp | 8 +-- core/header.h | 16 +++--- core/image.h | 1 + core/metadata/phase_encoding.cpp | 19 +++++-- core/metadata/phase_encoding.h | 27 +++++----- core/metadata/slice_encoding.cpp | 8 +-- src/dwi/tractography/GT/externalenergy.cpp | 10 ++-- src/dwi/tractography/GT/externalenergy.h | 56 +++++++++---------- src/dwi/tractography/GT/particlegrid.cpp | 36 +++++++++---- src/dwi/tractography/GT/particlegrid.h | 57 +++++++------------- src/dwi/tractography/algorithms/tensor_det.h | 2 +- src/dwi/tractography/tracking/shared.cpp | 5 +- src/dwi/tractography/tracking/shared.h | 3 +- 23 files changed, 225 insertions(+), 193 deletions(-) diff --git a/cmd/amp2sh.cpp b/cmd/amp2sh.cpp index e28d4995c1..cd137b7b19 100644 --- a/cmd/amp2sh.cpp +++ b/cmd/amp2sh.cpp @@ -201,8 +201,9 @@ class Amp2SH { MEMALIGN(Amp2SH) void run () { - auto amp = Image::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 bzeros, dwis; Eigen::MatrixXd dirs; @@ -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 dir_vector; for (auto line : split_lines (hit->second)) { auto v = parse_floats (line); @@ -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::create (argument[1], header); - + auto amp = header_in.get_image().with_direct_io (3); Amp2SHCommon common (sh2amp, bzeros, dwis, normalise); + auto SH = Image::create (argument[1], header_out); opt = get_options ("rician"); if (opt.size()) { diff --git a/cmd/dwi2adc.cpp b/cmd/dwi2adc.cpp index 49812e0950..f49d41b169 100644 --- a/cmd/dwi2adc.cpp +++ b/cmd/dwi2adc.cpp @@ -80,11 +80,11 @@ class DWI2ADC { MEMALIGN(DWI2ADC) void run () { - auto dwi = Header::open (argument[0]).get_image(); - 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)); @@ -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::create (argument[1], header); + auto dwi = header_in.get_image(); + auto adc = Image::create (argument[1], header_out); ThreadedLoop ("computing ADC values", dwi, 0, 3) .run (DWI2ADC (binv, dwi_axis), dwi, adc); diff --git a/cmd/dwi2mask.cpp b/cmd/dwi2mask.cpp index 466f62c5ad..c95451cc9d 100644 --- a/cmd/dwi2mask.cpp +++ b/cmd/dwi2mask.cpp @@ -66,13 +66,13 @@ OPTIONS void run () { - auto input = Image::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().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::scratch (dwi_brain_mask_filter, "brain mask"); dwi_brain_mask_filter (input, temp_mask); diff --git a/cmd/dwi2tensor.cpp b/cmd/dwi2tensor.cpp index cca2d52b87..8d77c0bfbd 100644 --- a/cmd/dwi2tensor.cpp +++ b/cmd/dwi2tensor.cpp @@ -202,14 +202,14 @@ inline Processor processor (const Eigen: void run () { - auto dwi = Header::open (argument[0]).get_image(); - auto grad = DWI::get_DW_scheme (dwi); + auto header_in = Header::open (argument[0]); + auto grad = DWI::get_DW_scheme (header_in); Image mask; auto opt = get_options ("mask"); if (opt.size()) { mask = Image::open (opt[0][0]); - check_dimensions (dwi, mask, 0, 3); + check_dimensions (header_in, mask, 0, 3); } bool ols = get_options ("ols").size(); @@ -217,37 +217,38 @@ void run () // 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 predict; opt = get_options ("predicted_signal"); if (opt.size()) - predict = Image::create (opt[0][0], header); + predict = Image::create (opt[0][0], header_out); - header.size(3) = 6; - auto dt = Image::create (argument[1], header); + header_out.size(3) = 6; + auto dt = Image::create (argument[1], header_out); Image b0; opt = get_options ("b0"); if (opt.size()) { - header.ndim() = 3; - b0 = Image::create (opt[0][0], header); + header_out.ndim() = 3; + b0 = Image::create (opt[0][0], header_out); } Image dkt; opt = get_options ("dkt"); if (opt.size()) { - header.ndim() = 4; - header.size(3) = 15; - dkt = Image::create (opt[0][0], header); + header_out.ndim() = 4; + header_out.size(3) = 15; + dkt = Image::create (opt[0][0], header_out); } - Eigen::MatrixXd b = -DWI::grad2bmatrix (grad, opt.size()>0); + Eigen::MatrixXd b = -DWI::grad2bmatrix (grad, dkt.valid()); + auto dwi = Header::open (argument[0]).get_image(); ThreadedLoop("computing tensors", dwi, 0, 3).run (processor (b, ols, iter, mask, b0, dkt, predict), dwi, dt); } diff --git a/cmd/dwiextract.cpp b/cmd/dwiextract.cpp index d60d947f1e..791eebd3b7 100644 --- a/cmd/dwiextract.cpp +++ b/cmd/dwiextract.cpp @@ -66,10 +66,10 @@ void usage () void run() { - auto input_image = Image::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 @@ -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"); @@ -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::create (argument[1], header); - DWI::export_grad_commandline (header); + auto input_image = Image::open (argument[0]); + auto output_image = Image::create (argument[1], header_out); + DWI::export_grad_commandline (header_out); auto input_volumes = Adapter::make (input_image, 3, volumes); threaded_copy_with_progress_message ("extracting volumes", input_volumes, output_image); diff --git a/cmd/tckglobal.cpp b/cmd/tckglobal.cpp index 7c4f4e560e..03a7352809 100644 --- a/cmd/tckglobal.cpp +++ b/cmd/tckglobal.cpp @@ -218,8 +218,7 @@ void run () using namespace DWI::Tractography::GT; // Inputs ----------------------------------------------------------------------------- - - auto dwi = Image::open(argument[0]).with_direct_io(3); + Header header_in = Header::open (argument[0]); Properties properties; properties.resp_WM = load_matrix (argument[1]); @@ -237,7 +236,7 @@ void run () opt = get_options("mask"); if (opt.size()) { mask = Image::open(opt[0][0]); - check_dimensions(dwi, mask, 0, 3); + check_dimensions(header_in, mask, 0, 3); } // Parameters ------------------------------------------------------------------------- @@ -304,9 +303,10 @@ void run () if (opt.size()) stats.open_stream(opt[0][0]); + auto dwi = header_in.get_image().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)); @@ -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::create (opt[0][0], header); + header_out.size(3) = Math::SH::NforL(properties.Lmax); + auto fODF = Image::create (opt[0][0], header_out); auto f = __copy_fod(properties.Lmax, properties.weight, !get_options("noapo").size()); ThreadedLoop(Eext->getTOD(), 0, 3).run(f, Eext->getTOD(), fODF); } @@ -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::create (opt[0][0], header); + header_out.size(3) = properties.resp_ISO.size(); + auto Fiso = Image::create (opt[0][0], header_out); threaded_copy(Eext->getFiso(), Fiso); } else { @@ -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::create (opt[0][0], header); + header_out.ndim() = 3; + auto EextI = Image::create (opt[0][0], header_out); threaded_copy(Eext->getEext(), EextI); } diff --git a/core/axes.h b/core/axes.h index f7d2b88cfd..12de59cbcd 100644 --- a/core/axes.h +++ b/core/axes.h @@ -19,6 +19,7 @@ #include +#include #include "types.h" @@ -31,15 +32,21 @@ namespace MR - using permutations_type = std::array; + using permutations_type = std::array; using flips_type = std::array; - 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(permutations.begin(), permutations.end()) == std::set({0, 1, 2}); + } permutations_type permutations; flips_type flips; }; diff --git a/core/dwi/gradient.cpp b/core/dwi/gradient.cpp index 8d93dc2704..95e55fffee 100644 --- a/core/dwi/gradient.cpp +++ b/core/dwi/gradient.cpp @@ -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); @@ -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); } @@ -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 diff --git a/core/file/nifti_utils.cpp b/core/file/nifti_utils.cpp index 9ddf624d06..ae71d7c272 100644 --- a/core/file/nifti_utils.cpp +++ b/core/file/nifti_utils.cpp @@ -588,7 +588,7 @@ namespace MR strides.resize (3); auto order = Stride::order (strides); Axes::Shuffle result; - result.permutations = {order[0], order[1], order[2]}; + result.permutations = {ssize_t(order[0]), ssize_t(order[1]), ssize_t(order[2])}; result.flips = {strides[order[0]] < 0, strides[order[1]] < 0, strides[order[2]] < 0}; return result; } @@ -600,7 +600,7 @@ namespace MR { const Axes::Shuffle shuffle = axes_on_write(H); axes = shuffle.permutations; - if (!shuffle) + if (shuffle.is_identity()) return H.transform(); const auto& M_in = H.transform().matrix(); diff --git a/core/filter/dwi_brain_mask.h b/core/filter/dwi_brain_mask.h index a6b2684b97..0b2d9f739e 100644 --- a/core/filter/dwi_brain_mask.h +++ b/core/filter/dwi_brain_mask.h @@ -46,8 +46,9 @@ namespace MR * * Typical usage: * \code + * auto header = Header::open (argument[0]); + * auto grad = DWI::get_DW_scheme (header); * auto input = Image::open (argument[0]); - * auto grad = DWI::get_DW_scheme (input); * Filter::DWIBrainMask dwi_brain_mask_filter (input, grad); * auto output = Image::create (argument[1], dwi_brain_mask_filter); * dwi_brain_mask_filter (input, output); diff --git a/core/header.cpp b/core/header.cpp index 6151a154fe..24408db1cb 100644 --- a/core/header.cpp +++ b/core/header.cpp @@ -632,7 +632,7 @@ namespace MR return; realignment_.shuffle_ = Axes::get_shuffle_to_make_RAS(transform()); - if (!realignment_) + if (realignment_.is_identity()) return; auto M (transform()); @@ -811,9 +811,9 @@ namespace MR - Header::Realignment::Realignment() : - applied_transform_(applied_transform_type::Identity()) { - orig_transform_.matrix().fill(std::numeric_limits::quiet_NaN()); + Header::Realignment::Realignment() { + applied_transform_.matrix().fill(0); + orig_transform_.matrix().fill(std::numeric_limits::signaling_NaN()); } diff --git a/core/header.h b/core/header.h index 6abaa98435..0822dba4b3 100644 --- a/core/header.h +++ b/core/header.h @@ -111,9 +111,11 @@ namespace MR Header (static_cast (original)) { } //! copy constructor from type of class other than Header - /*! This copies all relevant parameters over from \a original. */ + /*! This copies all relevant parameters over from \a original. + * Note that information about transform realignment on image load will not be available. + */ template ::value, void*>::type = nullptr> - Header (const HeaderType& original) : + explicit Header (const HeaderType& original) : transform_ (original.transform()), name_ (original.name()), keyval_ (original.keyval()), @@ -203,18 +205,20 @@ namespace MR // Class to store all information relating to internal transform realignment class Realignment { + MEMALIGN(Realignment) public: // From one image space to another image space; // linear component is permutations & flips only, // transformation is in voxel count, // therefore can store as integer + // TODO Calculate translations; turn into affine transform; verify using applied_transform_type = Eigen::Matrix; Realignment(); - Realignment (Header&); - operator bool() const { return bool(shuffle_); } - const std::array& permutations() const { return shuffle_.permutations; } + bool is_identity() const { return shuffle_.is_identity(); } + bool valid() const { return shuffle_.valid(); } + const Axes::permutations_type& permutations() const { return shuffle_.permutations; } size_t permutation (const size_t axis) const { assert(axis < 3); return shuffle_.permutations[axis]; } - const std::array& flips() const { return shuffle_.flips; } + const Axes::flips_type& flips() const { return shuffle_.flips; } bool flip (const size_t axis) const { assert(axis < 3); return shuffle_.flips[axis]; } const transform_type& orig_transform() const { return orig_transform_; } const applied_transform_type& applied_transform() const { return applied_transform_; } diff --git a/core/image.h b/core/image.h index 3dbb3aab6a..d628c0d786 100644 --- a/core/image.h +++ b/core/image.h @@ -61,6 +61,7 @@ namespace MR FORCE_INLINE const std::string& name() const { return buffer->name(); } FORCE_INLINE const transform_type& transform() const { return buffer->transform(); } + //FORCE_INLINE const Header::Realignment& realignment() const { return buffer->realignment(); } FORCE_INLINE size_t ndim () const { return buffer->ndim(); } FORCE_INLINE ssize_t size (size_t axis) const { return buffer->size (axis); } diff --git a/core/metadata/phase_encoding.cpp b/core/metadata/phase_encoding.cpp index 4794502ba0..d130f8cb7f 100644 --- a/core/metadata/phase_encoding.cpp +++ b/core/metadata/phase_encoding.cpp @@ -212,7 +212,7 @@ namespace MR { DEBUG("No phase encoding information found for transformation with load of image \"" + H.name() + "\""); return; } - if (!H.realignment()) { + if (H.realignment().is_identity()) { INFO("No transformation of phase encoding data for load of image \"" + H.name() + "\" required"); return; } @@ -223,7 +223,7 @@ namespace MR { scheme_type transform_for_image_load(const scheme_type& pe_scheme, const Header& H) { - if (!H.realignment()) + if (H.realignment().is_identity()) return pe_scheme; scheme_type result(pe_scheme.rows(), pe_scheme.cols()); for (ssize_t row = 0; row != pe_scheme.rows(); ++row) { @@ -251,7 +251,7 @@ namespace MR { if (pe_scheme.rows() == 0) return pe_scheme; Axes::Shuffle shuffle = File::NIfTI::axes_on_write(H); - if (!shuffle) { + if (shuffle.is_identity()) { INFO("No transformation of phase encoding data required for export to file:" " output image will be RAS"); return pe_scheme; @@ -358,6 +358,19 @@ namespace MR { + void save(const scheme_type& PE, const std::string& path) { + File::OFStream out(path); + for (ssize_t row = 0; row != PE.rows(); ++row) { + // Write phase-encode direction as integers; other information as floating-point + out << PE.template block<1, 3>(row, 0).template cast(); + if (PE.cols() > 3) + out << " " << PE.block(row, 3, 1, PE.cols() - 3); + out << "\n"; + } + } + + + } } } diff --git a/core/metadata/phase_encoding.h b/core/metadata/phase_encoding.h index 296984002e..e10f30ff4d 100644 --- a/core/metadata/phase_encoding.h +++ b/core/metadata/phase_encoding.h @@ -93,20 +93,17 @@ namespace MR { void transform_for_nifti_write(KeyValues& keyval, const Header& H); scheme_type transform_for_nifti_write(const scheme_type& pe_scheme, const Header& H); - namespace { - void _save(const scheme_type& PE, const std::string& path) { - File::OFStream out(path); - for (ssize_t row = 0; row != PE.rows(); ++row) { - // Write phase-encode direction as integers; other information as floating-point - out << PE.template block<1, 3>(row, 0).template cast(); - if (PE.cols() > 3) - out << " " << PE.block(row, 3, 1, PE.cols() - 3); - out << "\n"; - } - } - } // namespace + void save(const scheme_type& PE, const std::string& path); + + template + void save(const HeaderType &header, const std::string &path) { + const scheme_type scheme = get_scheme(header); + if (scheme.rows() == 0) + throw Exception ("No phase encoding scheme in header of image \"" + header.name() + "\" to save"); + save(scheme, header, path); + } - //! Save a phase-encoding scheme to file + //! Save a phase-encoding scheme associated with an image to file // Note that because the output table requires permutation / sign flipping // only if the output target image is a NIfTI, the output file name must have // already been set @@ -119,9 +116,9 @@ namespace MR { } if (Path::has_suffix(header.name(), {".mgh", ".mgz", ".nii", ".nii.gz", ".img"})) { - _save(transform_for_nifti_write(PE, header), path); + save(transform_for_nifti_write(PE, header), path); } else { - _save(PE, path); + save(PE, path); } } diff --git a/core/metadata/slice_encoding.cpp b/core/metadata/slice_encoding.cpp index e9109826d9..04235540ec 100644 --- a/core/metadata/slice_encoding.cpp +++ b/core/metadata/slice_encoding.cpp @@ -33,7 +33,7 @@ namespace MR { auto slice_encoding_it = keyval.find("SliceEncodingDirection"); auto slice_timing_it = keyval.find("SliceTiming"); if (!(slice_encoding_it == keyval.end() && slice_timing_it == keyval.end())) { - if (!header.realignment()) { + if (header.realignment().is_identity()) { INFO("No transformation of slice encoding direction for load of image \"" + header.name() + "\" required"); return; } @@ -71,7 +71,7 @@ namespace MR { return; const Axes::Shuffle shuffle = File::NIfTI::axes_on_write(H); - if (!shuffle) { + if (shuffle.is_identity()) { INFO("No need to transform slice encoding information for NIfTI image write:" " image is already RAS"); return; @@ -106,8 +106,8 @@ namespace MR { std::string resolve_slice_timing (const std::string &one, const std::string &two) { if (one == "variable" || two == "variable") return "variable"; - std::vector one_split = split(one, ","); - std::vector two_split = split(two, ","); + const vector one_split = split(one, ","); + const vector two_split = split(two, ","); if (one_split.size() != two_split.size()) { DEBUG("Slice timing vectors of inequal length"); return "invalid"; diff --git a/src/dwi/tractography/GT/externalenergy.cpp b/src/dwi/tractography/GT/externalenergy.cpp index 6014d56a27..012f92055a 100644 --- a/src/dwi/tractography/GT/externalenergy.cpp +++ b/src/dwi/tractography/GT/externalenergy.cpp @@ -28,17 +28,17 @@ namespace MR { namespace Tractography { namespace GT { - ExternalEnergyComputer::ExternalEnergyComputer(Stats& stat, const Image& dwimage, const Properties& props) + ExternalEnergyComputer::ExternalEnergyComputer(Stats& stat, Header& dwiheader, const Properties& props) : EnergyComputer(stat), - dwi(dwimage), - T(Transform(dwimage).scanner2voxel), + dwi(dwiheader.get_image().with_direct_io(3)), + T(Transform(dwiheader).scanner2voxel), lmax(props.Lmax), ncols(Math::SH::NforL(lmax)), nf(props.resp_ISO.size()), beta(props.beta), mu(props.ppot*M_sqrt4PI), dE(0.0) { DEBUG("Initialise computation of external energy."); // Create images -------------------------------------------------------------- - Header header (dwimage); + Header header (dwiheader); header.datatype() = DataType::Float32; header.size(3) = ncols; tod = Image::scratch(header, "TOD image"); @@ -54,7 +54,7 @@ namespace MR { eext = Image::scratch(header, "external energy"); // Set kernel matrices -------------------------------------------------------- - auto grad = DWI::get_DW_scheme(dwimage); + auto grad = DWI::get_DW_scheme(dwiheader); nrows = grad.rows(); DWI::Shells shells (grad); diff --git a/src/dwi/tractography/GT/externalenergy.h b/src/dwi/tractography/GT/externalenergy.h index 332c005d3a..9267c4a394 100644 --- a/src/dwi/tractography/GT/externalenergy.h +++ b/src/dwi/tractography/GT/externalenergy.h @@ -30,83 +30,83 @@ namespace MR { namespace DWI { namespace Tractography { namespace GT { - + class ExternalEnergyComputer : public EnergyComputer { MEMALIGN(ExternalEnergyComputer) public: - - ExternalEnergyComputer(Stats& stat, const Image& dwimage, const Properties& props); - - + + ExternalEnergyComputer(Stats& stat, Header& dwiheader, const Properties& props); + + Image& getTOD() { return tod; } Image& getFiso() { return fiso; } Image& getEext() { return eext; } - + void resetEnergy(); - + double stageAdd(const Point_t& pos, const Point_t& dir) { add(pos, dir, 1.0); return eval(); } - + double stageShift(const Particle* par, const Point_t& pos, const Point_t& dir) { add(par->getPosition(), par->getDirection(), -1.0); add(pos, dir, 1.0); return eval(); } - + double stageRemove(const Particle* par) { add(par->getPosition(), par->getDirection(), -1.0); return eval(); } - + void acceptChanges(); - + void clearChanges(); - + EnergyComputer* clone() const { return new ExternalEnergyComputer(*this); } - - - + + + protected: - + Image dwi; Image tod; Image fiso; Image eext; - + transform_type T; - - int lmax; + + int lmax; size_t nrows, ncols, nf; double beta, mu, dE; Eigen::MatrixXd K, Ak; Eigen::VectorXd y, t, d, fk; - + Math::ICLS::Problem nnls; - + vector changes_vox; vector changes_tod; vector changes_fiso; vector changes_eext; - - + + void add(const Point_t& pos, const Point_t& dir, const double factor = 1.); - + void add2vox(const Eigen::Vector3i& vox, const double w); - + double eval(); - + double calcEnergy(); - + inline double hanning(const double w) const { return (w <= (1.0-beta)/2) ? 0.0 : (w >= (1.0+beta)/2) ? 1.0 : (1 - std::cos(Math::pi * (w-(1.0-beta)/2)/beta )) / 2; } - + }; diff --git a/src/dwi/tractography/GT/particlegrid.cpp b/src/dwi/tractography/GT/particlegrid.cpp index fb2b8018eb..67019098e7 100644 --- a/src/dwi/tractography/GT/particlegrid.cpp +++ b/src/dwi/tractography/GT/particlegrid.cpp @@ -20,15 +20,33 @@ namespace MR { namespace DWI { namespace Tractography { namespace GT { - - + + + ParticleGrid::ParticleGrid(const Header& H) + { + DEBUG("Initialise particle grid."); + dims[0] = Math::ceil( H.size(0) * H.spacing(0) / (2.0*Particle::L) ); + dims[1] = Math::ceil( H.size(1) * H.spacing(1) / (2.0*Particle::L) ); + dims[2] = Math::ceil( H.size(2) * H.spacing(2) / (2.0*Particle::L) ); + grid.resize(dims[0]*dims[1]*dims[2]); + + // Initialise scanner-to-grid transform + Eigen::DiagonalMatrix newspacing (2.0*Particle::L, 2.0*Particle::L, 2.0*Particle::L); + Eigen::Vector3d shift (H.spacing(0)/2.0 - Particle::L, + H.spacing(1)/2.0 - Particle::L, + H.spacing(2)/2.0 - Particle::L); + T_s2g = H.transform() * newspacing; + T_s2g = T_s2g.inverse().translate(shift); + } + + void ParticleGrid::add(const Point_t &pos, const Point_t &dir) { Particle* p = pool.create(pos, dir); size_t gidx = pos2idx(pos); grid[gidx].push_back(p); } - + void ParticleGrid::shift(Particle *p, const Point_t& pos, const Point_t& dir) { size_t gidx0 = pos2idx(p->getPosition()); @@ -38,27 +56,27 @@ namespace MR { p->setDirection(dir); grid[gidx1].push_back(p); } - + void ParticleGrid::remove(Particle* p) { size_t gidx0 = pos2idx(p->getPosition()); grid[gidx0].erase(std::remove(grid[gidx0].begin(), grid[gidx0].end(), p), grid[gidx0].end()); pool.destroy(p); } - + void ParticleGrid::clear() { grid.clear(); pool.clear(); } - + const ParticleGrid::ParticleVectorType* ParticleGrid::at(const ssize_t x, const ssize_t y, const ssize_t z) const { if ((x < 0) || (size_t(x) >= dims[0]) || (y < 0) || (size_t(y) >= dims[1]) || (z < 0) || (size_t(z) >= dims[2])) // out of bounds return nullptr; return &grid[xyz2idx(x, y, z)]; } - + void ParticleGrid::exportTracks(Tractography::Writer &writer) { std::lock_guard lock (mutex); @@ -70,7 +88,7 @@ namespace MR { // Loop through all unvisited particles for (ParticleVectorType& gridvox : grid) { - for (Particle* par0 : gridvox) + for (Particle* par0 : gridvox) { par = par0; if (!par->isVisited()) @@ -114,7 +132,7 @@ namespace MR { } } } - + } } diff --git a/src/dwi/tractography/GT/particlegrid.h b/src/dwi/tractography/GT/particlegrid.h index d0e86bb8e5..fadd63bbaf 100644 --- a/src/dwi/tractography/GT/particlegrid.h +++ b/src/dwi/tractography/GT/particlegrid.h @@ -32,62 +32,45 @@ namespace MR { namespace DWI { namespace Tractography { namespace GT { - + /** * @brief The ParticleGrid class */ class ParticleGrid { MEMALIGN(ParticleGrid) public: - + using ParticleVectorType = vector; - - template - ParticleGrid(const HeaderType& image) - { - DEBUG("Initialise particle grid."); - dims[0] = Math::ceil( image.size(0) * image.spacing(0) / (2.0*Particle::L) ); - dims[1] = Math::ceil( image.size(1) * image.spacing(1) / (2.0*Particle::L) ); - dims[2] = Math::ceil( image.size(2) * image.spacing(2) / (2.0*Particle::L) ); - grid.resize(dims[0]*dims[1]*dims[2]); - - // Initialise scanner-to-grid transform - Eigen::DiagonalMatrix newspacing (2.0*Particle::L, 2.0*Particle::L, 2.0*Particle::L); - Eigen::Vector3d shift (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; - T_s2g = T_s2g.inverse().translate(shift); - } - + + ParticleGrid(const Header&); ParticleGrid(const ParticleGrid&) = delete; ParticleGrid& operator=(const ParticleGrid&) = delete; - + ~ParticleGrid() { clear(); } - + inline unsigned int getTotalCount() const { return pool.size(); } - + void add(const Point_t& pos, const Point_t& dir); - + void shift(Particle* p, const Point_t& pos, const Point_t& dir); - + void remove(Particle* p); - + void clear(); - + const ParticleVectorType* at(const ssize_t x, const ssize_t y, const ssize_t z) const; - + inline Particle* getRandom() { return pool.random(); } - + void exportTracks(Tractography::Writer& writer); - - + + protected: std::mutex mutex; ParticlePool pool; @@ -95,15 +78,15 @@ namespace MR { Math::RNG rng; transform_type T_s2g; size_t dims[3]; - - + + inline size_t pos2idx(const Point_t& pos) const { size_t x, y, z; pos2xyz(pos, x, y, z); return xyz2idx(x, y, z); } - + public: inline void pos2xyz(const Point_t& pos, size_t& x, size_t& y, size_t& z) const { @@ -113,13 +96,13 @@ namespace MR { y = Math::round(gpos[1]); z = Math::round(gpos[2]); } - + protected: inline size_t xyz2idx(const size_t x, const size_t y, const size_t z) const { return z + dims[2] * (y + dims[1] * x); } - + }; } diff --git a/src/dwi/tractography/algorithms/tensor_det.h b/src/dwi/tractography/algorithms/tensor_det.h index 44fc1781a0..870401d118 100644 --- a/src/dwi/tractography/algorithms/tensor_det.h +++ b/src/dwi/tractography/algorithms/tensor_det.h @@ -65,7 +65,7 @@ namespace MR properties["method"] = "TensorDet"; try { - auto grad = DWI::get_DW_scheme (source); + auto grad = DWI::get_DW_scheme (source_header); auto bmat_double = grad2bmatrix (grad); binv = Math::pinv (bmat_double).cast(); bmat = bmat_double.cast(); diff --git a/src/dwi/tractography/tracking/shared.cpp b/src/dwi/tractography/tracking/shared.cpp index 7109a9d2bd..01b8034a52 100644 --- a/src/dwi/tractography/tracking/shared.cpp +++ b/src/dwi/tractography/tracking/shared.cpp @@ -29,7 +29,8 @@ namespace MR SharedBase::SharedBase (const std::string& diff_path, Properties& property_set) : - source (Image::open (diff_path).with_direct_io (3)), + source_header (Header::open (diff_path)), + source (source_header.get_image().with_direct_io (3)), properties (property_set), init_dir ({ NaN, NaN, NaN }), min_num_points_preds (0), @@ -63,7 +64,7 @@ namespace MR properties.set (rk4, "rk4"); properties.set (stop_on_all_include, "stop_on_all_include"); - properties["source"] = source.name(); + properties["source"] = source_header.name(); max_num_seeds = Defaults::seed_to_select_ratio * max_num_tracks; properties.set (max_num_seeds, "max_num_seeds"); diff --git a/src/dwi/tractography/tracking/shared.h b/src/dwi/tractography/tracking/shared.h index e2bb586e76..911e4a44f4 100644 --- a/src/dwi/tractography/tracking/shared.h +++ b/src/dwi/tractography/tracking/shared.h @@ -56,6 +56,7 @@ namespace MR virtual ~SharedBase(); + Header source_header; Image source; Properties& properties; Eigen::Vector3f init_dir; @@ -76,7 +77,7 @@ namespace MR float vox () const { - return std::pow (source.spacing(0)*source.spacing(1)*source.spacing(2), float (1.0/3.0)); + return std::pow (source_header.spacing(0)*source_header.spacing(1)*source_header.spacing(2), float (1.0/3.0)); } void set_step_and_angle (const float stepsize, const float angle, bool is_higher_order);