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

Engineering
1

Hereby I would like to draw your attention to a class (template) that I just submitted: ENH: Added NeighborhoodRange for efficient modern C++ style iteration. It aims to support modern C++ style iteration on a pixel neighborhood, in an efficient way.

Features:

  • Can be used in C++11 range-based for loops, for example:
    for (PixelType pixel : neighborhoodRange) …
  • Can be used to easily construct an std container of pixels, for example:
    std::vector< PixelType > pixels(neighborhoodRange.begin(), neighborhoodRange.end())
  • Can be passed directly to std algorithms (std::for_each, std::copy, etc.)
  • Can also be used with std C++ < numeric> functions, e.g., std::inner_product
  • Supports bidirectional iteration (both ++it and --it)

Performance related properties:

  • No dynamic memory allocation (not even during construction)
  • No virtual functions (even its destructors are non-virtual)
  • No ‘mutable’ data members, making it easier to write thread-safe code.

itk::NeighborhoodRange would allow a significantly performance improvement of GaussianDerivativeImageFunction, as I observed from rerunning the test from
Removal of `virtual` keywords from ConstNeighborhoodIterator
On my Windows machine, I observed a reduction of runtime duration of more than 25%.

The patch includes unit tests (Core/Common/test/itkNeighborhoodRangeTest.cxx) which show various relevant use cases.

Your review and comments are appreciated: http://review.source.kitware.com/#/c/23326

4 Likes
2

Nice work!!

I know naming is hard, but why Range in the name? I am worried it could be confusing when the ranges TS lands in c++ (c++20?). https://github.com/ericniebler/range-v3

3

Thanks for asking, Pablo! My intention is to have NeighborhoodRange match the concept of a C++ range (according to the Ranges TS). Please let me know if you think my patch would need some adaptations, in order to be a real C++ range!

4

That’s great :smiley:

5

@phcerdan Triggered by your question whether ShapedNeighborhoodIterator-like functionality could also be supported, I’m thinking of a fundamental modification of the itk::NeighborhoodRange :thinking:. Instead of allowing the user to specify the neighborhood radius (typically Size<2> or Size<3>), what about allowing to specify directly the indices, relative to the center index? For example, a user might specify a cross shaped neigborhood by relative indices {{0, 0}, {-1, -1}, {1, -1}, {1, 1}, {-1, 1}}.

The range constructor would then look like this:

NeighborhoodRange(
  ImageType& image,
  const IndexType& centerIndex,
  const std::vector<IndexType>& relativeIndices) noexcept;

A helper function could be added to create an std::vector of relative indices for a radius, for a “traditional” rectangular neighborhood:

static std::vector<IndexType> GenerateRelativeIndicesForRadius(const RadiusSizeType&);

I tested this modification already on HoughTransform2DCircles filter->Update(), it performs at least as good as patch set 7. Maybe even better! :smiley:

However, it has a single drawback. :open_mouth: The std::vector of relative indices should be owned by the user, not by the NeighborhoodRange object. It would kill performance if NeighborhoodRange would have to store a copy of the std::vector. (Seriously!) Instead, NeighborhoodRange would just have a reference (or a pointer) to the std::vector.

So for example, GaussianDerivativeImageFunction would need to have an extra data member to store these relative indices. Fortunately, these indices can be reused between EvaluateAtIndex function calls. Even with simultaneous calls from multiple treads :slight_smile:

What do you think?

My modified version is at https://gist.github.com/N-Dekker/d91246687009ac8668281a5c46990c70

1 Like
6

I think it would be great! it would make easier implementing MathematicalTopology methods in ITK. For example the neighborhood of a pixel could be 26 (fully connected), 18, or 6 in 3D. 8 for fully connected in 2D or 4. Also the approach of having indices stored in a vector could reduce confusion (don’t having to dig in the iterator what is the order at visiting pixels) at creating look-up tables for speed computation of some properties of the neighborhood IsSimple, IsConnected, IsIsthmus, IsEndPoint,… depending on the topology/neighborhood.

I like it, the drawback of maintaining the vector of indices doesn’t sounds as bad for the full control it gives.

Maybe I would recommend that the vector is ordered, for consistency, a std::sort(relativeIndices) in the GenerateXXX helpers methods will do.

+1!

1 Like
7

@phcerdan Thanks for your reply! I just pushed a new version of NeighborhoodRange onto Gerrit, with added support for arbitrarily shaped neighborhoods: Patch Set 8, http://review.source.kitware.com/#/c/23326/8

Is it possible for you to try out the patch? If so, can you please do the same performance test again, from Removal of `virtual` keywords from ConstNeighborhoodIterator ?

8

Done, using the same script with Patch Set 8, 10 runs, results in seconds:
19, 20, 27, 25, 24, 24, 25, 25, 25, 25
mean = 23.9

With ITK master (so, with virtual neighborhood iterators)
40, 46, 45, 47, 53, 51, 51, 48, 50, 49
mean = 48.0

There are proper performance tests carried in the external module ITKPerformanceBenchmarking, that will be more accurate, but looks good :stuck_out_tongue:

1 Like
9

Take into account that current ITK iterators are more generic, they use ImageAccessor, so they work with every kind of image, and they are able to use any BoundaryCondition.

We should compare as well when the removal of virtual is done for all neighborhood iterators, patch.

1 Like
10

@phcerdan Thanks a lot for the encouraging performance results!

I will think about how to include ImageAccessor and BoundaryCondition without loosing performance. (@blowekamp already mentioned before that ImageAccessor is missing.)

Can you please think about those GenerateXXX helpers methods, to generate the shapes? With Patch Set 8, I just added an itk::GenerateOffsetsForRectangularNeighborhood(radius) function at the end of itkNeighborhoodRange.h I think that’s fine, as long as it’s just one function. But you want more of those GenerateXXX functions, right? Then I think it would be preferable to have them in a separate header file. How would you name such a file? “NeigborhoodOffsets.h”? “NeigborhoodShapes.h”? “NeigborhoodShapeGenerator.h”? And more importantly, how to name the functions therein? Which shapes would you want to support?

Note that I’m already happy about the interface towards NeighborhoodRange: it accepts any contiguous container of itk::Offset objects. That allows users to pass any arbitrary shape, even if it would not be provided “out of the box” by one of our GenerateXXX functions.

Please let me know what you think!

11

NeighborhoodOffsets.h sounds good to me!.

The name GenerateOffsetsForRectangularNeighborhood, might be misleading, because the offsets could be Cubic if Size is 3D. Or lines if Size<1>.
I would recommend the more formal GenerateHyperRectangleOffsets, or GenerateNRectangleOffsets?

GenerateNCrossOffsets(radius), for one pixel wide Greek cross (+), with the lengths of the cross taken from radius?

I would put there:

Generate3D26TopologyOffsets,
 equivalent to GenerateOffsetsForRectangularNeighborhood<Size<3>>({1,1,1}) 
with the index {0} removed.

Generate3D18TopologyOffsets
removing {0} and all indices with
 std::inner_product(indice.begin(), indice.end(), indice.begin(), 0) == 3, 
or putting other way, all indices in the corners of the cube removed.

Generate3D6TopologyOffsets,
 equivalent to GenerateNCrossOffsets<Size<3>>)({1,1,1}) with index {0} removed. 

Generate2D8Topology
 GenerateOffsetsForRectangularNeighborhood<Size<2>>({1,1}) with the index {0} removed.

Generate2D4Topology
equivalent to GenerateNCrossOffsets<Size<2>>)({1,1}) with index {0} removed.

Or maybe easier to just hard-code the relative indices, but I am concerned about the order of the indices:

For example, in the case of the Gaussian that you use.

    for (const TOutput kernelValue : operatorNeighborhood.GetBufferReference())
    {
      result += kernelValue * (*neighborhoodRangeIterator);
      ++neighborhoodRangeIterator;
    }

There, the indices of neighborhoodRangeIterator has the same order than the Buffer of the operatorNeighborhood, because you coded GenerateOffsetsForRectangularNeighborhood generating the offsets in the same order than the buffer. What about when a user puts indices at random order? Shouldn’t we provide a helper function to order them in the same way than Buffer is ordered if that is the standard?

Good stuff!

1 Like
12

Thanks for your suggestions, Pablo! I did not know the term “hyperrectangle” but found it at https://en.wikipedia.org/wiki/Hyperrectangle Yes, indeed it appears to be applicable here! I now consider GenerateHyperrectangularNeighborhoodOffsets. I would like the identifier to end with NeighborhoodOffsets, to match the header filename (NeighborhoodOffsets.h). And I chose hyperrectangular to express that it’s an adjective (according to English grammar) of the neighborhood. OK?

As an alternative, I’m thinking of GenerateBoxShapedNeighborhoodOffsets. Because ‘box’ is a term I’m more familiar with. :grinning: What do you think?

I haven’t thought of how to name the other shapes yet. My idea would be to just offer the hyperrectangular (box shaped) neigborhood offsets with the first commit, and add other shapes later. Or have other people (including you) adding more shapes later. Still I very much appreciate the suggestions you gave.

Would the helper function you’re suggesting look like this?

 template < typename TOffsets> void SortNeighborhoodOffsets(TOffsets&);

I’m not yet sure if it would be helpful, because if the user puts the offsets at random order, will the user not put the kernel values in random order either?!? I think the user should just be aware that the order of the offsets should match the order of the kernel values.

Maybe a function to check if the offsets are sorted could be useful:

 template < typename TOffsets> bool AreNeighborhoodOffsetsSorted(const TOffsets&);

But of course, that still does not tell whether the kernel values are in the right order. What do you think?

13

Both sounds ok to me, naming… is hard! :sweat_smile:
I would prefer the first one, but XXXShaped is flexible, and box might work as synonym for hyperrectangle.

About the order… If we are multiplying two NeighborhoodRange, they share the same relativeIndices, and that’s it. But in the GaussianDerivative, you are mixing a itkNeighborhood, with your relativeOffsets from GenerateRectangularOffsets.
I would have to remember or dig into itkNeighborhood.h to check what is the order of the Buffer, that I guess you already knew, or already checked.

      while ( it != dogNeighborhood.End() )
        {
        pt[0] = dogNeighborhood.GetOffset(i)[direction] * directionSpacing;
        ( *it ) = m_GaussianDerivativeFunction->Evaluate(pt);
        ++i;
        ++it;
        }

      m_OperatorArray[direction] = dogNeighborhood;

    const OperatorNeighborhoodType& operatorNeighborhood = m_OperatorArray[direction];

    const NeighborhoodRange<const InputImageType> neighborhoodRange(*image, index, m_NeighborhoodOffsets[direction]);
    auto neighborhoodRangeIterator = neighborhoodRange.cbegin();

    for (const TOutput kernelValue : operatorNeighborhood.GetBufferReference())
    {
      result += kernelValue * (*neighborhoodRangeIterator);
      ++neighborhoodRangeIterator;
    }
    gradient[direction] = result;

So documenting GenerateRectangleOffsets saying that the order is the same than a regular Neighborhood iterator, or its Buffer, that they can be mixed, and explicitly write down the order in its documentation is a good idea. What do you think?

Who knows if at some point somebody want to create two Ranges with different order of the indices to do a mixing!

I guess the only mix possible between NRanges is with NIterators and NOperators which all use the same Buffer order, so forget about ordering any other generator, but document this! :stuck_out_tongue:

14

OK, let’s name the function GenerateHyperrectangularNeighborhoodOffsets :smiley: (Update, April 16: the function is now included with itkNeighborhoodOffsets.h, Patch Set 9) Thanks, Pablo!

I certainly agree it’s a good idea to double-check that such a function generates the offsets in the right order. If you want to check it, please do!

There’s a trick aspect to the ordering of the offsets. For most cases, it’s preferable to have them in an order, in line with the memory locations of the pixels they refer to. For example, the following offsets are in the preferable order:

std::vector<itk::Offset<2>> offsets = { {1,1},{2,1},{3,1}, {1,2},{2,2},{3,2} };

However, I just tried to sort them, using std::sort:

std::sort(offsets.begin(), offsets.end());

Result: { {1,1},{1,2}, {2,1},{2,2}, {3,1},{3,2} }

So instead of preserving the order (which was the right order), std::sort sorts it in a way that is in general not good for our use cases! This has to do with the way itk::Offset has overloaded operator<.

For the time being, I would like to postpone something like a SortNeighborhoodOffsets function to “version 2”, and just be careful to make sure that GenerateHyperrectangularNeighborhoodOffsets does the right thing, in the right order. (Especially because the ImageAccessor and BoundaryCondition issues have a higher priority!)

15

Update:I just removed the “work-in-progress tag” (WIP) from the patch, as I think it’s ready for review now:

ENH: Added NeighborhoodRange for efficient modern C++ style iteration - Patch Set 10

It has no code changes compared to the previous patch set. Only the commit message was updated.

I really hope this patch can be included with the first release of ITK 5, as I think it would be a great new feature.

1 Like
16

Thanks for all your feedback so far! Patch Set 11 is there! I hope I have addressed all your comments properly. However, there’s one issue, raised by Bradley, that it does not solve: So far, the proposed itk::NeighborhoodRange< Image> class does not use the NeighborhoodAccessorFunctor of the image. Is that really necessary? Could we still postpone this issue to the next commit, once the initial version is accepted? Especially because the patch set is quite large already…

I did have a first glance at those accessors, and it appeared to me that they would make the implementation of NeighborhoodRange significantly more complicated. :hushed: Especially iterator::operator*(), which now returns a (possibly non-const) reference to the current pixel:

reference operator*() const noexcept
{
  const OffsetType currentOffset = (*m_CurrentOffset);
  auto* result = m_ImageBufferPointer + GetClampedIndex(currentOffset[0] + m_Location[0], m_ImageSize[0]);

  for (ImageDimensionType i = 1; i < ImageDimension; ++i)
  {
    result += GetClampedIndex(currentOffset[i] + m_Location[i], m_ImageSize[i]) * m_OffsetTable[i];
  }
  return *result;
}

I got the impression that especially non-const access to the pixel would be hard to implement, using those accessors, because they don’t offer a non-const reference to the pixel.

If the main reason for using accessors is just to support itk::VectorImage, I think VectorImage support can be achieved in a different way as well, just by specializing the relevant parts of NeighborhoodRange for VectorImage. What do you think?

17

I do like the idea of making new ITK iterators that are more compatible with STL, and modern C++11 iteration syntax. However, I think this abandons too many established ITK design requirements and conventions.

The design in ITK also support the Image Adaptors: https://itk.org/Doxygen/html/group__ImageAdaptors.html
As well as other data buffer layouts. This is part the the generic design of ITK and the image classes to make them as flexible as possible.

Because of this requirement for our image data structure the standard interface for ITK iterators with Set and Get was necessary.

The term “Neighborhood” in ITK is already heavily used with itk::Neighborhood, itk:NeighborhoodOperator and itk::NeighborhoodIterator. This new class is not compatible with these classes and methods. I don’t think it should use the term Neighborhood.

-grumpy old ITK developer

18

So what to do with it?

I’ll suggest that this really should be an Insight Journal Article and in a remote module (called ModernIterators?) and not in ITK proper yet. This is a new, and good idea here but it will take some time for people to use it and determine if it is generally suitable for ITK proper.

19

In practice, since a remote module cannot depend on another remote module, it is going to be difficult to test this iterator out with a remote module.

Good point about the term Neighborhood – it should not be in the name if does not have the meaning of the existing Neighboorhood classes. The name should also follow the ITK class naming convention so it can easily be identified and understood, i.e. it needs to end with Iterator. If it only supports itk::Image, then it should have Image in the name. Perhaps ImageRangeIterator?

The limitation of not supporting Image Adaptors could be noted in the class. In practice, ITK does not work with image data buffer layouts other than the block layout. While it may be sufficient to not support Image Adaptors, itk::VectorImage is commonly used, and there should be tests to confirm that the iterator works with this class.

And, documentation is critical, in the form of an Insight Journal article and documentation in the ITK Software Guide.

20

This new style seems a bit experimental to me. Perhaps it should be in a sub-namespace of itk to help indicate that then. As being experimental, I am not sure the ITK Software Guide is the place for things that have not stabilized and may need to be changed in the future.