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! :thinking: 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 it’s probably less useful than I thought :disappointed_relieved:

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:


#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
  /** Constructs a hyperrectangular shape whose size is specified by the radius */
  explicit HyperrectangularImageNeighborhoodShape(
    const Size<VImageDimension>& radius) ITK_NOEXCEPT

  /** Calculates the number of offsets needed for this shape. */
  constexpr std::size_t CalculateNumberOfOffsets() const noexcept
    return std::accumulate(m_Radius.begin(), m_Radius.end(), std::size_t{ 1 },
      [](const std::size_t accumulatedProduct, const std::size_t radiusValue)
      return ((2 * radiusValue) + 1) * accumulatedProduct;
  Size<VImageDimension> m_Radius;
} //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;


cmake_minimum_required(VERSION 3.1)
find_package(ITK COMPONENTS

add_executable(ComputeOffsets ComputeOffsets.cpp)
target_link_libraries(ComputeOffsets PUBLIC ${ITK_LIBRARIES})

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;

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

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

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
template <unsigned VImageDimension>
struct HyperrectangularImageNeighborhoodShape
  const Size<VImageDimension> m_Radius;

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

    for (const auto radiusValue : m_Radius.m_InternalArray)
      result *= 2 * radiusValue + 1;
    return result;

int main()
  constexpr Size<> radius{ { 3, 4 } };
  constexpr HyperrectangularImageNeighborhoodShape<2> shape{ radius };
  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{ {} };

Or click here:

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

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

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 :smiley: 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 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! :grinning:

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

By the way, there is a Factorial in 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: 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 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 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 (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 (


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:


Haha, nice! Great slides :smiley:

1 Like