Skip to content

Instantly share code, notes, and snippets.

Show Gist options
  • Save N-Dekker/0f17c0473b29494faea93244bbb43393 to your computer and use it in GitHub Desktop.
Save N-Dekker/0f17c0473b29494faea93244bbb43393 to your computer and use it in GitHub Desktop.
Early version of TopologicalConnectivityImageNeighborhoodShape (consider rename to simply "ConnectedImageNeighborhoodShape"!)
/*=========================================================================
*
* 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
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment