AACubeTree.h

Go to the documentation of this file.

#ifndef GfxTL__AACUBETREE_HEADER__
#define GfxTL__AACUBETREE_HEADER__
#include <deque>
#include <iostream>
 
#include <GfxTL/AABox.h>
#include <GfxTL/AACube.h>
#include <GfxTL/BaseTree.h>
#include <GfxTL/NullClass.h>
#include <GfxTL/VectorKernel.h>
 
namespace GfxTL
{
    template <unsigned int DimT, class BaseT>
    class AACubeTreeCell : public BaseT
    {
    public:
        enum
        {
            Dim = DimT,
            NChildren = 1 << DimT
        };
 
        typedef BaseT BaseType;
        typedef AACubeTreeCell<DimT, BaseT> ThisType;
 
        AACubeTreeCell()
        {
            for (size_t i = 0; i < 1 << DimT; ++i)
            {
                m_children[i] = NULL;
            }
        }
 
        ~AACubeTreeCell()
        {
            for (size_t i = 0; i < 1 << DimT; ++i)
                if (m_children[i] > (ThisType*)1)
                {
                    delete m_children[i];
                }
        }
 
        ThisType&
        operator[](unsigned int i)
        {
            return *m_children[i];
        }
 
        const ThisType&
        operator[](unsigned int i) const
        {
            return *m_children[i];
        }
 
        void
        Child(unsigned int i, ThisType* cell)
        {
            m_children[i] = cell;
        }
 
        ThisType** const
        Children()
        {
            return m_children;
        }
 
        const ThisType* const* const
        Children() const
        {
            return m_children;
        }
 
    private:
        ThisType* m_children[1 << DimT];
    };
 
    template <class DataStrategyT>
    struct BaseAACubeTreeStrategy
    {
        typedef typename DataStrategyT::value_type value_type;
 
        class CellData : public DataStrategyT::CellData
        {
        };
 
        template <class BaseT>
        class StrategyBase : public DataStrategyT::template StrategyBase<BaseT>
        {
        public:
            typedef typename DataStrategyT::template StrategyBase<BaseT> BaseType;
            typedef typename BaseType::CellType CellType;
 
            enum
            {
                Dim = CellType::Dim
            };
 
            typedef typename ScalarTypeDeferer<value_type>::ScalarType ScalarType;
            typedef AACube<VectorXD<Dim, ScalarType>> CubeType;
            typedef CubeType CellCubeType;
            typedef GfxTL::VectorXD<Dim, ScalarType> CellCenterType;
 
            const CubeType&
            RootCube() const
            {
                return m_rootCube;
            }
 
        protected:
            template <class BuildInformationT>
            bool
            ShouldSubdivide(const CellType& cell, BuildInformationT& bi) const
            {
                return false;
            }
 
            template <class BuildInformationT>
            void
            InitRootBuildInformation(BuildInformationT* bi)
            {
                BaseType::InitRootBuildInformation(bi);
                m_rootCube = bi->Cube();
            }
 
            template <class BuildInformationT>
            void
            InitGlobalBuildInformation(const CellType& root, const BuildInformationT& bi)
            {
            }
 
            template <class BuildInformationT>
            void
            EnterGlobalBuildInformation(const CellType& cell, BuildInformationT* bi) const
            {
            }
 
            template <class BuildInformationT>
            void
            LeaveGlobalBuildInformation(const CellType& cell, const BuildInformationT& bi) const
            {
            }
 
            template <class BuildInformationT>
            void
            InitLeaf(CellType* cell, const BuildInformationT& bi)
            {
            }
 
            template <class BuildInformationT>
            void
            InitSubdivided(const BuildInformationT& bi, CellType* cell)
            {
            }
 
            template <class TraversalBaseT>
            class CellCubeTraversalInformation : public TraversalBaseT
            {
            public:
                typedef AACube<VectorXD<Dim, ScalarType>> CubeType;
 
                CubeType&
                Cube()
                {
                    return m_cube;
                }
 
                const CubeType&
                Cube() const
                {
                    return m_cube;
                }
 
            private:
                CubeType m_cube;
            };
 
            template <class TraversalBaseT>
            class CellCenterTraversalInformation : public TraversalBaseT
            {
            public:
                typedef typename TraversalBaseT::CubeType CubeType;
 
                CubeType&
                Cube()
                {
                    return m_cube;
                }
 
                const CubeType&
                Cube() const
                {
                    return m_cube;
                }
 
            private:
                CubeType m_cube;
            };
 
            template <class TraversalBaseT>
            class TraversalInformationBase :
                public DataStrategyT::template StrategyBase<BaseT>::template TraversalInformation<
                    TraversalBaseT>
            {
            };
 
            template <class TraversalInformationT>
            void
            InitRootTraversalInformation(const CellType& root, TraversalInformationT* ti) const
            {
                BaseType::InitRootTraversalInformation(root, ti);
                InitCubeRootTraversalInformation(root, ti);
                InitCenterRootTraversalInformation(root, ti);
            }
 
            template <class TraversalInformationT>
            void
            InitTraversalInformation(const CellType& parent,
                                     const TraversalInformationT& pTi,
                                     unsigned int childIdx,
                                     TraversalInformationT* ti) const
            {
                BaseType::InitTraversalInformation(parent, pTi, childIdx, ti);
                InitCubeTraversalInformation(parent, pTi, childIdx, ti);
                InitCenterTraversalInformation(parent, pTi, childIdx, ti);
            }
 
            template <class TraversalBaseT>
            void
            CellCube(const CellType& cell,
                     const CellCubeTraversalInformation<TraversalBaseT>& ti,
                     CellCubeType* cube) const
            {
                *cube = ti.Cube();
            }
 
            template <class TraversalBaseT>
            void
            CellCenter(const CellType& cell,
                       const CellCenterTraversalInformation<TraversalBaseT>& ti,
                       CellCenterType* center) const
            {
                ti.Cube().Center(center);
            }
 
        private:
            template <class TraversalInformationT>
            void
            InitCubeRootTraversalInformation(const CellType& root, TraversalInformationT* ti) const
            {
            }
 
            template <class TraversalBaseT>
            void
            InitCubeRootTraversalInformation(const CellType& root,
                                             CellCubeTraversalInformation<TraversalBaseT>* ti) const
            {
                ti->Cube() = m_rootCube;
            }
 
            template <class TraversalInformationT>
            void
            InitCubeTraversalInformation(const CellType& parent,
                                         const TraversalInformationT& pTi,
                                         unsigned int childIdx,
                                         TraversalInformationT* ti) const
            {
            }
 
            template <class TraversalBaseT>
            void
            InitCubeTraversalInformation(const CellType& parent,
                                         const CellCubeTraversalInformation<TraversalBaseT>& pTi,
                                         unsigned int childIdx,
                                         CellCubeTraversalInformation<TraversalBaseT>* ti) const
            {
                ti->Cube().Width(pTi.Cube().Width() / ScalarType(2));
                for (unsigned int i = 0; i < Dim; ++i)
                    if (childIdx & (1 << (Dim - 1 - i)))
                        ti->Cube().Min()[i] = pTi.Cube().Min()[i] + ti->Cube().Width();
                    else
                    {
                        ti->Cube().Min()[i] = pTi.Cube().Min()[i];
                    }
            }
 
            template <class TraversalInformationT>
            void
            InitCenterRootTraversalInformation(const CellType& root,
                                               TraversalInformationT* ti) const
            {
            }
 
            template <class TraversalBaseT>
            void
            InitCenterRootTraversalInformation(
                const CellType& root,
                CellCenterTraversalInformation<TraversalBaseT>* ti) const
            {
                ti->Cube() = m_rootCube;
            }
 
            template <class TraversalInformationT>
            void
            InitCenterTraversalInformation(const CellType& parent,
                                           const TraversalInformationT& pTi,
                                           unsigned int childIdx,
                                           TraversalInformationT* ti) const
            {
            }
 
            template <class TraversalBaseT>
            void
            InitCenterTraversalInformation(
                const CellType& parent,
                const CellCenterTraversalInformation<TraversalBaseT>& pTi,
                unsigned int childIdx,
                CellCenterTraversalInformation<TraversalBaseT>* ti) const
            {
                ti->Cube().Width(pTi.Cube().Width() / ScalarType(2));
                for (unsigned int i = 0; i < Dim; ++i)
                    if (childIdx & (1 << (Dim - 1 - i)))
                        ti->Cube().Min()[i] = pTi.Cube().Min()[i] + ti->Cube().Width();
                    else
                    {
                        ti->Cube().Min()[i] = pTi.Cube().Min()[i];
                    }
            }
 
        private:
            CubeType m_rootCube;
        };
    };
 
    template <unsigned int DimT,
              class StrategiesT,
              template <unsigned int> class VectorKernelT = VectorKernelD>
    class AACubeTree :
        public StrategiesT::template StrategyBase<
            typename VectorKernelT<DimT>::template VectorKernelType<
                BaseTree<AACubeTreeCell<DimT, typename StrategiesT::CellData>>>>
    {
    public:
        typedef StrategiesT StrategiesType;
        typedef AACubeTreeCell<DimT, typename StrategiesT::CellData> CellType;
        typedef typename StrategiesT::template StrategyBase<
            typename VectorKernelT<DimT>::template VectorKernelType<
                BaseTree<AACubeTreeCell<DimT, typename StrategiesT::CellData>>>>
            BaseType;
        typedef AACubeTree<DimT, StrategiesT, VectorKernelT> ThisType;
        typedef typename BaseType::value_type value_type;
        typedef typename ScalarTypeDeferer<value_type>::ScalarType ScalarType;
        typedef typename BaseType::CellRange CellRange;
 
        void
        Build()
        {
            AABox<VectorXD<DimT, ScalarType>> bbox;
            bbox.Bound(BaseType::begin(), BaseType::end());
            ScalarType width = bbox.Max()[0] - bbox.Min()[0];
            for (unsigned int i = 1; i < DimT; ++i)
            {
                width = std::max(width, bbox.Max()[i] - bbox.Min()[i]);
            }
            AACube<VectorXD<DimT, ScalarType>> bcube(bbox.Min(), width);
            Build(bcube);
        }
 
        void
        Build(const AACube<VectorXD<DimT, ScalarType>>& bcube)
        {
            typedef std::pair<CellType*, BuildInformation> Pair;
            BaseType::Clear();
            BaseType::Root() = new CellType;
            std::deque<Pair> stack(1);
            // init build information directly on stack to avoid
            // copying.
            stack.back().first = BaseType::Root();
            this->InitRootBuildInformation(bcube, &stack.back().second);
            this->InitRoot(stack.back().second, BaseType::Root());
            BaseType::InitGlobalBuildInformation(*BaseType::Root(), stack.back().second);
            while (stack.size())
            {
                Pair& p = stack.back();
                if (p.second.CreateChild() == 1 << DimT)
                {
                    BaseType::LeaveGlobalBuildInformation(*p.first, p.second);
                    stack.pop_back();
                    continue;
                }
                if (this->IsLeaf(*p.first))
                {
                    if (!this->ShouldSubdivide(*p.first, p.second))
                    {
                        BaseType::InitLeaf(p.first, p.second);
                        stack.pop_back();
                        continue;
                    }
                    Subdivide(p.second, p.first);
                    if (this->IsLeaf(*p.first)) // couldn't subdivide?
                    {
                        BaseType::InitLeaf(p.first, p.second);
                        stack.pop_back();
                        continue;
                    }
                    BaseType::InitSubdivided(p.second, p.first);
                }
                else
                {
                    BaseType::LeaveGlobalBuildInformation(*p.first, p.second);
                }
                while (p.second.CreateChild() < (1 << DimT) &&
                       !this->ExistChild(*p.first, p.second.CreateChild()))
                {
                    ++p.second.CreateChild();
                }
                if (p.second.CreateChild() == (1 << DimT))
                {
                    stack.pop_back();
                    continue;
                }
                BaseType::EnterGlobalBuildInformation(*p.first, &p.second);
                stack.resize(stack.size() + 1); // create new entry
                stack.back().first = &(*p.first)[p.second.CreateChild()];
                this->InitBuildInformation(
                    *p.first, p.second, p.second.CreateChild(), &stack.back().second);
                this->InitCell(*p.first,
                               p.second,
                               p.second.CreateChild(),
                               stack.back().second,
                               &(*p.first)[p.second.CreateChild()]);
                do
                {
                    ++p.second.CreateChild();
                } while (p.second.CreateChild() < (1 << DimT) &&
                         !this->ExistChild(*p.first, p.second.CreateChild()));
            }
        }
 
        template <class PointT>
        const CellType*
        NodeContainingPoint(const PointT& point,
                            size_t maxLevel,
                            size_t minSize,
                            CellRange* range) const
        {
            if (!BaseType::Root())
            {
                range->first = 0;
                range->second = 0;
                return NULL;
            }
            typedef typename BaseType::template CellLevelTraversalInformation<
                typename BaseType::template CellCenterTraversalInformation<
                    typename BaseType::template TraversalInformationBase<NullClass>>>
                TraversalInfoType;
            TraversalInfoType rti;
            BaseType::InitRootTraversalInformation(*BaseType::Root(), &rti);
            BaseType::GetCellRange(*BaseType::Root(), rti, range);
            return NodeContainingPoint(point, maxLevel, minSize, *BaseType::Root(), rti, range);
        }
 
        void
        Serialize(std::ostream* out) const
        {
            typedef std::pair<const CellType*, SerializeInformation> Pair;
            std::deque<Pair> stack(1);
            // init build information directly on stack to avoid
            // copying.
            stack.back().first = BaseType::Root();
            InitRootSerializeInformation(&stack.back().second);
            while (stack.size())
            {
                Pair& p = stack.back();
                if (p.second.CreateChild() == 1 << DimT)
                {
                    stack.pop_back();
                    continue;
                }
                if (!ShouldSubdivide(*p.first, p.second))
                {
                    stack.pop_back();
                    continue;
                }
                if (!p.second.Written())
                {
                    const size_t byteCount =
                        ((1 << BaseType::m_dim) / 8) + (((1 << BaseType::m_dim) % 8) ? 1 : 0);
                    char b[byteCount];
                    for (unsigned int i = 0; i < (1 << DimT); ++i)
                        if (ExistChild(*p.first, i))
                        {
                            b[i >> 3] |= 1 << (i - ((i >> 3) << 3));
                        }
                        else
                        {
                            b[i >> 3] &= ~(1 << (i - ((i >> 3) << 3)));
                        }
                    out->write(b, byteCount);
                    p.second.Written(true);
                }
                else
                {
                    LeaveGlobalBuildInformation(*p.first, p.second);
                }
                while (p.second.CreateChild() < (1 << DimT) &&
                       !ExistChild(*p.first, p.second.CreateChild()))
                {
                    ++p.second.CreateChild();
                }
                if (p.second.CreateChild() == (1 << DimT))
                {
                    continue;
                }
                EnterGlobalBuildInformation(*p.first, &p.second);
                stack.resize(stack.size() + 1); // create new entry
                stack.back().first = &(*p.first)[p.second.CreateChild()];
                InitSerializeInformation(
                    *p.first, p.second, p.second.CreateChild(), &stack.back().second);
                do
                {
                    ++p.second.CreateChild();
                } while (p.second.CreateChild() < (1 << DimT) &&
                         !ExistChild(*p.first, p.second.CreateChild()));
            }
        }
 
        void
        SerializeBreadthFirst(std::ostream* out)
        {
            typedef std::pair<const CellType*, SerializeInformation> Pair;
            std::deque<Pair> queue(1);
            // init build information directly on stack to avoid
            // copying.
            queue.back().first = BaseType::Root();
            InitRootSerializeInformation(&queue.back().second);
            while (queue.size())
            {
                Pair& p = queue.front();
                if (p.second.CreateChild() == 1 << DimT)
                {
                    queue.pop_front();
                    continue;
                }
                if (!ShouldSubdivide(*p.first, p.second))
                {
                    queue.pop_front();
                    continue;
                }
                if (!p.second.Written())
                {
                    const size_t byteCount =
                        ((1 << BaseType::m_dim) / 8) + (((1 << BaseType::m_dim) % 8) ? 1 : 0);
                    char b[byteCount];
                    for (unsigned int i = 0; i < (1 << DimT); ++i)
                        if (ExistChild(*p.first, i))
                        {
                            b[i >> 3] |= 1 << (i - ((i >> 3) << 3));
                        }
                        else
                        {
                            b[i >> 3] &= ~(1 << (i - ((i >> 3) << 3)));
                        }
                    out->write(b, byteCount);
                    p.second.Written(true);
                }
                while (p.second.CreateChild() < (1 << DimT) &&
                       !ExistChild(*p.first, p.second.CreateChild()))
                {
                    ++p.second.CreateChild();
                }
                if (p.second.CreateChild() == (1 << DimT))
                {
                    continue;
                }
                queue.resize(queue.size() + 1); // create new entry
                queue.back().first = &(*p.first)[p.second.CreateChild()];
                InitSerializeInformation(
                    *p.first, p.second, p.second.CreateChild(), &queue.back().second);
                do
                {
                    ++p.second.CreateChild();
                } while (p.second.CreateChild() < (1 << DimT) &&
                         !ExistChild(*p.first, p.second.CreateChild()));
            }
        }
 
        void
        Serialize(std::istream* in)
        {
            VectorXD<DimT, ScalarType> min;
            min.Zero();
            AACube<VectorXD<DimT, ScalarType>> bcube(min, 1);
            Serialize(bcube, in);
        }
 
        void
        Serialize(const AACube<VectorXD<DimT, ScalarType>>& bcube, std::istream* in)
        {
            typedef std::pair<CellType*, BuildInformation> Pair;
            BaseType::Clear();
            BaseType::Root() = new CellType();
            std::deque<Pair> stack(1);
            // init build information directly on stack to avoid
            // copying.
            stack.back().first = BaseType::Root();
            InitRootBuildInformation(bcube, &stack.back().second);
            InitRoot(stack.back().second, BaseType::Root());
            while (stack.size())
            {
                Pair& p = stack.back();
                if (p.second.CreateChild() == 1 << DimT)
                {
                    LeaveGlobalBuildInformation(*p.first, p.second);
                    stack.pop_back();
                    continue;
                }
                if (IsLeaf(*p.first))
                {
                    if (!ShouldSubdivide(*p.first, p.second))
                    {
                        stack.pop_back();
                        continue;
                    }
                    // create the children
                    const size_t byteCount = ((1 << DimT) / 8) + (((1 << DimT) % 8) ? 1 : 0);
                    char b[byteCount];
                    in->read(b, byteCount);
                    for (size_t i = 0; i < (1 << DimT); ++i)
                    {
                        const size_t cmpB = 1 << (i - ((i >> 3) << 3));
                        if (b[i >> 3] & cmpB)
                        {
                            p.first->Children()[i] = new CellType();
                        }
                        else
                        {
                            p.first->Children()[i] = (CellType*)1;
                        }
                    }
                    Distribute(p.second, p.first);
                    if (IsLeaf(*p.first)) // couldn't subdivide?
                    {
                        stack.pop_back();
                        continue;
                    }
                }
                else
                {
                    LeaveGlobalBuildInformation(*p.first, p.second);
                }
                while (p.second.CreateChild() < (1 << DimT) &&
                       !ExistChild(*p.first, p.second.CreateChild()))
                {
                    ++p.second.CreateChild();
                }
                if (p.second.CreateChild() == (1 << DimT))
                {
                    continue;
                }
                EnterGlobalBuildInformation(*p.first, &p.second);
                stack.resize(stack.size() + 1); // create new entry
                stack.back().first = &(*p.first)[p.second.CreateChild()];
                InitBuildInformation(
                    *p.first, p.second, p.second.CreateChild(), &stack.back().second);
                InitCell(*p.first,
                         p.second,
                         p.second.CreateChild(),
                         stack.back().second,
                         &(*p.first)[p.second.CreateChild()]);
                do
                {
                    ++p.second.CreateChild();
                } while (p.second.CreateChild() < (1 << DimT) &&
                         !ExistChild(*p.first, p.second.CreateChild()));
            }
        }
 
        void
        SerializeBreadthFirst(std::istream* in)
        {
            VectorXD<DimT, ScalarType> min;
            min.Zero();
            AACube<VectorXD<DimT, ScalarType>> bcube(min, 1);
            SerializeBreadthFirst(bcube, in);
        }
 
        void
        SerializeBreadthFirst(const AACube<VectorXD<DimT, ScalarType>>& bcube, std::istream* in)
        {
            typedef std::pair<CellType*, BuildInformation> Pair;
            BaseType::Clear();
            BaseType::Root() = new CellType();
            std::deque<Pair> queue(1);
            // init build information directly on stack to avoid
            // copying.
            queue.back().first = BaseType::Root();
            InitRootBuildInformation(bcube, &queue.back().second);
            InitRoot(queue.back().second, BaseType::Root());
            while (queue.size())
            {
                Pair& p = queue.front();
                if (p.second.CreateChild() == 1 << DimT)
                {
                    queue.pop_front();
                    continue;
                }
                if (IsLeaf(*p.first))
                {
                    if (!ShouldSubdivide(*p.first, p.second))
                    {
                        queue.pop_front();
                        continue;
                    }
                    // create the children
                    const size_t byteCount = ((1 << DimT) / 8) + (((1 << DimT) % 8) ? 1 : 0);
                    char b[byteCount];
                    in->read(b, byteCount);
                    for (size_t i = 0; i < (1 << DimT); ++i)
                    {
                        const size_t cmpB = 1 << (i - ((i >> 3) << 3));
                        if (b[i >> 3] & cmpB)
                        {
                            p.first->Children()[i] = new CellType();
                        }
                        else
                        {
                            p.first->Children()[i] = (CellType*)1;
                        }
                    }
                    Distribute(p.second, p.first);
                    if (IsLeaf(*p.first)) // couldn't subdivide?
                    {
                        queue.pop_front();
                        continue;
                    }
                }
                while (p.second.CreateChild() < (1 << DimT) &&
                       !ExistChild(*p.first, p.second.CreateChild()))
                {
                    ++p.second.CreateChild();
                }
                if (p.second.CreateChild() == (1 << DimT))
                {
                    continue;
                }
                queue.resize(queue.size() + 1); // create new entry
                queue.back().first = &(*p.first)[p.second.CreateChild()];
                InitBuildInformation(
                    *p.first, p.second, p.second.CreateChild(), &queue.back().second);
                InitCell(*p.first,
                         p.second,
                         p.second.CreateChild(),
                         queue.back().second,
                         &(*p.first)[p.second.CreateChild()]);
                do
                {
                    ++p.second.CreateChild();
                } while (p.second.CreateChild() < (1 << DimT) &&
                         !ExistChild(*p.first, p.second.CreateChild()));
            }
        }
 
    protected:
        class BuildInformation : public BaseType::BuildInformation
        {
        public:
            unsigned int&
            CreateChild()
            {
                return m_createChild;
            }
 
            const unsigned int
            CreateChild() const
            {
                return m_createChild;
            }
 
            const AACube<VectorXD<DimT, ScalarType>>&
            Cube() const
            {
                return m_cube;
            }
 
            AACube<VectorXD<DimT, ScalarType>>&
            Cube()
            {
                return m_cube;
            }
 
        private:
            unsigned int m_createChild;
            AACube<VectorXD<DimT, ScalarType>> m_cube;
        };
 
        class SerializeInformation : public BaseType::BuildInformation
        {
        public:
            unsigned int&
            CreateChild()
            {
                return m_createChild;
            }
 
            const unsigned int
            CreateChild() const
            {
                return m_createChild;
            }
 
            bool
            Written() const
            {
                return m_written;
            }
 
            void
            Written(bool written)
            {
                m_written = written;
            }
 
        private:
            unsigned int m_createChild;
            bool m_written;
        };
 
        class AxisSplitter
        {
        public:
            AxisSplitter()
            {
            }
 
            AxisSplitter(unsigned int axis, ScalarType value) : m_axis(axis), m_value(value)
            {
            }
 
            unsigned int&
            Axis()
            {
                return m_axis;
            }
 
            ScalarType&
            Value()
            {
                return m_value;
            }
 
            template <class VectorT>
            bool
            operator()(const VectorT& v) const
            {
                return v[m_axis] <= m_value;
            }
 
        private:
            unsigned int m_axis;
            ScalarType m_value;
        };
 
        void
        InitRootBuildInformation(const AACube<VectorXD<DimT, ScalarType>>& bcube,
                                 BuildInformation* bi)
        {
            bi->CreateChild() = 0;
            bi->Cube() = bcube;
            BaseType::InitRootBuildInformation(bi);
        }
 
        void
        InitBuildInformation(const CellType& parent,
                             const BuildInformation& parentInfo,
                             unsigned int childIdx,
                             BuildInformation* bi)
        {
            bi->CreateChild() = 0;
            bi->Cube().Width(parentInfo.Cube().Width() / ScalarType(2));
            for (unsigned int i = 0; i < DimT; ++i)
                if (childIdx & (1 << (DimT - 1 - i)))
                    bi->Cube().Min()[i] = parentInfo.Cube().Min()[i] + bi->Cube().Width();
                else
                {
                    bi->Cube().Min()[i] = parentInfo.Cube().Min()[i];
                }
            BaseType::InitBuildInformation(parent, parentInfo, childIdx, bi);
        }
 
        void
        InitRootSerializeInformation(SerializeInformation* bi) const
        {
            bi->CreateChild() = 0;
            bi->Written(false);
            BaseType::InitRootBuildInformation(bi);
        }
 
        void
        InitSerializeInformation(const CellType& parent,
                                 const SerializeInformation& parentInfo,
                                 unsigned int childIdx,
                                 SerializeInformation* bi) const
        {
            bi->CreateChild() = 0;
            bi->Written(false);
            BaseType::InitBuildInformation(parent, parentInfo, childIdx, bi);
        }
 
        void
        Subdivide(BuildInformation& bi, CellType* cell)
        {
            // get the array of splitters
            VectorXD<DimT, ScalarType> center;
            bi.Cube().Center(&center);
            AxisSplitter splitters[DimT];
            for (unsigned int i = 0; i < DimT; ++i)
            {
                splitters[i].Axis() = i;
                splitters[i].Value() = center[i];
            }
            BaseType::SplitData(splitters, DimT, *cell, bi, cell->Children());
        }
 
        void
        Distribute(BuildInformation& bi, CellType* cell)
        {
            const size_t byteCount = ((1 << BaseType::m_dim) / 8) + (((1 << DimT) % 8) ? 1 : 0);
            char b[byteCount];
            // create missing cells
            for (unsigned int i = 0; i < (1 << DimT); ++i)
                if (!ExistChild(*cell, i))
                {
                    cell->Child(i, new CellType);
                    b[i >> 3] |= 1 << (i - ((i >> 3) << 3));
                }
                else
                {
                    b[i >> 3] &= ~(1 << (i - ((i >> 3) << 3)));
                }
            // get the array of splitters
            VectorXD<DimT, ScalarType> center;
            bi.Cube().Center(&center);
            AxisSplitter splitters[DimT];
            for (unsigned int i = 0; i < DimT; ++i)
            {
                splitters[i].Axis() = i;
                splitters[i].Value() = center[i];
            }
            SplitData(splitters, DimT, bi, cell->Children());
            for (unsigned int i = 0; i < (1 << DimT); ++i)
            {
                if (cell->Children()[i]->Size() && (b[i >> 3] & (1 << (i - ((i >> 3) << 3)))))
                {
                    std::cout << "Bug in distribute!" << std::endl;
                }
                if (!cell->Children()[i]->Size() && (b[i >> 3] & (1 << (i - ((i >> 3) << 3)))))
                {
                    delete cell->Children()[i];
                    cell->Children()[i] = (CellType*)1;
                }
            }
        }
 
        template <class PointT, class TraversalInformationT>
        const CellType*
        NodeContainingPoint(const PointT& point,
                            size_t maxLevel,
                            size_t minSize,
                            const CellType& cell,
                            const TraversalInformationT& ti,
                            CellRange* range) const
        {
            if (this->IsLeaf(cell))
            {
                this->GetCellRange(cell, ti, range);
                return &cell;
            }
            if (this->CellLevel(cell, ti) == maxLevel)
            {
                this->GetCellRange(cell, ti, range);
                return &cell;
            }
            // find the child containing the point
            //typename BaseType::CellCenterType center;
            //CellCenter(cell, ti, &center);
            size_t childIdx = 0;
            for (unsigned int i = 0; i < DimT; ++i)
            {
                if (point[i] > cell.Center()[i]) //center[i])
                {
                    childIdx |= 1 << (DimT - i - 1);
                }
            }
            if (this->ExistChild(cell, childIdx) && cell[childIdx].Size() >= minSize)
            {
                TraversalInformationT cti;
                BaseType::InitTraversalInformation(cell, ti, childIdx, &cti);
                return NodeContainingPoint(point, maxLevel, minSize, cell[childIdx], cti, range);
            }
            this->GetCellRange(cell, ti, range);
            return &cell;
        }
    };
}; // namespace GfxTL
 
#endif