Early version of TopologicalConnectivityImageNeighborhoodShape (consider rename to simply "ConnectedImageNeighborhoodShape"!)

itkTopologicalConnectivityImageNeighborhoodShape.h

/*=========================================================================
*
*  Copyright Insight Software Consortium
*
*  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
*
*         http://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 itkTopologicalConnectivityImageNeighborhoodShape_h
#define itkTopologicalConnectivityImageNeighborhoodShape_h

// Author: Niels Dekker, LKEB, Leiden University Medical Center, 2018

#include <cassert>
#include <vector>

#include "itkOffset.h"
#include "itkSize.h"

namespace itk
{
namespace Experimental
{

/**
 * \class TopologicalConnectivityImageNeighborhoodShape
 * Topological connectivity image-neighborhood shape.
 * Eases creating a sequence of offsets for ShapedImageNeighborhoodRange.
 * Can also be used for ShapedNeighborhoodIterator.
 *
 * \see ShapedNeighborhoodIterator
 * \see ShapedImageNeighborhoodRange
 * \ingroup ImageIterators
 * \ingroup ITKCommon
 */
template <unsigned int VImageDimension>
class TopologicalConnectivityImageNeighborhoodShape
{
public:
  static constexpr decltype(VImageDimension) ImageDimension = VImageDimension;

  /** Constructs a topological connectivity shape. Its offsets contain only the
  offset values -1, 0, and 1. The parameter maximumNumberOfNonZeroOffsetValues
  specifies the maximum number of non-zero offset values (might be called Manhattan distance).  */
  constexpr explicit TopologicalConnectivityImageNeighborhoodShape(
    const std::size_t maximumNumberOfNonZeroOffsetValues) ITK_NOEXCEPT
  :
    m_MaximumNumberOfNonZeroOffsetValues{ maximumNumberOfNonZeroOffsetValues },
    m_NumberOfOffsets
    {
      ((maximumNumberOfNonZeroOffsetValues == 0) || (ImageDimension == 0)) ? 0 :
      // Special cases for maximumNumberOfNonZeroOffsetValues, just to allow larger ImageDimension values.
      ((maximumNumberOfNonZeroOffsetValues == 1) ? (2 * ImageDimension) :
       (maximumNumberOfNonZeroOffsetValues >= ImageDimension) ? (CalculatePowerOfThree(ImageDimension) - 1) :
      CalculateSumOfNumberOfHypercubesOnBoundaryOfCube(ImageDimension - 1, (ImageDimension - maximumNumberOfNonZeroOffsetValues))) 
    }
  {
  }


  /** Calculates the number of offsets needed for this shape. */
  constexpr std::size_t GetNumberOfOffsets() const ITK_NOEXCEPT
  {
    return m_NumberOfOffsets;
  }


  /** Fills the specified buffer with the offsets for a neighborhood of this shape. */
  void FillOffsets(
    Offset<ImageDimension>* const offsets) const ITK_NOEXCEPT
  {
    if (m_NumberOfOffsets > 0)
    {
      assert(offsets != nullptr);
      Offset<ImageDimension> offset;
      std::fill_n(offset.begin(), ImageDimension, -1);

      std::size_t result = 0;
      std::size_t i = 0;

      while ( i < m_NumberOfOffsets )
      {
        const auto numberOfNonZeroOffsetValues = ImageDimension - static_cast<std::size_t>(std::count(offset.begin(), offset.end(), 0));

        if ((numberOfNonZeroOffsetValues > 0) && (numberOfNonZeroOffsetValues <= m_MaximumNumberOfNonZeroOffsetValues))
        {
          offsets[i] = offset;
          ++i;
        }

        // Go to the next offset:
        for (decltype(VImageDimension) direction = 0; direction < ImageDimension; ++direction)
        {
          auto& offsetValue = offset[direction];

          ++offsetValue;

          if (offsetValue <= 1)
          {
            break;
          }
          offsetValue = -1;
        }
      };
    }
  }

private:
  // The maximum number of non-zero offset values (might be called maximum
  // Manhattan distance).
  std::size_t m_MaximumNumberOfNonZeroOffsetValues;

  // The number of offsets needed to represent this shape.
  std::size_t m_NumberOfOffsets;

  // Calculate a * b
  static constexpr std::uintmax_t CalculateProduct(
    const std::uintmax_t a,
    const std::uintmax_t b)
  {
    return (((a * b) / a) == b) && (((a * b) / b) == a) ? (a * b) :
      (assert(!"CalculateProduct overflow!"), 0);
  }

  // Calculate a + b
  static constexpr std::uintmax_t CalculateSum(
    const std::uintmax_t a,
    const std::uintmax_t b)
  {
    return ((a + b) >= a) && ((a + b) >= b) ? (a + b) :
      (assert(!"CalculateSum overflow!"), 0);
  }

  // From Walter, June 23 2017:
  // https://stackoverflow.com/questions/44718971/calculate-binomial-coffeficient-very-reliable/44719165#44719165
  static constexpr std::uintmax_t CalculateBinomialCoefficient(const std::uintmax_t n, const std::uintmax_t k) noexcept
  {
    return
      (k > n) ? 0 :          // out of range
      (k == 0 || k == n) ? 1 :          // edge
      (k == 1 || k == n - 1) ? n :          // first
      (k + k < n) ?              // recursive:
      CalculateProduct(CalculateBinomialCoefficient(n - 1, k - 1), n) / k :       //  path to k=1   is faster
      CalculateProduct(CalculateBinomialCoefficient(n - 1, k), n) / (n - k);      //  path to k=n-1 is faster
  }
  // Calculate 3 ^ n
  static constexpr std::uintmax_t CalculatePowerOfThree(const std::size_t n)
  {
    return (n == 0) ? 1 :
      CalculateProduct(3, CalculatePowerOfThree(n - 1));
  }

  // Calculates the number of m-dimensional hypercubes on the boundary of an n-cube.
  // https://en.wikipedia.org/wiki/Hypercube#Elements 
  // Harold Scott MacDonald Coxeter, Regular polytopes, 1947, 1973, p.120
  static constexpr std::size_t CalculateNumberOfHypercubesOnBoundaryOfCube(const std::size_t m, const std::size_t n)
  {
    return CalculateProduct((std::uintmax_t{ 1 } << (n - m)),
      CalculateBinomialCoefficient(n, m));
  }


  // Iterates recursively from i = ImageDimension-1 down to m.
  static constexpr std::size_t CalculateSumOfNumberOfHypercubesOnBoundaryOfCube(const std::size_t i, const std::size_t m)
  {
    return assert(i >= m), CalculateSum(CalculateNumberOfHypercubesOnBoundaryOfCube(i, ImageDimension),
      ((i <= m)? 0 : CalculateSumOfNumberOfHypercubesOnBoundaryOfCube(i - 1, m)));
  }

};

} // namespace Experimental
} // namespace itk

#endif