From 609319760ae5c8c780c4b9be156a4a865ab79439 Mon Sep 17 00:00:00 2001 From: Roman Grothausmann Date: Thu, 7 Feb 2019 16:05:03 +0100 Subject: [PATCH] ENH: Stream ResampleImageFilter for linear transforms Method to compute InputRequestedRegion in itk::ResampleImageFilter 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) Co-authored-by: Matt McCormick --- .../Core/Common/include/itkImageAlgorithm.h | 5 +- .../Core/Common/include/itkImageAlgorithm.hxx | 13 +++-- Modules/Core/Common/test/itkImageTest.cxx | 6 ++- .../include/itkResampleImageFilter.hxx | 54 ++++++++++++++++--- .../ImageGrid/include/itkWarpImageFilter.hxx | 6 ++- .../Filtering/ImageGrid/test/CMakeLists.txt | 2 - .../test/itkResampleImageTest2Streaming.cxx | 17 +++--- .../ImageGrid/test/itkResampleImageTest4.cxx | 6 +-- .../ImageGrid/test/itkResampleImageTest7.cxx | 46 +++++----------- 9 files changed, 97 insertions(+), 58 deletions(-) diff --git a/Modules/Core/Common/include/itkImageAlgorithm.h b/Modules/Core/Common/include/itkImageAlgorithm.h index 09159368840..a3257936ae0 100644 --- a/Modules/Core/Common/include/itkImageAlgorithm.h +++ b/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 + template static typename OutputImageType::RegionType EnlargeRegionOverBox(const typename InputImageType::RegionType & inputRegion, const InputImageType* inputImage, - const OutputImageType* outputImage); + const OutputImageType* outputImage, + const TransformType* transform = nullptr); private: diff --git a/Modules/Core/Common/include/itkImageAlgorithm.hxx b/Modules/Core/Common/include/itkImageAlgorithm.hxx index 0117a04fb15..9ba6e835f0e 100644 --- a/Modules/Core/Common/include/itkImageAlgorithm.hxx +++ b/Modules/Core/Common/include/itkImageAlgorithm.hxx @@ -167,11 +167,12 @@ void ImageAlgorithm::DispatchedCopy( const InputImageType *inImage, } } -template +template 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 diff --git a/Modules/Core/Common/test/itkImageTest.cxx b/Modules/Core/Common/test/itkImageTest.cxx index 8517592ca7f..a0ccd278d85 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,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) || diff --git a/Modules/Filtering/ImageGrid/include/itkResampleImageFilter.hxx b/Modules/Filtering/ImageGrid/include/itkResampleImageFilter.hxx index 67fd2e386e3..558ca603e90 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. @@ -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, diff --git a/Modules/Filtering/ImageGrid/include/itkWarpImageFilter.hxx b/Modules/Filtering/ImageGrid/include/itkWarpImageFilter.hxx index dbe889baf15..70d57ae1c4f 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,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() ) diff --git a/Modules/Filtering/ImageGrid/test/CMakeLists.txt b/Modules/Filtering/ImageGrid/test/CMakeLists.txt index a0c39ea3d41..024a4b5223b 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 diff --git a/Modules/Filtering/ImageGrid/test/itkResampleImageTest2Streaming.cxx b/Modules/Filtering/ImageGrid/test/itkResampleImageTest2Streaming.cxx index 3636a9ab0c0..92bc68d66eb 100644 --- a/Modules/Filtering/ImageGrid/test/itkResampleImageTest2Streaming.cxx +++ b/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; - } diff --git a/Modules/Filtering/ImageGrid/test/itkResampleImageTest4.cxx b/Modules/Filtering/ImageGrid/test/itkResampleImageTest4.cxx index aee8a4897f9..05e13593ed8 100644 --- a/Modules/Filtering/ImageGrid/test/itkResampleImageTest4.cxx +++ b/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( 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; diff --git a/Modules/Filtering/ImageGrid/test/itkResampleImageTest7.cxx b/Modules/Filtering/ImageGrid/test/itkResampleImageTest7.cxx index b0fdd26e7e2..cb7befb62af 100644 --- a/Modules/Filtering/ImageGrid/test/itkResampleImageTest7.cxx +++ b/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; - 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