Streaming for itkResampleImageFilter in case of linear transforms

Modules/Core/Common/include/itkImageAlgorithm.h

@@ -110,11 +110,12 @@ struct ImageAlgorithm
    * the physical space covered by the input
    * region of the input image
    */
-  template<typename InputImageType, typename OutputImageType>
+  template<typename InputImageType, typename OutputImageType, typename TransformType>
   static typename OutputImageType::RegionType
   EnlargeRegionOverBox(const typename InputImageType::RegionType & inputRegion,
                        const InputImageType* inputImage,
-                       const OutputImageType* outputImage);
+                       const OutputImageType* outputImage,
+                       const TransformType* transform = nullptr);
 
 private:
 

Modules/Core/Common/include/itkImageAlgorithm.hxx

@@ -167,11 +167,12 @@ void ImageAlgorithm::DispatchedCopy( const InputImageType *inImage,
     }
 }
 
-template<typename InputImageType, typename OutputImageType>
+template<typename InputImageType, typename OutputImageType, typename TransformType>
 typename OutputImageType::RegionType
 ImageAlgorithm::EnlargeRegionOverBox(const typename InputImageType::RegionType & inputRegion,
                                      const InputImageType* inputImage,
-                                     const OutputImageType* outputImage)
+                                     const OutputImageType* outputImage,
+                                     const TransformType* transform)
 {
   typename OutputImageType::RegionType outputRegion;
 
@@ -217,8 +218,12 @@ ImageAlgorithm::EnlargeRegionOverBox(const typename InputImageType::RegionType &
 
     using PointType = Point< SpacePrecisionType, OutputImageType::ImageDimension >;
     PointType point;
-    inputImage->TransformContinuousIndexToPhysicalPoint(currentCornerIndex, point);
-    outputImage->TransformPhysicalPointToContinuousIndex(point, corners[count]);
+    inputImage->TransformContinuousIndexToPhysicalPoint( currentCornerIndex, point );
+    if( transform != nullptr )
+      {
+      point = transform->TransformPoint( point );
+      }
+    outputImage->TransformPhysicalPointToContinuousIndex( point, corners[count] );
     }
 
   // Compute a rectangular region from the vector of corner indexes

Modules/Core/Common/test/itkImageTest.cxx

@@ -20,6 +20,7 @@
 #include "itkImage.h"
 #include "itkFixedArray.h"
 #include "itkImageAlgorithm.h"
+#include "itkTransform.h"
 
 int itkImageTest(int, char* [] )
 {
@@ -122,9 +123,12 @@ int itkImageTest(int, char* [] )
   regionRef.SetSize(sizeRef);
   imageRef->SetRegions(regionRef);
 
+  typedef itk::Transform< double, Image::ImageDimension, Image::ImageDimension > TransformType;
+
   Image::RegionType boxRegion = itk::ImageAlgorithm::EnlargeRegionOverBox(image->GetLargestPossibleRegion(),
                                                                           image.GetPointer(),
-                                                                          imageRef.GetPointer());
+                                                                          imageRef.GetPointer(),
+                                                                          static_cast< TransformType *>(nullptr) );
   Image::IndexType correctIndex; correctIndex.Fill(0);
   Image::SizeType correctSize; correctSize.Fill(3);
   if( !(boxRegion.GetIndex() == correctIndex) ||

Modules/Filtering/ImageGrid/include/itkResampleImageFilter.hxx

@@ -26,6 +26,7 @@
 #include "itkImageScanlineIterator.h"
 #include "itkSpecialCoordinatesImage.h"
 #include "itkDefaultConvertPixelTraits.h"
+#include "itkImageAlgorithm.h"
 
 #include <type_traits>  // For is_same.
 
@@ -506,17 +507,58 @@ void
 ResampleImageFilter< TInputImage, TOutputImage, TInterpolatorPrecisionType, TTransformPrecisionType >
 ::GenerateInputRequestedRegion()
 {
-  // Call the superclass' implementation of this method
-  Superclass::GenerateInputRequestedRegion();
-
   // Get pointers to the input and output
-  InputImagePointer inputPtr  = const_cast< TInputImage * >( this->GetInput() );
+  InputImageType * input  = const_cast< InputImageType * >( this->GetInput() );
+
+  // Check whether the input or the output is a
+  // SpecialCoordinatesImage. If either are, then we cannot use the
+  // fast path since index mapping will definitely not be linear.
+  typedef SpecialCoordinatesImage< PixelType, ImageDimension >           OutputSpecialCoordinatesImageType;
+  typedef SpecialCoordinatesImage< InputPixelType, InputImageDimension > InputSpecialCoordinatesImageType;
+
+  const bool isSpecialCoordinatesImage = ( dynamic_cast< const InputSpecialCoordinatesImageType * >( this->GetInput() )
+       || dynamic_cast< const OutputSpecialCoordinatesImageType * >( this->GetOutput() ) );
+
+  const OutputImageType *output = this->GetOutput();
+  // Get the input transform
+  const TransformType *transform = this->GetTransform();
+
+  // Check whether we can use upstream streaming for resampling. Upstream streaming
+  // can be used if the transformation is linear. Transform respond
+  // to the IsLinear() call.
+  if ( !isSpecialCoordinatesImage && transform->GetTransformCategory() == TransformType::Linear )
+    {
+    typename TInputImage::RegionType inputRequestedRegion;
+    inputRequestedRegion = ImageAlgorithm::EnlargeRegionOverBox(output->GetRequestedRegion(),
+        output,
+        input,
+        transform);
+
+    const typename TInputImage::RegionType inputLargestRegion( input->GetLargestPossibleRegion() );
+    if( inputLargestRegion.IsInside( inputRequestedRegion.GetIndex() ) || inputLargestRegion.IsInside( inputRequestedRegion.GetUpperIndex() ) )
+      {
+      // Input requested region is partially outside the largest possible region.
+      //   or
+      // Input requested region is completely inside the largest possible region.
+      input->SetRequestedRegion( inputRequestedRegion );
+      }
+    else if( inputRequestedRegion.IsInside( inputLargestRegion ) )
+      {
+      // Input requested region completely surrounds the largest possible region.
+      input->SetRequestedRegion( inputLargestRegion );
+      }
+    else
+      {
+      // Input requested region is completely outside the largest possible region. Do not set the requested region in this case.
+      }
+    return;
+    }
 
-  // Determining the actual input region is non-trivial, especially
+  // Otherwise, determining the actual input region is non-trivial, especially
   // when we cannot assume anything about the transform being used.
   // So we do the easy thing and request the entire input image.
   //
-  inputPtr->SetRequestedRegionToLargestPossibleRegion();
+  input->SetRequestedRegionToLargestPossibleRegion();
 }
 
 template< typename TInputImage,

Modules/Filtering/ImageGrid/include/itkWarpImageFilter.hxx

@@ -27,6 +27,8 @@
 #include "itkProgressReporter.h"
 #include "itkContinuousIndex.h"
 #include "itkMath.h"
+#include "itkTransform.h"
+
 namespace itk
 {
 template< typename TInputImage, typename TOutputImage, typename TDisplacementField >
@@ -410,10 +412,12 @@ WarpImageFilter< TInputImage, TOutputImage, TDisplacementField >
     else
       {
       using DisplacementRegionType = typename TDisplacementField::RegionType;
+      using TransformType = itk::Transform< SpacePrecisionType, OutputImageType::ImageDimension, OutputImageType::ImageDimension >;
 
       DisplacementRegionType fieldRequestedRegion = ImageAlgorithm::EnlargeRegionOverBox(outputPtr->GetRequestedRegion(),
                                                                                          outputPtr,
-                                                                                         fieldPtr);
+                                                                                         fieldPtr,
+                                                                                         static_cast< TransformType *>(nullptr));
       fieldPtr->SetRequestedRegion( fieldRequestedRegion );
       }
     if ( !fieldPtr->VerifyRequestedRegion() )

Modules/Filtering/ImageGrid/test/CMakeLists.txt

@@ -276,7 +276,6 @@ itk_add_test(NAME itkResampleImageTest2Streaming
                           ${ITK_TEST_OUTPUT_DIR}/ResampleImageTest2bStreaming.mha
                           ${ITK_TEST_OUTPUT_DIR}/ResampleImageTest2cStreaming.mha
                           ${ITK_TEST_OUTPUT_DIR}/ResampleImageTest2dStreaming.mha)
-set_tests_properties(itkResampleImageTest2Streaming PROPERTIES WILL_FAIL TRUE) # should fail as long as itkResampleFilter does not support streaming for linear transforms GH PR #82
 itk_add_test(NAME itkResampleImageTest3
       COMMAND ITKImageGridTestDriver
     --compare DATA{Baseline/ResampleImageTest3.png}
@@ -297,7 +296,6 @@ itk_add_test(NAME itkResampleImageTest6
     itkResampleImageTest6 10 ${ITK_TEST_OUTPUT_DIR}/ResampleImageTest6.png)
 itk_add_test(NAME itkResampleImageTest7
       COMMAND ITKImageGridTestDriver itkResampleImageTest7)
-set_tests_properties(itkResampleImageTest7 PROPERTIES WILL_FAIL TRUE) # should fail as long as itkResampleFilter does not support streaming for linear transforms GH PR #82
 itk_add_test(NAME itkResamplePhasedArray3DSpecialCoordinatesImageTest
       COMMAND ITKImageGridTestDriver itkResamplePhasedArray3DSpecialCoordinatesImageTest)
 itk_add_test(NAME itkPushPopTileImageFilterTest

Modules/Filtering/ImageGrid/test/itkResampleImageTest2Streaming.cxx

@@ -50,8 +50,11 @@ class NonlinearAffineTransform:
   /** Run-time type information (and related methods). */
   itkTypeMacro(NonlinearAffineTransform, AffineTransform);
 
-  /** Override this. See test below. */
-    bool IsLinear() const override { return false; }
+  /** Override this so not linear. See test below. */
+  typename itk::TransformBaseTemplate< TCoordRepType >::TransformCategoryType GetTransformCategory() const override
+    {
+    return itk::TransformBaseTemplate< TCoordRepType >::UnknownTransformCategory;
+    }
 };
 }
 
@@ -149,11 +152,12 @@ int itkResampleImageTest2Streaming( int argc, char * argv [] )
   // This will use ResampleImageFilter::LinearThreadedGenerateData().
   std::cout << "Test with normal AffineTransform." << std::endl;
   writer1->SetNumberOfStreamDivisions(8); //split into 8 pieces for streaming.
+  monitor->ClearPipelineSavedInformation();
   TRY_EXPECT_NO_EXCEPTION( writer1->Update() );
 
-  // this verifies that the pipeline was executed as expected along
-  // with correct region propagation and output information
-  if (!monitor->VerifyAllInputCanStream(8))
+  // Note: We will only request the input 4 times because that last sampled
+  // chunk is completely outside the input and not requested
+  if( !monitor->VerifyInputFilterExecutedStreaming(4) )
     {
     std::cerr << "Streaming failed to execute as expected!" << std::endl;
     std::cerr << monitor;
@@ -178,7 +182,7 @@ int itkResampleImageTest2Streaming( int argc, char * argv [] )
   TRY_EXPECT_NO_EXCEPTION( writer2->Update() );
 
   // check that streaming is not possible for non-linear case
-  if (!monitor->VerifyAllInputCanStream(1))
+  if ( monitor->VerifyInputFilterExecutedStreaming(8) )
     {
     std::cerr << "Streaming succeeded for non-linear transform which should not be the case!" << std::endl;
     std::cerr << monitor;
@@ -205,5 +209,4 @@ int itkResampleImageTest2Streaming( int argc, char * argv [] )
 
   std::cout << "Test passed." << std::endl;
   return EXIT_SUCCESS;
-
 }

Modules/Filtering/ImageGrid/test/itkResampleImageTest4.cxx

@@ -56,8 +56,7 @@ int itkResampleImageTest4(int argc, char * argv [] )
   ImageRegionType region;
   region.SetSize ( size );
   region.SetIndex( index );
-  image->SetLargestPossibleRegion( region );
-  image->SetBufferedRegion( region );
+  image->SetRegions( region );
   image->Allocate();
 
   auto newDims = static_cast<unsigned int>( 64*scaling );
@@ -120,7 +119,8 @@ int itkResampleImageTest4(int argc, char * argv [] )
   // Run the resampling filter
   itk::TimeProbe clock;
   clock.Start();
-  resample->Update();
+  std::cout << "Input: " << image << std::endl;
+  resample->UpdateLargestPossibleRegion();
   clock.Stop();
 
   std::cout << "Resampling from " << size << " to " << osize << " took " << clock.GetMean() << " s" << std::endl;

Modules/Filtering/ImageGrid/test/itkResampleImageTest7.cxx

@@ -20,7 +20,6 @@
 
 #include "itkAffineTransform.h"
 #include "itkResampleImageFilter.h"
-#include "itkPipelineMonitorImageFilter.h"
 #include "itkStreamingImageFilter.h"
 #include "itkTestingMacros.h"
 #include "itkMath.h"
@@ -54,8 +53,7 @@ int itkResampleImageTest7( int , char *[] )
   ImageRegionType region;
   region.SetSize ( size );
   region.SetIndex( index );
-  image->SetLargestPossibleRegion( region );
-  image->SetBufferedRegion( region );
+  image->SetRegions( region );
   image->Allocate();
 
   // Fill image with a ramp
@@ -81,15 +79,10 @@ int itkResampleImageTest7( int , char *[] )
     itk::ResampleImageFilter< ImageType, ImageType >::New();
 
   EXERCISE_BASIC_OBJECT_METHODS( resample, ResampleImageFilter, ImageToImageFilter );
-
-  using MonitorFilter = itk::PipelineMonitorImageFilter<ImageType>;
-  MonitorFilter::Pointer monitor = MonitorFilter::New();
-  monitor->SetInput( image );
-
   resample->SetInterpolator(interp);
 
-  resample->SetInput( monitor->GetOutput() );
-  TEST_SET_GET_VALUE( monitor->GetOutput(), resample->GetInput() );
+  resample->SetInput( image );
+  TEST_SET_GET_VALUE( image.GetPointer(), resample->GetInput() );
 
   resample->SetSize( size );
   TEST_SET_GET_VALUE( size, resample->GetSize() );
@@ -125,36 +118,25 @@ int itkResampleImageTest7( int , char *[] )
   // Run the resampling filter without streaming, i.e. 1 StreamDivisions
   numStreamDiv= 1; // do not split, i.e. do not stream
   streamer->SetNumberOfStreamDivisions(numStreamDiv);
-  TRY_EXPECT_NO_EXCEPTION( streamer->Update() );
+  TRY_EXPECT_NO_EXCEPTION( streamer->UpdateLargestPossibleRegion() );
 
-  if (!monitor->VerifyAllInputCanStream(numStreamDiv))
-    {
-    std::cerr << "Avoiding streaming failed to execute as expected!" << std::endl;
-    std::cerr << monitor;
-    std::cerr << "Test failed." << std::endl;
-    return EXIT_FAILURE;
-    }
-  monitor->ClearPipelineSavedInformation();
-
-  ImagePointerType outputNoSDI= streamer->GetOutput(); // save output for later comparison
+  ImagePointerType outputNoSDI = streamer->GetOutput(); // save output for later comparison
   outputNoSDI->DisconnectPipeline(); // disconnect to create new output
 
   // Run the resampling filter with streaming
+  image->Modified();
   numStreamDiv= 8; // split into numStream pieces for streaming.
-  image->Modified(); // enforce re-execution even though nothing of filter changed
   streamer->SetNumberOfStreamDivisions(numStreamDiv);
-  TRY_EXPECT_NO_EXCEPTION( streamer->Update() );
+  TRY_EXPECT_NO_EXCEPTION( streamer->UpdateLargestPossibleRegion() );
 
-  if (!monitor->VerifyAllInputCanStream(numStreamDiv))
-    {
-    std::cerr << "Streaming failed to execute as expected!" << std::endl;
-    std::cerr << monitor;
-    std::cerr << "Test failed." << std::endl;
-    return EXIT_FAILURE;
-    }
-  monitor->ClearPipelineSavedInformation();
+  // Verify that we only requested a smaller chunk when streaming
+  const ImageRegionType finalRequestedRegion( image->GetRequestedRegion() );
+  TEST_SET_GET_VALUE( 0, finalRequestedRegion.GetIndex(0) );
+  TEST_SET_GET_VALUE( 49, finalRequestedRegion.GetIndex(1) );
+  TEST_SET_GET_VALUE( 59, finalRequestedRegion.GetSize(0) );
+  TEST_SET_GET_VALUE( 10, finalRequestedRegion.GetSize(1) );
 
-  ImagePointerType outputSDI= streamer->GetOutput();
+  ImagePointerType outputSDI = streamer->GetOutput();
   outputSDI->DisconnectPipeline();
 
   itk::ImageRegionIterator<ImageType>