DWI metadata handling fixes #2917
Conversation
- New class Axes::Shuffle, which encapsulates both permutations and flips so that these can in combination be returned by functions. - Move some handling of slice encoding direction transformation to core/metadata/slice_encoding.h. - Move handling of minor differences in slice timing to core/metadata/slice_encoding.h.
- For many functions relating to phase and slice encoding, operate on key-value dictionaries rather than Headers; this allows manipulation where that information is associated with external files rather than image hedaers. - DWI::load_bvecs_bvals(): Simplify implementation given availability of original image header on disk. - For JSON import, do not prioritise image header contents over JSON contents; instead use Header::merge_keyval() function to explicitly flag conflicts. - For JSON load, ensure that any operations relating to internal header transform realignment are applied only to those contents loaded from the JSON, not anything already stored in the header. - Always call Header::realign_transform() in order to capture on-disk information; just exit prematurely if do_realign_transform is false. - Homogenise function interfaces between phase and slice encoding manipulations. - If writing a JSON alongside a NIfTI, suppress "dw_scheme" from its contents.
![github-actions[bot]](https://avatars.githubusercontent.com/in/15368?s=60&v=4){kind=small}
| 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; |
There was a problem hiding this comment.
warning: narrowing conversion from 'value_type' (aka 'unsigned long') to signed type 'Index' (aka 'long') is implementation-defined [bugprone-narrowing-conversions]
result.flips[result.permutations[0]] = T(0, result.permutations[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; |
There was a problem hiding this comment.
warning: narrowing conversion from 'value_type' (aka 'unsigned long') to signed type 'Index' (aka 'long') is implementation-defined [bugprone-narrowing-conversions]
result.flips[result.permutations[1]] = T(1, result.permutations[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; |
There was a problem hiding this comment.
warning: narrowing conversion from 'value_type' (aka 'unsigned long') to signed type 'Index' (aka 'long') is implementation-defined [bugprone-narrowing-conversions]
result.flips[result.permutations[2]] = T(2, result.permutations[2]) < 0.0;
^| return (permutations[0] != 0 || permutations[1] != 1 || permutations[2] != 2 || // | ||
| flips[0] || flips[1] || flips[2]); | ||
| } | ||
| permutations_type permutations; |
There was a problem hiding this comment.
warning: member variable 'permutations' has public visibility [misc-non-private-member-variables-in-classes]
permutations_type permutations;
^| flips[0] || flips[1] || flips[2]); | ||
| } | ||
| permutations_type permutations; | ||
| flips_type flips; |
There was a problem hiding this comment.
warning: member variable 'flips' has public visibility [misc-non-private-member-variables-in-classes]
flips_type flips;
^| stride(i) = -stride(i); | ||
| } | ||
|
|
||
| // copy back 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)]}; |
There was a problem hiding this comment.
warning: variable 'a' of type 'Axis[3]' can be declared 'const' [misc-const-correctness]
| Axis a[] = {axes_[realignment_.permutation(0)], axes_[realignment_.permutation(1)], axes_[realignment_.permutation(2)]}; | |
| Axis const a[] = {axes_[realignment_.permutation(0)], axes_[realignment_.permutation(1)], axes_[realignment_.permutation(2)]}; |
| } | ||
|
|
||
| Header | ||
| concatenate(const std::vector<Header> &headers, const size_t axis_to_concat, const bool permit_datatype_mismatch) { | ||
| Header concatenate(const std::vector<Header> &headers, |
There was a problem hiding this comment.
warning: function 'concatenate' has cognitive complexity of 61 (threshold 25) [readability-function-cognitive-complexity]
Header concatenate(const std::vector<Header> &headers,
^Additional context
core/header.cpp:648: nesting level increased to 1
auto datatype_test = [&](const bool condition) {
^core/header.cpp:649: +2, including nesting penalty of 1, nesting level increased to 2
if (condition && !permit_datatype_mismatch) {
^core/header.cpp:649: +1
if (condition && !permit_datatype_mismatch) {
^core/header.cpp:656: nesting level increased to 1
auto concat_scheme = [](Eigen::MatrixXd &existing, const Eigen::MatrixXd &extra) {
^core/header.cpp:657: +2, including nesting penalty of 1, nesting level increased to 2
if (!existing.rows())
^core/header.cpp:659: +2, including nesting penalty of 1, nesting level increased to 2
if (!extra.rows() || (extra.cols() != existing.cols())) {
^core/header.cpp:659: +1
if (!extra.rows() || (extra.cols() != existing.cols())) {
^core/header.cpp:667: +1, including nesting penalty of 0, nesting level increased to 1
if (headers.empty())
^core/header.cpp:671: +1, including nesting penalty of 0, nesting level increased to 1
for (const auto &H : headers) {
^core/header.cpp:672: +2, including nesting penalty of 1, nesting level increased to 2
if (axis_to_concat > H.ndim() + 1) {
^core/header.cpp:677: +2, including nesting penalty of 1, nesting level increased to 2
for (this_max_nonunity_dim = H.ndim() - 1; this_max_nonunity_dim >= 0 && H.size(this_max_nonunity_dim) <= 1;
^core/header.cpp:677: +1
for (this_max_nonunity_dim = H.ndim() - 1; this_max_nonunity_dim >= 0 && H.size(this_max_nonunity_dim) <= 1;
^core/header.cpp:685: +1, including nesting penalty of 0, nesting level increased to 1
if (axis_to_concat >= result.ndim()) {
^core/header.cpp:692: +1, including nesting penalty of 0, nesting level increased to 1
for (size_t axis = 0; axis != result.ndim(); ++axis) {
^core/header.cpp:693: +2, including nesting penalty of 1, nesting level increased to 2
if (axis != axis_to_concat && result.size(axis) <= 1) {
^core/header.cpp:693: +1
if (axis != axis_to_concat && result.size(axis) <= 1) {
^core/header.cpp:694: +3, including nesting penalty of 2, nesting level increased to 3
for (const auto &H : headers) {
^core/header.cpp:695: +4, including nesting penalty of 3, nesting level increased to 4
if (H.ndim() > axis) {
^core/header.cpp:705: +1, including nesting penalty of 0, nesting level increased to 1
if (axis_to_concat == 3) {
^core/header.cpp:708: +2, including nesting penalty of 1, nesting level increased to 2
} catch (Exception &) {
^core/header.cpp:712: +2, including nesting penalty of 1, nesting level increased to 2
} catch (Exception &) {
^core/header.cpp:716: +1, including nesting penalty of 0, nesting level increased to 1
for (size_t i = 1; i != headers.size(); ++i) {
^core/header.cpp:720: +2, including nesting penalty of 1, nesting level increased to 2
for (size_t axis = 0; axis <= global_max_nonunity_dim; ++axis) {
^core/header.cpp:721: +3, including nesting penalty of 2, nesting level increased to 3
if (axis != axis_to_concat && axis < H.ndim() && H.size(axis) != result.size(axis)) {
^core/header.cpp:721: +1
if (axis != axis_to_concat && axis < H.ndim() && H.size(axis) != result.size(axis)) {
^core/header.cpp:729: +2, including nesting penalty of 1, nesting level increased to 2
result.size(axis_to_concat) += H.ndim() <= axis_to_concat ? 1 : H.size(axis_to_concat);
^core/header.cpp:732: +2, including nesting penalty of 1, nesting level increased to 2
if (axis_to_concat == 3) {
^core/header.cpp:736: +3, including nesting penalty of 2, nesting level increased to 3
} catch (Exception &) {
^core/header.cpp:742: +3, including nesting penalty of 2, nesting level increased to 3
} catch (Exception &) {
^core/header.cpp:752: +2, including nesting penalty of 1, nesting level increased to 2
if (datatype_test(!result.datatype().is_complex() && H.datatype().is_complex()))
^core/header.cpp:752: +1
if (datatype_test(!result.datatype().is_complex() && H.datatype().is_complex()))
^core/header.cpp:754: +2, including nesting penalty of 1, nesting level increased to 2
if (datatype_test(result.datatype().is_integer() && !result.datatype().is_signed() && H.datatype().is_signed()))
^core/header.cpp:754: +1
if (datatype_test(result.datatype().is_integer() && !result.datatype().is_signed() && H.datatype().is_signed()))
^core/header.cpp:756: +2, including nesting penalty of 1, nesting level increased to 2
if (datatype_test(result.datatype().is_integer() && H.datatype().is_floating_point()))
^core/header.cpp:756: +1
if (datatype_test(result.datatype().is_integer() && H.datatype().is_floating_point()))
^core/header.cpp:758: +2, including nesting penalty of 1, nesting level increased to 2
if (datatype_test(result.datatype().bytes() < H.datatype().bytes()))
^core/header.cpp:762: +1, including nesting penalty of 0, nesting level increased to 1
if (axis_to_concat == 3) {
^|
|
||
| Header::Realignment::Realignment() : | ||
| applied_transform_(applied_transform_type::Identity()), | ||
| orig_keyval_() { |
There was a problem hiding this comment.
warning: initializer for member 'orig_keyval_' is redundant [readability-redundant-member-init]
| orig_keyval_() { | |
| { |
| @@ -70,8 +69,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_) {} | |||
There was a problem hiding this comment.
warning: move constructor initializes class member by calling a copy constructor [performance-move-constructor-init]
realignment_(H.realignment_) {}
^Additional context
core/header.h:202: copy constructor being called
class Realignment {
^core/header.h:202: candidate move constructor here
class Realignment {
^| @@ -397,8 +412,7 @@ | |||
| void realign_transform(); | |||
| /*! store information about how image was | |||
| * realigned via realign_transform(). */ | |||
| std::array<size_t, 3> realign_perm_; | |||
| std::array<bool, 3> realign_flip_; | |||
| Realignment realignment_; | |||
There was a problem hiding this comment.
warning: member variable 'realignment_' has protected visibility [cppcoreguidelines-non-private-member-variables-in-classes]
Realignment realignment_;
^
| 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; |
There was a problem hiding this comment.
warning: constness of 'ondisk_transform' prevents automatic move [performance-no-automatic-move]
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)); |
There was a problem hiding this comment.
warning: narrowing conversion from 'ssize_t' (aka 'long') to 'default_type' (aka 'double') [bugprone-narrowing-conversions]
coord_switch(0, 3) = (header.size(header.realignment().permutation(0)) - 1) * header.spacing(header.realignment().permutation(0));
^| } 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()) + "\""); | ||
| } |
There was a problem hiding this comment.
warning: do not use 'else' after 'return' [readability-else-after-return]
| } 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()) + "\""); | |
| } | |
| } 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()) + "\""); | |
| } |
| 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 + "\""); |
There was a problem hiding this comment.
warning: do not use 'else' after 'return' [readability-else-after-return]
| 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 + "\""); | |
| 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 + "\""); |
| // clang-format on | ||
|
|
||
| void check(const scheme_type &PE) { | ||
| if (!PE.rows()) |
There was a problem hiding this comment.
warning: implicit conversion 'Index' (aka 'long') -> bool [readability-implicit-bool-conversion]
| if (!PE.rows()) | |
| if (PE.rows() == 0) |
| 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; |
There was a problem hiding this comment.
warning: variable 'time_match' of type 'bool' can be declared 'const' [misc-const-correctness]
| bool time_match = abs(PE(PE_row, 3) - config(config_row, 3)) < 1e-3; | |
| bool const time_match = abs(PE(PE_row, 3) - config(config_row, 3)) < 1e-3; |
| 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; |
There was a problem hiding this comment.
warning: narrowing conversion from 'ssize_t' (aka 'long') to signed type 'Scalar' (aka 'int') is implementation-defined [bugprone-narrowing-conversions]
indices[PE_row] = config_row + 1;
^| // 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(); |
There was a problem hiding this comment.
warning: narrowing conversion from 'Index' (aka 'long') to signed type 'Scalar' (aka 'int') is implementation-defined [bugprone-narrowing-conversions]
indices[PE_row] = config.rows();
^|
|
||
| void export_commandline(const Header &header) { | ||
| auto check = [&](const scheme_type &m) -> const scheme_type & { | ||
| if (!m.rows()) |
There was a problem hiding this comment.
warning: implicit conversion 'Index' (aka 'long') -> bool [readability-implicit-bool-conversion]
| if (!m.rows()) | |
| if (m.rows() == 0) |
| 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) { |
There was a problem hiding this comment.
warning: declaration uses identifier '__save', which is a reserved identifier [bugprone-reserved-identifier]
| void __save(const scheme_type &PE, const std::string &path) { | |
| void _save(const scheme_type &PE, const std::string &path) { |
core/metadata/phase_encoding.h:119:
- __save(PE, path);
+ _save(PE, path);
Closes #2477.
Work toward #2526.
This one has taken a lot more of my time than I would have liked. But it needed to be properly resolved, not just to fix issues for those disabling internal transform realignment but also because downstream work on BIDS BEP016 Diffusion models is dependent on these format definitions and conversion operations being robust.
The core mistake here was:
I've done a bunch of refactoring to store internally more exhaustive information about what (if anything) was performed at the internal header transform realignment phase. In addition to trivialising transformations between on-disk image header and scanner-space, this will also facilitate tracking exactly what differs between on-disk and interpreted versions of an image; not just transform & strides but potentially also metadata that are defined with respect to image axes.
According to my battery of tests, this all yields the correct results:
(DICOM, NIfTI / JSON / bvec / bval, MRtrix w. external metadata / MRtrix w. header metadata)
mrconvertordcm2niix(NIfTI / JSON / bvec / bval, MRtrix w. external metadata / MRtrix w. header metadata)
Questions:
Move
core/axes.*?From memory these used to be much larger; what's left there now is almost exclusively about transform realignment of spatial axes, which isn't really captured by the file name / namespace.
Move much of
dwi/gradient.*intocore/metadata/?It made sense for me to move
phase_encoding.*out ofcore/given I needed something similar for slice encoding and those are naturally grouped. But a gradient table is also "metadata that can be stored in an image header or external file(s) that necessitates advanced handling".Back-port to
master.Ultimately these are bug fixes, so should ideally be in
3.0.5.But maybe I'll try to restrict the scope of it by just making requisite changes and preserving a bit of code duplication that the refactoring removed / was done to hopefully facilitate mrconvert: fslgrad and RealignTransform: 0 #2477.
transformconvert flirt_import:File::NIfTI::axes_on_write()meshconvert -transform fs2real:File::NIfTI::axes_on_write()Check handling of user-specified phase encoding information
Particularly for
dwifslpreproc -pe_dir, there needs to be greater care around interpretation depending on whether the user's configuration does or does not include internal transform realignmentFix diffusion scheme docs page description of
bvecformatAppears to currently omit the left-handed vs. right-handed distinction, and also uses "x, y, z" to refer to image axes.