ENH: streaming for itkResampleImageFilter in case of linear transforms

Squashed commit of the following: commit 2cdc0ee Author: Roman Grothausmann <[email protected]> Date: Tue Feb 5 12:56:06 2019 +0100 Revert "ENH: new test are expected to fail as long as:": This reverts commit e5d421f. because #392 (replaced #84) was merged in bc4e197 commit d62cbd4 Author: Roman Grothausmann <[email protected]> Date: Fri Sep 28 11:55:28 2018 +0200 BUG: use ITK_NULLPTR instead of identityTransform commit 9c0b04c Author: Roman Grothausmann <[email protected]> Date: Thu Sep 27 12:45:50 2018 +0200 ENH: add transform to ImageAlgorithm::EnlargeRegionOverBox commit 5b8eaf5 Author: Roman Grothausmann <[email protected]> Date: Tue Sep 25 14:07:50 2018 +0200 ENH: use ImageAlgorithm::EnlargeRegionOverBox instead commit d6c44ea Author: Roman Grothausmann <[email protected]> Date: Thu Aug 30 13:05:31 2018 +0200 ENH: Method to compute InputRequestedRegion in itkResampleImageFilter based on suggestions from: https://discourse.itk.org/t/why-resampleimagefilter-is-slow/1217/14 https://itk.org/pipermail/insight-users/2015-April/051877.html and code from itkImageAlgorithm of v4.13.1 (based on commit 8510db2)
InsightSoftwareConsortium · Feb 5, 2019 · fb647ca · fb647ca
1 parent bc4e197
commit fb647ca
Show file tree

Hide file tree

Showing 6 changed files with 52 additions and 7 deletions.
diff --git a/Modules/Core/Common/include/itkImageAlgorithm.h b/Modules/Core/Common/include/itkImageAlgorithm.h
@@ -110,9 +110,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);
 

diff --git a/Modules/Core/Common/include/itkImageAlgorithm.hxx b/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 &
 
     using PointType = Point< SpacePrecisionType, OutputImageType::ImageDimension >;
     PointType point;
-    inputImage->TransformContinuousIndexToPhysicalPoint(currentCornerIndex, point);
+    if ( transformPtr == ITK_NULLPTR)
+      {
+      inputImage->TransformContinuousIndexToPhysicalPoint(currentCornerIndex, point);
+      }
+    else
+      {
+      point = transformPtr->TransformPoint(currentCornerIndex);
+      }
     outputImage->TransformPhysicalPointToContinuousIndex(point, corners[count]);
     }
 

diff --git a/Modules/Core/Common/test/itkImageTest.cxx b/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,7 +123,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);

diff --git a/Modules/Filtering/ImageGrid/include/itkResampleImageFilter.hxx b/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.
 
@@ -510,9 +511,38 @@ ResampleImageFilter< TInputImage, TOutputImage, TInterpolatorPrecisionType, TTra
   Superclass::GenerateInputRequestedRegion();
 
   // 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.
   //

diff --git a/Modules/Filtering/ImageGrid/include/itkWarpImageFilter.hxx b/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,8 +412,10 @@ 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(),
+                                                                                         static_cast< TransformType *>(ITK_NULLPTR),
                                                                                          outputPtr,
                                                                                          fieldPtr);
       fieldPtr->SetRequestedRegion( fieldRequestedRegion );

diff --git a/Modules/Filtering/ImageGrid/test/CMakeLists.txt b/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