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

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 Like

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

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 ?

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

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

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

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:

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

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.

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

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

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);

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

Thanks for your suggestions, Pablo! I did not know the term “hyperrectangle” but found it at 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?

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);

      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);
    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:

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!)

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

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?

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

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.

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.

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.

Yes, maybe itk::Experimental.

Yes, an Insight Journal article could be first, then some of its contents could be applied in the Software Guide.

@blowekamp, @matt.mccormick Thank you both for this fundamental discussion!

For now I would just like to add some remarks:

My current project at LKEB (Leiden University Medical Center) uses ITK’s implementation of Hough Transform for circle detection. If the range class I’m proposing here cannot be used to improve the performance of itk::HoughTransform2DCirclesImageFilter, I find it harder to justify the hours I’m spending on the range class, here at my work. I hope you understand! So this raises the question: if the range class is placed in an itk::Experimental namespace, can it still be used by other non-Experimental ITK classes?

If it’s really a show-stopper that the current NeighborhoodRange (Patch Set 12) does not use the NeighborhoodAccessor, I think it can be fixed, actually. :smiley: I do have it working right now, locally at my computer.

Regarding the name of the class, I did carefully choose “NeighborhoodRange”. The term “range” matches the concept of a modern C++ range (Ranges TS) The term “Neighborhood” is essential here as it only iterates over a neighborhood of pixels, not the entire image. And as far as I can see, it has the same meaning within my proposed class as it has with the “traditional” ITK neighborhood iterators. We could name it “ModernNeighborhoodRange”, but then we would have to rename it again in a few years, when it’s not that “modern” anymore. Anyway, I’m not opposed to “ModernNeighborhoodRange”, although I like “NeighborhoodRange” better.

BTW, I think NeighborhoodOffsets.h (which is part of my proposed patch) could be adapted to make use of itk::NeighborhoodAllocator, instead of std::vector, if that would bring the patch closer to the existing ITK Neighborhood family… would you like that?

1 Like

Yes, understood.

Most of the work for a class comes after the initial merge, e.g. making it work across all platforms, following up with bugs, etc. This is why it is important to have good testing and documentation, etc.


It at least needs Iterator if it is an iterator and Image if it only applies to images to follow ITK naming conventions.

This should not be a problem if it is not part of the non-Experimental API.

1 Like

:smiley: Patch Set 13 uses the NeighborhoodAccessor of the image. Please check it out!

Hereby I would like to clarify: The proposed itk::NeigborhoodRange class is not an iterator! It has two iterator types, though:


So the following line declares a range, not an iterator:

itk::NeighborhoodRange<ImageType> range{ *image, location, offsets };

But of course, you can get iterators from the range, as follows:

itk::NeighborhoodRange<ImageType>::iterator iterator1 = range.begin();
itk::NeighborhoodRange<ImageType>::iterator iterator2 = range.end();
itk::NeighborhoodRange<ImageType>::const_iterator iterator3 = range.cbegin();
itk::NeighborhoodRange<ImageType>::const_iterator iterator4 = range.cend();

So personally I still think the name of the class is correct. But certainly, placing the class in an Experimental namespace sounds like a good idea. I guess that would make it easier to change the name of the class after the first release, just in case we might regret the initial class name… :slight_smile:

1 Like