From 903ca572f99bacbde8cebe76dc87ed1c5d6c6968 Mon Sep 17 00:00:00 2001 From: Bradley Lowekamp Date: Mon, 1 Dec 2025 12:04:34 -0500 Subject: [PATCH 1/2] WIP: Comparing performance in CastImageFilter Comparing perfomance with using the cast image filter to convert between different vector pixel type. There are two option to compare: - Using ImageRegionRange appears to have performnace benifits, compared to the ImageScanelineIterators. - Using NumericTraits::GetLength is nearly a constant /compile time expression which is also performant. The challenge is choosing input vs output pixel type when only one has a compile time size. --- .../include/itkCastImageFilter.hxx | 32 +++++++++++- .../test/itkCastImageFilterTest.cxx | 50 +++++++++++++++++++ 2 files changed, 81 insertions(+), 1 deletion(-) diff --git a/Modules/Filtering/ImageFilterBase/include/itkCastImageFilter.hxx b/Modules/Filtering/ImageFilterBase/include/itkCastImageFilter.hxx index 16753ef1a3d..91a5cf469ea 100644 --- a/Modules/Filtering/ImageFilterBase/include/itkCastImageFilter.hxx +++ b/Modules/Filtering/ImageFilterBase/include/itkCastImageFilter.hxx @@ -20,6 +20,7 @@ #include "itkProgressReporter.h" #include "itkImageAlgorithm.h" +#include "itkImageRegionRange.h" namespace itk { @@ -139,18 +140,21 @@ CastImageFilter::DynamicThreadedGenerateDataDispatche this->CallCopyOutputRegionToInputRegion(inputRegionForThread, outputRegionForThread); - const unsigned int componentsPerPixel = outputPtr->GetNumberOfComponentsPerPixel(); +#if 0 // Define the iterators ImageScanlineConstIterator inputIt(inputPtr, inputRegionForThread); ImageScanlineIterator outputIt(outputPtr, outputRegionForThread); + + const unsigned int componentsPerPixel = inputPtr->GetNumberOfComponentsPerPixel(); while (!inputIt.IsAtEnd()) { while (!inputIt.IsAtEndOfLine()) { const InputPixelType & inputPixel = inputIt.Get(); OutputPixelType value{ outputIt.Get() }; + for (unsigned int k = 0; k < componentsPerPixel; ++k) { value[k] = static_cast(inputPixel[k]); @@ -163,6 +167,32 @@ CastImageFilter::DynamicThreadedGenerateDataDispatche inputIt.NextLine(); outputIt.NextLine(); } +#else + + auto inputRange = ImageRegionRange(*inputPtr, inputRegionForThread); + auto outputRange = ImageRegionRange(*outputPtr, outputRegionForThread); + + auto inputIt = inputRange.begin(); + auto outputIt = outputRange.begin(); + const auto inputEnd = inputRange.end(); + + + OutputPixelType outputPixel{ *outputIt }; + const unsigned int componentsPerPixel = NumericTraits::GetLength(outputPixel); + while (inputIt != inputEnd) + { + InputPixelType inputPixel = *inputIt; + for (unsigned int k = 0; k < componentsPerPixel; ++k) + { + outputPixel[k] = static_cast(inputPixel[k]); + } + *outputIt = outputPixel; + + ++inputIt; + ++outputIt; + } + +#endif } } // end namespace itk diff --git a/Modules/Filtering/ImageFilterBase/test/itkCastImageFilterTest.cxx b/Modules/Filtering/ImageFilterBase/test/itkCastImageFilterTest.cxx index fb6ec61a3f0..79963c71f4c 100644 --- a/Modules/Filtering/ImageFilterBase/test/itkCastImageFilterTest.cxx +++ b/Modules/Filtering/ImageFilterBase/test/itkCastImageFilterTest.cxx @@ -27,6 +27,7 @@ #include "itkFloatingPointExceptions.h" #include "itkRGBPixel.h" #include "itkRGBAPixel.h" +#include "itkTimeProbe.h" namespace itk {} @@ -444,6 +445,47 @@ TestVectorImageCast3() } +template +void +TimeImageCast() +{ + std::cout << "Timing cast from " << GetCastTypeName() << " to " + << GetCastTypeName() << " ... "; + + // Create a 512^3 image + auto image = TInputImage::New(); + + typename TInputImage::SizeType size; + size.Fill(512); + typename TInputImage::RegionType region{ size }; + image->SetRegions(region); + image->SetNumberOfComponentsPerPixel(3); + image->AllocateInitialized(); + + + using CastImageFilterType = itk::CastImageFilter; + auto castImageFilter = CastImageFilterType::New(); + castImageFilter->SetInput(image); + castImageFilter->SetNumberOfWorkUnits(1); // Single-threaded for consistent timing + + // Warm-up run + castImageFilter->Update(); + + // Timed run + itk::TimeProbe timer; + timer.Start(); + castImageFilter->Modified(); + castImageFilter->Update(); + timer.Stop(); + + const double totalPixels = size.CalculateProductOfElements(); + const double timeInSeconds = timer.GetMean(); + const double megapixelsPerSecond = (totalPixels / 1e6) / timeInSeconds; + + std::cout << megapixelsPerSecond << " Mpixels/sec (" << timeInSeconds << " seconds)" << std::endl; +} + + int itkCastImageFilterTest(int, char *[]) { @@ -477,6 +519,14 @@ itkCastImageFilterTest(int, char *[]) success &= TestVectorImageCast2(); success &= TestVectorImageCast3(); + + std::cout << std::endl; + std::cout << "Performance timing tests (512^3 images):" << std::endl; + TimeImageCast, 3>, itk::Image, 3>>(); + TimeImageCast, 3>, itk::VectorImage>(); + + TimeImageCast, itk::Image, 3>>(); + std::cout << std::endl; if (!success) { From a7204038ed2d4ab1386e654e27f858c1e5262150 Mon Sep 17 00:00:00 2001 From: Bradley Lowekamp Date: Mon, 1 Dec 2025 14:38:42 -0500 Subject: [PATCH 2/2] WIP: Add performance test to compare copy methods --- .../ImageFilterBase/test/CMakeLists.txt | 8 + .../test/itkImageIterationPerformanceTest.cxx | 416 ++++++++++++++++++ 2 files changed, 424 insertions(+) create mode 100644 Modules/Filtering/ImageFilterBase/test/itkImageIterationPerformanceTest.cxx diff --git a/Modules/Filtering/ImageFilterBase/test/CMakeLists.txt b/Modules/Filtering/ImageFilterBase/test/CMakeLists.txt index d0b79ae781f..366c85921c7 100644 --- a/Modules/Filtering/ImageFilterBase/test/CMakeLists.txt +++ b/Modules/Filtering/ImageFilterBase/test/CMakeLists.txt @@ -6,6 +6,7 @@ set( itkVectorNeighborhoodOperatorImageFilterTest.cxx itkMaskNeighborhoodOperatorImageFilterTest.cxx itkCastImageFilterTest.cxx + itkImageIterationPerformanceTest.cxx ) # Disable optimization on the tests below to avoid possible @@ -59,6 +60,13 @@ itk_add_test( ITKImageFilterBaseTestDriver itkCastImageFilterTest ) +itk_add_test( + NAME itkImageIterationPerformanceTest + COMMAND + ITKImageFilterBaseTestDriver + itkImageIterationPerformanceTest + 256 +) set(ITKImageFilterBaseGTests itkGeneratorImageFilterGTest.cxx) creategoogletestdriver(ITKImageFilterBase "${ITKImageFilterBase-Test_LIBRARIES}" "${ITKImageFilterBaseGTests}") diff --git a/Modules/Filtering/ImageFilterBase/test/itkImageIterationPerformanceTest.cxx b/Modules/Filtering/ImageFilterBase/test/itkImageIterationPerformanceTest.cxx new file mode 100644 index 00000000000..42be4f51aa0 --- /dev/null +++ b/Modules/Filtering/ImageFilterBase/test/itkImageIterationPerformanceTest.cxx @@ -0,0 +1,416 @@ +/*========================================================================= + * + * Copyright NumFOCUS + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * https://www.apache.org/licenses/LICENSE-2.0.txt + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + * + *=========================================================================*/ + +#include +#include "itkImage.h" +#include "itkVectorImage.h" +#include "itkVector.h" +#include "itkRGBPixel.h" +#include "itkImageRegionRange.h" +#include "itkImageScanlineIterator.h" +#include "itkImageScanlineConstIterator.h" +#include "itkImageRegionIteratorWithIndex.h" +#include "itkImageRegionConstIterator.h" +#include "itkTimeProbe.h" +#include "itkNumericTraits.h" +#include "itkTestingMacros.h" +#include + + +// Helper function to initialize an image with random values +template +typename TImage::Pointer +CreateAndInitializeImage(const typename TImage::SizeType & size, unsigned int numberOfComponentsPerPixel = 0) +{ + auto image = TImage::New(); + typename TImage::RegionType region{ size }; + image->SetRegions(region); + if (numberOfComponentsPerPixel > 0) + { + image->SetNumberOfComponentsPerPixel(numberOfComponentsPerPixel); + } + image->Allocate(); + + // Initialize with simple pattern (pixel index-based) + using PixelType = typename TImage::PixelType; + unsigned int count = 0; + + const unsigned int length = image->GetNumberOfComponentsPerPixel(); + + itk::ImageRegionIteratorWithIndex it(image, region); + for (it.GoToBegin(); !it.IsAtEnd(); ++it) + { + PixelType pixel{ it.Get() }; + for (unsigned int k = 0; k < length; ++k) + { + pixel[k] = static_cast::ValueType>(count + k); + } + it.Set(pixel); + ++count; + } + + return image; +} + + +// Method 1: ImageScanlineIterator approach +template +void +CopyScanlineIterator(const TInputImage * inputPtr, TOutputImage * outputPtr) +{ + using InputPixelType = typename TInputImage::PixelType; + using OutputPixelType = typename TOutputImage::PixelType; + using ImageScanlineConstIterator = itk::ImageScanlineConstIterator; + using ImageScanlineIterator = itk::ImageScanlineIterator; + + const typename TOutputImage::RegionType outputRegion = outputPtr->GetRequestedRegion(); + typename TInputImage::RegionType inputRegion = outputRegion; + + ImageScanlineConstIterator inputIt(inputPtr, inputRegion); + ImageScanlineIterator outputIt(outputPtr, outputRegion); + + const unsigned int componentsPerPixel = inputPtr->GetNumberOfComponentsPerPixel(); + while (!inputIt.IsAtEnd()) + { + while (!inputIt.IsAtEndOfLine()) + { + const InputPixelType & inputPixel = inputIt.Get(); + OutputPixelType value(outputIt.Get()); + for (unsigned int k = 0; k < componentsPerPixel; ++k) + { + value[k] = static_cast(inputPixel[k]); + } + outputIt.Set(value); + + ++inputIt; + ++outputIt; + } + inputIt.NextLine(); + outputIt.NextLine(); + } +} + + +// Method 1b: ImageScanlineIterator approach using NumericTraits::GetLength() +template +void +CopyScanlineIteratorNumericTraits(const TInputImage * inputPtr, TOutputImage * outputPtr) +{ + using InputPixelType = typename TInputImage::PixelType; + using OutputPixelType = typename TOutputImage::PixelType; + using ImageScanlineConstIterator = itk::ImageScanlineConstIterator; + using ImageScanlineIterator = itk::ImageScanlineIterator; + + const typename TOutputImage::RegionType outputRegion = outputPtr->GetRequestedRegion(); + typename TInputImage::RegionType inputRegion = outputRegion; + + ImageScanlineConstIterator inputIt(inputPtr, inputRegion); + ImageScanlineIterator outputIt(outputPtr, outputRegion); + + unsigned int componentsPerPixel = itk::NumericTraits::GetLength(outputIt.Get()); + while (!inputIt.IsAtEnd()) + { + while (!inputIt.IsAtEndOfLine()) + { + const InputPixelType & inputPixel = inputIt.Get(); + + OutputPixelType value{ outputIt.Get() }; + for (unsigned int k = 0; k < componentsPerPixel; ++k) + { + value[k] = static_cast(inputPixel[k]); + } + outputIt.Set(value); + + ++inputIt; + ++outputIt; + } + inputIt.NextLine(); + outputIt.NextLine(); + } +} + + +// Method 2: ImageRegionRange approach +template +void +CopyImageRegionRange(const TInputImage * inputPtr, TOutputImage * outputPtr) +{ + using InputPixelType = typename TInputImage::PixelType; + using OutputPixelType = typename TOutputImage::PixelType; + + const typename TOutputImage::RegionType outputRegion = outputPtr->GetRequestedRegion(); + typename TInputImage::RegionType inputRegion = outputRegion; + + auto inputRange = itk::ImageRegionRange(*inputPtr, inputRegion); + auto outputRange = itk::ImageRegionRange(*outputPtr, outputRegion); + + auto inputIt = inputRange.begin(); + auto outputIt = outputRange.begin(); + const auto inputEnd = inputRange.end(); + + const unsigned int componentsPerPixel = inputPtr->GetNumberOfComponentsPerPixel(); + while (inputIt != inputEnd) + { + const InputPixelType & inputPixel = *inputIt; + OutputPixelType outputPixel{ *outputIt }; + for (unsigned int k = 0; k < componentsPerPixel; ++k) + { + outputPixel[k] = static_cast(inputPixel[k]); + } + *outputIt = outputPixel; + + ++inputIt; + ++outputIt; + } +} + + +// Method 2b: ImageRegionRange approach using NumericTraits::GetLength() +template +void +CopyImageRegionRangeNumericTraits(const TInputImage * inputPtr, TOutputImage * outputPtr) +{ + using InputPixelType = typename TInputImage::PixelType; + using OutputPixelType = typename TOutputImage::PixelType; + + const typename TOutputImage::RegionType outputRegion = outputPtr->GetRequestedRegion(); + typename TInputImage::RegionType inputRegion = outputRegion; + + auto inputRange = itk::ImageRegionRange(*inputPtr, inputRegion); + auto outputRange = itk::ImageRegionRange(*outputPtr, outputRegion); + + auto inputIt = inputRange.begin(); + auto outputIt = outputRange.begin(); + const auto inputEnd = inputRange.end(); + + const unsigned int componentsPerPixel = itk::NumericTraits::GetLength(*outputIt); + while (inputIt != inputEnd) + { + const InputPixelType & inputPixel = *inputIt; + OutputPixelType outputPixel{ *outputIt }; + for (unsigned int k = 0; k < componentsPerPixel; ++k) + { + outputPixel[k] = static_cast(inputPixel[k]); + } + *outputIt = outputPixel; + ++inputIt; + ++outputIt; + } +} + + +// Method 2c: ImageRegionRange approach - BUGGY VERSION demonstrating incorrect pixel reuse +// NOTE: This method has a BUG - outputPixel is initialized once and reused, causing incorrect results +// when the pixel type holds a reference to the image buffer (modifies first pixel repeatedly) +template +void +CopyImageRegionRangeNumericTraitsBuggy(const TInputImage * inputPtr, TOutputImage * outputPtr) +{ + using InputPixelType = typename TInputImage::PixelType; + using OutputPixelType = typename TOutputImage::PixelType; + + const typename TOutputImage::RegionType outputRegion = outputPtr->GetRequestedRegion(); + typename TInputImage::RegionType inputRegion = outputRegion; + + auto inputRange = itk::ImageRegionRange(*inputPtr, inputRegion); + auto outputRange = itk::ImageRegionRange(*outputPtr, outputRegion); + + auto inputIt = inputRange.begin(); + auto outputIt = outputRange.begin(); + const auto inputEnd = inputRange.end(); + + const unsigned int componentsPerPixel = itk::NumericTraits::GetLength(*outputIt); + OutputPixelType outputPixel{ *outputIt }; + while (inputIt != inputEnd) + { + const InputPixelType & inputPixel = *inputIt; + for (unsigned int k = 0; k < componentsPerPixel; ++k) + { + outputPixel[k] = static_cast(inputPixel[k]); + } + *outputIt = outputPixel; + ++inputIt; + ++outputIt; + } +} + + +// Helper function to time a single method +template +void +TimeMethod(const std::string & methodName, + TCopyFunction copyFunc, + const TInputImage * inputImage, + typename TOutputImage::Pointer & outputImage, + double totalPixels, + typename TOutputImage::Pointer referenceImage = nullptr) +{ + // Allocate output image + outputImage = TOutputImage::New(); + outputImage->SetRegions(inputImage->GetLargestPossibleRegion()); + outputImage->SetNumberOfComponentsPerPixel(inputImage->GetNumberOfComponentsPerPixel()); + outputImage->Allocate(); + + // Warm-up run + copyFunc(inputImage, outputImage.GetPointer()); + + // Timed run + itk::TimeProbe timer; + timer.Start(); + copyFunc(inputImage, outputImage.GetPointer()); + timer.Stop(); + + const double time = timer.GetMean(); + const double mpixels = (totalPixels / 1e6) / time; + std::cout << std::left << std::setw(23) << methodName << mpixels << " Mpixels/sec (" << time << " seconds)" + << std::endl; + + // Verify correctness if reference provided + if (referenceImage) + { + bool match = true; + itk::ImageRegionConstIterator itRef(referenceImage, referenceImage->GetLargestPossibleRegion()); + itk::ImageRegionConstIterator it(outputImage, outputImage->GetLargestPossibleRegion()); + for (itRef.GoToBegin(), it.GoToBegin(); !itRef.IsAtEnd(); ++itRef, ++it) + { + if (itRef.Get() != it.Get()) + { + match = false; + break; + } + } + if (!match) + { + std::cout << " WARNING: Results differ from reference!" << std::endl; + } + } +} + + +// Performance testing function +template +void +TimeIterationMethods(const typename TInputImage::SizeType & size, const std::string & description) +{ + std::cout << "\n=== " << description << " ===" << std::endl; + + // Create and initialize input image + auto inputImage = CreateAndInitializeImage(size, 3); + + const double totalPixels = size.CalculateProductOfElements(); + + typename TOutputImage::Pointer referenceImage; + typename TOutputImage::Pointer outputImage; + + // Test Method 1: Scanline Iterator - serves as reference + TimeMethod("Scanline Iterator:", + CopyScanlineIterator, + inputImage.GetPointer(), + referenceImage, + totalPixels); + + // Test Method 2: ImageRegionRange + TimeMethod("ImageRegionRange:", + CopyImageRegionRange, + inputImage.GetPointer(), + outputImage, + totalPixels, + referenceImage); + + // Test Method 1b: Scanline Iterator with NumericTraits + TimeMethod("Scanline NumericTraits:", + CopyScanlineIteratorNumericTraits, + inputImage.GetPointer(), + outputImage, + totalPixels, + referenceImage); + + // Test Method 2b: ImageRegionRange with NumericTraits + TimeMethod("Range NumericTraits:", + CopyImageRegionRangeNumericTraits, + inputImage.GetPointer(), + outputImage, + totalPixels, + referenceImage); + + // Test Method 2c: ImageRegionRange with NumericTraits - buggy version + TimeMethod("Range NT (BUGGY):", + CopyImageRegionRangeNumericTraitsBuggy, + inputImage.GetPointer(), + outputImage, + totalPixels, + referenceImage); + + std::cout << "All methods validated." << std::endl; +} + + +int +itkImageIterationPerformanceTest(int argc, char * argv[]) +{ + std::cout << "CTEST_FULL_OUTPUT" << std::endl; + + // Test different image sizes and types + constexpr unsigned int Dimension = 3; + + // Get image size from command line or use default + unsigned int size = 64; + if (argc == 2) + { + size = std::stoi(argv[1]); + } + else + { + if (argc < 2) + { + std::cerr << "Missing parameters." << std::endl; + std::cerr << "Usage: " << itkNameOfTestExecutableMacro(argv) << " imageSize" << std::endl; + return EXIT_FAILURE; + } + } + + const auto imageSize = itk::Size::Filled(size); + + std::ostringstream oss; + oss << "Image Size: " << imageSize; + + std::string description = oss.str(); + + // Test 1: Image to Image + std::cout << "\n--- Test Suite 1: Image> to Image> ---" << std::endl; + TimeIterationMethods, Dimension>, itk::Image, Dimension>>( + imageSize, description); + + // Test 2: VectorImage to Image + std::cout << "\n--- Test Suite 2: VectorImage to Image> ---" << std::endl; + TimeIterationMethods, itk::Image, Dimension>>(imageSize, + description); + + // Test 3: Image to VectorImage + std::cout << "\n--- Test Suite 3: Image> to VectorImage ---" << std::endl; + TimeIterationMethods, Dimension>, itk::VectorImage>(imageSize, + description); + + // Test 4: VectorImage to VectorImage + std::cout << "\n--- Test Suite 4: VectorImage to VectorImage ---" << std::endl; + TimeIterationMethods, itk::VectorImage>(imageSize, description); + std::cout << "\n=====================================" << std::endl; + std::cout << "All tests completed successfully." << std::endl; + + return EXIT_SUCCESS; +}