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

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