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

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.

phcerdan · May 7, 2018, 7:25pm

I would like it indeed, but not sure if I will find time until next week

Niels_Dekker · May 8, 2018, 10:11am

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!

Niels_Dekker · May 15, 2018, 8:24pm

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

Niels_Dekker · May 17, 2018, 8:52am

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 Iterators++, Part 3 – Eric Niebler For more details, please check his blog!

Kind regards, Niels

Niels_Dekker · June 6, 2018, 8:59pm

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

Niels_Dekker · July 24, 2018, 1:25pm

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)

Niels_Dekker · July 27, 2018, 11:19am

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.

phcerdan · July 27, 2018, 2:43pm

Haha, nice! Great slides

matt.mccormick · July 27, 2018, 4:57pm

Absolutely fantastic slides @Niels_Dekker!!

What do you think about an image region range class? Is this feasible?