Compares performance of "old-school" and new range-based neighborhood iteration using http://review.source.kitware.com/#/c/23795/7

CompareOldAndRangeBasedConstantBoundary.cxx

/*
Compares performance of "old-school" neighborhood iteration, using
itk::ConstantBoundaryCondition, and range-based neighborhood iteration,
using the patch "WIP: Added policy for constant NeighborhoodRange values
outside image", http://review.source.kitware.com/#/c/23795/7.

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 "itkConstantBoundaryCondition.h"
#include "itkConstantBoundaryImageNeighborhoodPixelAccessPolicy.h"
#include "itkConstNeighborhoodIterator.h"
#include "itkImage.h"
#include "itkImageNeighborhoodOffsets.h"
#include "itkShapedImageNeighborhoodRange.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;
  }


  template<unsigned VDimension>
  void GoToNextPixelIndex(itk::Index<VDimension>& index, const itk::Size<VDimension>& size)
  {
    for (unsigned i = 0; i < VDimension; ++i)
    {
      auto& indexValue = index[i];

      ++indexValue;

      if (static_cast<itk::SizeValueType>(indexValue) < size[i])
      {
        break;
      }
      indexValue = 0;
    }
  }

  unsigned TestRangeBasedIteration(const ImageType& image, const SizeType& radius)
  {
    using namespace itk::Experimental;
    using RangeType = ShapedImageNeighborhoodRange<const ImageType,
      ConstantBoundaryImageNeighborhoodPixelAccessPolicy<ImageType>>;

    const auto region = image.GetBufferedRegion();
    const auto imageSize = region.GetSize();
    const auto offsets = GenerateRectangularImageNeighborhoodOffsets(radius);

    RangeType range{ image, IndexType{}, offsets };
    IndexType index{};
    unsigned sumOfNeighbors = 0;

    for (auto i = region.GetNumberOfPixels(); i > 0; --i)
    {
      range.SetLocation(index);

      for (const PixelType neighbor : range)
      {
        sumOfNeighbors += neighbor;
      }
      GoToNextPixelIndex(index, imageSize);
    }
    return sumOfNeighbors;
  }


  unsigned TestOldSchoolIteration(const ImageType& image, const SizeType& radius)
  {
    using NeighborhoodIteratorType = itk::ConstNeighborhoodIterator<ImageType,
      itk::ConstantBoundaryCondition<ImageType>>;

    const auto region = image.GetBufferedRegion();
    const auto imageSize = region.GetSize();

    NeighborhoodIteratorType neighborhoodIterator(radius, &image, region);
    const auto numberOfNeigbors = neighborhoodIterator.Size();
    unsigned sumOfNeighbors = 0;

    while (!neighborhoodIterator.IsAtEnd())
    {
      for (itk::SizeValueType i = 0; i < numberOfNeigbors; ++i)
      {
        sumOfNeighbors += neighborhoodIterator.GetPixel(i);
      }
      ++neighborhoodIterator;
    }
    return sumOfNeighbors;
  }


  template<typename TFunc>
  void ProfileFunc(TFunc func)
  {
    using namespace std::chrono;
    const auto timePointBefore = high_resolution_clock::now();
    const auto result = func();
    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;
  }
}


int main(int argc, char**)
{
#ifdef NDEBUG
  const ImageType::SizeType size = { { 26000 , 2600 } };
#else
  const ImageType::SizeType size = { { 8, 7 } };
#endif
  const auto image = CreateImageFilledWithSequenceOfNaturalNumbers<ImageType>(size[0], size[1]);
  const auto radiusValue = static_cast<unsigned>(argc);
  const SizeType radius = { {radiusValue, radiusValue} };

  for (int i = 0; i < 5; ++i)
  {
    std::cout << "Range-based neighborhood iteration:" << std::endl;
    ProfileFunc([=] { return TestRangeBasedIteration(*image, radius); });
    std::cout << "Old-school neighborhood iteration:" << std::endl;
    ProfileFunc([=] { return TestOldSchoolIteration(*image, radius); });
  }
}