BUG: Fix reading of displacement field vectors from NIFTI

[lassoan]

NIFTI image containing displacement field (intent = NIFTI_INTENT_DISPVECT) was read as scalar volume (all 3 vector components were set to the same value). Fixed by reading it the same way as the general-purpose vector image (intent = NIFTI_INTENT_VECTOR).

PR Checklist

• No API changes were made (or the changes have been approved)
• No major design changes were made (or the changes have been approved)
• Added test (or behavior not changed)
• Updated API documentation (or API not changed)
• Added license to new files (if any)
• Added Python wrapping to new files (if any) as described in ITK Software Guide Section 9.5
• Added ITK examples for all new major features (if any)

The changes already improve NIFTI compatibility, but we could potentially further improve correctness if we can answer these questions:

1. Should the displacement vectors be converted from RAS (native NIFTI) to LPS coordinate system (ITK convention)? The same way as the reader changes the image origin, spacing, axis directions to be in LPS coordinate system instead RAS.
2. NIFTI distinguishes between "vector" and "displacement vector" intent. Is there a way to keep this distinction in ITK? It would be nice if we read a displacement vector image and then we save it to file again to preserve the exact type.

[N-Dekker]

Thanks @lassoan Would it be an idea to add such information to the MetaDataDictionary of the image? For example:

  // Assume LPS if unspecified.
  EncapsulateMetaData(
    metaDataDictionary, "CoordinateSystem", SpatialOrientationEnums::ValidCoordinateOrientations::ITK_COORDINATE_ORIENTATION_RAS);
  EncapsulateMetaData(metaDataDictionary, "Intent", "DisplacementField");

[dzenanz]

Of course, the writer then needs to make of the information in MetaDataDictionary.

[lassoan]

Thanks for the suggestions!

MetaDataDictionary could work well for reading/writing the intent code. I'll add a commit to this PR to use these dictionary fields.

For the vector reorientation, I see two main options:

• Option A: ITK always converts to/from RAS if intent is "displacement vector" (does not reorient if intent is just generic "vector"). This would make all ITK-based tools handle "displacement vector" data sets as expected. Helper functions for converting between RAS/LPS coordinate system may be made available to users to manually deal with data sets that use generic "vector" intent for storing spatial data.
• Option B: Let the user reorient (regardless of "displacement vector" or "vector" intent), based on the MetadataDictionary field. This would be simpler for ITK but most likely users/applications would not be aware that flipping is necessary and they would just work incorrectly.

Would you have any preference?

[dzenanz]

Based on your description, I prefer A.

[N-Dekker]

I think I would prefer option A as well, unless some users would prefer to remain in RAS space...? If option A can be implemented, of course my suggested ("CoordinateSystem", ITK_COORDINATE_ORIENTATION_RAS) dictionary entry would not be necessary anymore. I don't know if an ("Intent", "DisplacementField") dictonary entry would then still be useful either. No problem.

[lassoan]

I was thinking about the automatic conversion to RAS, too. Slicer currently uses RAS internally, so after ITK reads a displacement field from a NRRD file, Slicer convert all the vectors to RAS. In case of NIFTI it would mean two conversions - RAS->LPS in ITK then LPS->RAS in Slicer. The time needed for conversion is usually not noticeable, so we could implement this simple, always-convert solution, and then if somebody runs into a performance issue then a bypass mechanism could be implemented.

[gdevenyi]

The ITK MINC reader needs this LPS/RAS conversion as it currently blindly assumes LPS, which is wrong, see #147

[dzenanz]

@lassoan I didn't deal much with DisplacementFields, but I assume that converting it to LPS coordinate system would be preferable. Or maybe have a setting in NiftiIO which controls it, and defaults to ConvertToLPS?

[dzenanz]

This PR looks good as is. The discussed improvements could be done in this PR, or a follow up PR.

[lassoan]

I didn't deal much with DisplacementFields, but I assume that converting it to LPS coordinate system would be preferable. Or maybe have a setting in NiftiIO which controls it, and defaults to ConvertToLPS

Yes - see the two options that I described above.

Probably Option A would be more consistent with what ITK does already (automatic conversion of image geometry to LPS).

[lassoan]

Could you recommend a simple and efficient way to rescale vector components with (-1,-1,1) factors to switch between LPS/RAS?

[dzenanz]

MultiplyImageFilter sounds as the right solution.

[lassoan]

Is there an example of using it within a filter?

The class already can do some rescaling but only for scalars. Maybe it could be just modified to allow applying different factors for each vector component?

[N-Dekker]

@lassoan Something like this?

using VectorType = itk::Vector<float, 3>;
using DisplacementFieldImageType = itk::Image<VectorType>;

void
ConvertCoordinates(DisplacementFieldImageType & displacementFieldImage)
{
  for (VectorType & displacementVector : itk::ImageBufferRange<DisplacementFieldImageType>(displacementFieldImage))
  {
    displacementVector[0] *= -1;
    displacementVector[1] *= -1;
  }
}

This range-based for-loop would iterate all the vectors in the displacement field image, and multiply the first and the second vector component "in place" by -1.

[dzenanz]

multiplyimagebyscalar name was given to differentiate it from multiplying pixels of two images. If you replace scalar types by vector types, it might work. If not, you might also need to provide definition for vector multiplication. But the filter should take care of parallelizing the operation. @N-Dekker's suggestion is not parallel.

[gdevenyi]

The ITK MINC reader has similar issues of blindly reading in as LPS without conversion, see #147. I would love to see a generic solution to this as its needed for that IO system as well. @vfonov

[lassoan]

Thanks for all the suggestions I'm still reviewing them to see which ones are the best applicable.

@gdevenyi This is indeed the same kind of issue as for MINC. It would be great to reuse the code that you implemented for MINC, but I find it hard (even if I copy-pasted the code) because the file formats are just too different and the IO classes are already way too complicated. Maybe by adding an adaptor layer outside could address this, but that would mean a much more complicated design, and I think it would only be needed for NIFTI and MINC, so I'm not sure if it is worth the effort.

[jhlegarreta]

Would adding a test or modifying an existing one to account for this be relevant/doable?

[lassoan]

I've added conversion of vectors between RAS/LPS coordinate systems. By default, conversion is enabled for displacement vector images and disabled for generic vector images. These defaults can be overridden by two member variables in the IO file.

I've also added the possibility of writing displacement vector images, by setting the "intent_code"="1006" in the metadata dictionary. I've kept the magic number stored as string in the metadata dictionary, as the NIFTI IO class already stored intent code like this and I wanted to preserve backward compatibility.

@dzenanz @N-Dekker it would be great if you could review.

[lassoan]

@N-Dekker Thank you very much for the review, I've addressed your comments and force-pushed. Please have a look.

[N-Dekker]

The library code looks fine to me, just two questions on the test!

[jhlegarreta]

Thanks for doing all this Andras 💯.

An in-line comment about printing m_ConvertRAS.

[lassoan]

I've force-pushed the requested changes (except adding more tests).

[jhlegarreta]

Sorry 🙏 for another review on the documentation Andras.

[N-Dekker]

Thanks for addressing my comments @lassoan. Your pull request looks fine to me now 👍

[jhlegarreta]

Valgrind is showing some memory defects related to the NIFTI IO class. Not sure which is the change that triggered them, as they have only started appearing today (Valgrind builds have been green for 3 days after this branch was merged):
https://open.cdash.org/viewDynamicAnalysisFile.php?id=9949304

Valgrind points to
https://github.com/InsightSoftwareConsortium/ITK/blob/master/Modules/IO/NIFTI/src/itkNiftiImageIO.cxx#L708

Which was not modified in this PR in reality.

The commits on Jan 26-28 do not seem to introduce related changes (unsure if #3852 has side-effects).

[jhlegarreta]

#3884 (comment) Looks like it was a hiccup from Valgrind: no defects have been notified since 3 days ago. The comment can thus be ignored. Sorry for the noise.

[lassoan]

Thanks for the update. I reviewed NIFTI IO class and noticed that the implementation was kind of fragile, because memory allocation happened very far from the place where the buffer is used; but I did not see any specific issue. So, I'm glad it turned out that there was no real failure.

[jhlegarreta]

#3884 (comment) Thanks for the effort Andras. Maybe something to watch then; if noticed again, maybe we can open an issue.

[blowekamp]

@jhlegarreta The Valgrind run I saw had many tests failing. This can indicate some type of abnormal execution of the tests. In such a case, the processes may be terminated prematurely ( due to some OS resource limitation encountered) where the memory will not be cleaned up correctly. This can trigger false positives of valgrind defects.

Therefore, if there are abnormal or inconsistent number of tests failing during valgrind then defects should not be considered actionable, just monitored to see if they are consistent.