From fb647ca094af89cd6f7708d803dc396473c4466b Mon Sep 17 00:00:00 2001 From: Roman Grothausmann Date: Tue, 5 Feb 2019 13:06:53 +0100 Subject: [PATCH] ENH: streaming for itkResampleImageFilter in case of linear transforms Squashed commit of the following: commit 2cdc0ee3f7993069907150c38f6c9b47758695ed Author: Roman Grothausmann Date: Tue Feb 5 12:56:06 2019 +0100 Revert "ENH: new test are expected to fail as long as:": This reverts commit e5d421f5d755cdfa208b1f1a6b052244716fe91e. because #392 (replaced #84) was merged in bc4e1978304df14672ff913c96f97a8abbf49160 commit d62cbd4baf18005722b3ff18f240af28c6b3a50a Author: Roman Grothausmann Date: Fri Sep 28 11:55:28 2018 +0200 BUG: use ITK_NULLPTR instead of identityTransform commit 9c0b04cd6a9628cd92b64a14ae7729a4eb745fc2 Author: Roman Grothausmann Date: Thu Sep 27 12:45:50 2018 +0200 ENH: add transform to ImageAlgorithm::EnlargeRegionOverBox commit 5b8eaf5fdf88b9b7f74dde87d587cffe1c389606 Author: Roman Grothausmann Date: Tue Sep 25 14:07:50 2018 +0200 ENH: use ImageAlgorithm::EnlargeRegionOverBox instead commit d6c44ea46606b68c5cc25e0e9904e1a5d774ea42 Author: Roman Grothausmann 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 8510db2df6606837126486c47aba58f0c278b801) --- .../Core/Common/include/itkImageAlgorithm.h | 3 +- .../Core/Common/include/itkImageAlgorithm.hxx | 12 +++++-- Modules/Core/Common/test/itkImageTest.cxx | 4 +++ .../include/itkResampleImageFilter.hxx | 34 +++++++++++++++++-- .../ImageGrid/include/itkWarpImageFilter.hxx | 4 +++ .../Filtering/ImageGrid/test/CMakeLists.txt | 2 -- 6 files changed, 52 insertions(+), 7 deletions(-) diff --git a/Modules/Core/Common/include/itkImageAlgorithm.h b/Modules/Core/Common/include/itkImageAlgorithm.h index 091593688409..848a63c106b6 100644 --- 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 + template 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 index 0117a04fb157..17a1c2558946 100644 --- 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 +template 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 index 8517592ca7f9..7a592fa99395 100644 --- 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 index 67fd2e386e3e..a28958c6dbd2 100644 --- 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 // 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 index dbe889baf155..2f0a913efdda 100644 --- 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 index a0c39ea3d414..024a4b5223b5 100644 --- 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