Measures performance of ConstNeighborhoodIterator::GetPixel(i,isInBounds), related to http://review.source.kitware.com/#/c/23826 by Bradley Lowekamp

ConstNeighborhoodIteratorGetPixelWithIsInBoundsPerf.cxx

/*
Measures performance of "old-school" ConstNeighborhoodIterator::GetPixel(i, isInBounds)
Related to the patch "PERF: Improve ConstNeighborhoodIterator::GetPixel(i, isInBounds )",
http://review.source.kitware.com/#/c/23826/ by Bradley Lowekamp

Niels Dekker, LKEB, Leiden University Medical Center, 2018

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
*/

#include "itkConstNeighborhoodIterator.h"
#include "itkImage.h"

#include <chrono>
#include <iostream>

namespace
{
  using PixelType = unsigned;
  using ImageType = itk::Image<PixelType>;
  using SizeType = ImageType::SizeType;
  using IndexType = ImageType::IndexType;

  template<typename TImage>
  typename TImage::Pointer CreateImage(const unsigned sizeX, const unsigned sizeY)
  {
    const auto image = TImage::New();
    const typename TImage::SizeType imageSize = { { sizeX , sizeY } };
    image->SetRegions(imageSize);
    image->Allocate();
    return image;
  }


  // Creates a test image, filled with a sequence of natural numbers, 1, 2, 3, ..., N.
  template<typename TImage>
  typename TImage::Pointer CreateImageFilledWithSequenceOfNaturalNumbers(const unsigned sizeX, const unsigned sizeY)
  {
    using PixelType = typename TImage::PixelType;
    const auto image = CreateImage<TImage>(sizeX, sizeY);

    const unsigned numberOfPixels = sizeX * sizeY;

    PixelType* const bufferPointer = image->GetBufferPointer();

    for (unsigned i = 0; i < numberOfPixels; ++i)
    {
      bufferPointer[i] = static_cast<typename TImage::PixelType>(i + 1);
    }
    return image;
  }


  unsigned TestOldSchoolIteration(const ImageType& image, const SizeType& radius)
  {
    const auto region = image.GetBufferedRegion();

    itk::ConstNeighborhoodIterator<ImageType> neighborhoodIterator(radius, &image, region);
    const auto numberOfNeigbors = neighborhoodIterator.Size();
    unsigned sumOfNeighbors = 0;
    bool isInBounds = true;

    while (!neighborhoodIterator.IsAtEnd())
    {
      for (itk::SizeValueType i = 0; i < numberOfNeigbors; ++i)
      {
        bool isPixelInBounds;
        sumOfNeighbors += neighborhoodIterator.GetPixel(i, isPixelInBounds);
        isInBounds = isInBounds && isPixelInBounds;
      }
      ++neighborhoodIterator;
    }
    // Do 'something' with isInBounds, so that it cannot be removed by the optimizer
    return sumOfNeighbors + (isInBounds ? 1 : 0);
  }
}


int main()
{
#ifdef NDEBUG
  const ImageType::SizeType size = { { 26000 , 2600 } };
#else
  const ImageType::SizeType size = { { 80, 70 } };
#endif
  const auto image = CreateImageFilledWithSequenceOfNaturalNumbers<ImageType>(size[0], size[1]);
  const SizeType radius = { { 1, 1 } };

  for (int i = 0; i < 5; ++i)
  {
    using namespace std::chrono;
    const auto timePointBefore = high_resolution_clock::now();
    const auto result = TestOldSchoolIteration(*image, radius);
    const auto timePointAfter = high_resolution_clock::now();
    const auto durationSeconds =
      duration_cast<duration<double>>(timePointAfter - timePointBefore);
    std::cout
      << " Duration: " << durationSeconds.count() << " seconds"
      << " Sum: " << result
      << std::endl;
  }
}