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

Niels_Dekker · April 23, 2018, 3:14pm

Update: This morning I submitted Patch Set 19, which can be reviewed, and possibly merged…

itkNeighborhoodRange.h has significant improvements to the reference types, NeighborhoodRange::iterator::reference (the return type of iterator::operator*()), and NeighborhoodRange::const_iterator::reference. Under the hood, these types are proxy types: Patch Set 19 implements them by an internal (private) template class, PixelProxy<VIsConst>. (VIsConst is true for const access to the pixel, and false for non-const access). Both const and non-const PixelProxy now have a private helper object of type PixelProxyHelper<VIsConst>, which communicates with the NeighborhoodAccessor of the image (calling NeighborhoodAccessor.Get(internalPixel) and NeighborhoodAccessor.Set(internalPixel, pixelValue)).

Tests were added to itkNeighborhoodRangeGTest.cxx to check that the reference types behave like real (built-in) C++ reference types (PixelType & and const PixelType &), and to check that a neighborhood iteration on an itk::VectorImage is also supported well.

Please have a look!

Niels_Dekker · April 24, 2018, 9:53am

@matt.mccormick I do appreciate that you care about proper class names. Now the following existing neighborhood related ITK classes deal specifically with images, and still don’t have Image in their identifier. Do you think they should?

NeighborhoodIterator
ConstNeighborhoodIterator
ConstNeighborhoodIteratorWithOnlyIndex

ShapedNeighborhoodIterator
ConstShapedNeighborhoodIterator

NeighborhoodAccessorFunctor
NeighborhoodInnerProduct

LevelSetNeighborhoodExtractor
LevelSetVelocityNeighborhoodExtractor

matt.mccormick · April 24, 2018, 1:45pm

Yes, I think they should have Image in their name. Since there is abundant client code that use these classes, it is difficult and disruptive to change the name now. Hence, all the fuss about the name

Niels_Dekker · April 24, 2018, 5:39pm

[Regarding NeighborhoodIterator, ConstNeighborhoodIterator, ConstNeighborhoodIteratorWithOnlyIndex, ShapedNeighborhoodIterator, ConstShapedNeighborhoodIterator, etc.]

OK, but still I don’t like the name ImageShapedNeighborhoodRange very much I would read it as “image shaped neighborhood range”. Which sounds to me like “heart shaped box”. The Nirvana song, you know? https://www.youtube.com/watch?v=n6P0SitRwy8 But while the box of Nirvana may be heart shaped, the neighborhood certainly isn’t image shaped. You get my point?

What do you think of ShapedImageNeighborhoodRange? It should be read as “shaped image-neighborhood range”: the range class for a shaped image-neighborhood.

matt.mccormick · April 24, 2018, 7:01pm

Good point – it could be perceived as image-shaped.

Agreed, that is even better!

Thanks for making the adjustments. It Smells Like Team Spirit.

Niels_Dekker · April 25, 2018, 2:18pm

I’m preparing another patch set which I’ll try to submit later today. It should include renaming the range class, as we discussed here. Hopefully this will be the final patch set before the merge to the master branch of ITK!

But there’s something I need to tell you… While developing the range class, I was informed by two C++ experts, Anthony Williams and Jonathan Wakely, that the iterators I’m proposing here do not entirely satisfy the iterator requirements as specified by the current C++ Standard! So what does that mean…?

itk::NeighborhoodRange::iterator supports everything you would expect from a bidirectonal iterator: ++it, --it, it++, it--, *it, auto it2{it1}, it2 = it1, it1 == it2, it1 != it2. However, it does not satisfy the following requirements (quoted from the current Working Draft of the C++ Standard):

From section [iterator.requirements.general]:

An iterator i for which the expression (*i).m is well-defined supports the expression i->m with the same semantics as (*i).m.

From section [forward.iterators]:

if X is a mutable iterator, reference is a reference to T; if X is a constant iterator, reference is a reference to const T

itk::NeighborhoodRange::iterator does not have an operator->(), so it does not support it->m. Moreover, itk::NeighborhoodRange::iterator::reference (the return type of iterator::operator*()) is not a real C++ reference type (PixelType &), instead, it’s a proxy class (PixelProxy).

Strictly speaking, that implies that itk::NeighborhoodRange::iterator is not guaranteed to work well as argument to a C++ Standard Library function that requires a standard compliant iterator.

Fortunately itkNeighborhoodRangeGTest.cxx tests such use cases, and they all pass successfully:

std::vector<PixelType>(range.cbegin(), range.cend());
std::reverse_copy(range.begin(), range.end(), stdVector.begin());
std::inner_product(range.begin(), range.end(), range.begin(), 0.0);
std::for_each(range.begin(), range.end(), [](const PixelType&){});

Moreover: some iterators in the C++ Standard Library don’t satisfy these particular requirement either! std::istreambuf_iterator does not have an operator->() (as reported by Jonathan Wakely) And std::vector<bool>::iterator also has a proxy as iterator::reference type.

Even more hopeful: There is a very detailed proposal to get full support for proxy iterators into the C++ Standard Library: Proxy Iterators for the Ranges Extensions by Eric Niebler.

Summary: itk::NeighborhoodRange::iterator (patch set 19) may not yet officially satisfy all iterator requirements from the current C++ Standard, but in practice, it just appears to work well, as argument to an std function. And of course, the other advantages are still there as well: a significant performance improvement, support for multi-threading and support for range-based for loops.

Niels_Dekker · April 26, 2018, 12:03pm

Very happy to see that itk::ShapedImageNeighborhoodRange - the class formerly known as “NeighborhoodRange” - has been merged onto the master! Thanks, Matt!!!

Thanks to you all: @matt.mccormick, @blowekamp, @dzenanz, @phcerdan, @jhlegarreta, @hjmjohnson, @LucH, and Kitware Robot, of course!

matt.mccormick · April 26, 2018, 2:05pm

Thank you, @Niels_Dekker and everyone who helped with the reviews! This is a huge contribution to ITK.

Niels_Dekker · April 28, 2018, 5:34pm

By the way, I found the review process for this patch quite exhausting. You see, it was only accepted after submitting the 21st patch set. While I felt that patch set 7 was already worth a commit. Patch set 7 already worked well enough for its initial use case, in my opinion: improving the performance of GaussianDerivativeImageFunction and HoughTransform2DCirclesImageFilter. OK, it only supported rectangular neighborhoods, but as a first commit, that seemed reasonable to me. Of course, adding custom shape support was a great enhancement, but it could have been added with a separate commit.

With patch set 12, custom shape support was fully implemented, and again, I was hoping that that would justify a commit. However, it was only after another major enhancement, using the NeighborhoodAccessor of the image, that the patch was accepted.

So basically I felt some pressure to keep adding more features to a single commit, which were holding back the merge. While I was afraid that these extra features would only make it harder to get everything reviewed carefully, as the implementation got more complicated, and that it would become too big to get it into ITK 5.0 in time.

I still feel it would have been better to have had three commits, instead of one, for this patch: commit 1 supporting rectangular neighborhoods, commit 2 adding custom shape support, and commit 3 using NeighborhoodAccessor. The intermediate commits might still have been useful to analyze the design decisions afterwards, from the git log. Hope you understand my point!

Anyway, thanks again for having ShapedImageNeighborhoodRange at the master branch! I’m proud that I was able to make this contribution

Kind regards, Niels

P.S. I’m now preparing to add a way to allow a custom border extrapolation (boundary conditions). How much time is still left before the next alpha?

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