Custom border extrapolation of ShapedImageNeighborhoodRange by ImageNeighborhoodPixelAccessPolicy

(Niels Dekker) #21

Thanks :grinning:

In the mean time I did figure out, roughly, how any number of run-time arguments could be supported for those pixel access policies that might need them, without any run-time penalty for those that don’t need them. By defining the range class as a variadic template class:

ShapedImageNeighborhoodRange<TImage, TPixelAccessPolicy, TPixelAccessParameters...>

An ITK user could then pass run-time arguments for those pixel access parameters to the constructor of the range. For example:

// For example: 1 run-time parameter, type = int, value = 255.
using RangeType = ShapedImageNeighborhoodRange<ImageType, PixelAccessPolicy, int>;
RangeType range { *image, location, offsets, 255 };

Internally, I think these optional pixel access parameters should then be wrapped in an std::tuple<TPixelAccessParameters...>. This tuple should then be stored as a private data member of the iterator object. Whenever a pixel is accessed, the iterator should pass the tuple to the pixel proxy, which should unwrap it again, and pass the arguments to the policy. Unwrapping/unrolling a tuple isn’t entirely trivial in C++11. It’s slightly easier in C++14 (using std::index_sequence), and possibly still easier in C++17. Anyway, here’s a link that I found helpful:

(Niels Dekker) #22

@phcerdan I just posted a possible implementation of ConstantBoundaryImageNeighborhoodPixelAccessPolicy, which has the same functionality of the policy classes we have been discussing here in July… but this one does allow specifying the constant value at run-time: WIP: Added policy for constant NeighborhoodRange values outside image

While I think it’s fully functional already, I just marked it WIP, because I don’t know if it has priority for the release of ITK 5.0. What do you think?

Eventually it could replace existing use cases of the old itk::ConstantBoundaryCondition, and possibly yield a performance improvement…

(Pablo Hernandez-Cerdan) #23

Hey Niels, great, good stuff :sparkles: I have an itch though, why does the Range class has to have a new template parameter for this? Didn’t you explore the options you mentioned in the past:

What problems did you find with those approaches? It seems to me they would generate a cleaner interface. Basically the user would only need to provide a template with the policy to the Range. Right know the Range class has to provide an optional parameter that would only be used for the case of the ConstantPolicy.

(Niels Dekker) #24

Thank you, @phcerdan, for having looked at WIP: Added policy for constant NeighborhoodRange values outside image. I hope I can clarify this patch a little more.

First of all, the patch still fully supports the original NeighborhoodRange interface: With the patch, it is still not required to specify a ‘pixel access parameter’. When you don’t specify the ‘pixel access parameter’, it is just set to a nullptr by default (and then it is simply ignored).

The policy introduced with the patch (ConstantBoundaryImageNeighborhoodPixelAccessPolicy) allows specifying the 'padding value' at run-time, through this pixel access parameter. My idea is that it would thereby have the same flexibility as the old itk::ConstantBoundaryCondition (as used in combination with the old itk::ConstNeighborhoodIterator).

Example use case:

// Retrieve padding-value at run-time, from user input: 
std::cout << "Please enter 'paddingValue':" << std::endl;
PixelType paddingValue;
std::cin >> paddingValue;

// Using the old itk::ConstNeighborhoodIterator
itk::ConstantBoundaryCondition<ImageType> boundaryCondition;

// The equivalent in ITK 5.0 + this patch.
// Using itk::Experimental::ShapedImageNeighborhoodRange:
using RangeType = ShapedImageNeighborhoodRange<ImageType,
  ConstantBoundaryImageNeighborhoodPixelAccessPolicy<ImageType>, PixelType>;
RangeType range{ *image, location, offsets, paddingValue };

The approaches we discussed in July did not support such a use case. They would just support specifying the padding value at compile-time. That’s why I still wanted to figure out how support for an optional run-time pixel access parameter could be added, without significant extra overhead for the non-parametric use cases!

Did I make myself more clear? :grinning:

(Pablo Hernandez-Cerdan) #25

Thanks Niels, crystal clear, I think I understood correctly the first time, but always appreciate the extra examples! :blush:

So, better than talking, going directly for the code (simplified):

template <typename TImage>
class ConstantBoundaryImageNeighborhoodPixelAccessPolicy final
//......... Appending to all the good stuff already there ...... //
  void SetConstant(const PixelType & inputValue) {m_Constant = inputValue; }

No new template parameter in Range:

template<typename TImage, typename TImageNeighborhoodPixelAccessPolicy >
class ShapedImageNeighborhoodRange final:
// ........ Appending .... //
TImageNeighborhoodPixelAccessPolicy & GetReferencePixelAccessPolicy() { return m_PixelAccessPolicy;}

The usage:

// The equivalent in ITK 5.0 + this patch +  this
// Using itk::Experimental::ShapedImageNeighborhoodRange:
using RangeType = ShapedImageNeighborhoodRange< ImageType,
ConstantBoundaryImageNeighborhoodPixelAccessPolicy<ImageType> >;

RangeType range{ *image, location, offsets}; // same constructor as before
range->GetReferencePixelAccessPolicy()->SetConstant(paddingValue); // we set the value at run time on the policy

Let me know what you think,

(Niels Dekker) #26

Hi Pablo!

I’m glad you appreciate the extra examples. I did test them, actually :grinning: Did you also test the (simplified) code from your repy? If not, I think it would be helpful if you would just try to make it work.

The current ShapedImageNeighborhoodRange class does not have a m_PixelAccessPolicy data member. Only the PixelProxy (which is returned when the user does *it on a Range iterator it) has a m_PixelAccessPolicy.

I fear an extra performance penalty when a PixelAccessPolicy data member would be added to the Range class. But please feel free to give it a try, if you want to!

P.S. In the meantime, I did replace std::nullptr_t by an empty struct, EmptyPixelAccessParameter, as a way to indicate that there is no pixel access parameter specified: Patch Set 4 Hope you like that better!

(Niels Dekker) #27

@matt.mccormick, @phcerdan I just did a performance comparison between “old-school” and new range-based neighborhood iteration, both doing zero-padding (constant boundary). Using Visual C++ 2017 64-bit Release, I observed that the range-based neighborhood iteration implemented by can be more than twice as fast as the “old-school” iteration :grinning:

The new range-based neighborhood iteration was tested by the following function. It calculates a number, sumOfNeighbors, just so that the optimizer cannot eliminate all the work.

unsigned TestRangeBasedIteration(const ImageType& image, const SizeType& radius)
  using namespace itk::Experimental;
  using RangeType = ShapedImageNeighborhoodRange<const 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)

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

And the old-school neighborhood iteration was tested by the following function:

unsigned TestOldSchoolIteration(const ImageType& image, const SizeType& radius)
  using NeighborhoodIteratorType = itk::ConstNeighborhoodIterator<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);
  return sumOfNeighbors;

The full test is at

Please check if you find it a fair comparison!

(Niels Dekker) #28

For the record, it was only after I wrote those words, that I found a way to significantly improve the performance of the “old-school” iteration (using Visual C++ 2017 Release): PERF: Made ConstNeighborhoodIterator::GetPixel(i) much faster

(Matt McCormick) #29

@Niels_Dekker Well done! :clap: :sparkles: To add more data :vulcan_salute: to the Visual C++ 2017 Release results, here are some metrics created on Linux with GCC 7.3.0 and ICC (the Intel Compiler) These results are created for this ultrasound application and your Gist.

They are very interesting! All we can conclude is that neither implementation is consistently faster across compilers! :astonished:

That said, the potential of these new image range classes to be used with STL algorithms is quite compelling because of potential future benefits:

Granted that we would have to figure out a good way to control data synchronization in the GPU case.

+1 to this API – setting the constant value separately from the range initialization is easier to follow, in my opinion. That said, this is not implemented, and maybe we could implement it in the future if it is really compelling. As an itk::Experimental namespaced feature, we should keep experimenting with nice things!

Let’s merge the patch, which @hjmjohnson has kindly migrated to GitHub! :octopus:

(Niels Dekker) #30

@matt.mccormick Thanks for your extensive performance tests! Very interesting! Clearly the new range-based iterators aren’t always faster (anymore) than the “old-school” NeighborhoodIterator classes, especially now that the old iterator classes have received some major PERF updates. :thinking: However, the new iterators will potentially be much faster than the old ones when the neighborhood location has to be reset frequently (for example, by means of SetLocation). As I observed with GaussianDerivativeImageFunction. Also I believe that the new range-based iterators could be significantly faster than the old ones on arbitrary neighborhood shapes (rather than just rectangular shapes). :grinning:

Sorry but I’m not really convinced. Allowing to modify the constant padding value after range initialization might make things harder to follow, in my opinion. Especially when the padding value is modified after retrieving iterators from the range (using begin() and end()). It would yield the question whether such a modification would still affect these iterators. It seems clearer to me to specify the constant padding value as part of the range initialization so that it would then clearly remain constant.

(There’s also the design issue that current ShapedImageNeighborhoodRange class does not have a policy data member, as I explained at Custom border extrapolation of ShapedImageNeighborhoodRange by ImageNeighborhoodPixelAccessPolicy)

Anyway, please go ahead if you would like to merge I haven’t yet got the GitHub/ITK entirely up and running but if @hjmjohnson would like to do the merge, please do!

(Pablo Hernandez-Cerdan) #31

The issue I see with the current API is that the optional parameter of a particular policy should be encoded in that policy, not in the Range that can use any policy. Right now Range<TPolicy, TOptionalParameterForConstantPolicy>

What would happen if another policy comes by with another parameter type, we will need to modify Range template parameter?

What do you think about this? Do you think another API would be possible?

I suggested the Policy as data member because it solves this mangling, it would be less clean to use, we will have to first construct the range and then modify the policy, but I don’t see where the performance will be reduced (there is only one way to know though).

In any case, this is good stuff! :tophat:

(Niels Dekker) #32

@matt.mccormick Is your ultrasound application just zero-padding outside the image borders, or does it use some specific (non-zero) constant value?

If it is just zero-padding, you could possibly use some ZeroBoundaryImageNeighborhoodPixelAccessPolicy (to be implemented), instead of the ConstantBoundaryImageNeighborhoodPixelAccessPolicy class that I’m proposing at PR 105.

Such a ZeroBoundaryImageNeighborhoodPixelAccessPolicy would not need the optional extra TOptionalPixelAccessParameter parameter from the NeighborhoodRange, so it could work independently of the NeighborhoodRange code changes from PR 105.

I’m saying this also because @phcerdan and I are still in discussion about the design of PR 105. A ZeroBoundaryImageNeighborhoodPixelAccessPolicy class could still be implemented based on the currently merged version of the NeighborhoodRange.

(Matt McCormick) #33

Yes, it is zero in this case.

Very cool! :cool: Thanks for implementing and contributing this! :clap:

In general, if we create these type-based optimized implementations, we are showcasing the power of C++.

(Niels Dekker) #34

I did not implement ZeroBoundaryImageNeighborhoodPixelAccessPolicy, I only implemented the more general ConstantBoundaryImageNeighborhoodPixelAccessPolicy, but that one is still under discussion. If you would like to, you could implement ZeroBoundaryImageNeighborhoodPixelAccessPolicy locally within your ultrasound application, by just copying the code from itkConstantBoundary*.h and simply removing its m_Constant data member, replacing it by a hard-code 0 (or Zero), within GetPixelValue. I would be interested to know if hard-coding the zero constant would improve the performance…