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

Not an expert, but looking at the formula and your code looks good to me! After this I am OK with the design Niels! :smiley:

By the way, there is a Factorial in https://itk.org/Doxygen/html/classitk_1_1GeometryUtilities.html that might be a good candidate to make it constrexpr and re-use it in your code.

Nice idea! I did not even know itk::GeometryUtilities::Factorial(long) yet!

I was about to say that a C++11 constexpr factorial implementation may not be ideal for calculations at run-time, because the C++11 constexpr can only be implemented by a recursive function call. For large numbers, that could potentially be slow at run-time, and use lots of stack space. So I expected the factorial function from GeometryUtilities to be implemented by a for-loop, instead of using recursion. But now I see GeometryUtilities::Factorial(long) is actually already implemented by a recursive function call!!!

Would you like to send send a patch, to make GeometryUtilities::Factorial(long) constexpr and move its implementation from the CXX file to the H file? (Or possibly add an overload for std::size_t?) Then it could certainly be used for TopologicalConnectivityImageNeighborhoodShape.

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

1 Like

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:

@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: http://review.source.kitware.com/#/c/23389/ 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.

1 Like

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 http://ericniebler.com/2015/03/03/iterators-plus-plus-part-3/#comment-175569 For more details, please check his blog!

Kind regards, Niels

2 Likes

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:

1 Like

The patch I just submitted should finally make ShapedImageNeighborhoodRange compatible with the Range Library from https://github.com/ericniebler/range-v3 (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 (http://review.source.kitware.com/#/c/23581)

3 Likes

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: https://nd.home.xs4all.nl/cpp/talks/NielsDekker-MeetUp-July2018-ITK-NeighborhoodRange.pptx

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

3 Likes

Haha, nice! Great slides :smiley:

1 Like

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

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

2 Likes

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?

1 Like

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.

1 Like

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: http://www.open-std.org/jtc1/sc22/wg21/docs/papers/2018/p0709r0.pdf

1 Like

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

1 Like

Nice…

2 Likes

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.

1 Like

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.

1 Like

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:

4 Likes