DigitalSetByOctree.h

 
#pragma once
 
// Inclusions
#include "DGtal/base/Common.h"
#include "DGtal/kernel/domains/HyperRectDomain.h"
 
//#include "DGtal/io/Display3D.h"
 
 
namespace DGtal
{
  template <class Space>
    class DigitalSetByOctree 
    {
      public:
        template<class>
        friend class SVOWriter;
        template<class>
        friend class SVOReader;
 
        using DimIndex = std::uint8_t;
        using CellIndex = typename Space::UnsignedInteger;
        using Size = size_t;
 
        // Only HyperRectDomain are supported 
        using Domain = HyperRectDomain<Space>;
        using Point  = typename Domain::Point;
        
        // Define a value type for this class to be somewhat compatible
        // with image-based code. 
        using Value = void;
 
        static constexpr DGtal::Dimension D = Space::dimension;
        static constexpr DGtal::Dimension dimension = Space::dimension;
        
        static constexpr CellIndex CELL_COUNT  = (1 << D);
        static constexpr CellIndex INVALID_IDX = std::numeric_limits<CellIndex>::max();
 
        // For each child node, store on which side of the octree split it is.
        static constexpr std::array<std::array<DimIndex, D>, CELL_COUNT> SIDES_FROM_INDEX = []()
        {
         std::array<std::array<DimIndex, D>, CELL_COUNT> sides_from_index{};
         for (CellIndex i = 0; i < CELL_COUNT; ++i) 
         {
           CellIndex coord = i;
           for (DimIndex d = 0; d < D; ++d) 
           {
             sides_from_index[i][d] = coord % 2;
             coord /= 2;
           }
         }
 
         return sides_from_index;
        }();

        struct Node 
        {
          Node() 
          {
            for (CellIndex i = 0; i < CELL_COUNT; ++i)
              children[i] = INVALID_IDX;
          }
 
          CellIndex children[CELL_COUNT];
        };
        
        struct ComputationCacheKey 
        {
          CellIndex parentLvl; 
          CellIndex parentIdx;
          std::vector<DimIndex> code;
 
          bool operator<(const ComputationCacheKey& other) const 
          {
            if (parentLvl != other.parentLvl) return parentLvl < other.parentLvl;
            if (parentIdx != other.parentIdx) return parentIdx < other.parentIdx;
            return code < other.code;
          }
        };

        struct TraversalMemory 
        {
          Domain domain; //< Domain represented by the node
          CellIndex lvl; //< Level of the node
          CellIndex idx; //< Index of the node 
          CellIndex currentChildIdx; //< Which child is currently explored.
 
          bool operator==(const TraversalMemory& other) const 
          {
            return lvl == other.lvl && idx == other.idx && currentChildIdx == other.currentChildIdx;
          }
        };

        struct OctreeIterator 
        {
          public:
            friend class DigitalSetByOctree;
 
            // Should be std::output_iterator_tag, but CDigitalSet
            // does not allow this; so we upgrade the category
            using iterator_category = std::forward_iterator_tag;
            using difference_type = std::ptrdiff_t;
            using value_type = Point;
            using reference = value_type&;
            using pointer = value_type*;

            OctreeIterator(const DigitalSetByOctree* container) 
            {
              myContainer = container;
            }

            OctreeIterator(const DigitalSetByOctree* container, 
                           TraversalMemory init) 
            {
              myContainer = container;
              myMemory.push_back(std::move(init));
 
              findNextLeaf();
            }

            OctreeIterator(const DigitalSetByOctree* container, 
                           std::vector<TraversalMemory>& memory) 
            {
              myContainer = container;
              myMemory = memory;
            }

            bool operator==(const OctreeIterator& other) const 
            {
              if (myContainer != other.myContainer) return false;
              if (myMemory.size() != other.myMemory.size()) return false;
 
              if (myMemory.size() != 0) 
              {
                return myMemory.back() == other.myMemory.back();
              }
              return true;
            }

            bool operator!=(const OctreeIterator& other) const 
            {
              return !(*this == other);
            }

            Point operator*() const 
            {
              const auto& sides = SIDES_FROM_INDEX[myMemory.back().currentChildIdx];
              return splitDomain(myMemory.back().domain, sides.data()).lowerBound();
            }

            OctreeIterator& operator++() 
            {
              findNextLeaf();
              return *this;
            }

            OctreeIterator operator++(int) 
            {
              auto it = *this;
              findNextLeaf();
              return it;
            }
          private:
            void findNextLeaf();
          private:
            const DigitalSetByOctree* myContainer; //< Pointer to the original octree
            std::vector<TraversalMemory> myMemory;  //< Current traversal information
        };
 
        using Iterator = OctreeIterator;
        using ConstIterator = OctreeIterator;

        DigitalSetByOctree(const Domain& d);

        const Domain& domain() const { return *myAdmissibleDomain; }

        CowPtr<Domain> domainPointer() const { return myAdmissibleDomain; }
      public:
        void insert(const Point& p);

        void insertNew(const Point& p) { insert(p); }

        bool operator()(const Point& p) const { return find(p) != end(); }

        Iterator find(const Point& p) const;

        size_t size() const { return mySize; }

        bool empty() const { return size() == 0; }

        size_t memoryFootprint() const 
        {
          size_t count = 0;
          for (CellIndex i = 0; i < myNodes.size(); ++i) 
            count += myNodes[i].capacity() * sizeof(Node);
          return count;
        }

        void clear() 
        {
          myNodes.clear();
          myState = State::OCTREE;
        }

        size_t erase(const Iterator& it);

        size_t erase(Iterator begin, Iterator end) 
        {
          size_t count = 0;
          for (; begin != end; ++begin) 
            count += erase(begin);
          return count;
        }

        size_t erase(const Point& p) 
        {
          return erase(find(p));
        }

        void shrinkToFit() 
        {
          for (CellIndex i = 0; i < myNodes.size(); ++i) 
            myNodes[i].shrink_to_fit();
        }

        DigitalSetByOctree& operator+=(const DigitalSetByOctree& other) 
        {
          if (myState == State::OCTREE) 
          {
            for (auto it = other.begin(); it != other.end(); ++it) 
              insert(*it);
          }
          return *this;
        }

        template<typename It>
        void computeComplement(It out) const 
        {
          DigitalSetByOctree complement(*myDomain);
          for (auto it = myDomain->begin(); it != myDomain->end(); ++it) 
          {
            if (!this->operator()(*it)) 
            {
              *out = *it;
              ++out;
            }
          }
        }

        template<class DSet>
        void assignFromComplement(const DSet& other) 
        {
          if (myState == State::OCTREE) 
          {
            clear();
            for (auto it = myDomain->begin(); it != myDomain->end(); ++it) 
            {
              if (!other(*it)) 
                  insert(*it);
            }
          }
        }
        
        void computeBoundingBox(Point& lb, Point& ub) const 
        {
          lb = myDomain->upperBound();
          ub = myDomain->lowerBound();
 
          for (auto it = begin(); it != end(); ++it) 
          {
            const auto p = *it;
            for (DimIndex d = 0; d < D; d++) 
            {
              lb[d] = std::min(lb[d], p[d]);
              ub[d] = std::max(ub[d], p[d]);
            }
          }
        }

        Iterator begin() const 
        {
          TraversalMemory mem;
          mem.domain = *myDomain;
          mem.lvl = 0;
          mem.idx = 0;
          mem.currentChildIdx = INVALID_IDX;
          
          return Iterator(this, mem);
        }

        Iterator begin() { return static_cast<const DigitalSetByOctree&>(*this).begin(); }

        Iterator end() { return Iterator(this); }

        Iterator end() const { return Iterator(this); }
 
    public:
        void convertToDAG();
 
        template<class Func>
        auto computeFunction(OctreeIterator start, OctreeIterator end, CellIndex range, const Func& f);

        void dumpOctree() const;
    private:
        static Domain splitDomain(const Domain& domain, const DimIndex* sides);
    
    private:
        enum class State 
        {
          OCTREE = 1, 
          DAG    = 2
        };
 
 
    private:
        // We store two domains:
        //  - myDomain is the domain used for computation, that 
        //    requires its size to be whole power of two for 
        //    computationnal purposes
        //  - myAdmissibleDomain is the domain of points that can be
        //    stored within the octree. It shares the same lower bound
        //    with domain, but its upper bound is lowered by one.
        CowPtr<Domain> myAdmissibleDomain;     // Domain of point that are allowed. 
        CowPtr<Domain> myDomain;       // Pointer to domain, as required by CDigitalSet.
        State myState = State::OCTREE; // Current state
        
        size_t mySize;                 // Current number of stored voxel.
        std::vector<std::vector<Node>> myNodes; // Nodes storage
    };
};
 
#include "DigitalSetByOctree.ih"