EstimateAxisAlignedBoundingBoxRegion (ITK5)

EstimateAxisAlignedBoundingBoxRegion.hxx

/*
Estimates the axis-aligned minimum bounding box of the buffered region of an image.

Roughly equivalent to ImageMaskSpatialObject::GetAxisAlignedBoundingBoxRegion() from
"itkImageMaskSpatialObject.hxx", at
github.com/InsightSoftwareConsortium/ITK/blob/91ab86f63eb5fc9799bb25fa49e76adf0d1b5251/Modules/Core/SpatialObjects/include

Related to:

"ENH: GTest for ImageMaskSpatialObject::GetAxisAlignedBoundingBoxRegion()"
https://github.com/InsightSoftwareConsortium/ITK/pull/681

"ENH: Add legacy interface to SpatialObjects"
https://github.com/InsightSoftwareConsortium/ITK/pull/676/

"ENH: Added ImageSpatialObject::GetAxisAlignedBoundingBoxRegion()"
https://github.com/InsightSoftwareConsortium/ITK/pull/645

"How to replace ImageMaskSpatialObject::GetAxisAlignedBoundingBoxRegion() calls?"
https://discourse.itk.org/t/how-to-replace-imagemaskspatialobject-getaxisalignedboundingboxregion-calls/1709

Niels Dekker, LKEB, Leiden University Medical Center, 2019

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
*/

#include <itkImage.h>
#include <itkImageRegion.h>
#include <itkImageRegionConstIterator.h>

namespace LKEB
{
  template<typename TImage>
  typename TImage::RegionType EstimateAxisAlignedBoundingBoxRegion(const TImage& image)
  {
    using namespace itk;

    using PixelType = typename TImage::PixelType;
    constexpr auto VImageDimension = TImage::ImageDimension;

    using RegionType = ImageRegion<VImageDimension>;
    using SizeType = Size<VImageDimension>;
    using IndexType = Index<VImageDimension>;

    const auto HasForegroundPixels = [&image](const RegionType& region)
    {
      for (ImageRegionConstIterator<TImage> it{ &image, region }; !it.IsAtEnd(); ++it)
      {
        constexpr auto zeroValue = NumericTraits<PixelType>::ZeroValue();

        if (it.Get() != zeroValue)
        {
          return true;
        }
      }
      return false;
    };

    const auto CreateRegion = [](const IndexType& minIndex, const IndexType& maxIndex)
    {
      SizeType regionSize;

      for (unsigned dim = 0; dim < SizeType::Dimension; ++dim)
      {
        regionSize[dim] = static_cast<SizeValueType>(maxIndex[dim] + 1 - minIndex[dim]);
      }
      return RegionType{ minIndex, regionSize };
    };

    const RegionType requestedRegion = image.GetRequestedRegion();

    if (requestedRegion.GetNumberOfPixels() == 0)
    {
      return {};
    }

    const SizeType imageSize = requestedRegion.GetSize();

    IndexType minIndex = requestedRegion.GetIndex();
    IndexType maxIndex = minIndex + imageSize;

    for (auto& maxIndexValue : maxIndex)
    {
      --maxIndexValue;
    }

    // Iterate from high to low (for significant performance reasons).
    for (int dim = VImageDimension - 1; dim >= 0; --dim)
    {
      auto subregion = CreateRegion(minIndex, maxIndex);
      subregion.SetSize(dim, 1);
      const auto initialMaxIndexValue = maxIndex[dim];

      // Estimate minIndex[dim]
      while(!HasForegroundPixels(subregion))
      {
        const auto indexValue = subregion.GetIndex(dim) + 1;

        if (indexValue > initialMaxIndexValue)
        {
          // The requested image region has only zero-valued pixels.
          return {};
        }
        subregion.SetIndex(dim, indexValue);
      }
      minIndex[dim] = subregion.GetIndex(dim);

      // Estimate maxIndex[dim]
      subregion.SetIndex(dim, initialMaxIndexValue);
      while (!HasForegroundPixels(subregion))
      {
        subregion.SetIndex(dim, subregion.GetIndex(dim) - 1);
      }
      maxIndex[dim] = subregion.GetIndex(dim);
    }
    return CreateRegion(minIndex, maxIndex);
  }
}

Nice!

One suggestion - line 74 makes a pass through the image, verifying that there is at least one non-zero voxel. That makes this function O(2N-1) in the worse case. Also, if the masked region is small, then the actual runtime is near O(2N-c). Perhaps it would be better to eliminate the test in line 74, and instead make the loops in lines 89 and 102 conditional on minIndex[i]<maxIndex[i] and vice versa?

Thanks for doing such a thorough code reading, @aylward I found it easier to just skip any image that are entirely black, to keep the loop conditions simple. But it may indeed be worth considering to eliminate the test in line 74, and adjust the loops, as you suggested. I would have to see the performance results.

At least, this implementation (revision 5) is already a huge performance improvement, compared to the old (legacy) ITK GetAxisAlignedBoundingBoxRegion()! I did a benchmark:

2D (1024x1024):
Entirely black image: > 7 x as fast
Black except for 1 white center pixel: > 4 x as fast
Entirely white image: > 5000 x as fast!!!

3D (256x256x256):
Entirely black image: > 80 x as fast
Black except for 1 white center voxel: > 30 x as fast
Entirely white image: > 30 x as fast

Quite promising, isn't it? :-) A problem with further optimization could be that it improved the performance of one use case at the cost of another use case. Should I care about the performance on an entirely black image, or is that an irrelevant use case?

@aylward Revision 6 includes your suggestion to skip verifying that there is at least one non-zero voxel :-)

Revision 7 has added a check if the requested region is empty, and somewhat simplified the code.