Itk: ConstShapedNeighborhoodIterator doesn't loop over all possibilities

[Leengit]

In a region of 12×12 pixels, an iterator is expected to hit all 144 possibilities. However it gets only 142 of them.

Steps to Reproduce

#include "itkImage.h"
#include "itkConstShapedNeighborhoodIterator.h"
#include "itkConstantBoundaryCondition.h"

int man(int argc, char * argv[])
{
    using PixelType = unsigned char;
    constexpr unsigned int Dimension = 2;
    using ImageType = itk::Image<PixelType, Dimension>;
    using ImagePointer = typename ImageType::Pointer;
    using RegionType = typename ImageType::RegionType;
    using IndexType = typename RegionType::IndexType;
    using SizeType = typename RegionType::SizeType;
    using OffsetType = typename RegionType::OffsetType;
    ImagePointer imageTest = ImageType::New();
    imageTest->SetRegions(RegionType{ IndexType{ 0, 0 }, SizeType{ 10, 10 } });
    imageTest->Allocate();
    using BoundaryConditionType = itk::ConstantBoundaryCondition<ImageType>;
    using SquareIterator = itk::ConstShapedNeighborhoodIterator<ImageType, BoundaryConditionType>;
    using RadiusType = typename SquareIterator::RadiusType;
    SquareIterator sqIt(RadiusType{ 1, 1 }, imageTest, RegionType{ IndexType{ -1, -1 }, SizeType{ 12, 12 } });
    sqIt.ActivateOffset(OffsetType{ 0, 0 });
    sqIt.ActivateOffset(OffsetType{ 0, 1 });
    sqIt.ActivateOffset(OffsetType{ 1, 0 });
    sqIt.ActivateOffset(OffsetType{ 1, 1 });
    int numberOfIterations{ 0 };
    for (sqIt.GoToBegin(); !sqIt.IsAtEnd(); ++sqIt)
    {
      ++numberOfIterations;
    }
    std::cout << "numberOfIterations = " << numberOfIterations << " but should be 144." << std::endl;
}

Expected behavior

Prints numberOfIterations = 144 but should be 144.

Actual behavior

Prints numberOfIterations = 142 but should be 144.

Reproducibility

100%

Versions

Master branch (currently 25dea2f73d974855fbdfb13f5f2d7f25e4ec8a4e).

Environment

Ubuntu 20.04. gcc 10.2.0.

Additional information.

I am not completely sure by any means, but maybe the problem is that the region for the square iterator includes values with x < 0. Because of this, the .end() value that is calculated as the left-most position of the row after the region can be misinterpreted as being near the right end of the previous row. If that's true, maybe it is also true that a computed WrapOffset is negative when it is implicitly assumed to always be positive.

[Leengit]

The above testing code will be added as Test 6 to itkContourExtractor2DImageFilterTest.cxx with Pull Request #2387. That pull request will fail Tests 5 and 6 until this issue is addressed.

[N-Dekker]

FWIW, The "range based" version would look something like this:

const std::array<OffsetType, 4> offsets = { { { 0, 0 }, { 0, 1 }, { 1, 0 }, { 1, 1 } } };
using Policy = itk::ConstantBoundaryImageNeighborhoodPixelAccessPolicy<ImageType>;
itk::ShapedImageNeighborhoodRange<ImageType, Policy> neighborhoodRange{ *imageTest, IndexType{}, offsets };

for (const auto index :
     itk::ImageRegionIndexRange<Dimension>(RegionType{ IndexType{ -1, -1 }, SizeType{ 12, 12 } }))
{
  neighborhoodRange.SetLocation(index);
  ++numberOfIterations;
}
std::cout << "numberOfIterations = " << numberOfIterations << " (144, indeed!)" << std::endl;

It does 144 iterations, at least according to my (local) observations.

[N-Dekker]

@Leengit I just added a print statement to your for loop, saying std::cout << sqIt.GetIndex(); It appears that the last two indices ([9, 10] and [10, 10]) aren't iterated by ConstShapedNeighborhoodIterator Interesting... 🤔

[Leengit]

@N-Dekker In the "range based" code example, how would I set the policy's constant? -- Never mind; it seems to be a fourth parameter to the itk::ShapedImageNeighborhoodRange<const InputImageType, itk::ConstantBoundaryImageNeighborhoodPixelAccessPolicy<InputImageType> > constructor.

[N-Dekker]

I'm not really sure, but it might be that ConstShapedNeighborhoodIterator never fully supported having the center pixel of the neighborhood outside the buffered region of the image. So then the question might be, should this feature (support of center pixels outside the image) still be added to this "old" iterator class?

FYI The initial reason for me to start developing the new "range based" ITK classes was to get a better performance. But honestly, in the mean time, the performance of the old iterator classes was also greatly improved, especially by making many of their member functions non-virtual. So now it actually depends on the use case, which of the two approaches is fastest. Of course, personally I would recommend using the "range based" classes by default 😃

[Leengit]

The above testing code will be added as Test 6 to itkContourExtractor2DImageFilterTest.cxx with Pull Request #2387. That pull request will fail Tests 5 and 6 until this issue is addressed.

The Test 6 has been commented out. When this issue is addressed that test should be reinstated, though not necessarily in that file. Test 5 no longer depends upon this issue being resolved; its implementation was changed to make use of a range rather than an iterator and that bypasses this issue.