grid.hpp

Go to the documentation of this file.

#ifndef _GRID_HPP
#define _GRID_HPP

#include <assert.h>
#include <list>
#include "geometry.hpp"

namespace Arak {

  template <class ItemType>
  class Grid {

  public:

    class Cell : public Geometry::Rectangle {

      friend class Grid;
      friend class Grid::LineCellIterator;
      friend class Grid::ConstLineCellIterator;
      friend class Grid::RectangleCellIterator;
      friend class Grid::ConstRectangleCellIterator;
      friend class Grid::TriangleCellIterator;
      friend class Grid::ConstTriangleCellIterator;
      friend class Grid::ConvexPolygonCellIterator;
      friend class Grid::ConstConvexPolygonCellIterator;
      
    public:

      typedef std::list<ItemType> ItemList;

    protected:

      typedef typename ItemList::iterator ItemIterator;

      int r;

      int c;

      ItemList items;
      
      ItemList unused;
      
      Cell* up;

      Cell* down;

      Cell* right;

      Cell* left;

    public:

      class Entry {

      private:

  Cell* cell;

  ItemIterator it;

      public:

  Entry() : cell(NULL), it() { }

  Entry(Cell* cell, ItemIterator it) : cell(cell), it(it) { }

  bool valid() const { return (cell != NULL); }

  void remove() {
    assert(valid());
    // Move the item to the unused list.
    cell->unused.splice(cell->unused.begin(), cell->items, it);
    cell = NULL;
  }

  const Cell& getCell() const { return *cell; }

      }; // end of class: Grid::Cell::Entry

      Cell(const Geometry::Rectangle& r, int row, int col) : 
  Geometry::Rectangle(r), r(row), c(col), items(),
  up(NULL), down(NULL), right(NULL), left(NULL) { }

      inline Entry add(const ItemType& e) {
  if (unused.empty()) {
    // We have run out of unused items.  Use a simple insertion.
    return Entry(this, items.insert(items.begin(), e));
  } else {
    // The unused list is not empty; reuse its first element.
    ItemIterator it = unused.begin();
    *it = e;
    items.splice(items.begin(), unused, it);
    return Entry(this, items.begin());
  }
      }

      inline const ItemList& getItemList() const {
  return items;
      }

      inline ItemList& getItemList() {
  return items;
      }

      int row() const { return r; }

      int col() const { return c; }

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

      Geometry::Point center() const {
  return CGAL::midpoint((*this)[0], (*this)[2]);
      }

    }; // end of class: Grid::Cell

    template <class GridType, class CellType>
    class LineCellIteratorBase {

    protected:
      
      Geometry::Line line;

      enum CompassDirection {
  NORTH,      
  NORTH_EAST, 
  EAST,       
  SOUTH_EAST, 
  SOUTH,      
  SOUTH_WEST, 
  WEST,       
  NORTH_WEST  
      };

      CompassDirection dir;

      CellType* curCell;

      CellType* lastCell;

      inline CompassDirection direction(const Geometry::Point& p,
          const Geometry::Point& q) {
  CGAL::Comparison_result xc = CGAL::compare_x(p, q);
  CGAL::Comparison_result yc = CGAL::compare_y(p, q);
  switch (xc) {
  case CGAL::SMALLER:
    switch (yc) {
    case CGAL::SMALLER:
      return NORTH_EAST; 
    case CGAL::LARGER:
      return SOUTH_EAST; 
    case CGAL::EQUAL:
    default:
      return EAST;
    }
  case CGAL::LARGER:
    switch (yc) {
    case CGAL::SMALLER:
      return NORTH_WEST;
    case CGAL::LARGER:
      return SOUTH_WEST;
    case CGAL::EQUAL:
    default:
      return WEST;
    }
  case CGAL::EQUAL:
  default:
    switch (yc) {
    case CGAL::SMALLER:
      return NORTH;
    case CGAL::LARGER:
      return SOUTH;
    case CGAL::EQUAL:
    default:
      assert(false);
    }
  }
      }

    public:

      struct SEGMENT { };

      struct RAY { };

      LineCellIteratorBase(GridType& g, 
         const Geometry::Point &p,
         const Geometry::Point &q,
         SEGMENT) : line(p, q) {
  dir = direction(p, q);
  if (g.boundary().has_on_unbounded_side(p) ||
      g.boundary().has_on_unbounded_side(q)) {
    Geometry::Segment s(p, q);
    Geometry::Point r;
    CGAL::Object intersection = CGAL::intersection(s, g.boundary());
    if (CGAL::assign(s, intersection)) {
      curCell = &g.getCell(Geometry::project(s.source(), g.boundary()));
      lastCell = &g.getCell(Geometry::project(s.target(), g.boundary()));
    } else if (CGAL::assign(r, intersection)) {
      curCell = lastCell = &g.getCell(Geometry::project(r, g.boundary()));
    } else {
      curCell = lastCell = NULL;
    }
  } else {
    curCell = &g.getCell(p);
    lastCell = &g.getCell(q);
  }
      }

      LineCellIteratorBase(GridType& g, 
         const Geometry::Point &p,
         const Geometry::Point &q,
         RAY) : line(p, q) {
  dir = direction(p, q);
  if (g.boundary().has_on_unbounded_side(p) ||
      g.boundary().has_on_unbounded_side(q)) {
    Geometry::Segment s(p, q);
    Geometry::Point r;
    CGAL::Object intersection = CGAL::intersection(s, g.boundary());
    if (CGAL::assign(s, intersection)) {
      curCell = &g.getCell(Geometry::project(s.source(), g.boundary()));
      lastCell = &g.getCell(Geometry::project(s.target(), g.boundary()));
    } else if (CGAL::assign(r, intersection)) {
      curCell = lastCell = &g.getCell(Geometry::project(r, g.boundary()));
    } else {
      curCell = lastCell = NULL;
    }
  } else {
    curCell = &g.getCell(p);
    lastCell = NULL;
  }
      }

      LineCellIteratorBase() 
  : line(1, 0, 0),
    curCell(NULL), 
    lastCell(NULL) { }

      inline bool 
      operator==(const LineCellIteratorBase<GridType,CellType>& other) const
      { return (curCell == other.curCell); }
      inline bool 
      operator!=(const LineCellIteratorBase<GridType,CellType>& other) const 
      { return !(*this == other); }

      inline LineCellIteratorBase<GridType,CellType>& operator++(int) {
  LineCellIteratorBase<GridType,CellType> tmp(*this);
  ++(*this);
  return tmp;
      }
      inline LineCellIteratorBase<GridType,CellType>& operator++() {
  using namespace Arak::Geometry;
  assert(curCell != NULL);
  if (curCell == lastCell) {
    curCell = NULL;
    return *this;
  }
  switch (dir) {
  case NORTH_EAST:
    switch (CGAL::compare_y_at_x(curCell->vertex(2), line)) {
    case CGAL::SMALLER:
      curCell = curCell->up; break;
    case CGAL::LARGER:
      curCell = curCell->right; break;
    case CGAL::EQUAL:
      curCell = curCell->right;
      if (curCell != NULL)
        curCell = curCell->up; 
      break;
    }
    break;
  case SOUTH_EAST:
    switch (CGAL::compare_y_at_x(curCell->vertex(1), line)) {
    case CGAL::LARGER:
      curCell = curCell->down; break;
    case CGAL::SMALLER:
      curCell = curCell->right; break;
    case CGAL::EQUAL:
      curCell = curCell->right;
      if (curCell != NULL)
        curCell = curCell->down; 
      break;
    }
    break;
  case SOUTH_WEST:
    switch (CGAL::compare_y_at_x(curCell->vertex(0), line)) {
    case CGAL::LARGER:
      curCell = curCell->down; break;
    case CGAL::SMALLER:
      curCell = curCell->left; break;
    case CGAL::EQUAL:
      curCell = curCell->left; 
      if (curCell != NULL)
        curCell = curCell->down; 
      break;
    }
    break;
  case NORTH_WEST:
    switch (CGAL::compare_y_at_x(curCell->vertex(3), line)) {
    case CGAL::SMALLER:
      curCell = curCell->up; break;
    case CGAL::LARGER:
      curCell = curCell->left; break;
    case CGAL::EQUAL:
      curCell = curCell->left; 
      if (curCell != NULL)
        curCell = curCell->up; 
      break;
    }
    break;
  case NORTH:
    curCell = curCell->up; 
    break;
  case SOUTH:
    curCell = curCell->down; 
    break;
  case EAST:
    curCell = curCell->right; 
    break;
  case WEST:
    curCell = curCell->left; 
    break;
  }
  return *this;
      }

      inline CellType* operator->() const {return curCell;}
      inline CellType& operator*() const { return *curCell;}

    }; // end of class: Grid::LineCellIteratorBase

    typedef LineCellIteratorBase<Grid, Cell> LineCellIterator;

    typedef LineCellIteratorBase<const Grid, const Cell> ConstLineCellIterator;

    template <class GridType, class CellType>
    class RectangleCellIteratorBase {

    protected:
      
      int i;

      int imax;

      int j;

      int jmin;

      int jmax;

      CellType* curCell;

      GridType* grid;

    public:

      RectangleCellIteratorBase(GridType& g, Geometry::Rectangle r)
  : curCell(NULL), grid(&g) { 
  if (CGAL::assign(r, CGAL::intersection(g.boundary(), r))) {
    assert(g.getCellCoords(r.min(), i, jmin));
    j = jmin;
    assert(g.getCellCoords(r.max(), imax, jmax));
    curCell = &g.getCell(i, j);
  } else {
    curCell = NULL;
  }
      }

      RectangleCellIteratorBase() : curCell(NULL), grid(NULL) { }

      inline bool operator==(const RectangleCellIteratorBase<GridType,CellType>& other) const
      { return (curCell == other.curCell); }
      inline bool operator!=(const RectangleCellIteratorBase<GridType,CellType>& other) const 
      { return !(*this == other); }

      inline RectangleCellIteratorBase<GridType,CellType>& operator++(int) {
  RectangleCellIteratorBase<GridType,CellType> tmp(*this);
  ++(*this);
  return tmp;
      }
      inline RectangleCellIteratorBase<GridType,CellType>& operator++() {
  using namespace Arak::Geometry;
  assert(curCell != NULL);
  if ((i == imax) && (j == jmax)) {
    curCell = NULL;
    return *this;
  } else if (j == jmax) {
    i++;
    j = jmin;
    curCell = &(grid->getCell(i, j));
  } else {
    j++;
    curCell = curCell->right;
  }
  return *this;
      }

      inline CellType* operator->() const { return curCell; }
      inline CellType& operator*() const { return *curCell; }

    }; // end of class: Grid::RectangleCellIteratorBase

    typedef RectangleCellIteratorBase<Grid, Cell> 
    RectangleCellIterator;

    typedef RectangleCellIteratorBase<const Grid, const Cell> 
    ConstRectangleCellIterator;

    template <class GridType, class CellType>
    class TriangleCellIteratorBase {

    protected:
      
      CellType* curCell;

      GridType* grid;

      int imin;

      int imax;

      std::vector<int> jmins;

      std::vector<int> jmaxs;

    public:

      TriangleCellIteratorBase(GridType& g, Geometry::Triangle t)
  : curCell(NULL), grid(&g) { 
  // Compute the smallest rectangle that covers the intersection
  // of the triangle and the grid's closure.
  Geometry::Rectangle r = Geometry::boundingRectangle(t[0], t[1], t[2]);
  if (!CGAL::assign(r, CGAL::intersection(g.boundary(), r))) {
    // The triangle does not intersect the grid's closure.
    curCell = NULL; 
    return;
  }
  // Compute the columns crossed by the triangle.
  int jmin, jmax;
  assert(g.getCellCoords(r.min(), imin, jmin));
  assert(g.getCellCoords(r.max(), imax, jmax));
  assert(imin <= imax);
  assert(jmin <= jmax);
  // Compute the values of jmins and jmaxs using line searches
  // along the boundary of the triangle.
  jmins.insert(jmins.end(), imax - imin + 1, 
         std::numeric_limits<int>::max());
  jmaxs.insert(jmaxs.end(), imax - imin + 1, 
         std::numeric_limits<int>::min());
  for (int v = 0; v < 3; v++) {
    ConstLineCellIterator it(g, t[v], t[v + 1], 
           ConstLineCellIterator::SEGMENT()), end;
    while (it != end) {
      const Cell& cell = *it;
      int i = cell.row();
      int j = cell.col();
      int k = i - imin;
      jmins[k] = std::min<int>(jmins[k], j);
      jmaxs[k] = std::max<int>(jmaxs[k], j);
      ++it;
    }
  }
  // Replace any unset limits using the bounding rectangle.
  std::replace(jmins.begin(), jmins.end(),
         std::numeric_limits<int>::max(), jmin);
  std::replace(jmaxs.begin(), jmaxs.end(),
         std::numeric_limits<int>::min(), jmax);
  // Initialize the iterator at the first cell.
  curCell = &g.getCell(imin, jmins[0]);
      }

      TriangleCellIteratorBase() : curCell(NULL), grid(NULL) { }

      inline bool operator==(const TriangleCellIteratorBase<GridType,CellType>& other) const
      { return (curCell == other.curCell); }
      inline bool operator!=(const TriangleCellIteratorBase<GridType,CellType>& other) const 
      { return !(*this == other); }

      inline TriangleCellIteratorBase<GridType,CellType>& operator++(int) {
  TriangleCellIterator tmp(*this);
  ++(*this);
  return tmp;
      }
      inline TriangleCellIteratorBase<GridType,CellType>& operator++() {
  using namespace Arak::Geometry;
  assert(curCell != NULL);
  int i = curCell->row();
  int j = curCell->col();
  int k = i - imin;
  if ((i == imax) && (j == jmaxs[k])) {
    // We are at the last cell; the iterator is finished.
    curCell = NULL;
    return *this;
  } else if (j == jmaxs[k])
    // We are at the last column but not the last row; wrap
    // around.
    curCell = &(grid->getCell(i + 1, jmins[k + 1]));
  else 
    // We have not reached the end of the current row; move
    // right.
    curCell = curCell->right;
  return *this;
      }

      inline CellType* operator->() const { return curCell; }
      inline CellType& operator*() const { return *curCell; }
      
    }; // end of class: Grid::TriangleCellIteratorBase

    typedef TriangleCellIteratorBase<Grid, Cell> 
    TriangleCellIterator;

    typedef TriangleCellIteratorBase<const Grid, const Cell> 
    ConstTriangleCellIterator;

    template <class GridType, class CellType>
    class ConvexPolygonCellIteratorBase {

    protected:
      
      CellType* curCell;

      GridType* grid;

      int imin;

      int imax;

      std::vector<int> jmins;

      std::vector<int> jmaxs;

    public:

      ConvexPolygonCellIteratorBase(GridType& g, const Geometry::Polygon& p)
  : curCell(NULL), grid(&g) { 
  // Compute the smallest rectangle that covers the intersection
  // of the triangle and the grid's closure.
  Geometry::Rectangle r(*(p.left_vertex()),
            *(p.right_vertex()),
            *(p.bottom_vertex()),
            *(p.top_vertex()));
  if (!CGAL::assign(r, CGAL::intersection(g.boundary(), r))) {
    // The polygon does not intersect the grid's closure.
    curCell = NULL; 
    return;
  }
  // Compute the columns crossed by the triangle.
  int jmin, jmax;
  assert(g.getCellCoords(r.min(), imin, jmin));
  assert(g.getCellCoords(r.max(), imax, jmax));
  assert(imin <= imax);
  assert(jmin <= jmax);
  // Compute the values of jmins and jmaxs using line searches
  // along the boundary of the polygon.
  jmins.insert(jmins.end(), imax - imin + 1, 
         std::numeric_limits<int>::max());
  jmaxs.insert(jmaxs.end(), imax - imin + 1, 
         std::numeric_limits<int>::min());
  const int size = p.size();
  for (int v = 0; v < size; v++) {
    ConstLineCellIterator it(g, p[v], p[(v + 1) % size], 
           ConstLineCellIterator::SEGMENT()), end;
    while (it != end) {
      const Cell& cell = *it;
      int i = cell.row();
      int j = cell.col();
      int k = i - imin;
      jmins[k] = std::min<int>(jmins[k], j);
      jmaxs[k] = std::max<int>(jmaxs[k], j);
      ++it;
    }
  }
  // Replace any unset limits using the bounding rectangle.
  std::replace(jmins.begin(), jmins.end(),
         std::numeric_limits<int>::max(), jmin);
  std::replace(jmaxs.begin(), jmaxs.end(),
         std::numeric_limits<int>::min(), jmax);
  // Initialize the iterator at the first cell.
  curCell = &g.getCell(imin, jmins[0]);
      }

      ConvexPolygonCellIteratorBase() : curCell(NULL), grid(NULL) { }

      inline bool operator==(const ConvexPolygonCellIteratorBase<GridType,CellType>& other) const
      { return (curCell == other.curCell); }
      inline bool operator!=(const ConvexPolygonCellIteratorBase<GridType,CellType>& other) const 
      { return !(*this == other); }

      inline ConvexPolygonCellIteratorBase<GridType,CellType>& operator++(int) {
  ConvexPolygonCellIterator tmp(*this);
  ++(*this);
  return tmp;
      }
      inline ConvexPolygonCellIteratorBase<GridType,CellType>& operator++() {
  using namespace Arak::Geometry;
  assert(curCell != NULL);
  int i = curCell->row();
  int j = curCell->col();
  int k = i - imin;
  if ((i == imax) && (j == jmaxs[k])) {
    // We are at the last cell; the iterator is finished.
    curCell = NULL;
    return *this;
  } else if (j == jmaxs[k])
    // We are at the last column but not the last row; wrap
    // around.
    curCell = &(grid->getCell(i + 1, jmins[k + 1]));
  else 
    // We have not reached the end of the current row; move
    // right.
    curCell = curCell->right;
  return *this;
      }

      inline CellType* operator->() const { return curCell; }
      inline CellType& operator*() const { return *curCell; }
      
    }; // end of class: Grid::ConvexPolygonCellIteratorBase

    typedef ConvexPolygonCellIteratorBase<Grid, Cell> 
    ConvexPolygonCellIterator;

    typedef ConvexPolygonCellIteratorBase<const Grid, const Cell> 
    ConstConvexPolygonCellIterator;

  protected:

    const Geometry::Rectangle bd;

    std::vector< std::vector<Cell*> > cells;

    const int rows;

    const int cols;

    Geometry::Kernel::FT cellWidth;
    
    Geometry::Kernel::FT cellHeight;

    std::vector<Geometry::Kernel::FT> xvals;

    std::vector<Geometry::Kernel::FT> yvals;

  public:
    
    Grid(const Geometry::Rectangle& boundary, int rows, int cols)
      : bd(boundary), rows(rows), cols(cols) {
      using namespace Arak::Geometry;
      // Compute the cell size.
      Kernel::FT xmin = bd.xmin();
      Kernel::FT xmax = bd.xmax();
      Kernel::FT ymin = bd.ymin();
      Kernel::FT ymax = bd.ymax();
      this->cellWidth = (xmax - xmin) / Kernel::FT(cols);
      this->cellHeight = (ymax - ymin) / Kernel::FT(rows);
      // Compute the grid lines.
      xvals.reserve(cols + 1);
      Kernel::FT x = xmin;
      for (int i = 0; i < cols; i++) {
  xvals.push_back(x);
  x += cellWidth;
      }
      xvals.push_back(xmax);
      yvals.reserve(rows + 1);
      Kernel::FT y = ymin;
      for (int i = 0; i < rows; i++) {
  yvals.push_back(y);
  y += cellHeight;
      }
      yvals.push_back(ymax);
      // Create the cells.
      cells.reserve(rows);
      for (int i = 0; i < rows; i++) {
  cells.push_back(std::vector<Cell*>());
  cells[i].reserve(cols);
  for (int j = 0; j < cols; j++) {
    Geometry::Rectangle r(xvals[j], yvals[i], 
        xvals[j + 1], yvals[i + 1]);
    cells[i][j] = new Cell(r, i, j);
    if (i > 0) {
      cells[i][j]->down = cells[i - 1][j];
      cells[i - 1][j]->up = cells[i][j];
    }
    if (j > 0) {
      cells[i][j]->left = cells[i][j - 1];
      cells[i][j - 1]->right = cells[i][j];
    }
  }
      }
    }
  
    ~Grid() {
      
    }

    inline int numRows() const { return rows; }

    inline int numCols() const { return cols; }

    const Geometry::Rectangle& boundary() const { return bd; }

    bool getCellCoords(const Geometry::Point& p, int& i, int& j) const {
      if (bd.has_on_unbounded_side(p)) return false;
      // First we use a simple computation to estimate the correct
      // cell.  This may be only approximate, though, due to floating
      // point error.  Then we shift the indexes as needed to compute
      // the exact cell.
      if (p.y() == bd.ymax()) 
  i = rows - 1;
      else {
  i = int(floor(CGAL::to_double((p.y() - bd.ymin()) / cellHeight)));
  while (p.y() < yvals[i]) i--;
  while (p.y() >= yvals[i + 1]) i++;
      }
      if (p.x() == bd.xmax()) 
  j = cols - 1;
      else {
  j = int(floor(CGAL::to_double((p.x() - bd.xmin()) / cellWidth)));
  while (p.x() < xvals[j]) j--;
  while (p.x() >= xvals[j + 1]) j++;
      }
      return true;
    }
    
    inline const Cell& getCell(const int i, const int j) const {
      return *(cells[i][j]);
    }

    inline Cell& getCell(const int i, const int j) {
      return *(cells[i][j]);
    }

    const Cell& getCell(const Geometry::Point& p) const {
      int i = -1, j = -1;
      assert(getCellCoords(p, i, j));
      return getCell(i, j);
    }

    Cell& getCell(const Geometry::Point& p) {
      int i = -1, j = -1;
      assert(getCellCoords(p, i, j));
      return *(cells[i][j]);
    }

  };  // end of class: Grid

}

#endif

---