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

matt.mccormick · April 29, 2018, 8:15pm

@Niels_Dekker thanks again for your contribution and your persistence!

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 .

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

Niels_Dekker · April 30, 2018, 11:53am

Hmmm… those warnings aren’t very helpful

https://open.cdash.org/viewBuildError.php?type=1&buildid=5357834 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 · April 30, 2018, 1:16pm

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

blowekamp · April 30, 2018, 1:55pm

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 · April 30, 2018, 2:37pm

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

Thanks, Matt, but that’s not the issue here. itk::Size has a default, VDimension = 2: ITK: itk::Size< VDimension > Struct Template Reference

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

blowekamp · April 30, 2018, 2:51pm

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 · April 30, 2018, 3:19pm

I would expect it to initialize radius.m_InternalArray by all zero’s 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: Aggregate initialization - cppreference.com So it does not call the default-constructor of itk::Size<>!!!

Hope that helps, Niels

Niels_Dekker · May 2, 2018, 6:21pm

@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 http://review.source.kitware.com/#/c/23389/

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?

phcerdan · May 3, 2018, 12:37am

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 · May 3, 2018, 10:47am

@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

phcerdan · May 4, 2018, 10:17am

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(
    const Size<VImageDimension>& radius) ITK_NOEXCEPT
  :
  m_Radius(radius)
  {
  }

  /** 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;
    });
  }
private:
  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;
}

CMakeLists.txt

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

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

phcerdan · May 4, 2018, 10:35am

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)

phcerdan · May 4, 2018, 12:13pm

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

Niels_Dekker · May 4, 2018, 2:36pm

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
{
  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: constexpr CalculateNumberOfOffsets() for ITK, C++ (clang) - rextester

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

Niels_Dekker · May 7, 2018, 2:17pm

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:

gist.github.com

https://gist.github.com/N-Dekker/0f17c0473b29494faea93244bbb43393

itkTopologicalConnectivityImageNeighborhoodShape.h

/*=========================================================================
*
*  Copyright Insight Software Consortium
*
*  Licensed under the Apache License, Version 2.0 (the "License");
*  you may not use this file except in compliance with the License.
*  You may obtain a copy of the License at
*
*         http://www.apache.org/licenses/LICENSE-2.0.txt
*

This file has been truncated. show original

phcerdan · May 7, 2018, 2:26pm

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 · May 7, 2018, 3:22pm

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!

phcerdan · May 7, 2018, 3:27pm

Not an expert, but looking at the formula and your code looks good to me! After this I am OK with the design Niels!

phcerdan · May 7, 2018, 3:35pm

By the way, there is a Factorial in ITK: itk::GeometryUtilities Class Reference that might be a good candidate to make it constrexpr and re-use it in your code.

Niels_Dekker · May 7, 2018, 4:16pm

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.