DomainSplitter.h

 
#pragma once
 
#include "DGtal/kernel/domains/CDomain.h"
#include <vector>
 
namespace DGtal
{
  template<typename Domain>
  struct SplitInfo
  {
    BOOST_CONCEPT_ASSERT(( concepts::CDomain< Domain > ));
 
    Domain domain; //< The actual split domain
    uint32_t hintVoxelCount = 0; //< An expected guess for the number of voxel
  };
 
  template<typename Domain>
  struct RegularDomainSplitter
  {
    BOOST_CONCEPT_ASSERT(( concepts::CDomain< Domain > ));
 
    //Output spllitted domain type
    typedef std::vector<SplitInfo<Domain>> SplitDomainsInfo;
 
    SplitDomainsInfo operator()(const Domain& d, uint32_t splitHint) const
    {
      // Find best match possible for even splitting
      const uint32_t splitCount = std::floor(std::log(splitHint) / std::log(Domain::dimension));
      const uint32_t totalSplits = std::pow(splitCount, Domain::dimension);
 
      if (splitCount == 0)
        return { SplitInfo{d, 0} };
 
      auto splitSize = (d.upperBound() - d.lowerBound()) / (int32_t)splitCount;
      SplitDomainsInfo result;
      result.reserve(totalSplits);
 
      for (uint32_t i = 0; i < totalSplits; ++i)
      {
        auto start = d.lowerBound();
        auto idx = i;
        for (uint32_t j = 0; j < Domain::dimension; ++j)
        {
          auto k = idx % splitCount;
          start[j] += k * splitSize[j] + k; // +k ensure no overlap between domains
          idx /= splitCount;
        }
 
        // Make correction to ensure it remains within the domain
        auto end = start + splitSize;
        for (uint32_t j = 0; j < Domain::dimension; ++j)
          end[j] = std::clamp(end[j], d.lowerBound()[j], d.upperBound()[j]);
 
        result.emplace_back(Domain(start, end), 0);
      }
 
      return result;
    };
  };
 
 
  template<typename Domain>
  struct AxisDomainSplitter
  {
 
    BOOST_CONCEPT_ASSERT(( concepts::CDomain< Domain > ));
 
    //Output spllitted domain type
    typedef std::vector<SplitInfo<Domain>> SplitDomainsInfo;
    typedef typename Domain::Dimension     Dimension;
 
    Dimension axis;
 
    explicit AxisDomainSplitter( Dimension dim = 0 )
      : axis( dim )
    {}
 
    SplitDomainsInfo operator()( const Domain& d, uint32_t splitHint ) const
    {
      return (*this)( d, splitHint, axis );
    }
 
    SplitDomainsInfo operator()(const Domain& d, uint32_t splitHint, Dimension dim) const
    {
      SplitDomainsInfo result;
      if (splitHint == 0)
        return result;
 
      auto lower = d.lowerBound();
      auto upper = d.upperBound();
      auto length = upper[dim] - lower[dim] + 1;
      uint32_t splitCount = splitHint;
      if (splitCount > length)
        splitCount = length;
 
      result.reserve(splitCount);
      auto base = length / splitCount;
      auto rem = length % splitCount;
 
      auto start = lower;
      for (uint32_t i = 0; i < splitCount; ++i)
      {
        auto size = base + (i < rem ? 1 : 0);
        auto end = upper;
        end[dim] = start[dim] + size - 1;
        result.emplace_back(Domain(start, end), 0);
        start[dim] = end[dim] + 1;
      }
 
      return result;
    };
  };
 
}