Octree.hpp

Go to the documentation of this file.

//  This file is part of ff3d - http://www.freefem.org/ff3d
//  Copyright (C) 2001, 2002, 2003 Stéphane Del Pino

//  This program is free software; you can redistribute it and/or modify
//  it under the terms of the GNU General Public License as published by
//  the Free Software Foundation; either version 2, or (at your option)
//  any later version.

//  This program is distributed in the hope that it will be useful,
//  but WITHOUT ANY WARRANTY; without even the implied warranty of
//  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
//  GNU General Public License for more details.

//  You should have received a copy of the GNU General Public License
//  along with this program; if not, write to the Free Software Foundation,
//  Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.  

//  $Id: Octree.hpp,v 1.10 2007/04/22 21:46:44 delpinux Exp $

#ifndef OCTREE_HPP
#define OCTREE_HPP

#include <TinyVector.hpp>
#include <StaticPow.hpp>

#include <ErrorHandler.hpp>

#include <stack>

template <typename Content,
          size_t GivenDimension>
class Octree
{
public:
  class iterator;

  enum {
    Dimension = GivenDimension
  };

private:
  class __SubTree;

  class __Node
  {
  public:
    enum Type {
      subtree,
      leaf
    };

    virtual const typename __Node::Type type() const = 0;

    __Node()
    {
      ;
    }
    
    __Node(const __Node& n)
    {
      ;
    }

    virtual ~__Node()
    {
      ;
    }
  };

  class __Leaf
    : public __Node
  {
  private:
    const TinyVector<Dimension, real_t> __position;
    Content __value;
  public:
    const typename __Node::Type type() const
    {
      return __Node::leaf;
    }

    void setValue(const Content& value)
    {
      __value = value;
    }

    const TinyVector<Dimension, real_t>& position() const
    {
      return __position;
    }

    Content& value()
    {
      return __value;
    }

    const Content& value() const
    {
      return __value;
    }

    __Leaf(const TinyVector<Dimension, real_t> position)
      :__position(position)
    {
      ;
    }

    __Leaf(const __Leaf& l)
      : __Node(l),
        __position(l.__position),
        __value(l.__value)
    {
      ;
    }

    ~__Leaf()
    {
      ;
    }
  };


  class __SubTree
    : public __Node
  {
  private:
    TinyVector<StaticPow<2,Dimension>::value,
               __Node*> __nodes;

  public:
    __Node* node(const size_t i)
    {
      return __nodes[i];
    }

    const typename __Node::Type type() const
    {
      return __Node::subtree;
    }

    typename Octree<Content, Dimension>::iterator&
    fuzzySearch(const TinyVector<Dimension, real_t>& X,
                const TinyVector<Dimension, real_t>& center,
                const real_t size,
                typename Octree<Content, Dimension>::iterator& cell)
    {
      int num=0;
      TinyVector<Dimension,int> I;

      for (size_t i=0; i<Dimension; ++i) {
        I[i] = (X[i]>center[i]);
        num += (I[i] << (Dimension-1-i));
      }

      cell.addAddress(this);
      cell.localNumber() = num;

      if (__nodes[num]==0) {
        cell.localNumber() = -1;        
        cell++;
      } else {
        double newS = 0.5*size;

        TinyVector<Dimension> newCenter = center;
        for (size_t i=0; i<Dimension; ++i) {
          newCenter[i] += (I[i]*2-1)*newS;
        }

        if ((*__nodes[num]).type() ==  __Node::subtree) {
          (static_cast<__SubTree&>(*__nodes[num])).fuzzySearch(X,newCenter,newS,cell);
        } else {
          cell.__leaf = static_cast<__Leaf*>(__nodes[num]);
        }
      }

      return cell;
    }

    void addLeaf(const TinyVector<Dimension, real_t>& center,
                 const real_t s,
                 __Leaf* l)
    {
      int num=0;
      TinyVector<Dimension,int> I;

      for (size_t i=0; i<Dimension; ++i) {
        I[i] = ((*l).position()[i]>center[i]);

        num += (I[i] << (Dimension-1-i));
      }

      if (__nodes[num]==0) {
        __nodes[num] = l;
      } else {
        double newS = 0.5*s;

        TinyVector<Dimension> newCenter = center;
        for (size_t i=0; i<Dimension; ++i) {
          newCenter[i] += (I[i]*2-1)*newS;
        }

        if ((*__nodes[num]).type() ==  __Node::subtree) {
          (static_cast<__SubTree&>(*__nodes[num])).addLeaf(newCenter,newS,l);
        } else {
          __Leaf* oldleaf=static_cast<__Leaf*>(__nodes[num]);
          __SubTree* s =new __SubTree();
          __nodes[num]=s;
          s->addLeaf(newCenter, newS, l);
          s->addLeaf(newCenter, newS, oldleaf);
        }
      }
    }

    iterator begin()
    {
      iterator i;
      i.addAddress(this);
      i++;
      return i;
    }

    __SubTree()
      : __nodes(0)
    {
      ;
    }

    __SubTree(const __SubTree& s)
      : __Node(s),
        __nodes(s)
    {
      ;
    }

    ~__SubTree()
    {
      ;
    }
  };

private:
  __Node* __node;               
  TinyVector<Dimension> __center;       
  real_t __size;                
public:
  class iterator
  {
  public:
    typedef std::pair <__SubTree*, int> LocalAddress;

    friend class Octree<Content, Dimension>;
    friend class __SubTree;

  private:
    __Leaf* __leaf;

    std::stack<LocalAddress> __address;

  public:

    size_t depth()
    {
      return __address.size();
    }

    __SubTree* currentFather()
    {
      if (__address.size() > 0) {
        return __address.top().first;
      } else {
        return 0;
      }
    }

    iterator& operator++(int)
    {
      __SubTree* father = currentFather();
      __leaf = 0;
      if (father != 0) {
        __address.top().second++;

        if (__address.top().second < 8) {
          __Node* n = father->node(__address.top().second);
          if (n != 0) {
            if(n->type() == __Node::leaf) {
              __leaf = static_cast<__Leaf*>(n);
            } else {
              __SubTree* t = static_cast<__SubTree*>(n);
              this->addAddress(t);
              (*this)++;
            }
          } else {
            (*this)++;
          }
        } else {
          __address.pop();
          (*this)++;
        }
      }
      return *this;
    }

    const int& localNumber() const
    {
      return __address.top().second;
    }

    int& localNumber()
    {
      return __address.top().second;
    }

    void addAddress(__SubTree* a)
    {
      LocalAddress l;
      l.first = a;
      l.second = -1;
      __address.push(l);
    }

    operator __Leaf*()
    {
      return __leaf;
    }

    iterator()
      : __leaf(0)
    {
      ;
    }

    iterator(const iterator& i)
      : __leaf(i.__leaf),
        __address(i.__address)
    {
      ;
    }

    ~iterator()
    {
      ;
    }
  };

  iterator begin() const
  {
    iterator i;
    if (__node != 0) {
      if ((*__node).type() == __Node::leaf) {
      } else {
        __SubTree& s = static_cast<__SubTree&>(*__node);
        i = s.begin();
      }
    }
    return i;
  }

  iterator end() const
  {
    return iterator();
  }


  iterator search(const TinyVector<Dimension, real_t>& X)
  {
    throw ErrorHandler(__FILE__,__LINE__,
                       "not implemented",
                       ErrorHandler::unexpected);
  }

  iterator fuzzySearch(const TinyVector<Dimension, real_t>& X) const
  {
    iterator i;
    if(__node == 0) {
      return this->end();
    } else if((*__node).type() != __Node::subtree) {
      i.__leaf = static_cast<__Leaf*>(__node);
      return i;
    } else {
      return static_cast<__SubTree&>(*__node).fuzzySearch(X, __center, __size, i);
    }

    return i;
  }

  void  setSize(const TinyVector<Dimension>& a, const TinyVector<Dimension>& b)
  {
    __center = a+b;
    __center *= 0.5;

    __size = std::abs(a[0]-__center[0]);
    for (size_t i=1; i<Dimension; ++i) {
      __size = std::max(__size,std::abs(a[i]-__center[i]));
    }
  }

  void add(const Content& c, const TinyVector<Dimension, real_t>& position)
  {
    if (__node == 0) {
      __node = new __Leaf(position);
      dynamic_cast<__Leaf&>(*__node).setValue(c);
    } else {
      if ((*__node).type() == __Node::leaf) {
        __Leaf* l = dynamic_cast<__Leaf*>(__node);
        if(__size == 0) {
          setSize(l->position(), position);
        }
        __SubTree* s = new __SubTree();
        __node = s;

        s->addLeaf(__center,__size,l);
        __Leaf* newLeaf = new __Leaf(position);
        newLeaf->setValue(c);
        s->addLeaf(__center,__size,newLeaf);

      } else {
        __Leaf* newLeaf = new __Leaf(position);
        newLeaf->setValue(c);
        (dynamic_cast<__SubTree&>(*__node)).addLeaf(__center,__size,newLeaf);
      }
    }
  }

  Octree()
    : __node(0),
      __center(0),
      __size(0)
  {
    ;
  }

  Octree(const TinyVector<Dimension, real_t>& a, const TinyVector<Dimension, real_t>& b)
    : __node(0)
  {
    this->setSize(a,b);
  }

  ~Octree()
  {
    ;
  }

};

#endif // OCTREE_HPP

---