geometry.hpp

Go to the documentation of this file.

#ifndef _GEOMETRY_HPP
#define _GEOMETRY_HPP

#include <algorithm>
#include "global.hpp"

/* A crucial decision when implementing algorithms that perform
 * geometric reasoning is what type of numeric representation to use.
 * Standard floating point has errors that often lead to
 * inconsistencies and bugs.  The CGAL library offers many choices,
 * but the crucial decisions seem to be whether the following are
 * supported:
 *
 *   - exact constructions (e.g., construct a point on this line)
 *   - exact predicates (e.g., is this point on that line?)
 * 
 * The code base has been written so that most portions rely on
 * neither of these.  Thus, we can use floating point representations
 * without (evident) problems.  The notable exceptions are
 *
 *   - sanity checking code makes extensive use of exact predicates, 
 *     and some isolated sanity checking code uses requires 
 *     constructions
 *   - using a triangulation to check for edge crossings
 *     (coloring.cpp) requires exact predicates
 *
 * Exact constructions may be used for debugging purposes, but they
 * lead to very slow programs.  By default, we use neither.  If an
 * inexplicable bug arises, the first action should be to turn on
 * exact geometric predicates.
 */

// #define EXACT_GEOM_PRED // Use exact predicates (~20% slower)
// #define EXACT_GEOM_CXN  // Use exact constructions (very slow)

/*
 * This file serves as an interface to the CGAL geometry library.
 */

#include <CGAL/basic.h>
#include <CGAL/Cartesian.h>
#include <CGAL/Simple_cartesian.h>
#include <CGAL/Gmpq.h>
#include <CGAL/Homogeneous.h>
#include <CGAL/Simple_homogeneous.h>
#include <CGAL/Filtered_kernel.h>
#include <CGAL/Iso_rectangle_2_Segment_2_intersection.h>
#include <CGAL/Segment_2_Segment_2_intersection.h>
#include <CGAL/Triangle_2.h>
#include <CGAL/Polygon_2.h>
#include <CGAL/Bbox_2.h>

namespace Arak {
  namespace Geometry {

    // This is a rational number represented by two arbitrary size
    // integers.  Computations are exact, but extremely slow.
    typedef CGAL::Cartesian<CGAL::Gmpq>          ExactCxnKernel;

    // The kernel below is inexact in its predicates and
    // constructions.
    typedef CGAL::Simple_cartesian<double>       InexactKernel;

    // This is a wrapper around a double representation that guarantees
    // exact predicates, but not exact constructions.  Exact predicates
    // are used in many parts of the code.
    typedef CGAL::Filtered_kernel<InexactKernel> ExactPredKernel;

#ifdef EXACT_GEOM_CXN
    #define EXACT_GEOM_PRED
    typedef ExactCxnKernel                       Kernel;
#elif defined(EXACT_GEOM_PRED)
    typedef ExactPredKernel                      Kernel;
#else
    typedef InexactKernel                        Kernel;
#endif
    typedef Kernel::Point_2                      Point;
    typedef Kernel::Line_2                       Line;
    typedef Kernel::Segment_2                    Segment;
    typedef Kernel::Direction_2                  Direction;
    typedef Kernel::Vector_2                     Vector;
    typedef Kernel::Ray_2                        Ray;
    typedef Kernel::Iso_rectangle_2              Rectangle;
    typedef Kernel::Triangle_2                   Triangle;
    typedef CGAL::Polygon_2<Kernel>                    Polygon;

    inline bool crosses(const Segment& s, const Point& a, const Point& b) {
      if ((s.source() == a) xor (s.target() == b)) 
  return false;
      else if ((s.source() == b) xor (s.target() == a)) 
  return false;
      else {
  // First do a quick bounding box check?
  Segment r(a, b);
  return CGAL::do_intersect(s, r);
      }
    }

    inline Rectangle boundingRectangle(const Point& u, 
               const Point& v, 
               const Point& w) {
      const Kernel::FT& xmin = 
  std::min(std::min(u.x(), v.x()), w.x());
      const Kernel::FT& xmax = 
  std::max(std::max(u.x(), v.x()), w.x());
      const Kernel::FT& ymin = 
  std::min(std::min(u.y(), v.y()), w.y());
      const Kernel::FT& ymax = 
  std::max(std::max(u.y(), v.y()), w.y());
      return Rectangle(xmin, ymin, xmax, ymax);
    }
    
    inline Rectangle boundingRectangle(const Point& t,
               const Point& u, 
               const Point& v, 
               const Point& w) {
      const Kernel::FT& xmin = 
  std::min(std::min(std::min(t.x(), u.x()), v.x()), w.x());
      const Kernel::FT& xmax = 
  std::max(std::max(std::max(t.x(), u.x()), v.x()), w.x());
      const Kernel::FT& ymin = 
  std::min(std::min(std::min(t.y(), u.y()), v.y()), w.y());
      const Kernel::FT& ymax = 
  std::max(std::max(std::max(t.y(), u.y()), v.y()), w.y());
      return Rectangle(xmin, ymin, xmax, ymax);
    }

    inline double logSineB(const Point& a, 
         const Point& b, 
         const Point& c) {
      Kernel::FT area = CGAL::abs(CGAL::area(a, b, c));
      Kernel::FT ab_len_sq = CGAL::squared_distance(a, b);
      Kernel::FT bc_len_sq = CGAL::squared_distance(b, c);
      return (ln(2.0 * CGAL::to_double(area))
        - ln(CGAL::to_double(ab_len_sq * bc_len_sq)) / 2.0);
    }

    inline Geometry::Point project(const Geometry::Point& p,
           const Geometry::Rectangle& r) {
      Geometry::Kernel::FT x = p.x();
      Geometry::Kernel::FT y = p.y();
      x = std::max<Geometry::Kernel::FT>(x, r.xmin());
      x = std::min<Geometry::Kernel::FT>(x, r.xmax());
      y = std::max<Geometry::Kernel::FT>(y, r.ymin());
      y = std::min<Geometry::Kernel::FT>(y, r.ymax());
      return Geometry::Point(x, y);
    }
    
    void n_gon(int n, double radius, Polygon& poly);

  }
}

#endif

---