BUG: Fix CastImageFilter for VariableLengthVector

[SimonRit]

The Get function of VectorImage iterators calls the VariableLengthVector constructor parameterized by a const pointer and a vector length (see DefaultVectorPixelAccessor::Get). The m_Data pointer of the constructed VariableLengthVector which was returned was then pointing to the returned pixel data. This led to an unexpected behavior in CastImageFilter::DynamicThreadedGenerateDataDispatched where modifications to the local variable

OutputPixelType value{ outputIt.Get() };

of CastImageFilter::DynamicThreadedGenerateDataDispatched were actually modifying the pixel of the image at the iterator position. The problem might occur elsewhere so it is safer to copy the data when a VariableLengthVector is constructed from a const pointer despite the potential loss in performances.

A CastImageFilter test has been added to illustrate the issue: when a VectorImage was casted to an Image< Vector <, the first pixel of the region was taking the value of the last pixel of the region. The test illustrated this by creating an image filled with 0 except the first pixel. The casted image was completely full of 0 before bug correction.

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)

Refer to the ITK Software Guide for
further development details if necessary.

[hjmjohnson]

How does this interact with "LetArrayManageMemory"? It seems like m_LetArrayManageMemory should always be "true" for this case.

[SimonRit]

Good point, I oversaw this. So this is not a good fix then... I suggest modifying the iterator Get function then.

[SimonRit]

I'm a bit lost and happy to get some advice. I hope the issue is clearly explained in the current commit log. The current proposal is not adequate as pointed out by @hjmjohnson's comment. I see two potential alternatives:

• modify DefaultVectorPixelAccessor::Get and create a local copy which would have a computional cost as pointed out in this comment.
• modify the CastImageFilter and replace this line

OutputPixelType value{ outputIt.Get() };

with a copy

OutputPixelType value;
value = outputIt.Get();

I'm concerned that people actually do not expect a difference in behavior between the two so I would lean towards the first suggestion. Any thought?

[SimonRit]

What I suggest is to revert @blowekamp commit 6a60327... What do you think @blowekamp?

[blowekamp]

I do not have time to look at this closely today. But this change has the potential to cause threaded-per-pixel memory allocations which can make filters runs 100s time slower.

[SimonRit]

I do not have time to look at this closely today. But this change has the potential to cause threaded-per-pixel memory allocations which can make filters runs 100s time slower.

I understand but the bug I encountered demonstrates that it is tricky and counter intuitive to use an iterator on a VariableVectorLength. I can limit this PR to patching the CastImageFilter if you want but I'm afraid the problem occurs in other places.

[SimonRit]

I do not have time to look at this closely today. But this change has the potential to cause threaded-per-pixel memory allocations which can make filters runs 100s time slower.

I understand but the bug I encountered demonstrates that it is tricky and counter intuitive to use an iterator on a VariableVectorLength. I can limit this PR to patching the CastImageFilter if you want but I'm afraid the problem occurs in other places.

One alternative proposition: create a private member DefaultVectorPixelAccessor::m_GetVector (of type ExternalType), resized and allocated by the DefaultVectorPixelAccessor constructor, which is filled and returned by DefaultVectorPixelAccessor::Get?

[dzenanz]

The compilers should have gotten better over the last 10 years, they should be better at applying return value optimization. We could try reverting this, and measure performance. ~100x performance difference would be painfully obvious. I doubt it will have such an impact.

[blowekamp]

Just a couple notes as I get into this.

1. In SimpleITK to convert VectorImage to/from Images of Vector Types, the underlying image buffer pointer is just copied so that its the same data. This has been an optimization that I have thought should go into the CastImageFilter when the "InPlace" options is enabled. Here is this SimpleITK code: https://github.com/SimpleITK/SimpleITK/blob/main/Code/Common/include/sitkImageConvert.hxx#L31-L70

2. SimpleITK uses VectorImages for multi-component images everywhere. Additionally, it has very good coverage of the filters and checks the results images. I am somewhat surprised of this bug, and I don't believe it is pervasive in the ITK filters based on the testing I have done in SimpleITK with VectorImages.

3. I would like so see a minimal reproducible test cast that does not use a filter. I may get to this as I dig into this issue.

[blowekamp]

Looking at the CastImageFilter:

ITK/Modules/Filtering/ImageFilterBase/include/itkCastImageFilter.hxx

Line 148 in a9e3442

OutputPixelType value{ outputIt.Get() };

That first line use to copy of the VariableLengthVector breaking the reference to the data buffer, but some how it changes to retain a reference to the data buffer.

The compilers should have gotten better over the last 10 years, they should be better at applying return value optimization. We could try reverting this, and measure performance. ~100x performance difference would be painfully obvious. I doubt it will have such an impact.

Per pixel memory allocation with 50+ cores was historically rather slow. Likely small memory concurrent mallocs have improved, but iterative per-pixel memory allocation is still a HUGE performance issue.

[dzenanz]

I have not disputed that, I merely suggested checking whether the compiler is now smart enough to generate same or similar assembly code for both variants of Get().

[blowekamp]

The change to add support to VariableLengthVector to the CastImageFilter occurred here: 87c3fb1

We can quickly address the issue with the filter with the following change:

@@ -145,7 +144,8 @@ CastImageFilter<TInputImage, TOutputImage>::DynamicThreadedGenerateDataDispatche
   ImageScanlineConstIterator inputIt(inputPtr, inputRegionForThread);
   ImageScanlineIterator      outputIt(outputPtr, outputRegionForThread);
 
-  OutputPixelType value{ outputIt.Get() };
+  const auto &outputValue{outputIt.Get()};
+  auto value{ outputValue };
   while (!inputIt.IsAtEnd())

corrected
Please let me know if you can add it to this PR, with the nice test or if I should create a separate PR.

If I have this correct conceptually, an ImageIterator's "Get" method should return a copy to the pixel, while the "Value" method should return a reference. In this case the PixelAccessor for the VectorImage is returning a (modifiable) reference and it is getting modified in the filter. This fix will be more involved, if needed.

[SimonRit]

• In SimpleITK to convert VectorImage to/from Images of Vector Types, the underlying image buffer pointer is just copied so that its the same data. This has been an optimization that I have thought should go into the CastImageFilter when the "InPlace" options is enabled. Here is this SimpleITK code: SimpleITK/SimpleITK@main/Code/Common/include/sitkImageConvert.hxx#L31-L70

I agree, it would be very nice to have this in CastImageFilter.

• SimpleITK uses VectorImages for multi-component images everywhere. Additionally, it has very good coverage of the filters and checks the results images. I am somewhat surprised of this bug, and I don't believe it is pervasive in the ITK filters based on the testing I have done in SimpleITK with VectorImages.

3. I would like so see a minimal reproducible test cast that does not use a filter. I may get to this as I dig into this issue.
Something like that:

#include <itkVectorImage.h>
#include <itkImageRegionIterator.h>

int
main(int argc, char * argv[])
{
  using ImageType = itk::VectorImage<float,2>;
  auto img = ImageType::New();
  img->SetVectorLength(3);
  ImageType::SizeType size = { { 2, 2 } };
  img->SetRegions(size);
  img->AllocateInitialized();
  itk::ImageRegionIterator<ImageType> it(img, img->GetLargestPossibleRegion());
  auto v { it.Get() };
  float c = 0.;
  while(!it.IsAtEnd())
  {
    v.Fill(c++);
    it.Set(v);
    ++it;
  }
  it.GoToBegin();
  while(!it.IsAtEnd())
  {
    std::cout << it.Get() << std::endl;
    ++it;
  }
  return EXIT_SUCCESS;
}

which returns

[3, 3, 3]
[1, 1, 1]
[2, 2, 2]
[3, 3, 3]

[SimonRit]

The change to add support to VariableLengthVector to the CastImageFilter occurred here: 87c3fb1

We can quickly address the issue with the filter with the following change:

@@ -145,7 +146,7 @@ CastImageFilter<TInputImage, TOutputImage>::DynamicThreadedGenerateDataDispatche
   ImageScanlineConstIterator inputIt(inputPtr, inputRegionForThread);
   ImageScanlineIterator      outputIt(outputPtr, outputRegionForThread);
 
-  OutputPixelType value{ outputIt.Get() };
+  itk::NumericTraits<typename OutputPixelType::ValueType> value[componentsPerPixel];

Please let me know if you can add it to this PR, with the nice test or if I should create a separate PR.

If I have this correct conceptually, an ImageIterator's "Get" method should return a copy to the pixel, while the "Value" method should return a reference. In this case the PixelAccessor for the VectorImage is returning a (modifiable) reference and it is getting modified in the filter. This fix will be more involved, if needed.

Done in the latest force-push. If you agree, let's merge this PR to fix CastImageFilter and then I'll open an issue with the minimal example above. BTW, GradientRecursiveGaussianImageFilter was failing with my previous "solution" so the filter might suffer from this iterator problem too.

[blowekamp]

Sorry I posted the wrong diff, this is what I had working:

+  const auto &outputValue{outputIt.Get()};
+  auto value{ outputValue };

It was not an elegant way to force a copy, and I began looking at other constructs.

[SimonRit]

Sorry I only saw this after pushing yet another version. Let me know which version you prefer, they're all fine to me.

[SimonRit]

@blowekamp Your suggestion was not compiling so I force-pushed a different version, see changes here.

[blowekamp]

@SimonRit I timed a couple options and this seems clean and perfrom the best:

@@ -145,12 +159,12 @@ CastImageFilter<TInputImage, TOutputImage>::DynamicThreadedGenerateDataDispatche
   ImageScanlineConstIterator inputIt(inputPtr, inputRegionForThread);
   ImageScanlineIterator      outputIt(outputPtr, outputRegionForThread);
 
-  OutputPixelType value{ outputIt.Get() };
   while (!inputIt.IsAtEnd())
   {
     while (!inputIt.IsAtEndOfLine())
     {
       const InputPixelType & inputPixel = inputIt.Get();
+      OutputPixelType value{outputIt.Get()};
       for (unsigned int k = 0; k < componentsPerPixel; ++k)
       {
         value[k] = static_cast<typename OutputPixelType::ValueType>(inputPixel[k]);

Sorry I didn't have a working diff to start.

[blowekamp]

I updated the branch with the recommendation.

[SimonRit]

I updated the branch with the recommendation.

Thanks for this. I approve my PR 😄 and I have opened #5668 to track the issue.

[blowekamp]

We should consider porting the fix to the 5.4 branch.