WIP: streaming for resample image filter in case of linear transforms

Modules/Core/Common/include/itkImageAlgorithm.h

@@ -138,9 +138,10 @@ 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 TransformType* transformPtr,
                        const InputImageType* inputImage,
                        const OutputImageType* outputImage);
 

Modules/Core/Common/include/itkImageAlgorithm.hxx

@@ -167,9 +167,10 @@ 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 TransformType* transformPtr,
                                      const InputImageType* inputImage,
                                      const OutputImageType* outputImage)
 {
@@ -217,7 +218,14 @@ ImageAlgorithm::EnlargeRegionOverBox(const typename InputImageType::RegionType &
 
     typedef Point< SpacePrecisionType, OutputImageType::ImageDimension > PointType;
     PointType point;
-    inputImage->TransformContinuousIndexToPhysicalPoint(currentCornerIndex, point);
+    if ( transformPtr == ITK_NULLPTR)
+      {
+      inputImage->TransformContinuousIndexToPhysicalPoint(currentCornerIndex, point);
+      }
+    else
+      {
+      point = transformPtr->TransformPoint(currentCornerIndex);
+      }
     outputImage->TransformPhysicalPointToContinuousIndex(point, corners[count]);
     }
 

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* [] )
 {
@@ -121,7 +122,10 @@ 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(),
+                                                                          static_cast< TransformType *>(ITK_NULLPTR),
                                                                           image.GetPointer(),
                                                                           imageRef.GetPointer());
   Image::IndexType correctIndex; correctIndex.Fill(0);

Modules/Filtering/ImageGrid/include/itkResampleImageFilter.hxx

@@ -26,6 +26,7 @@
 #include "itkImageScanlineIterator.h"
 #include "itkSpecialCoordinatesImage.h"
 #include "itkDefaultConvertPixelTraits.h"
+#include "itkImageAlgorithm.h"
 
 namespace itk
 {
@@ -492,10 +493,39 @@ ResampleImageFilter< TInputImage, TOutputImage, TInterpolatorPrecisionType, TTra
     }
 
   // Get pointers to the input and output
-  InputImagePointer inputPtr  =
-    const_cast< TInputImage * >( this->GetInput() );
+  InputImageType * inputPtr =
+    const_cast< InputImageType * >( this->GetInput() );
 
-  // Determining the actual input region is non-trivial, especially
+
+  // 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 *outputPtr = this->GetOutput();
+  // Get the input transform
+  const TransformType *transformPtr = 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 && transformPtr->GetTransformCategory() == TransformType::Linear )
+    {
+    typename TInputImage::RegionType inputRequestedRegion;
+    inputRequestedRegion = ImageAlgorithm::EnlargeRegionOverBox(outputPtr->GetRequestedRegion(),
+	transformPtr,
+	outputPtr,
+	inputPtr);
+
+    inputPtr->SetRequestedRegion( inputRequestedRegion );
+    return;
+    }
+
+  // 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.
   //

Modules/Filtering/ImageGrid/include/itkWarpImageFilter.hxx

@@ -27,6 +27,8 @@
 #include "itkProgressReporter.h"
 #include "itkContinuousIndex.h"
 #include "itkMath.h"
+#include "itkTransform.h"
+
 namespace itk
 {
 /**
@@ -438,8 +440,10 @@ WarpImageFilter< TInputImage, TOutputImage, TDisplacementField >
     else
       {
       typedef typename TDisplacementField::RegionType DisplacementRegionType;
+      typedef itk::Transform< SpacePrecisionType, OutputImageType::ImageDimension, OutputImageType::ImageDimension > TransformType;
 
       DisplacementRegionType fieldRequestedRegion = ImageAlgorithm::EnlargeRegionOverBox(outputPtr->GetRequestedRegion(),
+                                                                                         static_cast< TransformType *>(ITK_NULLPTR),
                                                                                          outputPtr,
                                                                                          fieldPtr);
       fieldPtr->SetRequestedRegion( fieldRequestedRegion );

Modules/Filtering/ImageGrid/test/Baseline/ResampleImageTest2.mha.sha512

@@ -0,0 +1 @@
+6dcda44ce36e39fb6a792e26c8cd93830bab8edbe9f33385dc59b24bd6b40c99bbadf9b6ada967b1d1b6624e986f1161960d86db63b2d6b71da58eec0a2b911c

Modules/Filtering/ImageGrid/test/Baseline/ResampleImageTest2NearestExtrapolate.mha.sha512

@@ -0,0 +1 @@
+bbc1a32a82e88334df11c79f54ba4260ea291f8316c882d6360e7546739eb92109a169cc07abe32250a4e12e49a4bb8bad22c45580ccb78c4ee482fa68675e4f

Modules/Filtering/ImageGrid/test/CMakeLists.txt

@@ -45,10 +45,12 @@ itkWrapPadImageTest.cxx
 itkMirrorPadImageTest.cxx
 itkResampleImageTest.cxx
 itkResampleImageTest2.cxx
+itkResampleImageTest2s.cxx
 itkResampleImageTest3.cxx
 itkResampleImageTest4.cxx
 itkResampleImageTest5.cxx
 itkResampleImageTest6.cxx
+itkResampleImageTest7.cxx
 itkResamplePhasedArray3DSpecialCoordinatesImageTest.cxx
 itkPushPopTileImageFilterTest.cxx
 itkShrinkImageStreamingTest.cxx
@@ -223,6 +225,22 @@ itk_add_test(NAME itkResampleImageTest2
                           ${ITK_TEST_OUTPUT_DIR}/ResampleImageTest2b.png
                           ${ITK_TEST_OUTPUT_DIR}/ResampleImageTest2c.png
                           ${ITK_TEST_OUTPUT_DIR}/ResampleImageTest2d.png)
+itk_add_test(NAME itkResampleImageTest2s
+      COMMAND ITKImageGridTestDriver
+    --compare DATA{Baseline/ResampleImageTest2.mha}
+              ${ITK_TEST_OUTPUT_DIR}/ResampleImageTest2a.mha
+    --compare DATA{Baseline/ResampleImageTest2.mha}
+              ${ITK_TEST_OUTPUT_DIR}/ResampleImageTest2b.mha
+    --compare DATA{Baseline/ResampleImageTest2NearestExtrapolate.mha}
+              ${ITK_TEST_OUTPUT_DIR}/ResampleImageTest2c.mha
+    --compare DATA{Baseline/ResampleImageTest2NearestExtrapolate.mha}
+              ${ITK_TEST_OUTPUT_DIR}/ResampleImageTest2d.mha
+    itkResampleImageTest2s DATA{${ITK_DATA_ROOT}/Input/cthead1.png}
+                          DATA{${ITK_DATA_ROOT}/Input/circle.png}
+                          ${ITK_TEST_OUTPUT_DIR}/ResampleImageTest2a.mha
+                          ${ITK_TEST_OUTPUT_DIR}/ResampleImageTest2b.mha
+                          ${ITK_TEST_OUTPUT_DIR}/ResampleImageTest2c.mha
+                          ${ITK_TEST_OUTPUT_DIR}/ResampleImageTest2d.mha)
 itk_add_test(NAME itkResampleImageTest3
       COMMAND ITKImageGridTestDriver
     --compare DATA{Baseline/ResampleImageTest3.png}
@@ -241,6 +259,8 @@ itk_add_test(NAME itkResampleImageTest6
     --compare DATA{Baseline/ResampleImageTest6.png}
               ${ITK_TEST_OUTPUT_DIR}/ResampleImageTest6.png
     itkResampleImageTest6 10 ${ITK_TEST_OUTPUT_DIR}/ResampleImageTest6.png)
+itk_add_test(NAME itkResampleImageTest7
+      COMMAND ITKImageGridTestDriver itkResampleImageTest7)
 itk_add_test(NAME itkResamplePhasedArray3DSpecialCoordinatesImageTest
       COMMAND ITKImageGridTestDriver itkResamplePhasedArray3DSpecialCoordinatesImageTest)
 itk_add_test(NAME itkPushPopTileImageFilterTest