Skip to content

ENH: bug #3241: Add qfac and qto_xyz to itkNiftiImageIO metadata#3242

Merged
dzenanz merged 1 commit intoInsightSoftwareConsortium:masterfrom
seanm:nifti-metadata
Mar 17, 2022
Merged

ENH: bug #3241: Add qfac and qto_xyz to itkNiftiImageIO metadata#3242
dzenanz merged 1 commit intoInsightSoftwareConsortium:masterfrom
seanm:nifti-metadata

Conversation

@seanm
Copy link
Copy Markdown
Contributor

@seanm seanm commented Mar 1, 2022

bug #3241: In f38b1dd I accidentally changed some behaviour. Specifically, the value of pixdim[0] in itk::NiftiImageIO metadata is different before and after the change.

Previous to the change, SetImageIOMetadataFromNIfTI() was populated from the nifti_1_header structure. After the change, it is populated from the nifti_image structure. The latter is created from the former with nifti_convert_nhdr2nim(). Most of the fields used in SetImageIOMetadataFromNIfTI() have 1-to-1 correspondences in these two structures. I expected all 8 elements of the pixdim array were the same in both structures, but in fact no; element 0 is special. nifti_convert_nhdr2nim() unconditionally sets pixdim[0] to 0.0 (by virtue of allocating the structure with calloc).

In my own app, I was using pixdim[0] (retrieved from itk::NiftiImageIO) as the last parameter to nifti_quatern_to_mat44(). I see now looking at niftilib itself that I should be passing qfac as the last parameter. But the itk::MetaDataObject doesn't include it, so this commit adds qfac to the metadata dictionary.

Additionally, we now include the whole qto_xyz matrix, so that ITK clients that call nifti_quatern_to_mat44() don't even need to. This plugs a leaky abstraction and allows my own app for example to not even know that ITK is implemented via niftilib, and never need to call down to it.

Note that the behaviour change of f38b1dd is not being reverted, because it's been four years now, no one else has noticed, and some people may be relying on this behaviour now.

Copy link
Copy Markdown
Member

@dzenanz dzenanz left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

LGTM

Comment thread Modules/IO/NIFTI/src/itkNiftiImageIO.cxx Outdated
Comment thread Modules/IO/NIFTI/src/itkNiftiImageIO.cxx Outdated
@thewtex thewtex requested a review from hjmjohnson March 3, 2022 01:42
@seanm seanm changed the title WIP: ENH: Add qfac and qto_xyz to itkNiftiImageIO metadata ENH: bug #3241: Add qfac and qto_xyz to itkNiftiImageIO metadata Mar 14, 2022
@seanm
Copy link
Copy Markdown
Contributor Author

seanm commented Mar 14, 2022

I've updated the commit message, and this is ready to merge IMHO.

@seanm
Copy link
Copy Markdown
Contributor Author

seanm commented Mar 14, 2022

But don't merge yet... I want to be 100% my own app works against this...

@seanm
Copy link
Copy Markdown
Contributor Author

seanm commented Mar 14, 2022

OK, good to go, I had a minor panic there for a second, but there was a typo in my own code :)

@dzenanz dzenanz requested a review from N-Dekker March 15, 2022 14:24
Comment on lines +1003 to +1018
qto_xyz.push_back(nim->qto_xyz.m[0][0]);
qto_xyz.push_back(nim->qto_xyz.m[0][1]);
qto_xyz.push_back(nim->qto_xyz.m[0][2]);
qto_xyz.push_back(nim->qto_xyz.m[0][3]);
qto_xyz.push_back(nim->qto_xyz.m[1][0]);
qto_xyz.push_back(nim->qto_xyz.m[1][1]);
qto_xyz.push_back(nim->qto_xyz.m[1][2]);
qto_xyz.push_back(nim->qto_xyz.m[1][3]);
qto_xyz.push_back(nim->qto_xyz.m[2][0]);
qto_xyz.push_back(nim->qto_xyz.m[2][1]);
qto_xyz.push_back(nim->qto_xyz.m[2][2]);
qto_xyz.push_back(nim->qto_xyz.m[2][3]);
qto_xyz.push_back(nim->qto_xyz.m[3][0]);
qto_xyz.push_back(nim->qto_xyz.m[3][1]);
qto_xyz.push_back(nim->qto_xyz.m[3][2]);
qto_xyz.push_back(nim->qto_xyz.m[3][3]);
Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

It appears that these 16 lines of code may be refactored to the following:

  for (const auto & row : nim->qto_xyz.m)
  {
    for (const float value : row)
    {
      qto_xyz.push_back(value);
    }
  }

Right?

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Or nested for r, for c loops.

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

For loops are more readable and less prone to errors. Let the compiler do the unrolling, if advantageous.

Copy link
Copy Markdown
Contributor

@Leengit Leengit Mar 16, 2022

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

If speed is somehow the issue then we might make use of the std::vector<>(begin, end) constructor, as in

const std::vector<float> qto_xyz(&nim->qto_xyz.m[0][0], &nim->qto_xyz.m[0][0] + 16);

Admittedly this is ugly and does assume something about the packing of the nim->qto_xyz.m object.

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think the nested r & c loops is best...

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Speed is not at all critical here. It's dwarfed by the file i/o to actually read the nifti file.

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Though not needed for the present case ...
It might be nice, it might be nice
To have

template <typename TUnderlying, int VRows, int VColumns>
constexpr inline std::vector<TUnderlying>
CArrayToStdVector(const TUnderlying (&arr)[VRows][VColumns])
{
  return { &arr[0][0], &arr[0][0] + VRows * VColumns };
}

on your side.

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Just made a pull request that would ease passing a C-style float [4][4] array directly to a itk::Matrix, please check: #3305 It would support the following:

using MatrixType = Matrix<float, 4, 4>;
EncapsulateMetaData<MatrixType>(thisDic, "qto_xyz", MatrixType(nim->qto_xyz.m));

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

template <typename TUnderlying, int VRows, int VColumns>
constexpr inline std::vector<TUnderlying>
CArrayToStdVector(const TUnderlying (&arr)[VRows][VColumns])
{
  return { &arr[0][0], &arr[0][0] + VRows * VColumns };
}

@Leengit It might be nice indeed, but it's still not obvious to me why a fixed-size two-dimensional array (matrix) would need to be flattened to a 1-dimensional vector, before storing it into the dictionary. It seems more intuitive to me to have an CArrayToStdArray(const TUnderlying (&arr)[VRows][VColumns]) function, returning an std::array of std::array elements. (std::array<std::array<TUnderlying, VColumns>, VRows>). Right?

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

returning an std::array of std::array elements

That would be nice.

qto_xyz.push_back(nim->qto_xyz.m[3][1]);
qto_xyz.push_back(nim->qto_xyz.m[3][2]);
qto_xyz.push_back(nim->qto_xyz.m[3][3]);
EncapsulateMetaData<std::vector<float>>(thisDic, "qto_xyz", qto_xyz);
Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Wouldn't it be preferable to pass nim->qto_xyz directly to EncapsulateMetaData, instead of copying its content to a std::vector<float>? As follows

EncapsulateMetaData<mat44>(thisDic, "qto_xyz", nim->qto_xyz);

I see, it not not yet compile, because there is no operator== overload for mat44. But that could be solved, of course.

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Encoding mat44 itself I don't like, because it's a niftilib data structure, which I don't want to leak outside ITK. Adding an operator== to it would also require changing niftilib, and it's C anyway.

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

mat44 is defined in the global namespace (obviously, because it is C), so technically it is allowed to add an operator== to "our" global namespace, for example as follows:

bool
operator==(const mat44 & lhs, const mat44 & rhs)
{
  return std::equal(
    std::begin(lhs.m), std::end(lhs.m), std::begin(rhs.m), [](const auto & leftRow, const auto & rightRow) {
      return std::equal(std::begin(leftRow), std::end(leftRow), std::begin(rightRow));
    });
}

Another option might be to adjust itk::MetaDataObject<MetaDataObjectType> to allow types without an operator==.

I realize that it is my fault that itk::MetaDataObject<MetaDataObjectType> requires operator== for the type to be stored: pull request #2246 commit 694bbc0 😮 It could be adjusted to use std::memcmp, for types that don't have an operator==. (@Leengit A SFINAE opportunity!!!)

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Frankly, this strikes me as an illegible complication, for little benefit.

Copy link
Copy Markdown
Contributor

@N-Dekker N-Dekker Mar 16, 2022

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Encoding mat44 itself I don't like, because it's a niftilib data structure

OK, but then, what would you think about encapsulating the matrix as a itk::Matrix<float, 4, 4>?

  using MatrixType = Matrix<float, 4, 4>;
  EncapsulateMetaData<MatrixType>(thisDic, "qto_xyz", MatrixType::InternalMatrixType(&nim->qto_xyz.m[0][0]));

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

That doesn't make me recoil with C++ whiplash. :) But it starts to be very different from the rest of the existing code in that function, which seems to prefer basic data structures.

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The other thing I like about this PR as-is is that's its easy to cherry pick just this change. For my own app, I'll be cherry-picking this into ITK 5.2 for a hot fix.

Copy link
Copy Markdown
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

OK, but do you agree that itk::Matrix<float, 4, 4> would be a proper data type to use when storing nim->qto_xyz into the itk::MetaDataDictionary?

Copy link
Copy Markdown
Member

@blowekamp blowekamp Mar 17, 2022

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The type used here should be one of the explicitly instantiated MetaDataObject types:
https://github.com/InsightSoftwareConsortium/ITK/blob/master/Modules/Core/Common/include/itkMetaDataObject.h#L298-L321

I believe only itk::Matrix<float,3,3> is there, itk::Matrix<float,4,4> would need to be added.

std::vector is already there.

…iftiImageIO metadata

bug InsightSoftwareConsortium#3241: In f38b1dd I accidentally changed some behaviour. Specifically, the value of pixdim[0] in itk::NiftiImageIO metadata is different before and after the change.

Previous to the change, SetImageIOMetadataFromNIfTI() was populated from the nifti_1_header structure. After the change, it is populated from the nifti_image structure. The latter is created from the former with nifti_convert_nhdr2nim(). Most of the fields used in SetImageIOMetadataFromNIfTI() have 1-to-1 correspondences in these two structures. I expected all 8 elements of the pixdim array were the same in both structures, but in fact no; element 0 is special. nifti_convert_nhdr2nim() unconditionally sets pixdim[0] to 0.0 (by virtue of allocating the structure with calloc).

In my own app, I was using pixdim[0] (retrieved from itk::NiftiImageIO) as the last parameter to nifti_quatern_to_mat44(). I see now looking at niftilib itself that I should be passing qfac as the last parameter.  But the itk::MetaDataObject doesn't include it, so this commit adds qfac to the metadata dictionary.

Additionally, we now include the whole qto_xyz matrix, so that ITK clients that call nifti_quatern_to_mat44() don't even need to.  This plugs a leaky abstraction and allows my own app for example to not even know that ITK is implemented via niftilib, and never need to call down to it.

Note that the behaviour change of f38b1dd is not being reverted, because it's been four years now, no one else has noticed, and some people may be relying on this behaviour now.
@github-actions github-actions Bot added the type:Enhancement Improvement of existing methods or implementation label Mar 16, 2022
Copy link
Copy Markdown
Member

@dzenanz dzenanz left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is good as-is, but if you want to improve it further feel free to modify.

Copy link
Copy Markdown
Contributor

@N-Dekker N-Dekker left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks @seanm I gave up my resistance 😺 While I still think it would be nicer to store nim->qto_xyz as a 4x4 float matrix (either mat44, or itk::Matrix, or an std::array of std::array's), Bradley (@blowekamp) has just explained that neither of them are currently supported as MetaDataObject type. You might still consider a vector-of-vectors, std::vector<std::vector<float>>, but that's up to you.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

area:IO Issues affecting the IO module type:Enhancement Improvement of existing methods or implementation

Projects

None yet

Development

Successfully merging this pull request may close these issues.

6 participants