diff --git a/cpp/cmd/amp2sh.cpp b/cpp/cmd/amp2sh.cpp index 20a65f4faf..09358d7f4d 100644 --- a/cpp/cmd/amp2sh.cpp +++ b/cpp/cmd/amp2sh.cpp @@ -21,7 +21,7 @@ #include "file/matrix.h" #include "image.h" #include "math/SH.h" -#include "phase_encoding.h" +#include "metadata/phase_encoding.h" #include "progressbar.h" using namespace MR; @@ -219,7 +219,7 @@ void run() { DWI::stash_DW_scheme(header, grad); } } - PhaseEncoding::clear_scheme(header); + Metadata::PhaseEncoding::clear_scheme(header.keyval()); auto sh2amp = DWI::compute_SH2amp_mapping(dirs, true, 8); diff --git a/cpp/cmd/dwi2adc.cpp b/cpp/cmd/dwi2adc.cpp index 268200fc1c..ef0dcd5671 100644 --- a/cpp/cmd/dwi2adc.cpp +++ b/cpp/cmd/dwi2adc.cpp @@ -19,7 +19,7 @@ #include "dwi/gradient.h" #include "image.h" #include "math/least_squares.h" -#include "phase_encoding.h" +#include "metadata/phase_encoding.h" #include "progressbar.h" using namespace MR; @@ -44,7 +44,7 @@ void usage () "adopts a 2-stage fitting strategy, in which the ADC is first estimated based on " "the DWI data with b > cutoff, and the other parameters are estimated subsequently. " "The output consists of 4 volumes, respectively S(0), D, f, and D'." - + + "Note that this command ignores the gradient orientation entirely. This approach is " "therefore only suited for mean DWI (trace-weighted) images or for DWI data collected " "with isotropically-distributed gradient directions."; @@ -60,12 +60,12 @@ void usage () + Argument ("bval").type_integer (0, 1000) + DWI::GradImportOptions(); - + REFERENCES + "Le Bihan, D.; Breton, E.; Lallemand, D.; Aubin, M.L.; Vignaud, J.; Laval-Jeantet, M. " "Separation of diffusion and perfusion in intravoxel incoherent motion MR imaging. " "Radiology, 1988, 168, 497–505." - + + "Jalnefjord, O.; Andersson, M.; Montelius; M.; Starck, G.; Elf, A.; Johanson, V.; Svensson, J.; Ljungberg, M. " "Comparison of methods for estimation of the intravoxel incoherent motion (IVIM) " "diffusion coefficient (D) and perfusion fraction (f). " @@ -159,7 +159,7 @@ void run() { Header header(dwi); header.datatype() = DataType::Float32; DWI::stash_DW_scheme(header, grad); - PhaseEncoding::clear_scheme(header); + Metadata::PhaseEncoding::clear_scheme(header.keyval()); header.ndim() = 4; header.size(3) = ivim ? 4 : 2; diff --git a/cpp/cmd/dwi2fod.cpp b/cpp/cmd/dwi2fod.cpp index a529e07f1f..d927ceb994 100644 --- a/cpp/cmd/dwi2fod.cpp +++ b/cpp/cmd/dwi2fod.cpp @@ -21,7 +21,7 @@ #include "header.h" #include "image.h" #include "math/SH.h" -#include "phase_encoding.h" +#include "metadata/phase_encoding.h" #include "dwi/sdeconv/csd.h" #include "dwi/sdeconv/msmt_csd.h" @@ -258,7 +258,7 @@ void run() { shared.init(); DWI::stash_DW_scheme(header_out, shared.grad); - PhaseEncoding::clear_scheme(header_out); + Metadata::PhaseEncoding::clear_scheme(header_out.keyval()); header_out.size(3) = shared.nSH(); auto fod = Image::create(argument[3], header_out); diff --git a/cpp/cmd/dwi2tensor.cpp b/cpp/cmd/dwi2tensor.cpp index d673a2eec9..4e8b00e7f6 100644 --- a/cpp/cmd/dwi2tensor.cpp +++ b/cpp/cmd/dwi2tensor.cpp @@ -22,7 +22,7 @@ #include "file/matrix.h" #include "image.h" #include "math/constrained_least_squares.h" -#include "phase_encoding.h" +#include "metadata/phase_encoding.h" #include "progressbar.h" using namespace MR; @@ -295,7 +295,8 @@ void run() { Header header(dwi); header.datatype() = DataType::Float32; header.ndim() = 4; - PhaseEncoding::clear_scheme(header); + DWI::stash_DW_scheme(header, grad); + Metadata::PhaseEncoding::clear_scheme(header.keyval()); Image predict; opt = get_options("predicted_signal"); diff --git a/cpp/cmd/dwiextract.cpp b/cpp/cmd/dwiextract.cpp index e077242afe..1130825916 100644 --- a/cpp/cmd/dwiextract.cpp +++ b/cpp/cmd/dwiextract.cpp @@ -19,7 +19,7 @@ #include "command.h" #include "dwi/gradient.h" #include "image.h" -#include "phase_encoding.h" +#include "metadata/phase_encoding.h" #include "progressbar.h" using namespace MR; @@ -65,8 +65,8 @@ void usage() { + DWI::GradImportOptions() + DWI::ShellsOption + DWI::GradExportOptions() - + PhaseEncoding::ImportOptions - + PhaseEncoding::SelectOptions + + Metadata::PhaseEncoding::ImportOptions + + Metadata::PhaseEncoding::SelectOptions + Stride::Options; } // clang-format off @@ -108,7 +108,7 @@ void run() { } auto opt = get_options("pe"); - const auto pe_scheme = PhaseEncoding::get_scheme(input_image); + const auto pe_scheme = Metadata::PhaseEncoding::get_scheme(input_image); if (!opt.empty()) { if (!pe_scheme.rows()) throw Exception("Cannot filter volumes by phase-encoding: No such information present"); @@ -154,7 +154,7 @@ void run() { 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]); - PhaseEncoding::set_scheme(header, new_scheme); + Metadata::PhaseEncoding::set_scheme(header.keyval(), new_scheme); } auto output_image = Image::create(argument[1], header); diff --git a/cpp/cmd/mrcalc.cpp b/cpp/cmd/mrcalc.cpp index 3513c708f9..98d9da44c0 100644 --- a/cpp/cmd/mrcalc.cpp +++ b/cpp/cmd/mrcalc.cpp @@ -320,7 +320,6 @@ UNARY_OP( #include "image.h" #include "math/rng.h" #include "memory.h" -#include "phase_encoding.h" using namespace MR; using namespace App; @@ -823,7 +822,7 @@ void get_header(const StackEntry &entry, Header &header) { header.spacing(n) = entry.image->spacing(n); } - header.merge_keyval(*entry.image); + header.merge_keyval(entry.image->keyval()); } class ThreadFunctor { diff --git a/cpp/cmd/mrcat.cpp b/cpp/cmd/mrcat.cpp index 8872cb4a50..08a80ad27d 100644 --- a/cpp/cmd/mrcat.cpp +++ b/cpp/cmd/mrcat.cpp @@ -18,7 +18,6 @@ #include "command.h" #include "dwi/gradient.h" #include "image.h" -#include "phase_encoding.h" #include "progressbar.h" using namespace MR; diff --git a/cpp/cmd/mrconvert.cpp b/cpp/cmd/mrconvert.cpp index f4f60bf762..14aa20372a 100644 --- a/cpp/cmd/mrconvert.cpp +++ b/cpp/cmd/mrconvert.cpp @@ -24,7 +24,7 @@ #include "file/ofstream.h" #include "header.h" #include "image.h" -#include "phase_encoding.h" +#include "metadata/phase_encoding.h" #include "transform.h" #include "types.h" @@ -244,8 +244,8 @@ void usage() { + DWI::bvalue_scaling_option + DWI::GradExportOptions() - + PhaseEncoding::ImportOptions - + PhaseEncoding::ExportOptions; + + Metadata::PhaseEncoding::ImportOptions + + Metadata::PhaseEncoding::ExportOptions; } // clang-format on @@ -271,7 +271,7 @@ void permute_DW_scheme(Header &H, const std::vector &axes) { } void permute_PE_scheme(Header &H, const std::vector &axes) { - auto in = PhaseEncoding::parse_scheme(H); + auto in = Metadata::PhaseEncoding::parse_scheme(H.keyval(), H); if (!in.rows()) return; @@ -285,16 +285,27 @@ void permute_PE_scheme(Header &H, const std::vector &axes) { for (int row = 0; row != in.rows(); ++row) out.block<1, 3>(row, 0) = in.block<1, 3>(row, 0) * permute; - PhaseEncoding::set_scheme(H, out); + Metadata::PhaseEncoding::set_scheme(H.keyval(), out); } void permute_slice_direction(Header &H, const std::vector &axes) { - auto it = H.keyval().find("SliceEncodingDirection"); - if (it == H.keyval().end()) + using Metadata::BIDS::axis_vector_type; + auto slice_encoding_it = H.keyval().find("SliceEncodingDirection"); + auto slice_timing_it = H.keyval().find("SliceTiming"); + if (slice_encoding_it == H.keyval().end() && slice_timing_it == H.keyval().end()) return; - 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); + if (slice_encoding_it == H.keyval().end()) { + const axis_vector_type orig_dir({0, 0, 1}); + const axis_vector_type new_dir(orig_dir[axes[0]], orig_dir[axes[1]], orig_dir[axes[2]]); + slice_encoding_it->second = Metadata::BIDS::vector2axisid(new_dir); + WARN("Header field \"SliceEncodingDirection\" inferred to be \"k\" in input image" + " and then transformed according to axis permutation" + " in order to preserve validity of existing header field \"SliceTiming\""); + return; + } + const axis_vector_type orig_dir = Metadata::BIDS::axisid2vector(slice_encoding_it->second); + const axis_vector_type new_dir(orig_dir[axes[0]], orig_dir[axes[1]], orig_dir[axes[2]]); + slice_encoding_it->second = Metadata::BIDS::vector2axisid(new_dir); } template inline std::vector set_header(Header &header, const ImageType &input) { @@ -352,7 +363,7 @@ void copy_permute(const InputType &in, Header &header_out, const std::string &ou const auto axes = set_header(header_out, in); auto out = Image::create(output_filename, header_out, add_to_command_history); DWI::export_grad_commandline(out); - PhaseEncoding::export_commandline(out); + Metadata::PhaseEncoding::export_commandline(out); auto perm = Adapter::make(in, axes); threaded_copy_with_progress(perm, out, 0, std::numeric_limits::max(), 2); } @@ -381,20 +392,17 @@ void run() { throw; e.display(2); } + auto opt = get_options("json_import"); + if (!opt.empty()) + File::JSON::load(header_in, opt[0][0]); + if (!get_options("import_pe_table").empty() || !get_options("import_pe_eddy").empty()) + Metadata::PhaseEncoding::set_scheme(header_in.keyval(), Metadata::PhaseEncoding::get_scheme(header_in)); Header header_out(header_in); header_out.datatype() = DataType::from_command_line(header_out.datatype()); - if (header_in.datatype().is_complex() && !header_out.datatype().is_complex()) WARN("requested datatype is real but input datatype is complex - imaginary component will be ignored"); - if (!get_options("import_pe_table").empty() || !get_options("import_pe_eddy").empty()) - PhaseEncoding::set_scheme(header_out, PhaseEncoding::get_scheme(header_in)); - - auto opt = get_options("json_import"); - if (!opt.empty()) - File::JSON::load(header_out, opt[0][0]); - opt = get_options("copy_properties"); if (!opt.empty()) { header_out.keyval().clear(); @@ -478,17 +486,17 @@ void run() { } Eigen::MatrixXd pe_scheme; try { - pe_scheme = PhaseEncoding::get_scheme(header_in); + pe_scheme = Metadata::PhaseEncoding::get_scheme(header_in); if (pe_scheme.rows()) { Eigen::MatrixXd extract_scheme(pos[3].size(), pe_scheme.cols()); for (size_t vol = 0; vol != pos[3].size(); ++vol) extract_scheme.row(vol) = pe_scheme.row(pos[3][vol]); - PhaseEncoding::set_scheme(header_out, extract_scheme); + Metadata::PhaseEncoding::set_scheme(header_out.keyval(), extract_scheme); } } catch (...) { - WARN("Phase encoding scheme of input file does not match number of image volumes; omitting information from " - "output image"); - PhaseEncoding::set_scheme(header_out, Eigen::MatrixXd()); + WARN("Phase encoding scheme of input file does not match number of image volumes;" + " omitting information from output image"); + Metadata::PhaseEncoding::clear_scheme(header_out.keyval()); } } } diff --git a/cpp/cmd/mrdegibbs.cpp b/cpp/cmd/mrdegibbs.cpp index 3fd59ae481..4e7496aa1e 100644 --- a/cpp/cmd/mrdegibbs.cpp +++ b/cpp/cmd/mrdegibbs.cpp @@ -19,6 +19,7 @@ #include "command.h" #include "degibbs/unring2d.h" #include "degibbs/unring3d.h" +#include "metadata/bids.h" using namespace MR; using namespace App; @@ -155,7 +156,8 @@ void run() { WARN("If data were acquired using multi-slice encoding, run in default 2D mode."); } else { try { - const Eigen::Vector3d slice_encoding_axis_onehot = Axes::id2dir(slice_encoding_it->second); + const Metadata::BIDS::axis_vector_type slice_encoding_axis_onehot = + Metadata::BIDS::axisid2vector(slice_encoding_it->second); std::vector auto_slice_axes = {0, 0}; if (slice_encoding_axis_onehot[0]) auto_slice_axes = {1, 2}; diff --git a/cpp/cmd/mrinfo.cpp b/cpp/cmd/mrinfo.cpp index 155811a9b6..9c3835ec67 100644 --- a/cpp/cmd/mrinfo.cpp +++ b/cpp/cmd/mrinfo.cpp @@ -22,7 +22,7 @@ #include "file/json_utils.h" #include "header.h" #include "image_io/pipe.h" -#include "phase_encoding.h" +#include "metadata/phase_encoding.h" #include "types.h" #include @@ -120,7 +120,7 @@ void usage() { + Option ("shell_sizes", "list the number of volumes in each shell") + Option ("shell_indices", "list the image volumes attributed to each b-value shell") - + PhaseEncoding::ExportOptions + + Metadata::PhaseEncoding::ExportOptions + Option ("petable", "print the phase encoding table") + OptionGroup ("Handling of piped images") @@ -253,7 +253,7 @@ void run() { ImageIO::Pipe::delete_piped_images = false; const bool export_grad = check_option_group(GradExportOptions); - const bool export_pe = check_option_group(PhaseEncoding::ExportOptions); + const bool export_pe = check_option_group(Metadata::PhaseEncoding::ExportOptions); if (export_grad && argument.size() > 1) throw Exception("can only export DW gradient table to file if a single input image is provided"); @@ -311,7 +311,7 @@ void run() { if (transform) print_transform(header); if (petable) - std::cout << PhaseEncoding::get_scheme(header) << "\n"; + std::cout << Metadata::PhaseEncoding::get_scheme(header) << "\n"; for (size_t n = 0; n < properties.size(); ++n) print_properties(header, properties[n][0]); @@ -330,7 +330,7 @@ void run() { } DWI::export_grad_commandline(header); - PhaseEncoding::export_commandline(header); + Metadata::PhaseEncoding::export_commandline(header); if (json_keyval) File::JSON::write(header, *json_keyval, (argument.size() > 1 ? std::string("") : std::string(argument[0]))); diff --git a/cpp/cmd/mrmath.cpp b/cpp/cmd/mrmath.cpp index f4b23c21b4..48469b3d47 100644 --- a/cpp/cmd/mrmath.cpp +++ b/cpp/cmd/mrmath.cpp @@ -24,8 +24,8 @@ #include "math/math.h" #include "math/median.h" #include "memory.h" +#include "metadata/phase_encoding.h" #include "misc/voxel2vector.h" -#include "phase_encoding.h" #include "progressbar.h" #include @@ -360,7 +360,7 @@ void run() { } catch (...) { } DWI::clear_DW_scheme(header_out); - PhaseEncoding::clear_scheme(header_out); + Metadata::PhaseEncoding::clear_scheme(header_out.keyval()); } header_out.datatype() = DataType::from_command_line(DataType::Float32); @@ -449,7 +449,7 @@ void run() { throw Exception("Image " + path + " has axis with non-unary dimension beyond first input image " + header.name()); } - header.merge_keyval(temp); + header.merge_keyval(temp.keyval()); } // Instantiate a kernel depending on the operation requested diff --git a/cpp/cmd/transformconvert.cpp b/cpp/cmd/transformconvert.cpp index 779ad93f4f..637fc4e160 100644 --- a/cpp/cmd/transformconvert.cpp +++ b/cpp/cmd/transformconvert.cpp @@ -14,6 +14,7 @@ * For more details, see http://www.mrtrix.org/. */ +#include "axes.h" #include "command.h" #include "file/key_value.h" #include "file/matrix.h" @@ -64,15 +65,15 @@ void usage() { // clang-format on transform_type get_flirt_transform(const Header &header) { - std::vector axes; - transform_type nifti_transform = File::NIfTI::adjust_transform(header, axes); - if (nifti_transform.matrix().topLeftCorner<3, 3>().determinant() < 0.0) - return nifti_transform; + const transform_type ondisk_transform = header.realignment().orig_transform(); + if (ondisk_transform.matrix().topLeftCorner<3, 3>().determinant() < 0.0) + return ondisk_transform; transform_type coord_switch; coord_switch.setIdentity(); coord_switch(0, 0) = -1.0f; - coord_switch(0, 3) = (header.size(axes[0]) - 1) * header.spacing(axes[0]); - return nifti_transform * coord_switch; + coord_switch(0, 3) = + (header.size(header.realignment().permutation(0)) - 1) * header.spacing(header.realignment().permutation(0)); + return ondisk_transform * coord_switch; } // transform_type parse_surfer_transform (const Header& header) { diff --git a/cpp/core/axes.cpp b/cpp/core/axes.cpp index 39fb023b6d..13765116ee 100644 --- a/cpp/core/axes.cpp +++ b/cpp/core/axes.cpp @@ -16,69 +16,21 @@ #include "axes.h" -#include "exception.h" -#include "mrtrix.h" - namespace MR::Axes { -std::string dir2id(const Eigen::Vector3d &axis) { - if (axis[0] == -1) { - assert(!axis[1]); - assert(!axis[2]); - return "i-"; - } else if (axis[0] == 1) { - assert(!axis[1]); - assert(!axis[2]); - return "i"; - } else if (axis[1] == -1) { - assert(!axis[0]); - assert(!axis[2]); - return "j-"; - } else if (axis[1] == 1) { - assert(!axis[0]); - assert(!axis[2]); - return "j"; - } else if (axis[2] == -1) { - assert(!axis[0]); - assert(!axis[1]); - return "k-"; - } else if (axis[2] == 1) { - assert(!axis[0]); - assert(!axis[1]); - return "k"; - } else { - throw Exception("Malformed axis direction: \"" + str(axis.transpose()) + "\""); - } -} - -Eigen::Vector3d id2dir(const std::string &id) { - if (id == "i-") - return {-1, 0, 0}; - else if (id == "i") - return {1, 0, 0}; - else if (id == "j-") - return {0, -1, 0}; - else if (id == "j") - return {0, 1, 0}; - else if (id == "k-") - return {0, 0, -1}; - else if (id == "k") - return {0, 0, 1}; - else - throw Exception("Malformed image axis identifier: \"" + id + "\""); -} - -void get_shuffle_to_make_axial(const transform_type &T, std::array &perm, std::array &flip) { - perm = closest(T.matrix().topLeftCorner<3, 3>()); +Shuffle get_shuffle_to_make_RAS(const transform_type &T) { + Shuffle result; + result.permutations = closest(T.matrix().topLeftCorner<3, 3>()); // Figure out whether any of the rows of the transform point in the // opposite direction to the MRtrix convention - flip[perm[0]] = T(0, perm[0]) < 0.0; - flip[perm[1]] = T(1, perm[1]) < 0.0; - flip[perm[2]] = T(2, perm[2]) < 0.0; + result.flips[result.permutations[0]] = T(0, result.permutations[0]) < 0.0; + result.flips[result.permutations[1]] = T(1, result.permutations[1]) < 0.0; + result.flips[result.permutations[2]] = T(2, result.permutations[2]) < 0.0; + return result; } -std::array closest(const Eigen::Matrix3d &M) { - std::array result{3, 3, 3}; +permutations_type closest(const Eigen::Matrix3d &M) { + permutations_type result{3, 3, 3}; // Find which row of the transform is closest to each scanner axis Eigen::Matrix3d::Index index(0); M.row(0).cwiseAbs().maxCoeff(&index); diff --git a/cpp/core/axes.h b/cpp/core/axes.h index eb8d5fe308..e3b2d620c1 100644 --- a/cpp/core/axes.h +++ b/cpp/core/axes.h @@ -16,25 +16,30 @@ #pragma once -#include +#include #include "types.h" namespace MR::Axes { -//! convert axis directions between formats -/*! these helper functions convert the definition of - * 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::Vector3d &); -Eigen::Vector3d id2dir(const std::string &); +using permutations_type = std::array; +using flips_type = std::array; +class Shuffle { +public: + Shuffle() : permutations({0, 1, 2}), flips({false, false, false}) {} + operator bool() const { + return (permutations[0] != 0 || permutations[1] != 1 || permutations[2] != 2 || // + flips[0] || flips[1] || flips[2]); + } + permutations_type permutations; + flips_type flips; +}; //! determine the axis permutations and flips necessary to make an image //! appear approximately axial -void get_shuffle_to_make_axial(const transform_type &T, std::array &perm, std::array &flip); +Shuffle get_shuffle_to_make_RAS(const transform_type &T); //! determine which vectors of a 3x3 transform are closest to the three axis indices -std::array closest(const Eigen::Matrix3d &M); +permutations_type closest(const Eigen::Matrix3d &M); } // namespace MR::Axes diff --git a/cpp/core/dwi/gradient.cpp b/cpp/core/dwi/gradient.cpp index de62290fa8..9b6433f0b7 100644 --- a/cpp/core/dwi/gradient.cpp +++ b/cpp/core/dwi/gradient.cpp @@ -139,27 +139,12 @@ Eigen::MatrixXd load_bvecs_bvals(const Header &header, const std::string &bvecs_ str(bvecs.cols()) + ", image: " + str(num_volumes) + ")"); // bvecs format actually assumes a LHS coordinate system even if image is - // stored using RHS - x axis is flipped to make linear 3x3 part of + // stored using RHS; first axis is flipped to make linear 3x3 part of // transform have negative determinant: - std::vector order; - auto adjusted_transform = File::NIfTI::adjust_transform(header, order); - if (adjusted_transform.linear().determinant() > 0.0) - bvecs.row(0) = -bvecs.row(0); - - // account for the fact that bvecs are specified wrt original image axes, - // which may have been re-ordered and/or inverted by MRtrix to match the - // expected anatomical frame of reference: - Eigen::MatrixXd G(bvecs.cols(), 3); - for (ssize_t n = 0; n < G.rows(); ++n) { - G(n, order[0]) = header.stride(order[0]) > 0 ? bvecs(0, n) : -bvecs(0, n); - G(n, order[1]) = header.stride(order[1]) > 0 ? bvecs(1, n) : -bvecs(1, n); - G(n, order[2]) = header.stride(order[2]) > 0 ? bvecs(2, n) : -bvecs(2, n); - } - - // rotate gradients into scanner coordinate system: - Eigen::MatrixXd grad(G.rows(), 4); - - grad.leftCols<3>().transpose() = header.transform().rotation() * G.transpose(); + if (header.realignment().orig_transform().linear().determinant() > 0.0) + bvecs.row(0) *= -1.0; + Eigen::MatrixXd grad(bvecs.cols(), 4); + grad.leftCols<3>().transpose() = header.realignment().orig_transform().linear() * bvecs; grad.col(3) = bvals.row(0); return grad; @@ -167,31 +152,16 @@ Eigen::MatrixXd load_bvecs_bvals(const Header &header, const std::string &bvecs_ void save_bvecs_bvals(const Header &header, const std::string &bvecs_path, const std::string &bvals_path) { const auto grad = parse_DW_scheme(header); - - // rotate vectors from scanner space to image space - Eigen::MatrixXd G = grad.leftCols<3>() * header.transform().rotation(); - - // deal with FSL requiring gradient directions to coincide with data strides - // also transpose matrices in preparation for file output - std::vector order; - auto adjusted_transform = File::NIfTI::adjust_transform(header, order); - Eigen::MatrixXd bvecs(3, grad.rows()); - Eigen::MatrixXd bvals(1, grad.rows()); - for (ssize_t n = 0; n < G.rows(); ++n) { - bvecs(0, n) = header.stride(order[0]) > 0 ? G(n, order[0]) : -G(n, order[0]); - bvecs(1, n) = header.stride(order[1]) > 0 ? G(n, order[1]) : -G(n, order[1]); - bvecs(2, n) = header.stride(order[2]) > 0 ? G(n, order[2]) : -G(n, order[2]); - bvals(0, n) = grad(n, 3); - } - + Axes::permutations_type order; + const auto adjusted_transform = File::NIfTI::adjust_transform(header, order); + Eigen::MatrixXd bvecs = adjusted_transform.inverse().linear() * grad.leftCols<3>().transpose(); // bvecs format actually assumes a LHS coordinate system even if image is - // stored using RHS - x axis is flipped to make linear 3x3 part of + // stored using RHS; first axis is flipped to make linear 3x3 part of // transform have negative determinant: if (adjusted_transform.linear().determinant() > 0.0) bvecs.row(0) = -bvecs.row(0); - File::Matrix::save_matrix(bvecs, bvecs_path, KeyValues(), false); - File::Matrix::save_matrix(bvals, bvals_path, KeyValues(), false); + File::Matrix::save_matrix(grad.col(3), bvals_path, KeyValues(), false); } void clear_DW_scheme(Header &header) { diff --git a/cpp/core/dwi/gradient.h b/cpp/core/dwi/gradient.h index ae6faf0bb5..de717be4f3 100644 --- a/cpp/core/dwi/gradient.h +++ b/cpp/core/dwi/gradient.h @@ -30,6 +30,8 @@ namespace App { class OptionGroup; } +// TODO Move some to metadata + namespace DWI { App::OptionGroup GradImportOptions(); diff --git a/cpp/core/file/dicom/mapper.cpp b/cpp/core/file/dicom/mapper.cpp index b9ca1ed2cf..2b5666cc9e 100644 --- a/cpp/core/file/dicom/mapper.cpp +++ b/cpp/core/file/dicom/mapper.cpp @@ -26,7 +26,7 @@ #include "image_io/default.h" #include "image_io/mosaic.h" #include "image_io/variable_scaling.h" -#include "phase_encoding.h" +#include "metadata/phase_encoding.h" namespace MR::File::Dicom { @@ -284,7 +284,7 @@ std::unique_ptr dicom_to_mapper(MR::Header &H, std::vectoris_boolean()) { @@ -103,56 +107,16 @@ KeyValues read(const nlohmann::json &json, const KeyValues &preexisting) { throw Exception("JSON entry \"" + i.key() + "\" contains mixture of elements and arrays"); } } - for (const auto &kv : preexisting) { - if (kv.first == "comments" && result.find("comments") != result.end()) { - add_line(result["comments"], kv.second); - } else { - // Will not overwrite existing entries - result.insert(kv); - } - } return result; } -void read(const nlohmann::json &json, Header &header, const bool realign) { - header.keyval() = read(json, header.keyval()); - const bool do_realign = realign && Header::do_realign_transform; - - // The corresponding header may have been rotated on image load prior to the JSON - // being loaded. If this is the case, any fields that indicate an image axis - // number / direction need to be correspondingly modified. - std::array perm; - std::array flip; - header.realignment(perm, flip); - if (perm[0] == 0 && perm[1] == 1 && perm[2] == 2 && !flip[0] && !flip[1] && !flip[2]) - return; - - auto pe_scheme = PhaseEncoding::get_scheme(header); - if (pe_scheme.rows()) { - if (do_realign) { - PhaseEncoding::set_scheme(header, PhaseEncoding::transform_for_image_load(pe_scheme, header)); - INFO("Phase encoding information read from JSON file modified to conform to prior MRtrix3 internal transform " - "realignment of image \"" + - header.name() + "\""); - } else { - INFO("Phase encoding information read from JSON file not modified"); - } - } - - auto slice_encoding_it = header.keyval().find("SliceEncodingDirection"); - if (slice_encoding_it != header.keyval().end()) { - if (do_realign) { - 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); - INFO("Slice encoding direction read from JSON file modified to conform to prior MRtrix3 internal transform " - "realignment of input image"); - } else { - INFO("Slice encoding information read from JSON file not modified"); - } - } +void read(const nlohmann::json &json, Header &header) { + KeyValues keyval = read(json); + // Reorientation based on image load should be applied + // exclusively to metadata loaded via JSON; not anything pre-existing + Metadata::PhaseEncoding::transform_for_image_load(keyval, header); + Metadata::SliceEncoding::transform_for_image_load(keyval, header); + header.merge_keyval(keyval); } namespace { @@ -251,36 +215,10 @@ void write(const Header &header, nlohmann::json &json, const std::string &image_ write(H_adj.keyval(), json); return; } - - std::vector order; - std::array flip; - File::NIfTI::axes_on_write(header, order, flip); - if (order[0] == 0 && order[1] == 1 && order[2] == 2 && !flip[0] && !flip[1] && !flip[2]) { - INFO( - "No need to transform orientation-based information written to JSON file to match image: image is already RAS"); - write(H_adj.keyval(), json); - return; - } - - auto pe_scheme = PhaseEncoding::get_scheme(header); - if (pe_scheme.rows()) { - // Assume that image being written to disk is going to have its transform adjusted, - // so modify the phase encoding scheme appropriately before writing to JSON - PhaseEncoding::set_scheme(H_adj, PhaseEncoding::transform_for_nifti_write(pe_scheme, header)); - INFO("Phase encoding information written to JSON file modified according to expected output NIfTI header transform " - "realignment"); - } - auto slice_encoding_it = H_adj.keyval().find("SliceEncodingDirection"); - if (slice_encoding_it != H_adj.keyval().end()) { - 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); - INFO("Slice encoding direction written to JSON file modified according to expected output NIfTI header transform " - "realignment"); - } - + if (App::get_options("export_grad_fsl").size()) + DWI::clear_DW_scheme(H_adj); + Metadata::PhaseEncoding::transform_for_nifti_write(H_adj.keyval(), H_adj); + Metadata::SliceEncoding::transform_for_nifti_write(H_adj.keyval(), H_adj); write(H_adj.keyval(), json); } diff --git a/cpp/core/file/json_utils.h b/cpp/core/file/json_utils.h index bf45eb2a10..be8ffd782f 100644 --- a/cpp/core/file/json_utils.h +++ b/cpp/core/file/json_utils.h @@ -27,8 +27,8 @@ namespace File::JSON { void load(Header &H, const std::string &path); void save(const Header &H, const std::string &json_path, const std::string &image_path); -KeyValues read(const nlohmann::json &json, const KeyValues &preexisting = KeyValues()); -void read(const nlohmann::json &json, Header &header, const bool realign); +KeyValues read(const nlohmann::json &json); +void read(const nlohmann::json &json, Header &header); void write(const KeyValues &keyval, nlohmann::json &json); void write(const Header &header, nlohmann::json &json, const std::string &image_path); diff --git a/cpp/core/file/mgh.h b/cpp/core/file/mgh.h index 4981d3e554..ff87cb2d4d 100644 --- a/cpp/core/file/mgh.h +++ b/cpp/core/file/mgh.h @@ -542,7 +542,7 @@ template void write_header(const Header &H, Output &out) { if (ndim > 4) throw Exception("MGH file format does not support images of more than 4 dimensions"); - std::vector axes; + Axes::permutations_type axes; auto M = File::NIfTI::adjust_transform(H, axes); store(1, out); // version diff --git a/cpp/core/file/nifti_utils.cpp b/cpp/core/file/nifti_utils.cpp index af7cac8d56..2cdffecc9d 100644 --- a/cpp/core/file/nifti_utils.cpp +++ b/cpp/core/file/nifti_utils.cpp @@ -358,7 +358,7 @@ template void store(NiftiHeader &NH, const Header &H, const bool is_BE = H.datatype().is_big_endian(); - std::vector axes; + Axes::permutations_type axes; auto M = File::NIfTI::adjust_transform(H, axes); memset(&NH, 0, sizeof(NH)); @@ -545,18 +545,21 @@ template void store(NiftiHeader &NH, const Header &H, const } } -void axes_on_write(const Header &H, std::vector &order, std::array &flip) { +Axes::Shuffle axes_on_write(const Header &H) { Stride::List strides = Stride::get(H); strides.resize(3); - order = Stride::order(strides); - flip = {strides[order[0]] < 0, strides[order[1]] < 0, strides[order[2]] < 0}; + auto order = Stride::order(strides); + Axes::Shuffle result; + result.permutations = {order[0], order[1], order[2]}; + result.flips = {strides[order[0]] < 0, strides[order[1]] < 0, strides[order[2]] < 0}; + return result; } -transform_type adjust_transform(const Header &H, std::vector &axes) { - std::array flip; - axes_on_write(H, axes, flip); +transform_type adjust_transform(const Header &H, Axes::permutations_type &axes) { + const Axes::Shuffle shuffle = axes_on_write(H); + axes = shuffle.permutations; - if (axes[0] == 0 && axes[1] == 1 && axes[2] == 2 && !flip[0] && !flip[1] && !flip[2]) + if (!shuffle) return H.transform(); const auto &M_in = H.transform().matrix(); @@ -568,7 +571,7 @@ transform_type adjust_transform(const Header &H, std::vector &axes) { auto translation = M_out.col(3); for (size_t i = 0; i < 3; ++i) { - if (flip[i]) { + if (shuffle.flips[i]) { auto length = default_type(H.size(axes[i]) - 1) * H.spacing(axes[i]); auto axis = M_out.col(i); axis = -axis; diff --git a/cpp/core/file/nifti_utils.h b/cpp/core/file/nifti_utils.h index d094b6539f..b57da8f541 100644 --- a/cpp/core/file/nifti_utils.h +++ b/cpp/core/file/nifti_utils.h @@ -29,8 +29,8 @@ class Header; namespace File::NIfTI { extern bool right_left_warning_issued; -void axes_on_write(const Header &H, std::vector &order, std::array &flip); -transform_type adjust_transform(const Header &H, std::vector &order); +Axes::Shuffle axes_on_write(const Header &H); +transform_type adjust_transform(const Header &H, Axes::permutations_type &order); bool check(int VERSION, Header &H, const size_t num_axes, const std::vector &suffixes); diff --git a/cpp/core/header.cpp b/cpp/core/header.cpp index 3aa1e0489d..45154e037e 100644 --- a/cpp/core/header.cpp +++ b/cpp/core/header.cpp @@ -26,8 +26,9 @@ #include "formats/list.h" #include "image_io/default.h" #include "image_io/scratch.h" +#include "metadata/phase_encoding.h" +#include "metadata/slice_encoding.h" #include "mrtrix.h" -#include "phase_encoding.h" #include "stride.h" #include "transform.h" @@ -63,38 +64,7 @@ void Header::check(const Header &H) const { throw Exception("scaling coefficients differ between image files for \"" + name() + "\""); } -namespace { -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, ","); - if (one_split.size() != two_split.size()) { - DEBUG("Slice timing vectors of inequal length"); - return "invalid"; - } - // Siemens CSA reports with 2.5ms precision = 0.0025s - // Allow slice times to vary by 1.5x this amount, but no more - for (size_t i = 0; i != one_split.size(); ++i) { - default_type f_one, f_two; - try { - f_one = to(one_split[i]); - f_two = to(two_split[i]); - } catch (Exception &e) { - DEBUG("Error converting slice timing vector to floating-point"); - return "invalid"; - } - const default_type diff = abs(f_two - f_one); - if (diff > 0.00375) { - DEBUG("Supra-threshold difference of " + str(diff) + "s in slice times"); - return "variable"; - } - } - return one; -} -} // namespace - -void Header::merge_keyval(const Header &H) { +void Header::merge_keyval(const KeyValues &in) { std::map new_keyval; std::set unique_comments; for (const auto &item : keyval()) { @@ -107,7 +77,7 @@ void Header::merge_keyval(const Header &H) { new_keyval.insert(item); } } - for (const auto &item : H.keyval()) { + for (const auto &item : in) { if (item.first == "comments") { const auto comments = split_lines(item.second); for (const auto &c : comments) { @@ -121,7 +91,7 @@ void Header::merge_keyval(const Header &H) { if (it == keyval().end() || it->second == item.second) new_keyval.insert(item); else if (item.first == "SliceTiming") - new_keyval["SliceTiming"] = resolve_slice_timing(item.second, it->second); + new_keyval["SliceTiming"] = Metadata::SliceEncoding::resolve_slice_timing(item.second, it->second); else new_keyval[item.first] = "variable"; } @@ -247,8 +217,7 @@ Header Header::open(const std::string &image_name) { } // End branching for [] notation H.sanitise(); - if (do_realign_transform) - H.realign_transform(); + H.realign_transform(); } catch (CancelException &e) { throw; } catch (Exception &E) { @@ -349,9 +318,9 @@ Header Header::create(const std::string &image_name, const Header &template_head DWI::clear_DW_scheme(H); } try { - pe_scheme = PhaseEncoding::parse_scheme(template_header); + pe_scheme = Metadata::PhaseEncoding::parse_scheme(template_header.keyval(), template_header); } catch (Exception &) { - PhaseEncoding::clear_scheme(H); + Metadata::PhaseEncoding::clear_scheme(H.keyval()); } if (split_4d_schemes) { try { @@ -362,11 +331,11 @@ Header Header::create(const std::string &image_name, const Header &template_head DWI::clear_DW_scheme(H); } try { - PhaseEncoding::check(pe_scheme, template_header); - PhaseEncoding::set_scheme(H, pe_scheme.row(0)); + Metadata::PhaseEncoding::check(pe_scheme, template_header); + Metadata::PhaseEncoding::set_scheme(H.keyval(), pe_scheme.row(0)); } catch (Exception &) { pe_scheme.resize(0, 0); - PhaseEncoding::clear_scheme(H); + Metadata::PhaseEncoding::clear_scheme(H.keyval()); } } @@ -400,7 +369,7 @@ Header Header::create(const std::string &image_name, const Header &template_head if (dw_scheme.rows()) DWI::set_DW_scheme(header, dw_scheme.row(counter)); if (pe_scheme.rows()) - PhaseEncoding::set_scheme(header, pe_scheme.row(counter)); + Metadata::PhaseEncoding::set_scheme(header.keyval(), pe_scheme.row(counter)); } std::shared_ptr io_handler((*format_handler)->create(header)); assert(io_handler); @@ -431,7 +400,7 @@ Header Header::create(const std::string &image_name, const Header &template_head if (split_4d_schemes) { DWI::set_DW_scheme(H, dw_scheme); - PhaseEncoding::set_scheme(H, pe_scheme); + Metadata::PhaseEncoding::set_scheme(H.keyval(), pe_scheme); } H.io->set_image_is_new(true); H.io->set_readwrite(true); @@ -610,12 +579,19 @@ void Header::sanitise_transform() { } void Header::realign_transform() { + realignment_.orig_transform_ = transform(); + realignment_.applied_transform_ = Realignment::applied_transform_type::Identity(); + realignment_.orig_strides_ = Stride::get(*this); + realignment_.orig_keyval_ = keyval(); + + if (!do_realign_transform) + return; + // find which row of the transform is closest to each scanner axis: - Axes::get_shuffle_to_make_axial(transform(), realign_perm_, realign_flip_); + realignment_.shuffle_ = Axes::get_shuffle_to_make_RAS(transform()); // check if image is already near-axial, return if true: - if (realign_perm_[0] == 0 && realign_perm_[1] == 1 && realign_perm_[2] == 2 && !realign_flip_[0] && - !realign_flip_[1] && !realign_flip_[2]) + if (!realignment_) return; auto M(transform()); @@ -623,22 +599,30 @@ void Header::realign_transform() { // modify translation vector: for (size_t i = 0; i < 3; ++i) { - if (realign_flip_[i]) { + if (realignment_.flip(i)) { const default_type length = (size(i) - 1) * spacing(i); auto axis = M.matrix().col(i); for (size_t n = 0; n < 3; ++n) { axis[n] = -axis[n]; translation[n] -= length * axis[n]; } + realignment_.applied_transform_.row(i) *= -1.0; } } // switch and/or invert rows if needed: for (size_t i = 0; i < 3; ++i) { - auto row = M.matrix().row(i).head<3>(); - row = Eigen::RowVector3d(row[realign_perm_[0]], row[realign_perm_[1]], row[realign_perm_[2]]); + auto row_transform = M.matrix().row(i).head<3>(); + row_transform = Eigen::RowVector3d(row_transform[realignment_.permutation(0)], + row_transform[realignment_.permutation(1)], + row_transform[realignment_.permutation(2)]); + + auto col_applied = realignment_.applied_transform_.matrix().col(i); + col_applied = Eigen::RowVector3i(col_applied[realignment_.permutation(0)], + col_applied[realignment_.permutation(1)], + col_applied[realignment_.permutation(2)]); - if (realign_flip_[i]) + if (realignment_.flip(i)) stride(i) = -stride(i); } @@ -646,42 +630,17 @@ void Header::realign_transform() { transform() = std::move(M); // switch axes to match: - Axis a[] = {axes_[realign_perm_[0]], axes_[realign_perm_[1]], axes_[realign_perm_[2]]}; + Axis a[] = {axes_[realignment_.permutation(0)], // + axes_[realignment_.permutation(1)], // + axes_[realignment_.permutation(2)]}; // axes_[0] = a[0]; axes_[1] = a[1]; axes_[2] = a[2]; INFO("Axes and transform of image \"" + name() + "\" altered to approximate RAS coordinate system"); - // If there's any phase encoding direction information present in the - // header, it's necessary here to update it according to the - // flips / permutations that have taken place - auto pe_scheme = PhaseEncoding::get_scheme(*this); - if (pe_scheme.rows()) { - for (ssize_t row = 0; row != pe_scheme.rows(); ++row) { - Eigen::VectorXd new_line(pe_scheme.row(row)); - for (ssize_t axis = 0; axis != 3; ++axis) { - new_line[axis] = pe_scheme(row, realign_perm_[axis]); - if (new_line[axis] && realign_flip_[realign_perm_[axis]]) - new_line[axis] = -new_line[axis]; - } - pe_scheme.row(row) = new_line; - } - PhaseEncoding::set_scheme(*this, pe_scheme); - INFO("Phase encoding scheme modified to conform to MRtrix3 internal header transform realignment"); - } - - // If there's any slice encoding direction information present in the - // header, that's also necessary to update here - auto slice_encoding_it = keyval().find("SliceEncodingDirection"); - if (slice_encoding_it != keyval().end()) { - 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); - INFO("Slice encoding direction has been modified to conform to MRtrix3 internal header transform realignment"); - } + Metadata::PhaseEncoding::transform_for_image_load(keyval(), *this); + Metadata::SliceEncoding::transform_for_image_load(keyval(), *this); } Header @@ -752,7 +711,7 @@ concatenate(const std::vector
&headers, const size_t axis_to_concat, con } catch (Exception &) { } try { - pe_scheme = PhaseEncoding::get_scheme(result); + pe_scheme = Metadata::PhaseEncoding::get_scheme(result); } catch (Exception &) { } } @@ -781,7 +740,7 @@ concatenate(const std::vector
&headers, const size_t axis_to_concat, con dw_scheme.resize(0, 0); } try { - const auto extra_pe = PhaseEncoding::get_scheme(H); + const auto extra_pe = Metadata::PhaseEncoding::get_scheme(H); concat_scheme(pe_scheme, extra_pe); } catch (Exception &) { pe_scheme.resize(0, 0); @@ -789,7 +748,7 @@ concatenate(const std::vector
&headers, const size_t axis_to_concat, con } // Resolve key-value pairs - result.merge_keyval(H); + result.merge_keyval(H.keyval()); // Resolve discrepancies in datatype; // also throw an exception if such mismatch is not permitted @@ -805,9 +764,13 @@ concatenate(const std::vector
&headers, const size_t axis_to_concat, con if (axis_to_concat == 3) { DWI::set_DW_scheme(result, dw_scheme); - PhaseEncoding::set_scheme(result, pe_scheme); + Metadata::PhaseEncoding::set_scheme(result.keyval(), pe_scheme); } return result; } +Header::Realignment::Realignment() : applied_transform_(applied_transform_type::Identity()), orig_keyval_() { + orig_transform_.matrix().fill(std::numeric_limits::quiet_NaN()); +} + } // namespace MR diff --git a/cpp/core/header.h b/cpp/core/header.h index 99ea597a3f..b182196eb7 100644 --- a/cpp/core/header.h +++ b/cpp/core/header.h @@ -20,6 +20,7 @@ #include #include "app.h" +#include "axes.h" #include "datatype.h" #include "debug.h" #include "file/mmap.h" @@ -53,12 +54,7 @@ class Header { }; Header() - : transform_(Eigen::Matrix::Constant(NaN)), - format_(nullptr), - offset_(0.0), - scale_(1.0), - realign_perm_{{0, 1, 2}}, - realign_flip_{{false, false, false}} {} + : transform_(Eigen::Matrix::Constant(NaN)), format_(nullptr), offset_(0.0), scale_(1.0) {} explicit Header(Header &&H) noexcept : axes_(std::move(H.axes_)), @@ -70,8 +66,7 @@ class Header { datatype_(std::move(H.datatype_)), offset_(H.offset_), scale_(H.scale_), - realign_perm_{{0, 1, 2}}, - realign_flip_{{false, false, false}} {} + realignment_(H.realignment_) {} Header &operator=(Header &&H) noexcept { axes_ = std::move(H.axes_); @@ -83,8 +78,7 @@ class Header { datatype_ = std::move(H.datatype_); offset_ = H.offset_; scale_ = H.scale_; - realign_perm_ = H.realign_perm_; - realign_flip_ = H.realign_flip_; + realignment_ = H.realignment_; return *this; } @@ -100,8 +94,7 @@ class Header { datatype_(H.datatype_), offset_(datatype().is_integer() ? H.offset_ : 0.0), scale_(datatype().is_integer() ? H.scale_ : 1.0), - realign_perm_(H.realign_perm_), - realign_flip_(H.realign_flip_) {} + realignment_(H.realignment_) {} //! copy constructor from type of class derived from Header /*! This invokes the standard Header(const Header&) copy-constructor. */ @@ -120,9 +113,7 @@ class Header { format_(nullptr), datatype_(DataType::from()), offset_(0.0), - scale_(1.0), - realign_perm_{{0, 1, 2}}, - realign_flip_{{false, false, false}} { + scale_(1.0) { axes_.resize(original.ndim()); for (size_t n = 0; n < original.ndim(); ++n) { size(n) = original.size(n); @@ -143,8 +134,7 @@ class Header { datatype_ = H.datatype_; offset_ = datatype().is_integer() ? H.offset_ : 0.0; scale_ = datatype().is_integer() ? H.scale_ : 1.0; - realign_perm_ = H.realign_perm_; - realign_flip_ = H.realign_flip_; + realignment_ = H.realignment_; io.reset(); return *this; } @@ -175,8 +165,7 @@ class Header { datatype_ = DataType::from(); offset_ = 0.0; scale_ = 1.0; - realign_perm_ = {{0, 1, 2}}; - realign_flip_ = {{false, false, false}}; + realignment_ = Realignment(); io.reset(); return *this; } @@ -207,11 +196,41 @@ class Header { //! get/set the 4x4 affine transformation matrix mapping image to world coordinates transform_type &transform() { return transform_; } + // Class to store all information relating to internal transform realignment + class 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 + using applied_transform_type = Eigen::Matrix; + Realignment(); + Realignment(Header &); + operator bool() const { return bool(shuffle_); } + const std::array &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; } + 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_; } + KeyValues &orig_keyval() { return orig_keyval_; } + + private: + Axes::Shuffle shuffle_; + transform_type orig_transform_; + Stride::List orig_strides_; + applied_transform_type applied_transform_; + KeyValues orig_keyval_; + friend class Header; + }; //! get information on how the transform was modified on image load - void realignment(std::array &perm, std::array &flip) const { - perm = realign_perm_; - flip = realign_flip_; - } + const Realignment &realignment() const { return realignment_; } class NDimProxy { public: @@ -359,8 +378,8 @@ class Header { const KeyValues &keyval() const { return keyval_; } //! get/set generic key/value text attributes KeyValues &keyval() { return keyval_; } - //! merge key/value entries from another header - void merge_keyval(const Header &H); + //! merge key/value entries from another dictionary + void merge_keyval(const KeyValues &); static Header open(const std::string &image_name); static Header @@ -397,8 +416,7 @@ class Header { void realign_transform(); /*! store information about how image was * realigned via realign_transform(). */ - std::array realign_perm_; - std::array realign_flip_; + Realignment realignment_; void sanitise_voxel_sizes(); void sanitise_transform(); diff --git a/cpp/core/metadata/bids.cpp b/cpp/core/metadata/bids.cpp new file mode 100644 index 0000000000..f5cbeec2eb --- /dev/null +++ b/cpp/core/metadata/bids.cpp @@ -0,0 +1,71 @@ +/* Copyright (c) 2008-2024 the MRtrix3 contributors. + * + * This Source Code Form is subject to the terms of the Mozilla Public + * License, v. 2.0. If a copy of the MPL was not distributed with this + * file, You can obtain one at http://mozilla.org/MPL/2.0/. + * + * Covered Software is provided under this License on an "as is" + * basis, without warranty of any kind, either expressed, implied, or + * statutory, including, without limitation, warranties that the + * Covered Software is free of defects, merchantable, fit for a + * particular purpose or non-infringing. + * See the Mozilla Public License v. 2.0 for more details. + * + * For more details, see http://www.mrtrix.org/. + */ + +#include "metadata/bids.h" + +#include "exception.h" +#include "mrtrix.h" + +namespace MR::Metadata::BIDS { + +std::string vector2axisid(const axis_vector_type &dir) { + if (dir[0] == -1) { + assert(!dir[1]); + assert(!dir[2]); + return "i-"; + } else if (dir[0] == 1) { + assert(!dir[1]); + assert(!dir[2]); + return "i"; + } else if (dir[1] == -1) { + assert(!dir[0]); + assert(!dir[2]); + return "j-"; + } else if (dir[1] == 1) { + assert(!dir[0]); + assert(!dir[2]); + return "j"; + } else if (dir[2] == -1) { + assert(!dir[0]); + assert(!dir[1]); + return "k-"; + } else if (dir[2] == 1) { + assert(!dir[0]); + assert(!dir[1]); + return "k"; + } else { + throw Exception("Malformed image axis vector: \"" + str(dir.transpose()) + "\""); + } +} + +axis_vector_type axisid2vector(const std::string &id) { + if (id == "i-") + return {-1, 0, 0}; + else if (id == "i") + return {1, 0, 0}; + else if (id == "j-") + return {0, -1, 0}; + else if (id == "j") + return {0, 1, 0}; + else if (id == "k-") + return {0, 0, -1}; + else if (id == "k") + return {0, 0, 1}; + else + throw Exception("Malformed image axis identifier: \"" + id + "\""); +} + +} // namespace MR::Metadata::BIDS diff --git a/cpp/core/metadata/bids.h b/cpp/core/metadata/bids.h new file mode 100644 index 0000000000..c07eac82c8 --- /dev/null +++ b/cpp/core/metadata/bids.h @@ -0,0 +1,33 @@ +/* Copyright (c) 2008-2024 the MRtrix3 contributors. + * + * This Source Code Form is subject to the terms of the Mozilla Public + * License, v. 2.0. If a copy of the MPL was not distributed with this + * file, You can obtain one at http://mozilla.org/MPL/2.0/. + * + * Covered Software is provided under this License on an "as is" + * basis, without warranty of any kind, either expressed, implied, or + * statutory, including, without limitation, warranties that the + * Covered Software is free of defects, merchantable, fit for a + * particular purpose or non-infringing. + * See the Mozilla Public License v. 2.0 for more details. + * + * For more details, see http://www.mrtrix.org/. + */ + +#pragma once + +#include +#include + +namespace MR::Metadata::BIDS { + +//! convert axis directions between formats +/*! these helper functions convert the definition of + * phase-encoding direction between a 3-vector (e.g. + * [0 1 0] ) and a BIDS NIfTI axis identifier (e.g. 'i-') + */ +using axis_vector_type = Eigen::Matrix; +std::string vector2axisid(const axis_vector_type &); +axis_vector_type axisid2vector(const std::string &); + +} // namespace MR::Metadata::BIDS diff --git a/cpp/core/metadata/phase_encoding.cpp b/cpp/core/metadata/phase_encoding.cpp new file mode 100644 index 0000000000..f77616083d --- /dev/null +++ b/cpp/core/metadata/phase_encoding.cpp @@ -0,0 +1,323 @@ +/* Copyright (c) 2008-2024 the MRtrix3 contributors. + * + * This Source Code Form is subject to the terms of the Mozilla Public + * License, v. 2.0. If a copy of the MPL was not distributed with this + * file, You can obtain one at http://mozilla.org/MPL/2.0/. + * + * Covered Software is provided under this License on an "as is" + * basis, without warranty of any kind, either expressed, implied, or + * statutory, including, without limitation, warranties that the + * Covered Software is free of defects, merchantable, fit for a + * particular purpose or non-infringing. + * See the Mozilla Public License v. 2.0 for more details. + * + * For more details, see http://www.mrtrix.org/. + */ + +#include "metadata/phase_encoding.h" + +namespace MR::Metadata::PhaseEncoding { + +// clang-format off +using namespace App; +const OptionGroup ImportOptions = + OptionGroup("Options for importing phase-encode tables") + + Option("import_pe_table", "import a phase-encoding table from file") + + Argument("file").type_file_in() + + Option("import_pe_eddy", "import phase-encoding information from an EDDY-style config / index file pair") + + Argument("config").type_file_in() + + Argument("indices").type_file_in(); + +const OptionGroup SelectOptions = + OptionGroup("Options for selecting volumes based on phase-encoding") + + Option("pe", + "select volumes with a particular phase encoding;" + " this can be three comma-separated values" + " (for i,j,k components of vector direction)" + " or four (direction & total readout time)") + + Argument("desc").type_sequence_float(); + +const OptionGroup ExportOptions = + OptionGroup("Options for exporting phase-encode tables") + + Option("export_pe_table", "export phase-encoding table to file") + + Argument("file").type_file_out() + + Option("export_pe_eddy", "export phase-encoding information to an EDDY-style config / index file pair") + + Argument("config").type_file_out() + + Argument("indices").type_file_out(); +// clang-format on + +void check(const scheme_type &PE) { + if (!PE.rows()) + throw Exception("No valid phase encoding table found"); + + if (PE.cols() < 3) + throw Exception("Phase-encoding matrix must have at least 3 columns"); + + for (ssize_t row = 0; row != PE.rows(); ++row) { + for (ssize_t axis = 0; axis != 3; ++axis) { + if (std::round(PE(row, axis)) != PE(row, axis)) + throw Exception("Phase-encoding matrix contains non-integral axis designation"); + } + } +} + +void check(const scheme_type &PE, const Header &header) { + check(PE); + const ssize_t num_volumes = (header.ndim() > 3) ? header.size(3) : 1; + if (num_volumes != PE.rows()) + throw Exception("Number of volumes in image \"" + header.name() + "\" (" + str(num_volumes) + + ") does not match that in phase encoding table (" + str(PE.rows()) + ")"); +} + +namespace { +void erase(KeyValues &keyval, const std::string &s) { + auto it = keyval.find(s); + if (it != keyval.end()) + keyval.erase(it); +}; +} // namespace +void set_scheme(KeyValues &keyval, const scheme_type &PE) { + if (!PE.rows()) { + erase(keyval, "pe_scheme"); + erase(keyval, "PhaseEncodingDirection"); + erase(keyval, "TotalReadoutTime"); + return; + } + std::string pe_scheme; + std::string first_line; + bool variation = false; + for (ssize_t row = 0; row < PE.rows(); ++row) { + std::string line = str(PE(row, 0)); + for (ssize_t col = 1; col < PE.cols(); ++col) + line += "," + str(PE(row, col), 3); + add_line(pe_scheme, line); + if (first_line.empty()) + first_line = line; + else if (line != first_line) + variation = true; + } + if (variation) { + keyval["pe_scheme"] = pe_scheme; + erase(keyval, "PhaseEncodingDirection"); + erase(keyval, "TotalReadoutTime"); + } else { + erase(keyval, "pe_scheme"); + const Metadata::BIDS::axis_vector_type dir{int(PE(0, 0)), int(PE(0, 1)), int(PE(0, 2))}; + keyval["PhaseEncodingDirection"] = Metadata::BIDS::vector2axisid(dir); + if (PE.cols() >= 4) + keyval["TotalReadoutTime"] = str(PE(0, 3), 3); + else + erase(keyval, "TotalReadoutTime"); + } +} + +void clear_scheme(KeyValues &keyval) { + auto erase = [&](const std::string &s) { + auto it = keyval.find(s); + if (it != keyval.end()) + keyval.erase(it); + }; + erase("pe_scheme"); + erase("PhaseEncodingDirection"); + erase("TotalReadoutTime"); +} + +scheme_type parse_scheme(const KeyValues &keyval, const Header &header) { + scheme_type PE; + const auto it = keyval.find("pe_scheme"); + if (it != keyval.end()) { + try { + PE = MR::parse_matrix(it->second); + } catch (Exception &e) { + throw Exception(e, "malformed PE scheme associated with image \"" + header.name() + "\""); + } + if (ssize_t(PE.rows()) != ((header.ndim() > 3) ? header.size(3) : 1)) + throw Exception("malformed PE scheme associated with image \"" + header.name() + "\":" + // + " number of rows does not equal number of volumes"); + } else { + const auto it_dir = keyval.find("PhaseEncodingDirection"); + if (it_dir != keyval.end()) { + const auto it_time = keyval.find("TotalReadoutTime"); + const size_t cols = it_time == keyval.end() ? 3 : 4; + Eigen::Matrix row(cols); + row.head(3) = BIDS::axisid2vector(it_dir->second).cast(); + if (it_time != keyval.end()) + row[3] = to(it_time->second); + PE.resize((header.ndim() > 3) ? header.size(3) : 1, cols); + PE.rowwise() = row.transpose(); + } + } + return PE; +} + +scheme_type get_scheme(const Header &header) { + DEBUG("searching for suitable phase encoding data..."); + using namespace App; + scheme_type result; + + try { + const auto opt_table = get_options("import_pe_table"); + if (!opt_table.empty()) + result = load(opt_table[0][0], header); + const auto opt_eddy = get_options("import_pe_eddy"); + if (!opt_eddy.empty()) { + if (!opt_table.empty()) + throw Exception("Phase encoding table can be provided" + " using either -import_pe_table or -import_pe_eddy option," + " but NOT both"); + result = load_eddy(opt_eddy[0][0], opt_eddy[0][1], header); + } + if (opt_table.empty() && opt_eddy.empty()) + result = parse_scheme(header.keyval(), header); + } catch (Exception &e) { + throw Exception(e, "error importing phase encoding table for image \"" + header.name() + "\""); + } + + if (!result.rows()) + return result; + + if (result.cols() < 3) + throw Exception("unexpected phase encoding table matrix dimensions"); + + INFO("found " + str(result.rows()) + "x" + str(result.cols()) + " phase encoding table"); + + return result; +} + +void transform_for_image_load(KeyValues &keyval, const Header &H) { + scheme_type pe_scheme = parse_scheme(keyval, H); + if (!pe_scheme.rows()) { + DEBUG("No phase encoding information found for transformation with load of image \"" + H.name() + "\""); + return; + } + if (!H.realignment()) { + INFO("No transformation of phase encoding data for load of image \"" + H.name() + "\" required"); + return; + } + set_scheme(keyval, transform_for_image_load(pe_scheme, H)); + INFO("Phase encoding data transformed to match RAS realignment of image \"" + H.name() + "\""); +} + +scheme_type transform_for_image_load(const scheme_type &pe_scheme, const Header &H) { + if (!H.realignment()) + return pe_scheme; + scheme_type result(pe_scheme.rows(), pe_scheme.cols()); + for (ssize_t row = 0; row != pe_scheme.rows(); ++row) { + Eigen::VectorXd new_line = pe_scheme.row(row); + new_line.head<3>() = (H.realignment().applied_transform() * new_line.head<3>().cast()).cast(); + result.row(row) = new_line; + } + return result; +} + +void transform_for_nifti_write(KeyValues &keyval, const Header &H) { + scheme_type pe_scheme = parse_scheme(keyval, H); + if (!pe_scheme.rows()) { + DEBUG("No phase encoding information found for transformation with save of NIfTI image \"" + H.name() + "\""); + return; + } + set_scheme(keyval, transform_for_nifti_write(pe_scheme, H)); +} + +scheme_type transform_for_nifti_write(const scheme_type &pe_scheme, const Header &H) { + if (!pe_scheme.rows()) + return pe_scheme; + Axes::Shuffle shuffle = File::NIfTI::axes_on_write(H); + if (!shuffle) { + INFO("No transformation of phase encoding data required for export to file:" + " output image will be RAS"); + return pe_scheme; + } + scheme_type result(pe_scheme.rows(), pe_scheme.cols()); + for (ssize_t row = 0; row != pe_scheme.rows(); ++row) { + Eigen::VectorXd new_line = pe_scheme.row(row); + for (ssize_t axis = 0; axis != 3; ++axis) + new_line[axis] = // + pe_scheme(row, shuffle.permutations[axis]) && shuffle.flips[axis] // + ? -pe_scheme(row, shuffle.permutations[axis]) // + : pe_scheme(row, shuffle.permutations[axis]); // + result.row(row) = new_line; + } + INFO("Phase encoding data transformed to match NIfTI / MGH image export prior to writing to file"); + return result; +} + +void scheme2eddy(const scheme_type &PE, Eigen::MatrixXd &config, Eigen::Array &indices) { + try { + check(PE); + } catch (Exception &e) { + throw Exception(e, "Cannot convert phase-encoding scheme to eddy format"); + } + if (PE.cols() != 4) + throw Exception("Phase-encoding matrix requires 4 columns to convert to eddy format"); + config.resize(0, 0); + indices = Eigen::Array::Constant(PE.rows(), PE.rows()); + for (ssize_t PE_row = 0; PE_row != PE.rows(); ++PE_row) { + for (ssize_t config_row = 0; config_row != config.rows(); ++config_row) { + bool dir_match = PE.template block<1, 3>(PE_row, 0).isApprox(config.block<1, 3>(config_row, 0)); + bool time_match = abs(PE(PE_row, 3) - config(config_row, 3)) < 1e-3; + if (dir_match && time_match) { + // FSL-style index file indexes from 1 + indices[PE_row] = config_row + 1; + break; + } + } + if (indices[PE_row] == PE.rows()) { + // No corresponding match found in config matrix; create a new entry + config.conservativeResize(config.rows() + 1, 4); + config.row(config.rows() - 1) = PE.row(PE_row); + indices[PE_row] = config.rows(); + } + } +} + +scheme_type eddy2scheme(const Eigen::MatrixXd &config, const Eigen::Array &indices) { + if (config.cols() != 4) + throw Exception("Expected 4 columns in EDDY-format phase-encoding config file"); + scheme_type result(indices.size(), 4); + for (ssize_t row = 0; row != indices.size(); ++row) { + if (indices[row] > config.rows()) + throw Exception("Malformed EDDY-style phase-encoding information:" + " index exceeds number of config entries"); + result.row(row) = config.row(indices[row] - 1); + } + return result; +} + +void export_commandline(const Header &header) { + auto check = [&](const scheme_type &m) -> const scheme_type & { + if (!m.rows()) + throw Exception("no phase-encoding information found within image \"" + header.name() + "\""); + return m; + }; + + auto scheme = parse_scheme(header.keyval(), header); + + auto opt = get_options("export_pe_table"); + if (!opt.empty()) + save(check(scheme), header, opt[0][0]); + + opt = get_options("export_pe_eddy"); + if (!opt.empty()) + save_eddy(check(scheme), header, opt[0][0], opt[0][1]); +} + +scheme_type load(const std::string &path, const Header &header) { + const scheme_type PE = File::Matrix::load_matrix(path); + check(PE, header); + // As with JSON import, need to query the header to discover if the + // strides / transform were modified on image load to make the image + // data appear approximately axial, in which case we need to apply the + // same transforms to the phase encoding data on load + return transform_for_image_load(PE, header); +} + +scheme_type load_eddy(const std::string &config_path, const std::string &index_path, const Header &header) { + const Eigen::MatrixXd config = File::Matrix::load_matrix(config_path); + const Eigen::Array indices = File::Matrix::load_vector(index_path); + const scheme_type PE = eddy2scheme(config, indices); + check(PE, header); + return transform_for_image_load(PE, header); +} + +} // namespace MR::Metadata::PhaseEncoding diff --git a/cpp/core/metadata/phase_encoding.h b/cpp/core/metadata/phase_encoding.h new file mode 100644 index 0000000000..510cae46f1 --- /dev/null +++ b/cpp/core/metadata/phase_encoding.h @@ -0,0 +1,145 @@ +/* Copyright (c) 2008-2024 the MRtrix3 contributors. + * + * This Source Code Form is subject to the terms of the Mozilla Public + * License, v. 2.0. If a copy of the MPL was not distributed with this + * file, You can obtain one at http://mozilla.org/MPL/2.0/. + * + * Covered Software is provided under this License on an "as is" + * basis, without warranty of any kind, either expressed, implied, or + * statutory, including, without limitation, warranties that the + * Covered Software is free of defects, merchantable, fit for a + * particular purpose or non-infringing. + * See the Mozilla Public License v. 2.0 for more details. + * + * For more details, see http://www.mrtrix.org/. + */ + +#pragma once + +#include +#include + +#include "app.h" +#include "axes.h" +#include "file/matrix.h" +#include "file/nifti_utils.h" +#include "file/ofstream.h" +#include "header.h" +#include "metadata/bids.h" +#include "types.h" + +namespace MR::Metadata::PhaseEncoding { + +extern const App::OptionGroup ImportOptions; +extern const App::OptionGroup SelectOptions; +extern const App::OptionGroup ExportOptions; + +using scheme_type = Eigen::MatrixXd; + +//! check that a phase-encoding table is valid +void check(const scheme_type &PE); + +//! check that the PE scheme matches the DWI data in \a header +void check(const scheme_type &PE, const Header &header); + +//! store the phase encoding matrix in a header +/*! this will store the phase encoding matrix into the + * Header::keyval() structure of \a header. + * - If the phase encoding direction and/or total readout + * time varies between volumes, then the information + * will be stored in field "pe_scheme"; if not, it + * will instead be stored in fields "PhaseEncodingDirection" + * and "TotalReadoutTime" + */ +void set_scheme(KeyValues &keyval, const scheme_type &PE); + +//! clear the phase encoding matrix from a key-value dictionary +/*! this will delete any trace of phase encoding information + * from the dictionary. + */ +void clear_scheme(KeyValues &keyval); + +//! parse the phase encoding matrix from a key-value dictionary +/*! extract the phase encoding matrix stored in \a the key-value dictionary + * if one is present. The key-value dictionary is not in all use cases + * the "keyval" member of the Header class. + */ +Eigen::MatrixXd parse_scheme(const KeyValues &, const Header &); + +//! get a phase encoding matrix +/*! get a valid phase-encoding matrix, either from files specified at + * the command-line that exclusively provide phase encoding information + * (ie. NOT from .json; that is handled elsewhere), + * or from the contents of the image header. + */ +Eigen::MatrixXd get_scheme(const Header &); + +//! Convert a phase-encoding scheme into the EDDY config / indices format +void scheme2eddy(const scheme_type &PE, Eigen::MatrixXd &config, Eigen::Array &indices); + +//! Convert phase-encoding infor from the EDDY config / indices format into a standard scheme +scheme_type eddy2scheme(const Eigen::MatrixXd &, const Eigen::Array &); + +//! Modifies a phase encoding scheme if being imported alongside a non-RAS image +// and internal header realignment is performed +void transform_for_image_load(KeyValues &keyval, const Header &H); +scheme_type transform_for_image_load(const scheme_type &pe_scheme, const Header &H); + +//! Modifies a phase encoding scheme if being exported alongside a NIfTI image +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 + +//! Save a phase-encoding scheme 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 +template void save(const scheme_type &PE, const HeaderType &header, const std::string &path) { + try { + check(PE, header); + } catch (Exception &e) { + throw Exception(e, "Cannot export phase-encoding table to file \"" + path + "\""); + } + + if (Path::has_suffix(header.name(), {".mgh", ".mgz", ".nii", ".nii.gz", ".img"})) { + __save(transform_for_nifti_write(PE, header), path); + } else { + __save(PE, path); + } +} + +//! Save a phase-encoding scheme to EDDY format config / index files +template +void save_eddy(const scheme_type &PE, + const HeaderType &header, + const std::string &config_path, + const std::string &index_path) { + Eigen::MatrixXd config; + Eigen::Array indices; + scheme2eddy(transform_for_nifti_write(PE, header), config, indices); + File::Matrix::save_matrix(config, config_path, KeyValues(), false); + File::Matrix::save_vector(indices, index_path, KeyValues(), false); +} + +//! Save the phase-encoding scheme from a header to file depending on command-line options +void export_commandline(const Header &); + +//! Load a phase-encoding scheme from a matrix text file +scheme_type load(const std::string &path, const Header &header); + +//! Load a phase-encoding scheme from an EDDY-format config / indices file pair +scheme_type load_eddy(const std::string &config_path, const std::string &index_path, const Header &header); + +} // namespace MR::Metadata::PhaseEncoding diff --git a/cpp/core/metadata/slice_encoding.cpp b/cpp/core/metadata/slice_encoding.cpp new file mode 100644 index 0000000000..91ef1d2c37 --- /dev/null +++ b/cpp/core/metadata/slice_encoding.cpp @@ -0,0 +1,127 @@ +/* Copyright (c) 2008-2024 the MRtrix3 contributors. + * + * This Source Code Form is subject to the terms of the Mozilla Public + * License, v. 2.0. If a copy of the MPL was not distributed with this + * file, You can obtain one at http://mozilla.org/MPL/2.0/. + * + * Covered Software is provided under this License on an "as is" + * basis, without warranty of any kind, either expressed, implied, or + * statutory, including, without limitation, warranties that the + * Covered Software is free of defects, merchantable, fit for a + * particular purpose or non-infringing. + * See the Mozilla Public License v. 2.0 for more details. + * + * For more details, see http://www.mrtrix.org/. + */ + +#include "metadata/slice_encoding.h" + +#include "axes.h" +#include "file/nifti_utils.h" +#include "header.h" +#include "metadata/bids.h" + +namespace MR::Metadata::SliceEncoding { + +void transform_for_image_load(KeyValues &keyval, const Header &header) { + // If there's any slice encoding direction information present in the + // header, that's also necessary to update here + 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()) { + INFO("No transformation of slice encoding direction for load of image \"" + header.name() + "\" required"); + return; + } + const Metadata::BIDS::axis_vector_type orig_dir(slice_encoding_it == keyval.end() // + ? Metadata::BIDS::axis_vector_type({0, 0, 1}) // + : Metadata::BIDS::axisid2vector(slice_encoding_it->second)); // + const Metadata::BIDS::axis_vector_type new_dir = header.realignment().applied_transform() * orig_dir; + if (slice_encoding_it != keyval.end()) { + slice_encoding_it->second = Metadata::BIDS::vector2axisid(new_dir); + INFO("Slice encoding direction has been modified" // + " to conform to MRtrix3 internal header transform realignment of image \"" // + + header.name() + "\""); // + } else if ((new_dir * -1).dot(orig_dir) == 1) { + auto slice_timing = parse_floats(slice_timing_it->second); + std::reverse(slice_timing.begin(), slice_timing.end()); + slice_timing_it->second = join(slice_timing, ","); + INFO("Slice timing vector reversed to conform to MRtrix3 internal transform realignment" + " of image \"" + + header.name() + "\""); + } else { + keyval["SliceEncodingDirection"] = Metadata::BIDS::vector2axisid(new_dir); + WARN("Slice encoding direction of image \"" // + + header.name() // + + "\" inferred to be \"k\" in order to preserve interpretation of existing \"SliceTiming\" field" // + " after MRtrix3 internal transform realignment"); // + } + } +} + +void transform_for_nifti_write(KeyValues &keyval, const Header &H) { + 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()) + return; + + const Axes::Shuffle shuffle = File::NIfTI::axes_on_write(H); + if (!shuffle) { + INFO("No need to transform slice encoding information for NIfTI image write:" + " image is already RAS"); + return; + } + + const Metadata::BIDS::axis_vector_type orig_dir(slice_encoding_it == keyval.end() // + ? Metadata::BIDS::axis_vector_type({0, 0, 1}) // + : Metadata::BIDS::axisid2vector(slice_encoding_it->second)); // + Metadata::BIDS::axis_vector_type new_dir; + for (size_t axis = 0; axis != 3; ++axis) + new_dir[axis] = shuffle.flips[axis] ? -orig_dir[shuffle.permutations[axis]] : orig_dir[shuffle.permutations[axis]]; + + if (slice_encoding_it != keyval.end()) { + slice_encoding_it->second = Metadata::BIDS::vector2axisid(new_dir); + INFO("Slice encoding direction modified according to output NIfTI strides"); + } else if ((new_dir * -1).dot(orig_dir) == 1) { + auto slice_timing = parse_floats(slice_timing_it->second); + std::reverse(slice_timing.begin(), slice_timing.end()); + slice_timing_it->second = join(slice_timing, ","); + INFO("Slice timing vector reversed to conform to output NIfTI strides"); + } else { + keyval["SliceEncodingDirection"] = Metadata::BIDS::vector2axisid(new_dir); + WARN("Slice encoding direction added to metadata" + " in order to preserve interpretation of existing \"SliceTiming\" field" + " following output NIfTI strides"); + } +} + +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, ","); + if (one_split.size() != two_split.size()) { + DEBUG("Slice timing vectors of inequal length"); + return "invalid"; + } + // Siemens CSA reports with 2.5ms precision = 0.0025s + // Allow slice times to vary by 1.5x this amount, but no more + for (size_t i = 0; i != one_split.size(); ++i) { + default_type f_one, f_two; + try { + f_one = to(one_split[i]); + f_two = to(two_split[i]); + } catch (Exception &e) { + DEBUG("Error converting slice timing vector to floating-point"); + return "invalid"; + } + const default_type diff = abs(f_two - f_one); + if (diff > 0.00375) { + DEBUG("Supra-threshold difference of " + str(diff) + "s in slice times"); + return "variable"; + } + } + return one; +} + +} // namespace MR::Metadata::SliceEncoding diff --git a/cpp/core/metadata/slice_encoding.h b/cpp/core/metadata/slice_encoding.h new file mode 100644 index 0000000000..91999e03a5 --- /dev/null +++ b/cpp/core/metadata/slice_encoding.h @@ -0,0 +1,35 @@ +/* Copyright (c) 2008-2024 the MRtrix3 contributors. + * + * This Source Code Form is subject to the terms of the Mozilla Public + * License, v. 2.0. If a copy of the MPL was not distributed with this + * file, You can obtain one at http://mozilla.org/MPL/2.0/. + * + * Covered Software is provided under this License on an "as is" + * basis, without warranty of any kind, either expressed, implied, or + * statutory, including, without limitation, warranties that the + * Covered Software is free of defects, merchantable, fit for a + * particular purpose or non-infringing. + * See the Mozilla Public License v. 2.0 for more details. + * + * For more details, see http://www.mrtrix.org/. + */ + +#pragma once + +#include + +#include "types.h" + +namespace MR { +class Header; +} // namespace MR + +namespace MR::Metadata::SliceEncoding { + +void transform_for_image_load(KeyValues &keyval, const Header &H); + +void transform_for_nifti_write(KeyValues &keyval, const Header &H); + +std::string resolve_slice_timing(const std::string &one, const std::string &two); + +} // namespace MR::Metadata::SliceEncoding diff --git a/cpp/core/phase_encoding.cpp b/cpp/core/phase_encoding.cpp deleted file mode 100644 index ac7f7dfec8..0000000000 --- a/cpp/core/phase_encoding.cpp +++ /dev/null @@ -1,153 +0,0 @@ -/* Copyright (c) 2008-2025 the MRtrix3 contributors. - * - * This Source Code Form is subject to the terms of the Mozilla Public - * License, v. 2.0. If a copy of the MPL was not distributed with this - * file, You can obtain one at http://mozilla.org/MPL/2.0/. - * - * Covered Software is provided under this License on an "as is" - * basis, without warranty of any kind, either expressed, implied, or - * statutory, including, without limitation, warranties that the - * Covered Software is free of defects, merchantable, fit for a - * particular purpose or non-infringing. - * See the Mozilla Public License v. 2.0 for more details. - * - * For more details, see http://www.mrtrix.org/. - */ - -#include "phase_encoding.h" - -namespace MR::PhaseEncoding { - -// clang-format off -using namespace App; -const OptionGroup ImportOptions = - OptionGroup("Options for importing phase-encode tables") - + Option("import_pe_table", "import a phase-encoding table from file") - + Argument("file").type_file_in() - + Option("import_pe_eddy", "import phase-encoding information from an EDDY-style config / index file pair") - + Argument("config").type_file_in() - + Argument("indices").type_file_in(); - -const OptionGroup SelectOptions = - OptionGroup("Options for selecting volumes based on phase-encoding") - + Option("pe", - "select volumes with a particular phase encoding;" - " this can be three comma-separated values" - " (for i,j,k components of vector direction)" - " or four (direction & total readout time)") - + Argument("desc").type_sequence_float(); - -const OptionGroup ExportOptions = - OptionGroup("Options for exporting phase-encode tables") - + Option("export_pe_table", "export phase-encoding table to file") - + Argument("file").type_file_out() - + Option("export_pe_eddy", "export phase-encoding information to an EDDY-style config / index file pair") - + Argument("config").type_file_out() - + Argument("indices").type_file_out(); -// clang-format on - -void clear_scheme(Header &header) { - auto erase = [&](const std::string &s) { - auto it = header.keyval().find(s); - if (it != header.keyval().end()) - header.keyval().erase(it); - }; - erase("pe_scheme"); - erase("PhaseEncodingDirection"); - erase("TotalReadoutTime"); -} - -Eigen::MatrixXd parse_scheme(const Header &header) { - Eigen::MatrixXd PE; - const auto it = header.keyval().find("pe_scheme"); - if (it != header.keyval().end()) { - try { - PE = MR::parse_matrix(it->second); - } catch (Exception &e) { - throw Exception(e, "malformed PE scheme in image \"" + header.name() + "\""); - } - if (ssize_t(PE.rows()) != ((header.ndim() > 3) ? header.size(3) : 1)) - throw Exception("malformed PE scheme in image \"" + header.name() + "\":" + // - " number of rows does not equal number of volumes"); - } else { - const auto it_dir = header.keyval().find("PhaseEncodingDirection"); - if (it_dir != header.keyval().end()) { - const auto it_time = header.keyval().find("TotalReadoutTime"); - const size_t cols = it_time == header.keyval().end() ? 3 : 4; - Eigen::Matrix row(cols); - row.head(3) = Axes::id2dir(it_dir->second); - if (it_time != header.keyval().end()) - row[3] = to(it_time->second); - PE.resize((header.ndim() > 3) ? header.size(3) : 1, cols); - PE.rowwise() = row.transpose(); - } - } - return PE; -} - -Eigen::MatrixXd get_scheme(const Header &header) { - DEBUG("searching for suitable phase encoding data..."); - using namespace App; - Eigen::MatrixXd result; - - try { - const auto opt_table = get_options("import_pe_table"); - if (!opt_table.empty()) - result = load(opt_table[0][0], header); - const auto opt_eddy = get_options("import_pe_eddy"); - if (!opt_eddy.empty()) { - if (!opt_table.empty()) - throw Exception("Phase encoding table can be provided" - " using either -import_pe_table or -import_pe_eddy option," - " but NOT both"); - result = load_eddy(opt_eddy[0][0], opt_eddy[0][1], header); - } - if (opt_table.empty() && opt_eddy.empty()) - result = parse_scheme(header); - } catch (Exception &e) { - throw Exception(e, "error importing phase encoding table for image \"" + header.name() + "\""); - } - - if (!result.rows()) - return result; - - if (result.cols() < 3) - throw Exception("unexpected phase encoding table matrix dimensions"); - - INFO("found " + str(result.rows()) + "x" + str(result.cols()) + " phase encoding table"); - - return result; -} - -Eigen::MatrixXd eddy2scheme(const Eigen::MatrixXd &config, const Eigen::Array &indices) { - if (config.cols() != 4) - throw Exception("Expected 4 columns in EDDY-format phase-encoding config file"); - Eigen::MatrixXd result(indices.size(), 4); - for (ssize_t row = 0; row != indices.size(); ++row) { - if (indices[row] > config.rows()) - throw Exception("Malformed EDDY-style phase-encoding information:" - " index exceeds number of config entries"); - result.row(row) = config.row(indices[row] - 1); - } - return result; -} - -void export_commandline(const Header &header) { - auto check = [&](const Eigen::MatrixXd &m) -> const Eigen::MatrixXd & { - if (!m.rows()) - throw Exception("no phase-encoding information found within image \"" + header.name() + "\""); - return m; - }; - - auto scheme = parse_scheme(header); - - auto opt = get_options("export_pe_table"); - if (!opt.empty()) - save(check(scheme), header, opt[0][0]); - - opt = get_options("export_pe_eddy"); - if (!opt.empty()) - save_eddy(check(scheme), header, opt[0][0], opt[0][1]); -} - -} // namespace MR::PhaseEncoding diff --git a/cpp/core/phase_encoding.h b/cpp/core/phase_encoding.h deleted file mode 100644 index eb12d79591..0000000000 --- a/cpp/core/phase_encoding.h +++ /dev/null @@ -1,281 +0,0 @@ -/* Copyright (c) 2008-2025 the MRtrix3 contributors. - * - * This Source Code Form is subject to the terms of the Mozilla Public - * License, v. 2.0. If a copy of the MPL was not distributed with this - * file, You can obtain one at http://mozilla.org/MPL/2.0/. - * - * Covered Software is provided under this License on an "as is" - * basis, without warranty of any kind, either expressed, implied, or - * statutory, including, without limitation, warranties that the - * Covered Software is free of defects, merchantable, fit for a - * particular purpose or non-infringing. - * See the Mozilla Public License v. 2.0 for more details. - * - * For more details, see http://www.mrtrix.org/. - */ - -#pragma once - -#include -#include - -#include "app.h" -#include "axes.h" -#include "file/matrix.h" -#include "file/nifti_utils.h" -#include "file/ofstream.h" -#include "header.h" - -namespace MR::PhaseEncoding { - -extern const App::OptionGroup ImportOptions; -extern const App::OptionGroup SelectOptions; -extern const App::OptionGroup ExportOptions; - -//! check that a phase-encoding table is valid -template void check(const MatrixType &PE) { - if (!PE.rows()) - throw Exception("No valid phase encoding table found"); - - if (PE.cols() < 3) - throw Exception("Phase-encoding matrix must have at least 3 columns"); - - for (ssize_t row = 0; row != PE.rows(); ++row) { - for (ssize_t axis = 0; axis != 3; ++axis) { - if (std::round(PE(row, axis)) != PE(row, axis)) - throw Exception("Phase-encoding matrix contains non-integral axis designation"); - } - } -} - -//! check that the PE scheme matches the DWI data in \a header -template void check(const MatrixType &PE, const HeaderType &header) { - check(PE); - const ssize_t num_volumes = (header.ndim() > 3) ? header.size(3) : 1; - if (num_volumes != PE.rows()) - throw Exception("Number of volumes in image \"" + header.name() + "\" (" + str(num_volumes) + - ") does not match that in phase encoding table (" + str(PE.rows()) + ")"); -} - -//! store the phase encoding matrix in a header -/*! this will store the phase encoding matrix into the - * Header::keyval() structure of \a header. - * - If the phase encoding direction and/or total readout - * time varies between volumes, then the information - * will be stored in field "pe_scheme"; if not, it - * will instead be stored in fields "PhaseEncodingDirection" - * and "TotalReadoutTime" - */ -template void set_scheme(Header &header, const MatrixType &PE) { - auto erase = [&](const std::string &s) { - auto it = header.keyval().find(s); - if (it != header.keyval().end()) - header.keyval().erase(it); - }; - if (!PE.rows()) { - erase("pe_scheme"); - erase("PhaseEncodingDirection"); - erase("TotalReadoutTime"); - return; - } - check(PE, header); - std::string pe_scheme; - std::string first_line; - bool variation = false; - for (ssize_t row = 0; row < PE.rows(); ++row) { - std::string line = str(PE(row, 0)); - for (ssize_t col = 1; col < PE.cols(); ++col) - line += "," + str(PE(row, col), 3); - add_line(pe_scheme, line); - if (first_line.empty()) - first_line = line; - else if (line != first_line) - variation = true; - } - if (variation) { - header.keyval()["pe_scheme"] = pe_scheme; - erase("PhaseEncodingDirection"); - erase("TotalReadoutTime"); - } else { - erase("pe_scheme"); - 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); - else - erase("TotalReadoutTime"); - } -} - -//! clear the phase encoding matrix from a header -/*! this will delete any trace of phase encoding information - * from the Header::keyval() structure of \a header. - */ -void clear_scheme(Header &header); - -//! parse the phase encoding matrix from a header -/*! extract the phase encoding matrix stored in the \a header if one - * is present. This is expected to be stored in the Header::keyval() - * structure, under the key 'pe_scheme'. Alternatively, if the phase - * encoding direction and bandwidth is fixed for all volumes in the - * series, this information may be stored using the keys - * 'PhaseEncodingDirection' and 'TotalReadoutTime'. - */ -Eigen::MatrixXd parse_scheme(const Header &); - -//! get a phase encoding matrix -/*! get a valid phase-encoding matrix, either from files specified at - * the command-line, or from the contents of the image header. - */ -Eigen::MatrixXd get_scheme(const Header &); - -//! Convert a phase-encoding scheme into the EDDY config / indices format -template -void scheme2eddy(const MatrixType &PE, Eigen::MatrixXd &config, Eigen::Array &indices) { - try { - check(PE); - } catch (Exception &e) { - throw Exception(e, "Cannot convert phase-encoding scheme to eddy format"); - } - if (PE.cols() != 4) - throw Exception("Phase-encoding matrix requires 4 columns to convert to eddy format"); - config.resize(0, 0); - indices = Eigen::Array::Constant(PE.rows(), PE.rows()); - for (ssize_t PE_row = 0; PE_row != PE.rows(); ++PE_row) { - for (ssize_t config_row = 0; config_row != config.rows(); ++config_row) { - bool dir_match = PE.template block<1, 3>(PE_row, 0).isApprox(config.block<1, 3>(config_row, 0)); - bool time_match = abs(PE(PE_row, 3) - config(config_row, 3)) < 1e-3; - if (dir_match && time_match) { - // FSL-style index file indexes from 1 - indices[PE_row] = config_row + 1; - break; - } - } - if (indices[PE_row] == PE.rows()) { - // No corresponding match found in config matrix; create a new entry - config.conservativeResize(config.rows() + 1, 4); - config.row(config.rows() - 1) = PE.row(PE_row); - indices[PE_row] = config.rows(); - } - } -} - -//! Convert phase-encoding infor from the EDDY config / indices format into a standard scheme -Eigen::MatrixXd eddy2scheme(const Eigen::MatrixXd &, const Eigen::Array &); - -//! Modifies a phase encoding scheme if being imported alongside a non-RAS image -template -Eigen::MatrixXd transform_for_image_load(const MatrixType &pe_scheme, const HeaderType &H) { - std::array perm; - std::array flip; - H.realignment(perm, flip); - if (perm[0] == 0 && perm[1] == 1 && perm[2] == 2 && !flip[0] && !flip[1] && !flip[2]) { - INFO("No transformation of external phase encoding data required to accompany image \"" + H.name() + "\""); - return pe_scheme; - } - Eigen::MatrixXd result(pe_scheme.rows(), pe_scheme.cols()); - for (ssize_t row = 0; row != pe_scheme.rows(); ++row) { - Eigen::VectorXd new_line = pe_scheme.row(row); - for (ssize_t axis = 0; axis != 3; ++axis) { - new_line[axis] = pe_scheme(row, perm[axis]); - if (new_line[axis] && flip[perm[axis]]) - new_line[axis] = -new_line[axis]; - } - result.row(row) = new_line; - } - INFO("External phase encoding data transformed to match RAS realignment of image \"" + H.name() + "\""); - return result; -} - -//! Modifies a phase encoding scheme if being exported alongside a NIfTI image -template -Eigen::MatrixXd transform_for_nifti_write(const MatrixType &pe_scheme, const HeaderType &H) { - std::vector order; - std::array flip; - File::NIfTI::axes_on_write(H, order, flip); - if (order[0] == 0 && order[1] == 1 && order[2] == 2 && !flip[0] && !flip[1] && !flip[2]) { - INFO("No transformation of phase encoding data required for export to file"); - return pe_scheme; - } - Eigen::Matrix result(pe_scheme.rows(), pe_scheme.cols()); - for (ssize_t row = 0; row != pe_scheme.rows(); ++row) { - Eigen::VectorXd new_line = pe_scheme.row(row); - for (ssize_t axis = 0; axis != 3; ++axis) - new_line[axis] = - pe_scheme(row, order[axis]) && flip[axis] ? -pe_scheme(row, order[axis]) : pe_scheme(row, order[axis]); - result.row(row) = new_line; - } - INFO("Phase encoding data transformed to match NIfTI / MGH image export prior to writing to file"); - return result; -} - -namespace { -template void __save(const MatrixType &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 - -//! Save a phase-encoding scheme 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 -template -void save(const MatrixType &PE, const HeaderType &header, const std::string &path) { - try { - check(PE, header); - } catch (Exception &e) { - throw Exception(e, "Cannot export phase-encoding table to file \"" + path + "\""); - } - - if (Path::has_suffix(header.name(), {".mgh", ".mgz", ".nii", ".nii.gz", ".img"})) { - __save(transform_for_nifti_write(PE, header), path); - } else { - __save(PE, path); - } -} - -//! Save a phase-encoding scheme to EDDY format config / index files -template -void save_eddy(const MatrixType &PE, - const HeaderType &header, - const std::string &config_path, - const std::string &index_path) { - Eigen::MatrixXd config; - Eigen::Array indices; - scheme2eddy(transform_for_nifti_write(PE, header), config, indices); - File::Matrix::save_matrix(config, config_path, KeyValues(), false); - File::Matrix::save_vector(indices, index_path, KeyValues(), false); -} - -//! Save the phase-encoding scheme from a header to file depending on command-line options -void export_commandline(const Header &); - -//! Load a phase-encoding scheme from a matrix text file -template Eigen::MatrixXd load(const std::string &path, const HeaderType &header) { - const Eigen::MatrixXd PE = File::Matrix::load_matrix(path); - check(PE, header); - // As with JSON import, need to query the header to discover if the - // strides / transform were modified on image load to make the image - // data appear approximately axial, in which case we need to apply the - // same transforms to the phase encoding data on load - return transform_for_image_load(PE, header); -} - -//! Load a phase-encoding scheme from an EDDY-format config / indices file pair -template -Eigen::MatrixXd load_eddy(const std::string &config_path, const std::string &index_path, const HeaderType &header) { - const Eigen::MatrixXd config = File::Matrix::load_matrix(config_path); - const Eigen::Array indices = File::Matrix::load_vector(index_path); - const Eigen::MatrixXd PE = eddy2scheme(config, indices); - check(PE, header); - return transform_for_image_load(PE, header); -} - -} // namespace MR::PhaseEncoding diff --git a/cpp/core/surface/filter/vertex_transform.cpp b/cpp/core/surface/filter/vertex_transform.cpp index ab80f71f33..4317156116 100644 --- a/cpp/core/surface/filter/vertex_transform.cpp +++ b/cpp/core/surface/filter/vertex_transform.cpp @@ -17,6 +17,7 @@ #include "surface/filter/vertex_transform.h" #include "file/nifti_utils.h" +#include "axes.h" #include "exception.h" namespace MR::Surface::Filter { @@ -82,14 +83,13 @@ void VertexTransform::operator()(const Mesh &in, Mesh &out) const { break; case transform_t::FS2REAL: - std::vector axes(3); - auto M = File::NIfTI::adjust_transform(header, axes); + auto M = header.realignment().orig_transform(); + const Axes::permutations_type &axes = header.realignment().permutations(); Eigen::Vector3d cras(3, 1); for (size_t i = 0; i < 3; i++) { cras[i] = M(i, 3); - for (size_t j = 0; j < 3; j++) { + for (size_t j = 0; j < 3; j++) cras[i] += 0.5 * header.size(axes[j]) * header.spacing(axes[j]) * M(i, j); - } } for (size_t i = 0; i != V; ++i) { vertices.push_back(in.vert(i) + cras);