diff --git a/Modules/Filtering/ImageGrid/test/itkResampleImageFilterGTest.cxx b/Modules/Filtering/ImageGrid/test/itkResampleImageFilterGTest.cxx index d56c051c6179..b44c7e900bf3 100644 --- a/Modules/Filtering/ImageGrid/test/itkResampleImageFilterGTest.cxx +++ b/Modules/Filtering/ImageGrid/test/itkResampleImageFilterGTest.cxx @@ -19,7 +19,13 @@ // The header file to be tested: #include "itkResampleImageFilter.h" +#include "itkAffineTransform.h" +#include "itkCastImageFilter.h" +#include "itkGaussianInterpolateImageFunction.h" #include "itkImage.h" +#include "itkImageRegionIterator.h" +#include "itkImageRegionIteratorWithIndex.h" +#include "itkStreamingImageFilter.h" // Google Test header file: #include @@ -172,3 +178,66 @@ TEST(ResampleImageFilter, ThrowsOnIncompleteConfiguration) { Expect_ResampleImageFilter_thows_on_incomplete_configuration(128.0); } + + +// Streaming a Gaussian-interpolated resample must reproduce the non-streamed result (issue #1011). +TEST(ResampleImageFilter, StreamingMatchesNonStreamedWithGaussianInterpolator) +{ + constexpr unsigned int Dimension{ 3 }; + using PixelType = float; + using ImageType = itk::Image; + + constexpr unsigned int kNumberOfStreamDivisions{ 8 }; + + const auto input = ImageType::New(); + const auto size = ImageType::SizeType::Filled(32); + input->SetRegions(size); + input->Allocate(); + for (itk::ImageRegionIteratorWithIndex it(input, input->GetLargestPossibleRegion()); !it.IsAtEnd(); ++it) + { + const auto idx = it.GetIndex(); + it.Set(static_cast(idx[0] + idx[1] + idx[2])); + } + + const auto transform = itk::AffineTransform::New(); + transform->Scale(0.9); + + using InterpolatorType = itk::GaussianInterpolateImageFunction; + const auto interpolator = InterpolatorType::New(); + auto sigma = itk::MakeFilled(1.0); + interpolator->SetSigma(sigma); + interpolator->SetAlpha(1.0); + + // CastImageFilter is a streaming-aware no-op; it produces a per-chunk BufferedRegion + // upstream of the resampler, which is required to exercise the interpolator bug. + const auto upstream = itk::CastImageFilter::New(); + upstream->SetInput(input); + + const auto resampler = itk::ResampleImageFilter::New(); + resampler->SetInput(upstream->GetOutput()); + resampler->SetTransform(transform); + resampler->SetInterpolator(interpolator); + resampler->SetSize(size); + + const auto streamer = itk::StreamingImageFilter::New(); + streamer->SetInput(resampler->GetOutput()); + + streamer->SetNumberOfStreamDivisions(1); + ASSERT_NO_THROW(streamer->UpdateLargestPossibleRegion()); + const ImageType::Pointer unstreamed = streamer->GetOutput(); + unstreamed->DisconnectPipeline(); + + input->Modified(); + streamer->SetNumberOfStreamDivisions(kNumberOfStreamDivisions); + ASSERT_NO_THROW(streamer->UpdateLargestPossibleRegion()); + const ImageType::Pointer streamed = streamer->GetOutput(); + streamed->DisconnectPipeline(); + + itk::ImageRegionIterator itU(unstreamed, unstreamed->GetLargestPossibleRegion()); + itk::ImageRegionIterator itS(streamed, streamed->GetLargestPossibleRegion()); + for (; !itU.IsAtEnd() && !itS.IsAtEnd(); ++itU, ++itS) + { + EXPECT_FLOAT_EQ(itU.Get(), itS.Get()); + } + EXPECT_EQ(itU.IsAtEnd(), itS.IsAtEnd()); +}