Not an expert, but looking at the formula and your code looks good to me! After this I am OK with the design Niels!
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
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
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!
@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! 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.
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:
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::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_movecustomization 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
For your information,
ShapedImageNeighborhoodRange now also supports reverse iteration over the neighbor pixels, by offering member functions
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
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)
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
Haha, nice! Great slides
Absolutely fantastic slides @Niels_Dekker!!
What do you think about an image region range class? Is this feasible?
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?
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.
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
A plain image iterator does not need extrapolation. Only neighborhood iterators need that.
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.
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.
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:
- ZeroFluxNeumannImageNeighborhoodPixelAccessPolicy (the default policy): replicating the nearest image boundary values
- ConstantBoundaryImageNeighborhoodPixelAccessPolicy: padding with a user-specified constant value
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:
- RectangularImageNeighborhoodShape for a (hyper) rectangular neighborhood
- ConnectedImageNeighborhoodShape for an N-connected neighborhood
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
Anyway, that’s it for now