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

Nice Niels, I am not sure about imposing those three functions for every shape. They might be good for the Hyperrectangular case with a variable radius where some calculations are needed, but for simpler cases… do we need a whole class?

If you want to have the hyperrectangular class in its own header we can use the existing NeighborhoodOffsets to just hard-code these offsets (I am not sure if this order is optimal for performance):

``````std::array<itk::Offset<2>, 4> Topology4NeighborhoodOffsets{{1,0},{-1,0},{0,1}{0,-1}};
std::array<itk::Offset<2>, 8> Topology8NeighborhoodOffsets{{1,0},{-1,0},{0,1}{0,-1}, {1,1},{1,-1},{-1,1},{-1,-1}};
std::array<itk::Offset<3>, 6> Topology6NeighborhoodOffsets{{1,0,0},{-1,0,0},{0,1,0}{0,-1,0}, {0,0,1},{0,0,-1}};
std::array<itk::Offset<3>, 18> Topology18NeighborhoodOffsets{{1,0,0},{-1,0,0},{0,1,0}{0,-1,0}, {0,0,1},{0,0,-1},
{1,1,0},{1,-1,0},{-1,1,0},{-1,-1,0}, {1,0,1},{1,0,-1},{-1,0,1},{-1,0,-1}, {0,1,1},{0,1,-1},{0,-1,1},{0,-1,-1}};
std::array<itk::Offset<3>, 26> Topology26NeighborhoodOffsets{{1,0,0},{-1,0,0},{0,1,0}{0,-1,0}, {0,0,1},{0,0,-1},
{1,1,0},{1,-1,0},{-1,1,0},{-1,-1,0}, {1,0,1},{1,0,-1},{-1,0,1},{-1,0,-1}, {0,1,1},{0,1,-1},{0,-1,1},{0,-1,-1},
{1,1,1},{1,1,-1},{1,-1,1},{-1,1,1}, {1,-1,-1},{-1,1,-1},{-1,-1,1},{-1,-1,-1}};
``````

These offsets are not designed for kernel multiplications as you do in your Gaussian application, but for setting properties of the pixel/voxel based on its neighborhood.

Also, if you are worried about performance, why not use `std::array` instead of `std::vector`? Is it possible in c++11 to constexpr `CalculateNumberOfOffsets` to be used as the dimension parameter of `std::array`?

1 Like

@phcerdan Thanks for your feedback. First of all, the change that I have in mind with WIP: Added class to create offsets for a hyperrectangular neighborhood shape keeps the interface of `ShapedImageNeighborhoodRange` as it is. In my opinion, the neighborhood range class should definitely keep allowing users to specify the shape by passing a pointer to a contiguous sequence of `itk::Offset` objects.

But I do think it could be useful for various shapes to offer all these three function: `CalculateNumberOfOffsets()`, `FillOffsets(offsets, numberOfOffsets)`, and `GenerateOffsets()`. In that case, I think it would be nice to have these three functions grouped together in a class, specific for that shape.

The size of the neighborhood is not always not always known at compile time. Certainly not for `GaussianDerivativeImageFunction`! (If you want to know the details: For this class, the neighborhood size depends on member variables m_Sigma and m_Extent.) As a rule of thumb, you might consider: use `std::array<OffsetType, N>` when the size is known at compile time, otherwise use `std::vector<OffsetType>`.

Good question! Yes, I would like `CalculateNumberOfOffsets()` to be constexpr. I’m pretty sure it’s possible in C++14. For C++11, I’d have to try… Please let me know if you already have it working, for the hyperrectangular case!

Could these also be implemented in a less “hard-coded” way? Possibly parameterized by their size?

Anyway, if there are no other shapes of interest to ITK users that could be implemented by the class design that I proposed at http://review.source.kitware.com/#/c/23389/ it’s probably less useful than I thought

I have tested if it would be possible to make `CalculateNumberOfOffsets` constexpr, but with no luck. It seems `std::accumulate` is not `constexpr` (neither in c++11, 14 or 17)

Here are my tests, just in case:

`ComputeOffsets.cpp`

``````#include <algorithm> // For transform.
#include <cassert>
#include <numeric>  // For accumulate.
#include <vector>

#include "itkOffset.h"
#include "itkSize.h"

namespace itk {
template <unsigned VImageDimension>
class HyperrectangularImageNeighborhoodShape
{
public:
/** Constructs a hyperrectangular shape whose size is specified by the radius */
explicit HyperrectangularImageNeighborhoodShape(
:
{
}

/** Calculates the number of offsets needed for this shape. */
constexpr std::size_t CalculateNumberOfOffsets() const noexcept
{
[](const std::size_t accumulatedProduct, const std::size_t radiusValue)
{
return ((2 * radiusValue) + 1) * accumulatedProduct;
});
}
private:
};
} //end ns

int main()
{
itk::HyperrectangularImageNeighborhoodShape<2> shape({1,1});
auto noffsets = shape.CalculateNumberOfOffsets();
std::array<itk::Offset<2>,shape.CalculateNumberOfOffsets()> hola;

std::cout << "Number of Offsets: " << noffsets << std::endl;
}
``````

`CMakeLists.txt`

``````cmake_minimum_required(VERSION 3.1)
project(ComputeOffsets)
find_package(ITK COMPONENTS
ITKImageFeature
REQUIRED
)
include(\${ITK_USE_FILE})

``````

configure with:
`cmake -DCMAKE_CXX_STANDARD=11 -DITK_DIR:PATH=//ITK/build ../src`

1 Like

Maybe with?

``````template<VDimension, TNumberOfOffsets>
struct NearestNeighborhoodOffsets{
std::array<itk::Offsets<VDimension>, TNumberOfOffsets> m_Offsets;
};

template<>
std::array<itk::Offsets<2, 8> NearestNeighborhoodOffsets<2,8>::m_Offsets {{}};
...
``````

(Not tested)

Don’t get me wrong I think it is still useful for complicated shapes! My only concern was imposing it to simpler cases as NearestNeighborOffsets for example

1 Like

Thanks very much for giving it a try, @phcerdan I hadn’t realized that `std::accumulate` is not constexpr. But even in the latest draft of the C++ Standard, it is not.

However, in the meantime I did manage to write a constexpr `CalculateNumberOfOffsets()` that does compile… at least on g++ 5.4.0, clang 3.8.0 and MSVC 2017:

``````#include <array>

// Minimal version of itk::Size<VDimension>
template <unsigned int VDimension = 2>
struct Size
{
std::size_t m_InternalArray[VDimension];
};

// Minimal version of proposed shape class from http://review.source.kitware.com/#/c/23389/
template <unsigned VImageDimension>
struct HyperrectangularImageNeighborhoodShape
{

constexpr std::size_t CalculateNumberOfOffsets() const noexcept
{
std::size_t result = 1;

{
result *= 2 * radiusValue + 1;
}
return result;
}
};

int main()
{
constexpr Size<> radius{ { 3, 4 } };
constexpr std::size_t numberOfOffsets = shape.CalculateNumberOfOffsets();
static_assert(numberOfOffsets == (2 * 3 + 1) * (2 * 4 + 1), "Yes!");

// Yes: Can be used as the dimension parameter of std::array!  :-)
struct OffsetType {};
std::array<OffsetType, numberOfOffsets> offsets{ {} };
}
``````

Of course, this cannot yet be used to implement CalculateNumberOfOffsets() for ITK 5.0, because C++11 does not support a for loop or a variable declaration inside a constexpr function. But any ITK user can actually write such code, when using a C++14 compiler. Isn’t that cool?

I’ll still look at your other posts here… hope you’re not impatient

In the meantime I did write a class, `TopologicalConnectivityImageNeighborhoodShape`, to produce those offsets. It ensures that the order is optimal for performance, and it has a constexpr `CalculateNumberOfOffsets()`. Moreover, it can in principle do any connectivity, and any number of dimensions

Here are the offsets that the shape class produces for 1-D, 2-D, and 3-D:

1-D 2-connected: {[-1],[1]}
2-D 4-connected: {[0, -1],[-1, 0],[1, 0],[0, 1]}
2-D 8-connected: {[-1, -1],[0, -1],[1, -1],[-1, 0],[1, 0],[-1, 1],[0, 1],[1, 1]}
3-D 6-connected: {[0, 0, -1],[0, -1, 0],[-1, 0, 0],[1, 0, 0],[0, 1, 0],[0, 0, 1]}
3-D 18-connected: {[0, -1, -1],[-1, 0, -1],[0, 0, -1],[1, 0, -1],[0, 1, -1],[-1, -1, 0],[0, -1, 0],[1, -1, 0],[-1, 0, 0],[1, 0, 0],[-1, 1, 0],[0, 1, 0],[1, 1, 0],[0, -1, 1],[-1, 0, 1],[0, 0, 1],[1, 0, 1],[0, 1, 1]}
3-D 26-connected: {[-1, -1, -1],[0, -1, -1],[1, -1, -1],[-1, 0, -1],[0, 0, -1],[1, 0, -1],[-1, 1, -1],[0, 1, -1],[1, 1, -1],[-1, -1, 0],[0, -1, 0],[1, -1, 0],[-1, 0, 0],[1, 0, 0],[-1, 1, 0],[0, 1, 0],[1, 1, 0],[-1, -1, 1],[0, -1, 1],[1, -1, 1],[-1, 0, 1],[0, 0, 1],[1, 0, 1],[-1, 1, 1],[0, 1, 1],[1, 1, 1]}

Would such a class be interesting to you?

Note: It’s still an early version, and I would only submit it to ITK if anyone would actually like to use it:

1 Like

Nice!! I am super interested! Could you copy those results in the documentation of the class? Most of the time the usage would be with those values, but the generecity of it is great indeed!. Good stuff, thanks for this.

1 Like

Thanks for your enthusiasm, Pablo Here is the code that I wrote to print the offsets:

``````Print(TopologicalConnectivityImageNeighborhoodShape<1>{ 1 });
Print(TopologicalConnectivityImageNeighborhoodShape<2>{ 1 });
Print(TopologicalConnectivityImageNeighborhoodShape<2>{ 2 });
Print(TopologicalConnectivityImageNeighborhoodShape<3>{ 1 });
Print(TopologicalConnectivityImageNeighborhoodShape<3>{ 2 });
Print(TopologicalConnectivityImageNeighborhoodShape<3>{ 3 });
``````

With the following `Print` function:

``````template <unsigned VDimension>
void Print(const TopologicalConnectivityImageNeighborhoodShape<VDimension>& shape)
{
std::cout << VDimension << "-D " << shape.CalculateNumberOfOffsets() << "-connected: ";
const std::vector<Offset<VDimension>> offsets = shape.GenerateOffsets();

std::cout << "{";
for (const Offset<VDimension>& offset : offsets)
{
std::cout << offset << ((offset == offsets.back()) ? "" : ",");
}
std::cout << "}" << std::endl;
}
``````

The hardest part of the implementation of TopologicalConnectivityImageNeighborhoodShape was how to write its constexpr `CalculateNumberOfOffsets()` member function. I’m not a mathematician! Fortunately I got a lot of help from my LKEB colleague Baldur van Lew, who suggested me to use the formula from H.S.M. Coxeter, at https://en.wikipedia.org/wiki/Hypercube#Elements Is that formula familiar to you? If so, maybe you can check the source file to see if I applied the formula correctly!

But before I would submit a topological/connectivity shape class to ITK, I would like to be sure that the general class design for shape classes is OK. As we discussed before: having a shape class, offering those three member functions, `CalculateNumberOfOffsets()`, `FillOffsets(offsets, numberOfOffsets)`, and `GenerateOffsets()`. It’s still open for discussion!

1 Like

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

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!

@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.

1 Like

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

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.

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.

3 Likes

Haha, nice! Great slides

1 Like