itk::NeighborhoodRange: a new class for efficient modern C++ style iteration

(Pablo Hernandez-Cerdan) #61

I would like it indeed, but not sure if I will find time until next week :roll_eyes:

(Niels Dekker) #62

Thanks Pablo! That would be cool! Take your time, no problem.

By the way, I would especially like to see a larger return type for itk::GeometryUtilities::Factorial. Maybe even std::uintmax_t (the largest available unsigned integer type). Because those factorials can become quite large! For example, 12! is equal to 479001600, that’s already the largest factorial that can fit in a 32-bit long integer. With a 64-bit return type, we can go up to 20! :star_struck:

(Niels Dekker) #63

@phcerdan Update: I just learned that it appears better not to call a Factorial function, when calculating the binomial coefficient, “n over m”, as I did before to calculate the number of offsets N needed for an N-connected neigborhood (TopologicalConnectivityImageNeighborhoodShape). My colleague Baldur van Lew suggested me to calculate the binomial coefficient as described by Walter (June 23 2017), at:

When using this clever binom function, the value N of an N-connected neighborhood can be calculated for images that have far more than 20 dimensions!!! Without integer overflow! :smiley: So although it might be nice to make the existing itk::GeometryUtilities::Factorial function constexpr, it’s certainly not necessary for the implementation this shape class.

By the way, the hyperrectangular shape class is still “under review”, it has its first “+1” from @dzenanz: See also: A shape class for ShapedImageNeighborhoodRange and ShapedNeighborhoodIterator As you know, I would like the hyperrectangular shape class to be accepted before proposing an N-connected shape class.

(Niels Dekker) #64

For your information,

The ShapedImageNeighborhoodRange iterators now have all the operators to support random access iteration, rather than just bidirectional iteration:

Which means that you can do all those random access operations on ShapedImageNeighborhoodRange iterators, as you would probably expect already: it1-it2, it+=n, it+n, n+it, it-=n, it-n, it[n], it1<it2, it1<=it2, it1>it, it1>=it2.

It also means that you could pass such iterators to a C++11 standard library (std) algorithm that requires Random Access Iterators, for example, std::sort and std::nth_element. However, please note that there is no 100% guarantee that passing such proxy iterators to an std algorithm actually works, with C++11. In practice, it appears to work just well, as tested at itkShapedImageNeighborhoodRangeGTest.cxx However, the C++ standard library does not yet officially support proxy iterators.

Eric Niebler extensively investigated the topic of proxy iterators, and how to add proxy iterator support to the C++ standard library. I asked him if we’re just lucky that the ShapedImageNeighborhoodRange iterators happen to work well, with std algorithms. Here is his reply:

I’ll say that you’re just getting lucky. Now, it’s more than just luck; STL implementors try to support proxy iterators, I’m told, because of bastard types like vector<bool>. But without the proper iter_swap and iter_move customization points I describe in my blog posts, the support can be half-baked at best. YMMV.

See For more details, please check his blog!

Kind regards, Niels

(Niels Dekker) #65

For your information, ShapedImageNeighborhoodRange now also supports reverse iteration over the neighbor pixels, by offering member functions rbegin(), rend(), crbegin(), and crend(). As Matt just merged:

I just found out how useful reverse iterators can be in a different context (using the reverse iterators from itk::Offset<N> when testing that a sequence of offsets is in colexicographic order). Hopefully, these reverse-image-neighborhood-iterators can be just as useful :slightly_smiling_face:

(Niels Dekker) #66

The patch I just submitted should finally make ShapedImageNeighborhoodRange compatible with the Range Library from (which is proposed to be added to the C++ standard library). So for example, when neighborhoodRange is an object of type ShapedImageNeighborhoodRange<ImageType>, the patch allows ITK users to write:

#include <range/v3/numeric/accumulate.hpp> // from the Range Library 
sum = ranges::accumulate(neighborhoodRange, 0.0);

Which is equivalent to:

#include <numeric>  // from the Standard C++ Library 
sum = std::accumulate(neighborhoodRange.begin(), neighborhoodRange.end(), 0.0);

(I’m not suggesting to use the Range Library within the implementation of ITK itself already, but of course, it is available to ITK users.)

The fix basically just adds a default-constructor to the iterator types of ShapedImageNeighborhoodRange<ImageType>. Thereby it fixes the Visual C++ 2015 compile errors (_error C2672: 'operator _surrogate_func’: no matching overloaded function found), that I previously saw when trying out the example, ranges::accumulate(neighborhoodRange, 0.0), using the Range Library for Visual Studio.

Please review the patch: BUG: Added default-constructor iterator ShapedImageNeighborhoodRange (

(Niels Dekker) #67

Yesterday evening I did a talk about the neighborhood range, at a C++ meetup in Amsterdam, which was very well received. Here are the slides:

@matt.mccormick , @blowekamp, @phcerdan there’s a picture of you on one of the slides. :slight_smile:

(Pablo Hernandez-Cerdan) #68

Haha, nice! Great slides :smiley:

(Matt McCormick) #69

Absolutely fantastic slides @Niels_Dekker!! :sparkler::heart::dizzy:

What do you think about an image region range class? Is this feasible?

(Niels Dekker) #70

Thanks Matt! I haven’t entirely studied itk::ImageRegionIterator but it seems to me that a modern C++ image region range class is feasible. An image region is always rectangular, right? Is it safe to assume that an image region is always entirely within its associated image?

(Dženan Zukić) #71

If an image region is not entirely withing image’s buffered region, the iterator’s constructor should raise an exception. I think that is the convention.

(Niels Dekker) #72

Thanks, @dzenanz I actually meant to ask whether (customizable?) border extrapolation/boundary conditions would be necessary for image region iteration.

I’m not sure if it’s a good convention to throw a C++ exception in Release mode, when the caller violates a precondition of the called function. Especially now that we have the noexcept keyword! But that’s a different discussion! Error handling is still a hot topic in C++. Herb Sutter just did a new proposal:

(Dženan Zukić) #73

A plain image iterator does not need extrapolation. Only neighborhood iterators need that.

(Pablo Hernandez-Cerdan) #74


(Matt McCormick) #75

Cool. The image region is rectangular. As @dzenanz noted, we can check at construction time to ensure the iterated region is within the buffered region of the image. The challenge, which makes the image region range not just a trivial wrapper, is that the iterated image region is a subset of the buffered image region.

(Dženan Zukić) #76

I don’t know whether it would be possible to stuff into the new iterator, but it would be good if it were amenable to auto-vectorization in a similar way to ImageScanlineIterator.

(Niels Dekker) #77

Summary of the current status of ShapedImageNeighborhoodRange

First of all, thanks to @matt.mccormick for merging pull request 212 onto the master branch of ITK. With this merge, users of ShapedImageNeighborhoodRange<TImage, TImageNeighborhoodPixelAccessPolicy> will now have a choice between two policies to extrapolate pixel values outside the image boundaries:

Other policies could be implemented as well, for example, a wrapping (circular) policy, like itk::PeriodicBoundaryCondition, or a mirror-reflecting (symmetric) policy. But I think this could wait until there’s an actual user request for such policies.

ShapedImageNeighborhoodRange supports any arbitrary user specified neighborhood shape, but two shapes are already directly available to the user:

A diamond shape (Von Neumann neighborhood) could still be added later, as well. But again, I think this could wait until there’s an actual user request.

Compared to the “old-school” itk::ConstNeighborhoodIterator, ShapedImageNeighborhoodRange is especially fast when the neighborhood location frequently has to be set to an arbitrary index. itk::GaussianDerivativeImageFunction::EvaluateAtIndex(index) became much faster, once it started using ShapedImageNeighborhoodRange internally, instead itk::ConstNeighborhoodIterator. It seems likely to me that the performance of other EvaluateAtIndex member functions, from other itk::ImageFunction-derived classes, could also be improved significantly by using ShapedImageNeighborhoodRange. But I don’t know how many people actually use those ImageFunction-derived classes…

Anyway, that’s it for now :grinning: