itkConnectedImageNeighborhoodShape.h

Go to the documentation of this file.

/*=========================================================================
 *
 *  Copyright NumFOCUS
 *
 *  Licensed under the Apache License, Version 2.0 (the "License");
 *  you may not use this file except in compliance with the License.
 *  You may obtain a copy of the License at
 *
 *         https://www.apache.org/licenses/LICENSE-2.0.txt
 *
 *  Unless required by applicable law or agreed to in writing, software
 *  distributed under the License is distributed on an "AS IS" BASIS,
 *  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 *  See the License for the specific language governing permissions and
 *  limitations under the License.
 *
 *=========================================================================*/
 
#ifndef itkConnectedImageNeighborhoodShape_h
#define itkConnectedImageNeighborhoodShape_h
 
#include "itkIntTypes.h" // For uintmax_t
#include "itkMath.h"
#include "itkOffset.h"
 
#include <array>
#include <cassert>
#include <limits>
 
namespace itk
{
 
template <unsigned int VImageDimension>
class ConnectedImageNeighborhoodShape
{
public:
  static constexpr unsigned int ImageDimension = VImageDimension;
  constexpr explicit ConnectedImageNeighborhoodShape(const size_t maximumCityblockDistance,
                                                     const bool   includeCenterPixel) noexcept
    : m_MaximumCityblockDistance{ maximumCityblockDistance }
    , m_IncludeCenterPixel{ includeCenterPixel }
    , m_NumberOfOffsets{ CalculateNumberOfOffsets(maximumCityblockDistance, includeCenterPixel) }
  {}
  constexpr size_t
  GetNumberOfOffsets() const noexcept
  {
    return m_NumberOfOffsets;
  }
 
 
  void
  FillOffsets(Offset<ImageDimension> * const offsets) const noexcept
  {
    if (m_NumberOfOffsets > 0)
    {
      assert(offsets != nullptr);
      Offset<ImageDimension> offset;
      std::fill_n(offset.begin(), ImageDimension, -1);
      size_t i = 0;
 
      while (i < m_NumberOfOffsets)
      {
        const size_t numberOfNonZeroOffsetValues =
          ImageDimension - static_cast<size_t>(std::count(offset.begin(), offset.end(), 0));
 
        if ((m_IncludeCenterPixel || (numberOfNonZeroOffsetValues > 0)) &&
            (numberOfNonZeroOffsetValues <= m_MaximumCityblockDistance))
        {
          offsets[i] = offset;
          ++i;
        }
 
        // Go to the next offset:
        for (unsigned int direction = 0; direction < ImageDimension; ++direction)
        {
          auto & offsetValue = offset[direction];
 
          ++offsetValue;
 
          if (offsetValue <= 1)
          {
            break;
          }
          offsetValue = -1;
        }
      }
    }
  }
 
private:
  // The maximum city-block distance (Manhattan distance) between the center
  // pixel and each connected neighbor pixel.
  size_t m_MaximumCityblockDistance;
 
  // Specifies whether or not the center pixel (offset zero) should be included
  // with the offsets for this shape.
  bool m_IncludeCenterPixel;
 
  // The number of offsets needed to represent this shape.
  size_t m_NumberOfOffsets;
 
 
  // Calculates a + b. Numeric overflow triggers a compilation error in
  // "constexpr context" and a debug assert failure at run-time.
  static constexpr uintmax_t
  CalculateSum(const uintmax_t a, const uintmax_t b) noexcept
  {
    return ((a + b) >= a) && ((a + b) >= b) ? (a + b) : (assert(!"CalculateSum overflow!"), 0);
  }
 
 
  // Calculates 2 ^ n. Numeric overflow triggers a compilation error in
  // "constexpr context" and a debug assert failure at run-time.
  static constexpr uintmax_t
  CalculatePowerOfTwo(const size_t n) noexcept
  {
    return (n < std::numeric_limits<uintmax_t>::digits) ? (uintmax_t{ 1 } << n)
                                                        : (assert(!"CalculatePowerOfTwo overflow!"), 0);
  }
 
 
  // Calculates the binomial coefficient, 'n' over 'k'.
  // Inspired by the 'binom' function from Walter, June 23, 2017:
  // https://stackoverflow.com/questions/44718971/calculate-binomial-coffeficient-very-reliable/44719165#44719165
  // Optimized for small values of 'k' (k <= n/2).
  static constexpr uintmax_t
  CalculateBinomialCoefficient(const uintmax_t n, const uintmax_t k) noexcept
  {
    return (k > n)    ? (assert(!"Out of range!"), 0)
           : (k == 0) ? 1
                      : Math::UnsignedProduct(n, CalculateBinomialCoefficient(n - 1, k - 1)) / k;
  }
 
 
  // Calculates the number of m-dimensional hypercubes on the boundary of an
  // n-cube, as described at https://en.wikipedia.org/wiki/Hypercube#Elements
  // (Which refers to H.S.M. Coxeter, Regular polytopes, 3rd ed., 1973, p.120.)
  static constexpr size_t
  CalculateNumberOfHypercubesOnBoundaryOfCube(const size_t m, const size_t n) noexcept
  {
    // Calculate 2^(n-m) * BinomialCoefficient(n, m)
    return Math::UnsignedProduct(CalculatePowerOfTwo(n - m),
                                 (((2 * m) < n) ?
                                                // Calculate either the binomial coefficient of (n, m) or (n, n - m).
                                                // Mathematically, both should yield the same number, but the
                                                // calculation is optimized for a smaller second argument.
                                    CalculateBinomialCoefficient(n, m)
                                                : CalculateBinomialCoefficient(n, n - m)));
  }
 
 
  // Iterates recursively from i = ImageDimension-1 down to m (inclusive).
  static constexpr size_t
  CalculateSumOfNumberOfHypercubesOnBoundaryOfCube(const size_t i, const size_t m) noexcept
  {
    return assert(i >= m), CalculateSum(CalculateNumberOfHypercubesOnBoundaryOfCube(i, ImageDimension),
                                        ((i <= m) ? 0 : CalculateSumOfNumberOfHypercubesOnBoundaryOfCube(i - 1, m)));
  }
 
 
  static constexpr size_t
  CalculateNumberOfConnectedNeighbors(const size_t maximumCityblockDistance) noexcept
  {
    return (((maximumCityblockDistance == 0) || (ImageDimension == 0))
              ? 0
              : ((maximumCityblockDistance >= ImageDimension)
                   ? (Math::UnsignedPower(3, ImageDimension) - 1)
                   : CalculateSumOfNumberOfHypercubesOnBoundaryOfCube(ImageDimension - 1,
                                                                      (ImageDimension - maximumCityblockDistance))));
  }
 
 
  static constexpr size_t
  CalculateNumberOfOffsets(const size_t maximumCityblockDistance, const bool includeCenterPixel) noexcept
  {
    return (includeCenterPixel ? 1 : 0) + CalculateNumberOfConnectedNeighbors(maximumCityblockDistance);
  }
};
template <unsigned int VImageDimension, size_t VMaximumCityblockDistance, bool VIncludeCenterPixel>
auto
GenerateConnectedImageNeighborhoodShapeOffsets() noexcept
{
  constexpr ConnectedImageNeighborhoodShape<VImageDimension> shape{ VMaximumCityblockDistance, VIncludeCenterPixel };
  std::array<Offset<VImageDimension>, shape.GetNumberOfOffsets()> offsets;
  shape.FillOffsets(offsets.data());
  return offsets;
}
} // namespace itk
 
#endif