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

(Matt McCormick) #41

@Niels_Dekker thanks again for your contribution and your persistence! :sweat_smile:

Yes, we can do better in terms of making incremental improvements. If additional features could optionally be added in future commits, we should create an issue in the issue tracker and move forward.

At the same time, changes to improve regression test coverage, style changes, API changes, and similar should be made before the patch is merged, otherwise they may (likely) never happen.

Merging faster is also possible when follow-up’s come quickly for issues identified on the Nightly Dashboard builds. If these are not fixed quickly, they mask issues from other patches that other community members are trying to develop.

Currently, there are four issues causing warnings from the ShapedImageNeighborhoodRange patch – I will take care of these, unless you get to them before me :wink: .

Boundary condition improvements sound excellent! We are targeting 5.0 Alpha 2 sometime over this week…

(Niels Dekker) #42

Hmmm… those warnings aren’t very helpful :disappointed: says:

[CTest: warning matched] /scratch/dashboards/Linux-x86_64-gcc4.8-m32/ITK/Modules/Core/Common/test/itkImageNeighborhoodOffsetsGTest.cxx:33:31:
warning: missing initializer for member 'itk::Size<>::m_InternalArray' [-Wmissing-field-initializers]
   const itk::Size<> radius = {};

That’s a perfect line of C++! Would GCC4.8 like one the following instead?

   const itk::Size<> radius = { {} };

Or the following:

   const itk::Size<> radius = { {0} };

Or the following modern C++11:

   const auto radius = itk::Size<>{};

Or this one:

   const auto radius = itk::Size<>();

At least, I wouldn’t like to see all elements spelled out:

   const itk::Size<> radius = { {0, 0} };  // Please no!  :-(

PS Another option would be to just remove the ‘=’ sign from the original line of code:

const itk::Size<> radius{};

Would that compile without a warning?

(Matt McCormick) #43

I assume we need to tell the compiler how long we want the itk::Size type to be, either in the template arguments, < Dimension > or in the initializer, { {0, 0} }

(Bradley Lowekamp) #44

Why not just do:

const itk::Size<> radius;

I’m confused about when the aggregate initializer is called vs the wonderful “modern” style initializers with you suggestions. And with the default initialization I’d have to lookup if the member C-array actually gets initialized or not when it’s a member variable.

I agree with Matt that being explicit is good too!

(Niels Dekker) #45
const itk::Size<> radius = {};

Thanks, Matt, but that’s not the issue here. itk::Size has a default, VDimension = 2:

For testing purposes (the issue is specifically about the tests, itkImageNeighborhoodOffsetsGTest and itkShapedImageNeighborhoodRangeGTest), I think it’s just fine to use the default (VDimension = 2).

However, even if the dimension of the itk::Size<VDimension> would be specified explicitly, I guess that wouldn’t stop the old GCC4.8 from producing the warning.

Looks fine, but I’m afraid that won’t compile. A “const” object needs to be initialized during its declaration. Again, the following looks like a proper C++11 way to me:

const itk::Size<> radius{};

I think that’s what Bjarne Stroustrup would prefer :slightly_smiling_face: (Bjarne designed the new C++11 {} uniform initialization syntax.)

As you already offered taking care of it, @matt.mccormick, please go ahead, if you like…

But please please have a look at my proposal to allow ITK developers and users to customize the border extrapolation (specifying the boundary condition), allowing to specify a custom policy class: ENH: Added ImageNeighborhoodPixelAccessPolicy for custom border extrapolation

(Bradley Lowekamp) #46

What do you expect the following code to initialize m_InternalArray to?

const itk::Size<> radius{};

I would expect this to call the default constructor which is declared as =default. I don’t believe for aggregate types the C-array gets initialized by default. So your proposed code suffers from a similar problem as my suggestion, it does not initialize m_InternalArray. Is my understanding not correct here?

(Niels Dekker) #47

I would expect it to initialize radius.m_InternalArray by all zero’s :slightly_smiling_face: And I think it does!

Fortunately the C++11 initialization semantics are slightly different. Because itk::Size<> is an aggregate (all its data members are public, and it has no virtual functions), the code itk::Size<> radius{} does aggregate-initialization, which properly initializes all data: So it does not call the default-constructor of itk::Size<>!!!

Hope that helps, Niels

(Niels Dekker) #48

@phcerdan As you specifically asked for shape support in the first place, can you maybe have a look at my neighborhood shape class? WIP: Added class to create offsets for a hyperrectangular neighborhood shape

It aims to replace the function GenerateHyperrectangularImageNeighborhoodOffsets(radius), that we discussed here before, by a class, HyperrectangularImageNeighborhoodShape<ImageDimension>, which has a GenerateOffsets() member function.

I found a dedicated shape class preferable to a bunch of functions, because actually there are three of them, not just one:

size_t CalculateNumberOfOffsets();
void FillOffsets(OffsetType* offsets, size_t numberOfOffsets);
std::vector<OffsetType> GenerateOffsets();

GenerateOffsets() is fine if you just want to have the offsets stored in an std::vector, but in some cases, you might prefer to store the offsets elsewhere (typically for performance reasons), for example in an std::array<OffsetType, N>, or an itk::NeighborhoodAllocator<OffsetType>. For such use cases, you could use CalculateNumberOfOffsets() and FillOffsets(offsets, numberOfOffsets).

GenerateOffsets() is implemented, simply by calling CalculateNumberOfOffsets() and FillOffsets(offsets, numberOfOffsets), to fill an std::vector of offsets.

Do you think this would be a proper design for other shape classes as well?

(Pablo Hernandez-Cerdan) #49

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?

(Niels Dekker) #50

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

(Pablo Hernandez-Cerdan) #51

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

(Pablo Hernandez-Cerdan) #52

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)

(Pablo Hernandez-Cerdan) #53

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:

(Niels Dekker) #54

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:

A shape class for ShapedImageNeighborhoodRange and ShapedNeighborhoodIterator
(Niels Dekker) #55

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:

A shape class for ShapedImageNeighborhoodRange and ShapedNeighborhoodIterator
(Pablo Hernandez-Cerdan) #56

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.

(Niels Dekker) #57

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:

(Pablo Hernandez-Cerdan) #58

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:

(Pablo Hernandez-Cerdan) #59

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.

(Niels Dekker) #60

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.