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

Engineering
61

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

1 Like
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:

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: 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
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 Iterators++, Part 3 ā€“ Eric Niebler For more details, please check his blog!

Kind regards, Niels

2 Likes
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:

1 Like
66

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
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: 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:

highres_473244168.jpeg3968Ɨ2976 622 KB

3 Likes
68

Haha, nice! Great slides :smiley:

1 Like
69

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

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

2 Likes
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?

1 Like
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.

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

1 Like
73

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

1 Like
74

Niceā€¦

2 Likes
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.

1 Like
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.

1 Like
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:

4 Likes