Quadtree.hpp

Go to the documentation of this file.

/* $Id$ */

#ifndef LGL_LGLQUADTREE_HPP
#define LGL_LGLQUADTREE_HPP

#include <iostream>

#include <list>
#include <vector>

#include <lgl/Decomp.hpp>
#include <lgl/Decomposer.hpp>
#include <lgl/Raster.hpp>

using namespace std;

namespace jafar { //namespace jafar
  namespace lgl { //namespace lgl

  // some declaration to speed up the compilation
  template<typename T, class RasterT> class Quadtree;
  template<typename T, class RasterT> class QuadNode;
  template<typename T, class RasterT> class QuadNodeLeaf;
  template<typename T, class RasterT> class QuadNodeInternal;

  template <
    typename T,
    class RasterT
    >
  class Quadtree : public TDecomposer<T,RasterT> { /*{{{ class Quadtree*/

  public:

    enum {  QUADTREE_TOP_LEFT_CELL=0,   QUADTREE_TOP_RIGHT_CELL,
            QUADTREE_BOTTOM_LEFT_CELL,  QUADTREE_BOTTOM_RIGHT_CELL };

    enum  UNWISE_SIZE_DIVISION_POLITIC {LEFTTOP, BOTTOMRIGHT, BALANCE};

    // Default dull constructor
    //Quadtree();

    Quadtree(RasterT* _raster, UNWISE_SIZE_DIVISION_POLITIC _uspolitic=LEFTTOP);

    ~Quadtree();

      QuadNodeInternal<T,RasterT>* getRootCell() const;

    bool getClusterForCell(const RasterCellIndex& cell, RasterRect& cluster) const;

    bool getClusterForCell(const RasterCellIndex& cell, RasterCellIndex& clusterRoot) const;

    bool getAdjacentClusters(const RasterCellIndex& root_cell, Decomp::SIDE side,std::list<RasterCellIndex>& clusters) const;

    T getData(const RasterRect& rect) const;

    T getData(const RasterCellIndex& cluster) const;


    std::list<QuadNodeLeaf<T,RasterT>*> getFullLeavesList() const;

   private:

    QuadNodeInternal<T,RasterT>* root;

    UNWISE_SIZE_DIVISION_POLITIC uspolitic;

    friend class QuadNode<T,RasterT>;
    friend class QuadNodeInternal<T,RasterT>;
    friend class QuadNodeLeaf<T,RasterT>;

    QuadNodeLeaf<T,RasterT>* getLeafForCell(const RasterCellIndex& cell) const;
    QuadNodeLeaf<T,RasterT>* getLeafForRect(const RasterRect& rect) const;

    void setNewLeaf(QuadNodeLeaf<T,RasterT>* leaf, bool multiCluster=false);


    std::list<QuadNodeLeaf<T,RasterT>*> allLeaves;

   }; /*}}} end of class Quadtree*/


   template <typename T, class RasterT>
   class QuadNode { /*{{{ class QuadNode*/

     public:

      QuadNode(Quadtree<T,RasterT>* _tree, QuadNodeInternal<T,RasterT>* _parent, RasterRect _r);

      virtual ~QuadNode();

      enum  NODE_TYPE {INTERNAL_NODE, LEAF_NODE, ROOT_NODE};

      bool contains(const RasterCellIndex& cell) const;
      bool contains(const RasterRect& _rect) const;

      QuadNodeInternal<T,RasterT>* getParent() const;

      Quadtree<T,RasterT>* getTree() const;

      void setParent(QuadNodeInternal<T,RasterT>* _parent);

      // abstract methods
      /*
      virtual int getxroot() const = 0;
      virtual int getyroot() const = 0;
      virtual int getxsize() const = 0;
      virtual int getysize() const = 0;

      virtual RasterRect getRect()const=0;
      */
      int getxroot() const { return rect.xtopleft; };
      int getyroot() const { return rect.ytopleft; };
      int getxsize() const { return rect.xsize(); }
      int getysize() const { return rect.ysize(); }
      RasterRect  getRect()   const { return rect; }


      virtual QuadNodeLeaf<T,RasterT>* getLeafForCell(const RasterCellIndex& cell) =0;

      virtual NODE_TYPE getType() const=0;

    protected:

      Quadtree<T,RasterT>* tree;

      QuadNodeInternal<T,RasterT>* parent;

      RasterRect rect;

  }; /* end of class QuadNode }}}*/

   template <typename T, class RasterT>
   class QuadNodeInternal : public QuadNode<T,RasterT> { /*{{{ class QuadNodeInternal*/

     private:

      bool divide( const RasterRect& _rect, bool xbalance=true, bool ybalance=true);

      //RasterRect rect;
     public:

      QuadNodeInternal(Quadtree<T,RasterT>* _tree, QuadNodeInternal<T,RasterT>* _parent, RasterRect rect, bool xbalance=true, bool ybalance=true);

      ~QuadNodeInternal();

      /*
      int getxroot() const;
      int getyroot() const;
      int getxsize() const;
      int getysize() const;
      RasterRect getRect() const;
      */


      typename QuadNode<T,RasterT>::NODE_TYPE getType() const;

      // get the cluster (tree leaf) for cell
      QuadNodeLeaf<T,RasterT>* getLeafForCell(const RasterCellIndex& cell);

      QuadNode<T,RasterT>* getChild(int index);

      // get the position (topleft,topright,bottomleft,bottoright) of a child
      int getChildPos(QuadNode<T,RasterT>* child);

      bool setChildren(QuadNode<T,RasterT>* S1, QuadNode<T,RasterT>* S2, QuadNode<T,RasterT>* S3, QuadNode<T,RasterT>* S4);

      bool setChild(QuadNode<T,RasterT>* child,int index);

     protected:

      //std::vector<QuadNode<T,RasterT>*> children;
      QuadNode<T,RasterT>* children[4];

      template<typename U, class URasterT>
      friend std::ostream& operator<< (std::ostream& out, const QuadNodeInternal<U,URasterT>& node);

   }; /* end of class QuadNodeInternal }}}*/
   
   template <typename T, class RasterT>
   class QuadNodeLeaf: public QuadNode<T,RasterT> { /*{{{ class QuadNodeLeaf */
     private:

      int nodePos;

      T data;

      //RasterRect rect;

      // adjacent clusters to this cluster represented by the quadtree leaf
      //std::vector<std::vector<QuadNodeLeaf*>> neighboors;

      //template<typename U, class URasterT>
      //friend std::ostream& operator<< (std::ostream& out, const QuadNodeLeaf<U,URasterT>& node);

    public:

      QuadNodeLeaf(Quadtree<T,RasterT>* _tree, QuadNodeInternal<T,RasterT>* _parent, int _nodePos, RasterRect _rect, T _data=RasterT::DEFAULT_DATA_VALUE());

      ~QuadNodeLeaf();

      /*
      int getxroot() const;
      int getyroot() const;
      int getxsize() const;
      int getysize() const;
      RasterRect getRect() const;
      */

      typename QuadNode<T,RasterT>::NODE_TYPE getType() const;


      // if the leaf contains only one cluster
      bool hasManyClusters() const;

      bool isAtomic() const;

      QuadNodeLeaf<T,RasterT>* getLeafForCell(const RasterCellIndex& cell);

      bool getClusterForCell(const RasterCellIndex& cell, RasterRect& cluster);

      // return the data of the cluster represented by the leaf
      T getData(const RasterCellIndex& clusterRoot);

  }; /* class QuadNodeLeaf }}}*/

   /**************************************************
    *                 IMPLEMENTATION                 *
    **************************************************/

  /*--------------------- Quadtree ---------------------- {{{*/
  template <typename T, class RasterT>
   Quadtree<T,RasterT>::Quadtree(RasterT* _raster, UNWISE_SIZE_DIVISION_POLITIC _uspolitic):
    TDecomposer<T,RasterT>(_raster),
    root(NULL),
    uspolitic(_uspolitic)
   { //{{{

    RasterRect rect = this->raster->getRasterRect();
    // check if the raster is homogene, if not decompose in quadtree
    if ((rect.xsize()>0)&&(rect.ysize()>0)) {
      if (this->raster->isHomogeneous(rect)==false) {
        root = new QuadNodeInternal<T,RasterT>(this,NULL,rect,true,true);
      } else {
        // We don't set any root cell
        root = NULL;
        this->decomp->setCluster(rect);
      }
    }
   } //}}}

  template <typename T, class RasterT>
  Quadtree<T,RasterT>::~Quadtree() {
    // Clear the leaves container
    allLeaves.clear();
    // Recursive deletion of all nodes
    delete root;
    // Clearing the clusters in the decomposition structure
    delete this->decomp;
  }

  template <typename T, class RasterT>
  QuadNodeInternal<T,RasterT>* Quadtree<T,RasterT>::getRootCell() const {
    return root;
  }

   template <typename T, class RasterT>
   std::list<QuadNodeLeaf<T,RasterT>*> Quadtree<T,RasterT>::getFullLeavesList() const {
    return allLeaves;
   }

   template <typename T, class RasterT>
   bool Quadtree<T,RasterT>::getClusterForCell(const RasterCellIndex& cell, RasterRect& cluster) const {
    QuadNodeLeaf<T,RasterT>* leaf = getLeafForCell(cell);

    if (leaf==NULL)
    {
      if (root == NULL)
      {
        // No root in tree means there were no need to decompose.
        // So directly returning the full Raster's rectangle.
        cluster = this->raster->getRasterRect();
        return true;
      } 
      // Else if root exist and no leaves are found there is a prob.
      return false;
    }

    // if leaf is not a cluster but a set of (1 or more) cluster of size (1,1)
    if (leaf->hasManyClusters())
      cluster.setRect(cell.i,cell.j,1,1);
    else
      cluster=leaf->getRect();

    return true;
   }

   template <typename T, class RasterT>
   bool Quadtree<T,RasterT>::getClusterForCell(const RasterCellIndex& cell, RasterCellIndex& clusterRoot) const {
    QuadNodeLeaf<T,RasterT>* leaf = getLeafForCell(cell);

    if (leaf==NULL)
      return false;

    // if leaf is not a cluster but a set of (1 or more) cluster of size (1,1)
    if (leaf->hasManyClusters())
      clusterRoot=cell;
    else
      clusterRoot.set(leaf->getRect().xroot(),leaf->getRect().yroot());
    return true;
   }

   template <typename T, class RasterT>
   bool Quadtree<T,RasterT>::getAdjacentClusters(const RasterCellIndex& root_cell, Decomp::SIDE side, std::list<RasterCellIndex>& clusters) const
   { //{{{

    DecompCluster* curCell;
    RasterCellIndex curCluster;
    int i,j;
    int clusterXSize = this->decomp->at(root_cell.i,root_cell.j)->xsize;
    int clusterYSize = this->decomp->at(root_cell.i,root_cell.j)->ysize;
    int curClusterXRoot, curClusterYRoot;
    bool curCellIsRootCluster;

    switch(side) {
    case Decomp::CLUSTER_TOP_SIDE:
      {
      int topSide_i=root_cell.i-1;
      if (topSide_i < this->raster->getxroot()) {
      break;
      }
      j=root_cell.j;
      while (j<root_cell.j+clusterYSize) {
        curCell=this->decomp->at(topSide_i,j);
        if (curCell != NULL) {
          curCellIsRootCluster=!curCell->isNegative();
          //curClusterXRoot=curCellIsRootCluster?topSide_i:(topSide_i+curCell->xsize+1);
          //curClusterYRoot=curCellIsRootCluster?j:(j+curCell->ysize+1);
          if (curCellIsRootCluster)
          {
            curClusterXRoot = topSide_i;
            curClusterYRoot = j;
            j += curCell->ysize;
          } else {
            curClusterXRoot = (topSide_i+curCell->xsize+1);
            curClusterYRoot = (j+curCell->ysize+1);
          }
          clusters.push_back(RasterCellIndex(curClusterXRoot,curClusterYRoot));
        }
        j++;
      }
      if (clusters.empty()) {
        getClusterForCell(RasterCellIndex(topSide_i,root_cell.j),curCluster);
        clusters.push_back(curCluster);
      }
      break;
      }
    case Decomp::CLUSTER_BOTTOM_SIDE:
      {
       int bottomSide_i=root_cell.i+clusterXSize;
       if (bottomSide_i<this->raster->getxroot()) {
        break;
       }
       j=root_cell.j;
       while (j<root_cell.j+clusterYSize) {
        curCell=this->decomp->at(bottomSide_i,j);
        if (curCell!=NULL) {
          curCellIsRootCluster=!curCell->isNegative();
          //curClusterXRoot=curCellIsRootCluster?bottomSide_i:(bottomSide_i+curCell->xsize+1);
          //curClusterYRoot=curCellIsRootCluster?j:(j+curCell->ysize+1);
          if (curCellIsRootCluster)
          {
            curClusterXRoot = bottomSide_i;
            curClusterYRoot = j;
           j += curCell->ysize;
          } else {
            curClusterXRoot = (bottomSide_i+curCell->xsize+1);
            curClusterYRoot = (j+curCell->ysize+1);
          }

          clusters.push_back(RasterCellIndex(curClusterXRoot,curClusterYRoot));
        }
        j++;
       }
       if (clusters.empty()) {
        getClusterForCell(RasterCellIndex(bottomSide_i,root_cell.j),curCluster);
        clusters.push_back(curCluster);
       }
       break;
      }
    case Decomp::CLUSTER_LEFT_SIDE:
    {
      int leftSide_j=root_cell.j-1;
      if (leftSide_j<this->raster->getyroot()) {
        break;
      }
      i=root_cell.i;
      while (i<root_cell.i+clusterXSize) {
        curCell=this->decomp->at(i,leftSide_j);
        if (curCell!=NULL) {
          curCellIsRootCluster=!curCell->isNegative();
          //curClusterXRoot=curCellIsRootCluster?i:(i+curCell->xsize+1);
          //curClusterYRoot=curCellIsRootCluster?leftSide_j:(leftSide_j+curCell->ysize+1);
          if (curCellIsRootCluster)
          {
            curClusterXRoot = i;
            curClusterYRoot = leftSide_j;
            i += curCell->xsize;
          } else {
            curClusterXRoot = (i+curCell->xsize+1);
            curClusterYRoot = (leftSide_j+curCell->ysize+1);
          }

          clusters.push_back(RasterCellIndex(curClusterXRoot,curClusterYRoot));
        }
        i++;
      }
      if (clusters.empty()) {
        getClusterForCell(RasterCellIndex(root_cell.i,leftSide_j),curCluster);
        clusters.push_back(curCluster);
      }
      break;
      }
    case Decomp::CLUSTER_RIGHT_SIDE:
    {
      int rightSide_j=root_cell.j+clusterYSize;
      if (rightSide_j<this->raster->getyroot()) {
        break;
      }
      i=root_cell.i;
      while (i<root_cell.i+clusterXSize) {
        curCell=this->decomp->at(i,rightSide_j);
        if (curCell!=NULL) {
          curCellIsRootCluster=!curCell->isNegative();
          //curClusterXRoot=curCellIsRootCluster?i:(i+curCell->xsize+1);
          //curClusterYRoot=curCellIsRootCluster?rightSide_j:(rightSide_j+curCell->ysize+1);
          if (curCellIsRootCluster)
          {
            curClusterXRoot = i;
            curClusterYRoot = rightSide_j;
            i += curCell->xsize;
          } else {
            curClusterXRoot = (i+curCell->xsize+1);
            curClusterYRoot = (rightSide_j+curCell->ysize+1);
          }

          clusters.push_back(RasterCellIndex(curClusterXRoot,curClusterYRoot));
        }
        i++;
      }
      if (clusters.empty()) {
        getClusterForCell(RasterCellIndex(root_cell.i,rightSide_j),curCluster);
        clusters.push_back(curCluster);
      }
      break;
    }
    default:
      return false;
    }
    return true;
   } //}}}

   template <typename T, class RasterT>
   T Quadtree<T,RasterT>::getData(const RasterRect& rect) const {
    // look for the leaf
    QuadNodeLeaf<T,RasterT>* leaf = getLeafForRect(rect);
    if (leaf==NULL)
      return RasterT::DEFAULT_DATA_VALUE();
    else
      return leaf->getData(RasterCellIndex(rect.xtopleft,rect.ytopleft));
   }

   template <typename T, class RasterT>
   T Quadtree<T,RasterT>::getData(const RasterCellIndex& cluster) const {
    // looking for the leaf that contains this cluster
    QuadNodeLeaf<T,RasterT>* leaf = getLeafForCell(cluster);
    if (leaf==NULL)
      return RasterT::DEFAULT_DATA_VALUE();
    else
      return leaf->getData(cluster);
   }

   template <typename T, class RasterT>
   QuadNodeLeaf<T,RasterT>* Quadtree<T,RasterT>::getLeafForCell(const RasterCellIndex& cell) const {
    if (root == NULL)
    { return NULL; }

    QuadNodeLeaf<T,RasterT>* cluster =  root->getLeafForCell(cell);
    if (cluster == NULL)
    {
      std::cout << "[Quadtree] cannot get cluster for cell " << cell.i << "," << cell.j << std::endl;
    }

    return cluster;
   }

   template <typename T, class RasterT>
   QuadNodeLeaf<T,RasterT>* Quadtree<T,RasterT>::getLeafForRect(const RasterRect& rect) const {
    if (root==NULL)
    { return NULL; }

    QuadNodeLeaf<T,RasterT>* cluster =  root->getLeafForCell(RasterCellIndex(rect.xtopleft,rect.ytopleft));
    if (cluster==NULL) {
      std::cout << "[Quadtree] cannot get cluster for rect " << rect << std::endl;
    }

    if (cluster->getRect().contains(rect)) {
      return cluster;
    } else {
      return NULL;
    }
   }

  template <typename T, class RasterT>
  void Quadtree<T,RasterT>::setNewLeaf(QuadNodeLeaf<T,RasterT>* leaf, bool multiCluster) {
    RasterRect leafRect = leaf->getRect();
    if (multiCluster) {
      // Case where the leaf is not a square and quadtree division can not be homogeneous anymore.
      for (int i = leafRect.xroot(); i <= leafRect.xend(); i++) {
        for (int j = leafRect.yroot();j <= leafRect.yend(); j++) {
          this->decomp->setLeafCluster(RasterRect(i,j,1,1));
        }
      }
    } else {
      // update the decomp structure in the nominal case
      this->decomp->setLeafCluster(leafRect);
    }

    allLeaves.push_front(leaf);
  }

  /*---}}}*/
  /*--------------------- QuadNode ---------------------- {{{*/

   template <typename T, class RasterT>
   QuadNode<T,RasterT>::QuadNode(Quadtree<T,RasterT>* _tree, QuadNodeInternal<T,RasterT>* _parent, RasterRect _rect):
    tree(_tree),
    parent(_parent),
    rect(_rect)
  {
  }
       
  template <typename T, class RasterT>
  QuadNode<T,RasterT>::~QuadNode() {}

  template <typename T, class RasterT>
  bool QuadNode<T,RasterT>::contains(const RasterCellIndex& cell) const {
    return rect.contains(cell);
  }

  template <typename T, class RasterT>
  bool QuadNode<T,RasterT>::contains(const RasterRect& _rect) const {
    return rect.contains(_rect);
  }

  template <typename T, class RasterT>
  QuadNodeInternal<T,RasterT>* QuadNode<T,RasterT>::getParent() const
  {
    return parent;
  }

  template <typename T, class RasterT>
  Quadtree<T,RasterT>* QuadNode<T,RasterT>::getTree() const {
    return tree;
  }

  template <typename T, class RasterT>
  void QuadNode<T,RasterT>::setParent(QuadNodeInternal<T,RasterT>* _parent) {
    parent=_parent;
  }

  /*---}}}*/
  /*--------------------- QuadNodeInternal ---------------------- {{{*/

   template <typename T, class RasterT>
   QuadNodeInternal<T,RasterT>::QuadNodeInternal(Quadtree<T,RasterT>* _tree, QuadNodeInternal<T,RasterT>* _parent, RasterRect _rect, bool xbalance, bool ybalance):
    QuadNode<T,RasterT>(_tree,_parent,_rect)
   { //{{{
    // resize the vector of children nodes
    //children.resize(4, NULL);
    children[0] = NULL;
    children[1] = NULL;
    children[2] = NULL;
    children[3] = NULL;

    // divide this internal node
    if (divide(this->rect, xbalance, ybalance)==false) {
      std::cout << "#EEE# [Quadtree] Error while dividing the node" << std::endl;
      return;
    }

   } //}}}

   template <typename T, class RasterT>
   QuadNodeInternal<T,RasterT>::~QuadNodeInternal() {
     for (int i =  Quadtree<T,RasterT>::QUADTREE_TOP_LEFT_CELL;
              i <= Quadtree<T,RasterT>::QUADTREE_BOTTOM_RIGHT_CELL;
          i++)
     {
       delete (children[i]);
     }
   }

  template <typename T, class RasterT>
  bool QuadNodeInternal<T,RasterT>::divide( const RasterRect& _rect, bool xbalance, bool ybalance) {//{{{

    // std::cout << "[Quadtree] Devide " << _rect << std::endl;

    int _xroot = _rect.xtopleft;
    int _yroot = _rect.ytopleft;
    int _xsize = _rect.xsize();
    int _ysize = _rect.ysize();

    if ((_xsize<=1)||(_ysize<=1)) {
      return false;
    }

    // divide to 4 children
    // x size is multiple of 2
    int topleft_xsize, topleft_ysize;
    if (_xsize%2==0) {
      topleft_xsize=_xsize/2;
    } else {
      if (((this->tree)->uspolitic ==  Quadtree<T,RasterT>::LEFTTOP)||(((this->tree)->uspolitic ==  Quadtree<T,RasterT>::BALANCE)&&(xbalance==true))) {
        topleft_xsize = (_xsize+1)/2;
      } else {
        topleft_xsize = (_xsize-1)/2;
      }
    }
    if (_ysize%2==0) {
      topleft_ysize=_ysize/2;
    } else {
      if ((this->tree->uspolitic ==  Quadtree<T,RasterT>::LEFTTOP)||((this->tree->uspolitic ==  Quadtree<T,RasterT>::BALANCE)&&(ybalance==true))) {
       topleft_ysize = (_ysize+1)/2;
      } else {
       topleft_ysize = (_ysize-1)/2;
      }
    }

    // create 4 children nodes
    // 0 - top left
    RasterRect rectTopLeft(_xroot, _yroot, topleft_xsize, topleft_ysize);
    if ((topleft_xsize<=1)||(topleft_ysize<=1)||(this->tree->raster->isHomogeneous(rectTopLeft))) {
      // create a leaf
       // std::cout << "[Quadtree] Create leaf " << rectTopLeft << std::endl;
      QuadNodeLeaf<T,RasterT>* topleft = new QuadNodeLeaf<T,RasterT>(this->tree,this,Quadtree<T,RasterT>::QUADTREE_TOP_LEFT_CELL,rectTopLeft);
       // std::cout << "[Quadtree] Create leaf created " << rectTopLeft << std::endl;
      setChild(topleft,Quadtree<T,RasterT>::QUADTREE_TOP_LEFT_CELL);
      // std::cout << "[Quadtree] Create leaf DONE" << rectTopLeft << std::endl;
    } else {
      // else create another QuadNodeInternal to continue the quadtree decomposition
      QuadNodeInternal<T,RasterT>* topleft = new QuadNodeInternal<T,RasterT>(this->tree,this,rectTopLeft,!xbalance,!ybalance);
      setChild(topleft,Quadtree<T,RasterT>::QUADTREE_TOP_LEFT_CELL);
    }

    // 1 - top right
    //RasterRect rectTopRight(_xroot,_yroot+topleft_ysize,topleft_xsize,_ysize-topleft_ysize);
    RasterRect rectTopRight(_xroot+topleft_xsize,_yroot, _xsize-topleft_xsize, topleft_ysize);

    if ((rectTopRight.xsize()<=1)||(rectTopRight.ysize()<=1)||(this->tree->raster->isHomogeneous(rectTopRight))) {
      QuadNodeLeaf<T,RasterT>* topright = new QuadNodeLeaf<T,RasterT>(this->tree,this,Quadtree<T,RasterT>::QUADTREE_TOP_RIGHT_CELL,rectTopRight);
      setChild(topright,Quadtree<T,RasterT>::QUADTREE_TOP_RIGHT_CELL);
    } else {
      QuadNodeInternal<T,RasterT>* topright = new QuadNodeInternal<T,RasterT>(this->tree,this,rectTopRight,!xbalance,!ybalance);
      setChild(topright,Quadtree<T,RasterT>::QUADTREE_TOP_RIGHT_CELL);
    }

    // 2 - bottom left
    //RasterRect rectBottomLeft(_xroot+topleft_xsize,_yroot,_xsize-topleft:_xsize,topleft_ysize);
    RasterRect rectBottomLeft(_xroot,_yroot+topleft_ysize, topleft_xsize, _ysize-topleft_ysize);

    if (( rectBottomLeft.xsize() <=1)||(rectBottomLeft.ysize() <= 1)||(this->tree->raster->isHomogeneous(rectBottomLeft))) {
      QuadNodeLeaf<T,RasterT>*  bottomleft = new QuadNodeLeaf<T,RasterT>(this->tree,this,Quadtree<T,RasterT>::QUADTREE_BOTTOM_LEFT_CELL,rectBottomLeft);
      setChild(bottomleft,Quadtree<T,RasterT>::QUADTREE_BOTTOM_LEFT_CELL);
    } else {
      QuadNodeInternal<T,RasterT>*  bottomleft = new QuadNodeInternal<T,RasterT>(this->tree,this,rectBottomLeft,!xbalance,!ybalance);
      setChild(bottomleft,Quadtree<T,RasterT>::QUADTREE_BOTTOM_LEFT_CELL);
    }

    // 3 - bottom right
    RasterRect rectBottomRight(_xroot+topleft_xsize,_yroot+topleft_ysize,_xsize-topleft_xsize,_ysize-topleft_ysize);

    if (( rectBottomRight.xsize()<=1)||(rectBottomRight.ysize()<=1)||(this->tree->raster->isHomogeneous(rectBottomRight))) {
      QuadNodeLeaf<T,RasterT>* bottomright = new QuadNodeLeaf<T,RasterT>(this->tree,this,Quadtree<T,RasterT>::QUADTREE_BOTTOM_RIGHT_CELL,rectBottomRight);
      setChild(bottomright,Quadtree<T,RasterT>::QUADTREE_BOTTOM_RIGHT_CELL);          
    } else {
      QuadNodeInternal<T,RasterT>* bottomright = new QuadNodeInternal<T,RasterT>(this->tree,this,rectBottomRight,!xbalance,!ybalance);
      setChild(bottomright,Quadtree<T,RasterT>::QUADTREE_BOTTOM_RIGHT_CELL);          
    }

    // std::cout << "[Quadtree] Devide end " << _rect << std::endl;

    return true;

   } //}}}
/*
   template <typename T, class RasterT>
   int QuadNodeInternal<T,RasterT>::getxroot() const {
    return children[Quadtree<T,RasterT>::QUADTREE_TOP_LEFT_CELL]->getxroot();
   }

   template <typename T, class RasterT>
   int QuadNodeInternal<T,RasterT>::getyroot() const {
    return children[Quadtree<T,RasterT>::QUADTREE_TOP_LEFT_CELL]->getyroot();
   }

   template <typename T, class RasterT>
   int QuadNodeInternal<T,RasterT>::getxsize() const {
    return children[Quadtree<T,RasterT>::QUADTREE_TOP_LEFT_CELL]->getxsize() + children[Quadtree<T,RasterT>::QUADTREE_TOP_RIGHT_CELL]->getxsize();
   }

   template <typename T, class RasterT>
   int QuadNodeInternal<T,RasterT>::getysize() const {
    return children[Quadtree<T,RasterT>::QUADTREE_TOP_LEFT_CELL]->getysize() + children[Quadtree<T,RasterT>::QUADTREE_BOTTOM_LEFT_CELL]->getysize();
   }

  template <typename T, class RasterT>
  RasterRect QuadNodeInternal<T,RasterT>::getRect() const {
    return RasterRect(getxroot(),getyroot(),getxsize(),getysize());
  }
*/

   template <typename T, class RasterT>
   typename QuadNode<T,RasterT>::NODE_TYPE QuadNodeInternal<T,RasterT>::getType() const {
    return  QuadNode<T,RasterT>::INTERNAL_NODE;
   }

   // Get the cluster (tree leaf) for cell
   template <typename T, class RasterT>
   QuadNodeLeaf<T,RasterT>* QuadNodeInternal<T,RasterT>::getLeafForCell(const RasterCellIndex& cell) {

    if (children[Quadtree<T,RasterT>::QUADTREE_TOP_LEFT_CELL]->contains(cell)) {
      return children[Quadtree<T,RasterT>::QUADTREE_TOP_LEFT_CELL]->getLeafForCell(cell);
    }

    if (children[Quadtree<T,RasterT>::QUADTREE_TOP_RIGHT_CELL]->contains(cell)) {
      return children[Quadtree<T,RasterT>::QUADTREE_TOP_RIGHT_CELL]->getLeafForCell(cell);
    }

    if (children[Quadtree<T,RasterT>::QUADTREE_BOTTOM_LEFT_CELL]->contains(cell)) {
      return children[Quadtree<T,RasterT>::QUADTREE_BOTTOM_LEFT_CELL]->getLeafForCell(cell);
    }

    if (children[Quadtree<T,RasterT>::QUADTREE_BOTTOM_RIGHT_CELL]->contains(cell)) {
      return children[Quadtree<T,RasterT>::QUADTREE_BOTTOM_RIGHT_CELL]->getLeafForCell(cell);
    }

    // At this point the only possibility is that we are at the ROOT node
    // and thus we must return the root node as a leaf... tough !
    // For now a bad MEMORY solution is to create a leaf-copy of the root cell
    // on the fly.. on the fly means it is freed no where !
    //return (new QuadNodeLeaf<T,RasterT>(this->tree, this, 0, this->rect, this->tree->getRaster()->getData(cell) ));
    // This would be bad so we'd rather deal with it at a higher level
    return NULL;
   }

   template <typename T, class RasterT>
   QuadNode<T,RasterT>* QuadNodeInternal<T,RasterT>::getChild(int index) {
    if ((index < 0) || (index > 4)) {
      std::cout << "[Quadtree] Child index is invalid" << std::endl;
      return NULL;
    }
    return children[index];
   }

   // get the position (topleft,topright,bottomleft,bottoright) of a child
   template <typename T, class RasterT>
   int QuadNodeInternal<T,RasterT>::getChildPos(QuadNode<T,RasterT>* child) {
    for (int i=0;i<4;i++) {
      if (children[i]==child)
       return i;
    }
    return -1;
   }

   template <typename T, class RasterT>
   bool QuadNodeInternal<T,RasterT>::setChildren(QuadNode<T,RasterT>* S1, QuadNode<T,RasterT>* S2, QuadNode<T,RasterT>* S3, QuadNode<T,RasterT>* S4) {
    // free memory
    delete children[Quadtree<T,RasterT>::QUADTREE_TOP_LEFT_CELL];
    delete children[Quadtree<T,RasterT>::QUADTREE_TOP_RIGHT_CELL];
    delete children[Quadtree<T,RasterT>::QUADTREE_BOTTOM_LEFT_CELL];
    delete children[Quadtree<T,RasterT>::QUADTREE_BOTTOM_RIGHT_CELL];

    children[Quadtree<T,RasterT>::QUADTREE_TOP_LEFT_CELL]=S1;
    children[Quadtree<T,RasterT>::QUADTREE_TOP_RIGHT_CELL]=S2;
    children[Quadtree<T,RasterT>::QUADTREE_BOTTOM_LEFT_CELL]=S3;
    children[Quadtree<T,RasterT>::QUADTREE_BOTTOM_RIGHT_CELL]=S4;
    return true;
   }

   template <typename T, class RasterT>
   bool QuadNodeInternal<T,RasterT>::setChild(QuadNode<T,RasterT>* child,int index) {
    if (child == NULL) {
      std::cout << "[Quadtree] Internal node error, empty children" << std::endl;
      return false;
    } else if ((index<0)||(index>4)) {
      std::cout << "[Quadtree] Child index is invalid" << std::endl;       
      delete child;
        child = NULL;
      return false;
    } else {
      // free mem at corresponding index
      delete children[index];
      // set child
      children[index]=child;
    }
    return true;
   }

  /*}}}*/
  /*--------------------- QuadNodeLeaf ---------------------- {{{*/

  template <typename T, class RasterT>
  QuadNodeLeaf<T,RasterT>::QuadNodeLeaf(Quadtree<T,RasterT>* _tree, QuadNodeInternal<T,RasterT>* _parent, int _nodePos, RasterRect _rect, T _data):
    QuadNode<T,RasterT>(_tree,_parent, _rect),
    nodePos(_nodePos),
    data(_data)
  {
    // update decomp structure
    // if the leaf contains more than 1 cluster
    if (hasManyClusters())  {
      // update the decomp structure in the tree
      this->getTree()->setNewLeaf(this,true);
    } else {
      // update the decomp structure in the tree
      this->getTree()->setNewLeaf(this);
    }

    // std::cout << "[Leaf] update decomp done " << _rect << std::endl;

    if ((data==RasterT::DEFAULT_DATA_VALUE()) && (!isAtomic()) && (!(hasManyClusters()))) {
      // looking data for the cluster (non atomic)
      data = this->tree->getRaster()->getData(this->rect);
    }

    // std::cout << "[Leaf] get data done " << _rect << std::endl;

  }

   template <typename T, class RasterT>
   QuadNodeLeaf<T,RasterT>::~QuadNodeLeaf() {

   }

/*
   template <typename T, class RasterT>
   int QuadNodeLeaf<T,RasterT>::getxroot() const {
    return rect.xtopleft;
   }

   template <typename T, class RasterT>
   int QuadNodeLeaf<T,RasterT>::getyroot() const {
    return rect.ytopleft;
   }

   template <typename T, class RasterT>
   int QuadNodeLeaf<T,RasterT>::getxsize() const {
    return rect.xsize();
   }

   template <typename T, class RasterT>
   int QuadNodeLeaf<T,RasterT>::getysize() const {
    return rect.ysize();
   }

   template <typename T, class RasterT>
   RasterRect QuadNodeLeaf<T,RasterT>::getRect() const {
    return RasterRect(getxroot(),getyroot(),getxsize(),getysize());
   }

*/
   template <typename T, class RasterT>
   typename QuadNode<T,RasterT>::NODE_TYPE QuadNodeLeaf<T,RasterT>::getType() const {
    return QuadNode<T,RasterT>::LEAF_NODE;
   }

   // if the leaf contains only one cluster
   template <typename T, class RasterT>
   bool QuadNodeLeaf<T,RasterT>::hasManyClusters() const {
    if (   ((this->rect.xsize()==1)&&(this->rect.ysize()>1))
       || ((this->rect.xsize()>1)&&(this->rect.ysize()==1)))
      return true;
    return false;
   }

   template <typename T, class RasterT>
   bool QuadNodeLeaf<T,RasterT>::isAtomic() const {
    return ((this->rect.xsize()==1)&&(this->rect.ysize()==1));
   }

  template <typename T, class RasterT>
  QuadNodeLeaf<T,RasterT>* QuadNodeLeaf<T,RasterT>::getLeafForCell(const RasterCellIndex& cell) {//{{{
    if (this->contains(cell)) {
      return this;
    } else {
      return NULL;
    }
  }//}}}

  template <typename T, class RasterT>
  bool QuadNodeLeaf<T,RasterT>::getClusterForCell(const RasterCellIndex& cell, RasterRect& cluster) { //{{{
    if (this->contains(cell)) {
      if (hasManyClusters()) {
        cluster.setRect(cell.i,cell.j,1,1);
      } else {
        cluster=this->rect;
      }
      return true;
    } else {
      return false;
    }
  }//}}}

   template <typename T, class RasterT>
   T QuadNodeLeaf<T,RasterT>::getData(const RasterCellIndex& clusterRoot) {//{{{
    // if the leaf has only one cluster => return the 0-index element inside the data vector
    // if not, find the element by the index computed from substraction of the wanted cluster root indices and the leaf cluster root indices
    // if the leaf cluster is one raster cell
    if (isAtomic()||hasManyClusters())
      return this->tree->getRaster()->getData(clusterRoot);
    else
      return data;
 } //}}}

   /*}}}*/
   /*------------------------- IO operators -----------------------{{{*/

    template <typename T, class RasterT>
    std::ostream& operator<< (std::ostream& out, const QuadNodeLeaf<T,RasterT>& node) {
      out << node.getRect() << std::endl;
      return out;
    }

    template <typename T, class RasterT>
    std::ostream& operator<< (std::ostream& out, const QuadNodeInternal<T,RasterT>& node) {
      out << node.getRect() << std::endl;
      return out;
    }

    template <typename T, typename RasterT>
    std::ostream& operator<< (std::ostream& out, const Quadtree<T,RasterT>& tree) { /*{{{*/

      std::list<QuadNodeLeaf<T,RasterT>*> leaves = tree.getFullLeavesList();

      // "typename" is required here to make the compiler understand that
      // "iterator" is not a static variable but really a type of STL list.
      typename std::list<QuadNodeLeaf<T,RasterT>* >::iterator  it;

      if (tree.getRootCell() != NULL)
      {
        out << "# Quadtree root node : " << tree.getRootCell() << std::endl;
        out << "# Quadtree root node : " << (tree.getRootCell())->getRect() << std::endl;
        out << "# Number of leaves : " << leaves.size() << std::endl;
        out << "# Xcenter    Ycenter   Xsize    Ysize" << std::endl;

        double xleft = 0.0;
        double ytop = 0.0;
        double xsize = 0.0;
        double ysize = 0.0;

        // Gnuplot header: some useful parameters
        out << "\tset key off" << std::endl;

        for (it = leaves.begin(); it != leaves.end(); ++it)
        {
          /* Prints each leaf's coordinates */
          xleft = (double)(*it)->getxroot();
          ytop  = -(double)(*it)->getyroot(); /* Minus to get into Raster origin on top left corner */
          xsize = (double)(*it)->getxsize();
          ysize = (double)(*it)->getysize();

          /* Converting to be gnuplot compatible.
           * Gnuplot needs : rectangle center & rectangle size
           * (or left down corner coord & top right corner coord)
           */
          xleft = xleft + xsize/2.0; /* xleft is to be considered as xcenter now */
          ytop  = ytop  - ysize/2.0; /* ytop  is to be considered as ycenter now */

          out << "\tset object rect center "<< xleft <<","<< ytop
              <<" size "<<xsize<<","<<ysize << " back fs empty border 3" << std::endl;
        }

        out << "\tplot [0:" << ((tree.getRootCell())->getRect()).xsize() << "] [-" << ((tree.getRootCell())->getRect()).ysize() << ":0] 0"
          << std::endl << "\tpause -1" << std::endl;
      } else {
        out << "# Quadtree has not root cell. Containing only 1 leaf." << std::endl;
        out << "# Number of leaves : " << leaves.size() << std::endl;
      }

      // Clearing the copied list
      leaves.clear();
      return out;
    } /*}}}*/

   /*--------------- end of IO Operators }}}*/

  } /* end of namespace lgl */
} /* end of namespace jafar */

#endif // _QUADTREE_HPP