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.

Author: blowekamp

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<TInputImage, TOutputImage>::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<typename OutputPixelType::ValueType>(inputPixel[k]);
@@ -163,6 +167,32 @@ CastImageFilter<TInputImage, TOutputImage>::DynamicThreadedGenerateDataDispatche
     inputIt.NextLine();
     outputIt.NextLine();
   }
+#else
+
+  auto inputRange = ImageRegionRange<const TInputImage>(*inputPtr, inputRegionForThread);
+  auto outputRange = ImageRegionRange<TOutputImage>(*outputPtr, outputRegionForThread);
+
+  auto       inputIt = inputRange.begin();
+  auto       outputIt = outputRange.begin();
+  const auto inputEnd = inputRange.end();
+
+
+  OutputPixelType    outputPixel{ *outputIt };
+  const unsigned int componentsPerPixel = NumericTraits<OutputPixelType>::GetLength(outputPixel);
+  while (inputIt != inputEnd)
+  {
+    InputPixelType inputPixel = *inputIt;
+    for (unsigned int k = 0; k < componentsPerPixel; ++k)
+    {
+      outputPixel[k] = static_cast<typename OutputPixelType::ValueType>(inputPixel[k]);
+    }
+    *outputIt = outputPixel;
+
+    ++inputIt;
+    ++outputIt;
+  }
+
+#endif
 }
 
 } // end namespace itk

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 <typename TInputImage, typename TOutputImage>
+void
+TimeImageCast()
+{
+  std::cout << "Timing cast from " << GetCastTypeName<typename TInputImage::PixelType>() << " to "
+            << GetCastTypeName<typename TOutputImage::PixelType>() << " ... ";
+
+  // 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<TInputImage, TOutputImage>;
+  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<itk::Image<itk::Vector<float, 3>, 3>, itk::Image<itk::RGBPixel<double>, 3>>();
+  TimeImageCast<itk::Image<itk::Vector<float, 3>, 3>, itk::VectorImage<double, 3>>();
+
+  TimeImageCast<itk::VectorImage<float, 3>, itk::Image<itk::RGBPixel<double>, 3>>();
+
   std::cout << std::endl;
   if (!success)
   {