cone_tracing.hpp

Go to the documentation of this file.

#ifndef _CONE_TRACING_HPP
#define _CONE_TRACING_HPP

#include "handle.hpp"
#include "geometry.hpp"
#include "coloring.hpp"
#include "grid.hpp"

namespace Arak {

  inline Geometry::Triangle outerApprox(const Geometry::Point& vertex,
          const Geometry::Direction& min,
          const Geometry::Direction& max,
          const Geometry::Kernel::FT& radius) {
    // Compute the other two vertices of the largest inner triangular
    // approximation to the cone.
    Geometry::Point p = vertex + min.vector();
    double length = sqrt(CGAL::to_double(CGAL::squared_distance(p, vertex)));
    p = vertex + (p - vertex) * (CGAL::to_double(radius) / length);
    Geometry::Point q = vertex + max.vector();
    length = sqrt(CGAL::to_double(CGAL::squared_distance(q, vertex)));
    q = vertex + (q - vertex) * (CGAL::to_double(radius) / length);
    // Squared hypotenuse length of inner triangle.
    const Geometry::Kernel::FT dsq = CGAL::squared_distance(vertex, p);
    // Squared base length of the inner triangle
    const Geometry::Kernel::FT bsq = CGAL::squared_distance(p, q);
    // Squared height of the inner triangle
    const double hsq = CGAL::to_double(dsq) - CGAL::to_double(bsq) / 4.0;
    // The inner and outer triangles are similar.  So the ratio
    // between d and h (the heights of the outer and inner
    // triangles) is the same as the ratio between the hypotenuse
    // length of the outer triangle and d, the hypotenuse length of
    // the inner triangle.
    double scale = sqrt(CGAL::to_double(dsq / hsq));
    return Geometry::Triangle(vertex, 
            vertex + (p - vertex) * scale, 
            vertex + (q - vertex) * scale);
  }

  inline Geometry::Point intersection(const Geometry::Point& p,
              const Geometry::Direction& d,
              const Geometry::Segment& s) {
    const Geometry::Line line1 = s.supporting_line();
    const Geometry::Line line2(p, d);
    Geometry::Point contact;
    assert(CGAL::assign(contact, CGAL::intersection(line1, line2)));
    return contact;
  }

  // Forward declaration.
  class ScanEdge;

  // A handle to a scan edge.
  typedef PointerHandle<ScanEdge> ScanEdgeHandle;

  struct ScanVertex {

    ScanEdgeHandle edge;

    bool isMin;

    ScanVertex(ScanEdgeHandle edge, bool isMin) : edge(edge), isMin(isMin) { }

  };

  struct ScanEdge {

    Coloring::IntEdgeHandle edge;

    Geometry::Direction dmin;

    Geometry::Point pmin;

    Geometry::Direction dmax;

    Geometry::Point pmax;

    bool pruned;

    std::list<ScanEdgeHandle>::iterator it;

    void init(const Geometry::Point& vertex,
        const Geometry::Direction& min,
        const Geometry::Direction& max,
        Coloring::IntEdgeHandle e) {
      this->pruned = true;
      this->edge = e;
      // Get the source, target, and midpoint of the edge.
      const Geometry::Point& source = e->source();
      const Geometry::Point& target = e->target();
      Geometry::Point midpoint = CGAL::midpoint(source, target);
      // Compute the directions to these points from the vertex.
      Geometry::Segment toSource(vertex, source);
      Geometry::Segment toTarget(vertex, target);
      Geometry::Segment toMidpoint(vertex, midpoint);
      Geometry::Direction ds = toSource.direction();
      Geometry::Direction dt = toTarget.direction();
      Geometry::Direction dm = toMidpoint.direction();
      // Determine the (counter-clockwise) order of these directions.
      bool sourceIsMin = dm.counterclockwise_in_between(ds, dt);
      if (sourceIsMin) {
  dmin = ds;
  dmax = dt;
      } else {
  dmin = dt;
  dmax = ds;
      }
      // Compute the minimum visible point on the segment.
      if (min.counterclockwise_in_between(dmin, dmax)) {
  dmin = min;
  pmin = intersection(vertex, min, e->segment());
      } else
  pmin = (sourceIsMin ? source : target);
      // Compute the maximum visible point on the segment.
      if (max.counterclockwise_in_between(dmin, dmax)) {
  dmax = max;
  pmax = intersection(vertex, max, e->segment());
      } else
  pmax = (sourceIsMin ? target : source);
    }

    inline ScanVertex min() {
      return ScanVertex(ScanEdgeHandle(this), true);
    };

    inline ScanVertex max() {
      return ScanVertex(ScanEdgeHandle(this), false);
    };
  };

  inline CGAL::Comparison_result compare(const Geometry::Direction& a,
           const Geometry::Direction& b,
           const Geometry::Direction& min) {
    if (a == b) return CGAL::EQUAL;
    if (a == min) return CGAL::SMALLER;
    if (b == min) return CGAL::LARGER;
    if (a.counterclockwise_in_between(min, b)) return CGAL::SMALLER;
    return CGAL::LARGER;
  }

  class ScanVertexComparator {

  private:

    Geometry::Direction min;

  public:

    ScanVertexComparator(const Geometry::Direction& min) : min(min) { }

    inline bool operator() (const ScanVertex& a,
          const ScanVertex& b) const {
      const Geometry::Direction& adir = 
  (a.isMin ? a.edge->dmin : a.edge->dmax);
      const Geometry::Direction& bdir = 
  (b.isMin ? b.edge->dmin : b.edge->dmax);
      return (compare(adir, bdir, min) == CGAL::SMALLER);
    }
  };

  inline Geometry::Point intersection(const Geometry::Point& p,
              const Geometry::Direction& d,
              const ScanEdgeHandle se) {
    if (d == se->dmin) return se->pmin;
    if (d == se->dmax) return se->pmax;
    return intersection(p, d, se->edge->segment());
  }

  inline bool cones_intersect(const Geometry::Direction& amin,
            const Geometry::Direction& amax,
            const Geometry::Direction& bmin,
            const Geometry::Direction& bmax) {
    return ((amin == bmin) ||
      (amax == bmax) ||
      bmin.counterclockwise_in_between(amin, amax) ||
      bmax.counterclockwise_in_between(amin, amax) ||
      amin.counterclockwise_in_between(bmin, bmax) ||
      amax.counterclockwise_in_between(bmin, bmax));
  }

  inline bool cone_subset(const Geometry::Direction& amin,
        const Geometry::Direction& amax,
        const Geometry::Direction& bmin,
        const Geometry::Direction& bmax) {
    return (((amin == bmin) || 
       amin.counterclockwise_in_between(bmin, bmax))
      &&
      ((amax == bmax) || 
       amax.counterclockwise_in_between(bmin, bmax)));
  }

  inline bool closer(const Geometry::Point& p,
         const Geometry::Direction& dir,
         ScanEdgeHandle a,
         ScanEdgeHandle b) {
    Geometry::Point ap = intersection(p, dir, a);
    Geometry::Point bp = intersection(p, dir, b);
    if (ap == bp) {
      bool a_at_min = (ap == a->pmin);
      bool b_at_min = (bp == b->pmin);
      if (a_at_min) {
  if (b_at_min) {
    // a_at_min && b_at_min
    return (CGAL::orientation(b->pmin, b->pmax, a->pmax) ==
      CGAL::orientation(b->pmin, b->pmax, p));
  } else {
    // a_at_min && !b_at_min
    return false;
  }
      } else {
  if (b_at_min) {
    // !a_at_min && b_at_min
    return false;
  } else {
    // !a_at_min && !b_at_min
    return (CGAL::orientation(b->pmin, b->pmax, a->pmax) ==
      CGAL::orientation(b->pmin, b->pmax, p));
  }
      }
    } else
      return (CGAL::compare_distance_to_point(p, ap, bp) == CGAL::SMALLER);
  }

  struct ScanEvent {

    Coloring::IntEdgeHandle edge;

    Geometry::Direction dmin;

    Geometry::Direction dmax;

    Geometry::Point pmin;

    Geometry::Point pmax;
  };

  template <class OutputIterator>
  void coneTrace(const Coloring& coloring,
     const Geometry::Point& vertex,
     const Geometry::Direction& min,
     const Geometry::Direction& max,
     const Geometry::Kernel::FT& radius,
     const Geometry::Triangle& filter,
     OutputIterator out) {
    typedef typename std::set<Coloring::IntEdgeHandle> EdgeHandleSet;
    static EdgeHandleSet edges; // Warning: this is not thread-safe
    edges.clear();
    const Coloring::InteriorEdgeIndex& index = 
      coloring.getInteriorEdgeIndex();
#ifdef RECTANGULAR_FILTER
    Geometry::Rectangle rFilter = boundingRectangle(filter[0],
                filter[1],
                filter[2]);
    Coloring::InteriorEdgeIndex::ConstRectangleCellIterator 
      it(index, rFilter), end;
#else
    Coloring::InteriorEdgeIndex::ConstTriangleCellIterator 
      it(index, filter), end;
#endif
    while (it != end) {
      const Coloring::InteriorEdgeIndex::Cell::ItemList& obsList = 
  it->getItemList();
      edges.insert(obsList.begin(), obsList.end());
      ++it;
    }
    // If there are no edges, we are done.
    if (edges.empty()) return;
    // Create scan edges and vertices.  The scan edge objects are
    // allocated statically for efficiency.
    const size_t MAX_NUM_SCAN_EDGES = 1024;
    static ScanEdge _se[MAX_NUM_SCAN_EDGES];
    assert(edges.size() <= MAX_NUM_SCAN_EDGES);
    // ScanEdge* _se = new ScanEdge[edges.size()];
    typedef typename std::vector<ScanVertex> ScanVertexSequence;
    static ScanVertexSequence scanVertices; // Warning: not thread-safe
    scanVertices.clear();
    int i = 0;
    for (EdgeHandleSet::iterator it = edges.begin(); it != edges.end(); it++) {
      _se[i].init(vertex, min, max, *it);
      ScanEdgeHandle seh(_se[i]);
      i++;
      // If the scan edge does not intersect the scan cone, ignore it.
      if (!cones_intersect(min, max, seh->dmin, seh->dmax)) continue;
      // Store the associated scan vertices.
      scanVertices.push_back(seh->max());
      scanVertices.push_back(seh->min());
    }
    // If there are no scan vertices, we are done.
    if (scanVertices.empty()) return;
    // Sort the scan vertices by direction.
    ScanVertexComparator comp(min);
    std::sort(scanVertices.begin(), scanVertices.end(), comp);
    /* Scan through the sorted vertices.  When a minimum vertex is
     * encountered, insert its edge into the list in order of its
     * depth along the current angle.  When a maximum vertex is
     * encountered, remove its edge from the list.  Whenever a change
     * occurs, append an event to the list of events.
     */
    // The list of scan edges, sorted by depth.
    typedef std::list<ScanEdgeHandle> ScanEdgeHandleList;
    ScanEdgeHandleList front_to_back;
    // Insert a sentinel in the list which means 'no current edge'.
    ScanEdgeHandle front;
    front_to_back.push_back(front);
    // The current scan event.
    ScanEvent curEvent;
    curEvent.edge = Coloring::IntEdgeHandle();
    // Scan through the vertices.
    for (ScanVertexSequence::iterator svi = scanVertices.begin(); 
   svi != scanVertices.end(); svi++) {
      // Get a handle to the edge that is currently in front.
      front = front_to_back.front();
      // Get the direction to the current vertex.
      const Geometry::Direction& dir = 
  (svi->isMin ? svi->edge->dmin : svi->edge->dmax);
      // Check to see if the current scan vertex is a minimum vertex
      // or a maximum vertex.
      if (svi->isMin) {
  // The current vertex is a minimum vertex.  Insert its scan
  // edge into sorted position in the list.
  for (ScanEdgeHandleList::iterator sehi = front_to_back.begin();
       sehi != front_to_back.end(); ++sehi) {
    ScanEdgeHandle seh = *sehi;
    if (!seh.valid() || closer(vertex, dir, svi->edge, seh)) {
      svi->edge->it = front_to_back.insert(sehi, svi->edge);
      svi->edge->pruned = false;
      // Prune any deeper edges that are dominated by this one.
      if (seh.valid()) {
        ++sehi;
        while ((seh = *sehi).valid()) {
    if (cone_subset(seh->dmin, seh->dmax,
        svi->edge->dmin, svi->edge->dmax)) {
      seh->pruned = true;
      sehi = front_to_back.erase(sehi);
    } else
      ++sehi;
        }
      }
      break;
    } else if (cone_subset(svi->edge->dmin, svi->edge->dmax, 
         seh->dmin, seh->dmax)) {
      // This edge is dominated by the previous edge.  Do not
      // insert it in the list.
      svi->edge->pruned = true;
      break;
    }
  }
      } else {
  // The current vertex is a maximum vertex.  Remove its scan
  // edge from the list.
  if (!(svi->edge->pruned))
    front_to_back.erase(svi->edge->it);
      }
      // If the edge in front has changed, update the current event.
      if (front != front_to_back.front()) {
  // Check to see if there is an event in progress that must be
  // reported.
  if (curEvent.edge.valid() && (curEvent.dmin != dir)) {
    // Report the current event.
    curEvent.dmax = dir;
    curEvent.pmax = intersection(vertex, dir, front);
    out = curEvent;
    ++out;
  }
  // Update the current event.
  if (front_to_back.front().valid()) {
    curEvent.edge = front_to_back.front()->edge;
    curEvent.dmin = dir;
    curEvent.pmin = intersection(vertex, dir, front_to_back.front());
  } else
    curEvent.edge = Coloring::IntEdgeHandle();
      }
    }
  }

  template <class OutputIterator>
  void coneTrace(const Coloring& coloring,
     const Geometry::Point& vertex,
     const Geometry::Direction& min,
     const Geometry::Direction& max,
     const Geometry::Kernel::FT& radius,
     OutputIterator out) {
    Geometry::Triangle outer = outerApprox(vertex, min, max, radius);
    coneTrace(coloring, vertex, min, max, radius, outer, out);
  }

} // End of namespace: Arak

#endif

---