Skip to content

Instantly share code, notes, and snippets.

@N-Dekker
Last active November 7, 2018 17:19
Show Gist options
  • Save N-Dekker/eb94ba817b2a16e75d5e7c3fe54ec66f to your computer and use it in GitHub Desktop.
Save N-Dekker/eb94ba817b2a16e75d5e7c3fe54ec66f to your computer and use it in GitHub Desktop.
Compares performance of "old-school" and new range-based neighborhood iteration using http://review.source.kitware.com/#/c/23795/7
/*
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); });
}
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment