[Boost-commit] svn:boost r80263 - in trunk: boost/polygon boost/polygon/detail libs/polygon/example libs/polygon/test

Subject: [Boost-commit] svn:boost r80263 - in trunk: boost/polygon boost/polygon/detail libs/polygon/example libs/polygon/test
From: sydorchuk.andriy_at_[hidden]
Date: 2012-08-27 14:49:54

Author: asydorchuk
Date: 2012-08-27 14:49:52 EDT (Mon, 27 Aug 2012)
New Revision: 80263
URL: http://svn.boost.org/trac/boost/changeset/80263

Log:
Polygon: refactoring & redesigning voronoi diagram structure; removing voronoi utils from the main package.
Added:
   trunk/libs/polygon/example/voronoi_visual_utils.hpp   (contents, props changed)
Removed:
   trunk/boost/polygon/voronoi_utils.hpp
Text files modified: 
   trunk/boost/polygon/detail/voronoi_structures.hpp   |    88 ++++++-                                 
   trunk/boost/polygon/voronoi_builder.hpp             |    22 +                                       
   trunk/boost/polygon/voronoi_diagram.hpp             |   411 +++++++++++++++++++++----------------   
   trunk/libs/polygon/example/voronoi_visualizer.cpp   |   427 +++++++++++++++++++++++++++------------ 
   trunk/libs/polygon/test/voronoi_builder_test.cpp    |   286 +++++++++++++-------------              
   trunk/libs/polygon/test/voronoi_structures_test.cpp |    51 +++-                                    
   trunk/libs/polygon/test/voronoi_test_helper.hpp     |   125 ++++------                              
   7 files changed, 836 insertions(+), 574 deletions(-)

Modified: trunk/boost/polygon/detail/voronoi_structures.hpp
==============================================================================
--- trunk/boost/polygon/detail/voronoi_structures.hpp	(original)
+++ trunk/boost/polygon/detail/voronoi_structures.hpp	2012-08-27 14:49:52 EDT (Mon, 27 Aug 2012)
@@ -61,6 +61,34 @@
   coordinate_type y_;
 };
 
+// Represents topology type of the voronoi site.
+enum GeometryCategory {
+  GEOMETRY_CATEGORY_POINT = 0x0,
+  GEOMETRY_CATEGORY_SEGMENT = 0x1,
+};
+
+// Represents category of the input source that forms Voronoi cell.
+enum SourceCategory {
+  // Point subtypes.
+  SOURCE_CATEGORY_SINGLE_POINT = 0x0,
+  SOURCE_CATEGORY_SEGMENT_START_POINT = 0x1,
+  SOURCE_CATEGORY_SEGMENT_END_POINT = 0x2,
+
+  // Segment subtypes.
+  SOURCE_CATEGORY_INITIAL_SEGMENT = 0x8,
+  SOURCE_CATEGORY_REVERSE_SEGMENT = 0x9,
+
+  SOURCE_CATEGORY_GEOMETRY_SHIFT = 0x3,
+  SOURCE_CATEGORY_BITMASK = 0x1F,
+};
+
+bool belongs(
+    const SourceCategory& source_category,
+    const GeometryCategory& geometry_category) {
+  return (static_cast<std::size_t>(source_category) >> SOURCE_CATEGORY_GEOMETRY_SHIFT) ==
+      static_cast<std::size_t>(geometry_category);
+}
+
 // Site event type.
 // Occurs when the sweepline sweeps over one of the initial sites:
 //   1) point site
@@ -91,39 +119,42 @@
   site_event() :
       point0_(0, 0),
       point1_(0, 0),
-      sorted_index_(0) {}
+      sorted_index_(0),
+      flags_(0) {}
 
   site_event(coordinate_type x, coordinate_type y) :
       point0_(x, y),
       point1_(x, y),
-      sorted_index_(0) {}
+      sorted_index_(0),
+      flags_(0) {}
 
   site_event(const point_type &point) :
       point0_(point),
       point1_(point),
-      sorted_index_(0) {}
+      sorted_index_(0),
+      flags_(0) {}
 
   site_event(coordinate_type x1, coordinate_type y1,
              coordinate_type x2, coordinate_type y2):
       point0_(x1, y1),
       point1_(x2, y2),
-      sorted_index_(0) {}
+      sorted_index_(0),
+      flags_(0) {}
 
   site_event(const point_type &point1, const point_type &point2) :
       point0_(point1),
       point1_(point2),
-      sorted_index_(0) {}
+      sorted_index_(0),
+      flags_(0) {}
 
   bool operator==(const site_event &that) const {
     return (this->point0_ == that.point0_) &&
-           (this->point1_ == that.point1_) &&
-           (this->sorted_index_ == that.sorted_index_);
+           (this->point1_ == that.point1_);
   }
 
   bool operator!=(const site_event &that) const {
     return (this->point0_ != that.point0_) ||
-           (this->point1_ != that.point1_) ||
-           (this->sorted_index_ != that.sorted_index_);
+           (this->point1_ != that.point1_);
   }
 
   coordinate_type x(bool oriented = false) const {
@@ -170,46 +201,61 @@
     return is_inverse() ? point0_ : point1_;
   }
 
+  std::size_t sorted_index() const {
+    return sorted_index_;
+  }
+
   site_event& sorted_index(std::size_t index) {
-    sorted_index_ = (index << 1) + (sorted_index_ & 1);
+    sorted_index_ = index;
     return *this;
   }
 
+  std::size_t initial_index() const {
+    return initial_index_;
+  }
+
   site_event& initial_index(std::size_t index) {
     initial_index_ = index;
     return *this;
   }
 
+  bool is_inverse() const {
+    return (flags_ & IS_INVERSE) ? true : false;
+  }
+
   site_event& inverse() {
-    sorted_index_ ^= 1;
+    flags_ ^= IS_INVERSE;
     return *this;
   }
 
-  std::size_t sorted_index() const {
-    return sorted_index_ >> 1;
+  SourceCategory source_category() const {
+    return static_cast<SourceCategory>(flags_ & SOURCE_CATEGORY_BITMASK);
   }
 
-  std::size_t initial_index() const {
-    return initial_index_;
+  site_event& source_category(SourceCategory source_category) {
+    flags_ |= source_category;
+    return *this;
   }
 
   bool is_point() const {
-    return point0_.x() == point1_.x() && point0_.y() == point1_.y();
+    return (point0_.x() == point1_.x()) && (point0_.y() == point1_.y());
   }
 
   bool is_segment() const {
-    return point0_.x() != point1_.x() || point0_.y() != point1_.y();
-  }
-
-  bool is_inverse() const {
-    return static_cast<bool>(sorted_index_ & 1);
+    return (point0_.x() != point1_.x()) || (point0_.y() != point1_.y());
   }
 
 private:
+  enum Bits {
+    SOURCE_CATEGORY_BITMASK = SOURCE_CATEGORY_BITMASK,
+    IS_INVERSE = 0x20,  // 32
+  };
+
   point_type point0_;
   point_type point1_;
   std::size_t sorted_index_;
   std::size_t initial_index_;
+  std::size_t flags_;
 };
 
 // Circle event type.

Modified: trunk/boost/polygon/voronoi_builder.hpp
==============================================================================
--- trunk/boost/polygon/voronoi_builder.hpp	(original)
+++ trunk/boost/polygon/voronoi_builder.hpp	2012-08-27 14:49:52 EDT (Mon, 27 Aug 2012)
@@ -51,6 +51,7 @@
   std::size_t insert_point(const int_type& x, const int_type& y) {
     site_events_.push_back(site_event_type(x, y));
     site_events_.back().initial_index(index_);
+    site_events_.back().source_category(detail::SOURCE_CATEGORY_SINGLE_POINT);
     return index_++;
   }
 
@@ -61,16 +62,29 @@
   std::size_t insert_segment(
       const int_type& x1, const int_type& y1,
       const int_type& x2, const int_type& y2) {
+    // Set up start point site.
     point_type p1(x1, y1);
-    point_type p2(x2, y2);
     site_events_.push_back(site_event_type(p1));
     site_events_.back().initial_index(index_);
+    site_events_.back().source_category(
+        detail::SOURCE_CATEGORY_SEGMENT_START_POINT);
+
+    // Set up end point site.
+    point_type p2(x2, y2);
     site_events_.push_back(site_event_type(p2));
     site_events_.back().initial_index(index_);
+    site_events_.back().source_category(
+        detail::SOURCE_CATEGORY_SEGMENT_END_POINT);
+
+    // Set up segment site.
     if (point_comparison_(p1, p2)) {
       site_events_.push_back(site_event_type(p1, p2));
+      site_events_.back().source_category(
+          detail::SOURCE_CATEGORY_INITIAL_SEGMENT);
     } else {
       site_events_.push_back(site_event_type(p2, p1));
+      site_events_.back().source_category(
+          detail::SOURCE_CATEGORY_REVERSE_SEGMENT);
     }
     site_events_.back().initial_index(index_);
     return index_++;
@@ -93,8 +107,7 @@
       } else if (site_event_iterator_ == site_events_.end()) {
         process_circle_event(output);
       } else {
-        if (event_comparison(*site_event_iterator_,
-                             circle_events_.top().first)) {
+        if (event_comparison(*site_event_iterator_, circle_events_.top().first)) {
           process_site_event(output);
         } else {
           process_circle_event(output);
@@ -114,7 +127,6 @@
   void clear() {
     index_ = 0;
     site_events_.clear();
-    beach_line_.clear();
   }
 
 private:
@@ -169,8 +181,6 @@
 
   template <typename OUTPUT>
   void init_beach_line(OUTPUT *output) {
-    if (!beach_line_.empty())
-      beach_line_.clear();
     if (site_events_.empty())
       return;
     if (site_events_.size() == 1) {

Modified: trunk/boost/polygon/voronoi_diagram.hpp
==============================================================================
--- trunk/boost/polygon/voronoi_diagram.hpp	(original)
+++ trunk/boost/polygon/voronoi_diagram.hpp	2012-08-27 14:49:52 EDT (Mon, 27 Aug 2012)
@@ -24,96 +24,118 @@
 
 // Represents Voronoi cell.
 // Data members:
-//   1) pointer to the incident edge
-//   2) site inside cell
+//   1) index of the source within the initial input set
+//   2) pointer to the incident edge
 //   3) mutable color member
 // Cell may contain point or segment site inside.
 template <typename T>
 class voronoi_cell {
 public:
   typedef T coordinate_type;
-  typedef detail::point_2d<coordinate_type> point_type;
   typedef std::size_t color_type;
   typedef voronoi_edge<coordinate_type> voronoi_edge_type;
+  typedef std::size_t source_index_type;
+  typedef detail::SourceCategory source_category_type;
 
-  voronoi_cell(const point_type &p1, voronoi_edge_type *edge) :
-      point0_(p1),
-      point1_(p1),
-      color_(0),
-      incident_edge_(edge) {}
-
-  voronoi_cell(const point_type &p1,
-               const point_type &p2,
-               voronoi_edge_type *edge) :
-      point0_(p1),
-      point1_(p2),
-      color_(0),
-      incident_edge_(edge) {}
+  voronoi_cell(const source_index_type& source_index,
+               const source_category_type& source_category,
+               voronoi_edge_type* edge) :
+      source_index_(source_index),
+      incident_edge_(edge),
+      color_(source_category) {}
 
   // Returns true if the cell contains point site, false else.
-  bool contains_point() const { return point0_ == point1_; }
+  bool contains_point() const {
+    source_category_type source_category = this->source_category();
+    return detail::belongs(source_category, detail::GEOMETRY_CATEGORY_POINT);
+  }
 
   // Returns true if the cell contains segment site, false else.
-  bool contains_segment() const { return point0_ != point1_; }
+  bool contains_segment() const {
+    source_category_type source_category = this->source_category();
+    return detail::belongs(source_category, detail::GEOMETRY_CATEGORY_SEGMENT);
+  }
+
+  source_index_type source_index() const {
+    return source_index_;
+  }
+
+  source_category_type source_category() const {
+    return static_cast<source_category_type>(
+        color_ & detail::SOURCE_CATEGORY_BITMASK);
+  }
 
   // Degenerate cells don't have any incident edges.
   bool is_degenerate() const { return incident_edge_ == NULL; }
 
-  // Returns site point in case cell contains point site,
-  // the first endpoint of the segment site else.
-  const point_type &point0() const { return point0_; }
-
-  // Returns site point in case cell contains point site,
-  // the second endpoint of the segment site else.
-  const point_type &point1() const { return point1_; }
-
-  color_type color() const { return color_; }
-  void color(const color_type& color) const { color_ = color; }
-
-  voronoi_edge_type *incident_edge() { return incident_edge_; }
-  const voronoi_edge_type *incident_edge() const { return incident_edge_; }
-  void incident_edge(voronoi_edge_type *e) { incident_edge_ = e; }
+  voronoi_edge_type* incident_edge() { return incident_edge_; }
+  const voronoi_edge_type* incident_edge() const { return incident_edge_; }
+  void incident_edge(voronoi_edge_type* e) { incident_edge_ = e; }
+
+  color_type color() const { return color_ >> BITS_SHIFT; }
+  void color(const color_type& color) const {
+    color_ &= BITS_MASK;
+    color_ |= color << BITS_SHIFT;
+  }
 
 private:
-  point_type point0_;
-  point_type point1_;
+  // 5 color bits are reserved.
+  enum Bits {
+    BITS_SHIFT = 0x5,
+    BITS_MASK = 0x1F,
+  };
+
+  source_index_type source_index_;
+  voronoi_edge_type* incident_edge_;
   mutable color_type color_;
-  voronoi_edge_type *incident_edge_;
 };
 
 // Represents Voronoi vertex.
 // Data members:
-//   1) vertex point itself
+//   1) vertex coordinates
 //   2) pointer to the incident edge
 //   3) mutable color member
 template <typename T>
 class voronoi_vertex {
 public:
   typedef T coordinate_type;
-  typedef detail::point_2d<T> point_type;
   typedef std::size_t color_type;
   typedef voronoi_edge<coordinate_type> voronoi_edge_type;
 
-  voronoi_vertex(const point_type &vertex, voronoi_edge_type *edge) :
-      vertex_(vertex),
-      color_(0),
-      incident_edge_(edge) {}
+  voronoi_vertex(const coordinate_type& x,
+                 const coordinate_type& y,
+                 voronoi_edge_type* edge) :
+      x_(x),
+      y_(y),
+      incident_edge_(edge),
+      color_(0) {}
 
-  const point_type &vertex() const { return vertex_; }
+  const coordinate_type& x() const { return x_; }
+  const coordinate_type& y() const { return y_; }
 
   bool is_degenerate() const { return incident_edge_ == NULL; }
 
-  color_type color() const { return color_; }
-  void color(const color_type& color) const { color_ = color; }
-
-  voronoi_edge_type *incident_edge() { return incident_edge_; }
-  const voronoi_edge_type *incident_edge() const { return incident_edge_; }
-  void incident_edge(voronoi_edge_type *e) { incident_edge_ = e; }
+  voronoi_edge_type* incident_edge() { return incident_edge_; }
+  const voronoi_edge_type* incident_edge() const { return incident_edge_; }
+  void incident_edge(voronoi_edge_type* e) { incident_edge_ = e; }
+
+  color_type color() const { return color_ >> BITS_SHIFT; }
+  void color(const color_type& color) const {
+    color_ &= BITS_MASK;
+    color_ |= color << BITS_SHIFT;
+  }
 
 private:
-  point_type vertex_;
+  // 5 color bits are reserved.
+  enum Bits {
+    BITS_SHIFT = 0x5,
+    BITS_MASK = 0x1F,
+  };
+
+  coordinate_type x_;
+  coordinate_type y_;
+  voronoi_edge_type* incident_edge_;
   mutable color_type color_;
-  voronoi_edge_type *incident_edge_;
 };
 
 // Half-edge data structure. Represents Voronoi edge.
@@ -134,53 +156,58 @@
   typedef voronoi_edge<coordinate_type> voronoi_edge_type;
   typedef std::size_t color_type;
 
-  voronoi_edge() :
+  voronoi_edge(bool is_linear, bool is_primary) :
       cell_(NULL),
       vertex_(NULL),
       twin_(NULL),
       next_(NULL),
       prev_(NULL),
-      color_(0) {}
-
-  voronoi_cell_type *cell() { return cell_; }
-  const voronoi_cell_type *cell() const { return cell_; }
-  void cell(voronoi_cell_type *c) { cell_ = c; }
-
-  voronoi_vertex_type *vertex0() { return vertex_; }
-  const voronoi_vertex_type *vertex0() const { return vertex_; }
-  void vertex0(voronoi_vertex_type *v) { vertex_ = v; }
-
-  voronoi_vertex_type *vertex1() { return twin_->vertex0(); }
-  const voronoi_vertex_type *vertex1() const { return twin_->vertex0(); }
-  void vertex1(voronoi_vertex_type *v) { twin_->vertex0(v); }
-
-  voronoi_edge_type *twin() { return twin_; }
-  const voronoi_edge_type *twin() const { return twin_; }
-  void twin(voronoi_edge_type *e) { twin_ = e; }
-
-  voronoi_edge_type *next() { return next_; }
-  const voronoi_edge_type *next() const { return next_; }
-  void next(voronoi_edge_type *e) { next_ = e; }
-
-  voronoi_edge_type *prev() { return prev_; }
-  const voronoi_edge_type *prev() const { return prev_; }
-  void prev(voronoi_edge_type *e) { prev_ = e; }
+      color_(0) {
+    if (is_linear)
+      color_ |= BIT_IS_LINEAR;
+    if (is_primary)
+      color_ |= BIT_IS_PRIMARY;
+  }
+
+  voronoi_cell_type* cell() { return cell_; }
+  const voronoi_cell_type* cell() const { return cell_; }
+  void cell(voronoi_cell_type* c) { cell_ = c; }
+
+  voronoi_vertex_type* vertex0() { return vertex_; }
+  const voronoi_vertex_type* vertex0() const { return vertex_; }
+  void vertex0(voronoi_vertex_type* v) { vertex_ = v; }
+
+  voronoi_vertex_type* vertex1() { return twin_->vertex0(); }
+  const voronoi_vertex_type* vertex1() const { return twin_->vertex0(); }
+  void vertex1(voronoi_vertex_type* v) { twin_->vertex0(v); }
+
+  voronoi_edge_type* twin() { return twin_; }
+  const voronoi_edge_type* twin() const { return twin_; }
+  void twin(voronoi_edge_type* e) { twin_ = e; }
+
+  voronoi_edge_type* next() { return next_; }
+  const voronoi_edge_type* next() const { return next_; }
+  void next(voronoi_edge_type* e) { next_ = e; }
+
+  voronoi_edge_type* prev() { return prev_; }
+  const voronoi_edge_type* prev() const { return prev_; }
+  void prev(voronoi_edge_type* e) { prev_ = e; }
 
   // Returns a pointer to the rotation next edge
   // over the starting point of the half-edge.
-  voronoi_edge_type *rot_next() {
+  voronoi_edge_type* rot_next() {
     return (vertex_) ? prev_->twin() : NULL;
   }
-  const voronoi_edge_type *rot_next() const {
+  const voronoi_edge_type* rot_next() const {
     return (vertex_) ? prev_->twin() : NULL;
   }
 
   // Returns a pointer to the rotation prev edge
   // over the starting point of the half-edge.
-  voronoi_edge_type *rot_prev() {
+  voronoi_edge_type* rot_prev() {
     return (vertex_) ? twin_->next() : NULL;
   }
-  const voronoi_edge_type *rot_prev() const {
+  const voronoi_edge_type* rot_prev() const {
     return (vertex_) ? twin_->next() : NULL;
   }
 
@@ -188,67 +215,68 @@
   // Returns false if the edge is infinite (ray, line).
   bool is_finite() const { return vertex0() && vertex1(); }
 
+  // Returns true if the edge is infinite (ray, line).
+  // Returns false if the edge is finite (segment, parabolic arc).
+  bool is_infinite() const { return !vertex0() || !vertex1(); }
+
   // Returns true if the edge is linear (segment, ray, line).
   // Returns false if the edge is curved (parabolic arc).
   bool is_linear() const {
-    if (!is_primary())
-      return true;
-    return !(cell()->contains_segment() ^ twin()->cell()->contains_segment());
+    return (color_ & BIT_IS_LINEAR) ? true : false;
   }
 
   // Returns true if the edge is curved (parabolic arc).
   // Returns false if the edge is linear (segment, ray, line).
   bool is_curved() const {
-    if (!is_primary())
-      return false;
-    return (cell()->contains_segment() ^ twin()->cell()->contains_segment());
+    return (color_ & BIT_IS_LINEAR) ? false : true;
   }
 
   // Returns false if edge goes through the endpoint of the segment.
   // Returns true else.
   bool is_primary() const {
-    bool flag1 = cell_->contains_segment();
-    bool flag2 = twin_->cell()->contains_segment();
-    if (flag1 && !flag2) {
-      return cell_->point0() != twin_->cell()->point0() &&
-             cell_->point1() != twin_->cell()->point0();
-    }
-    if (!flag1 && flag2) {
-      return twin_->cell()->point0() != cell_->point0() &&
-             twin_->cell()->point1() != cell_->point0();
-    }
-    return true;
+    return (color_ & BIT_IS_PRIMARY) ? true : false;
+  }
+
+  // Returns true if edge goes through the endpoint of the segment.
+  // Returns false else.
+  bool is_secondary() const {
+    return (color_ & BIT_IS_PRIMARY) ? false : true;
   }
 
-  color_type color() const { return color_; }
-  void color(const color_type& color) const { color_ = color; }
+  color_type color() const { return color_ >> BITS_SHIFT; }
+  void color(const color_type& color) const {
+    color_ &= BITS_MASK;
+    color_ |= color << BITS_SHIFT;
+  }
 
 private:
-  voronoi_cell_type *cell_;
-  voronoi_vertex_type *vertex_;
-  voronoi_edge_type *twin_;
-  voronoi_edge_type *next_;
-  voronoi_edge_type *prev_;
+  // 5 color bits are reserved.
+  enum Bits {
+    BIT_IS_LINEAR = 0x1,  // linear is opposite to curved
+    BIT_IS_PRIMARY = 0x2,  // primary is opposite to secondary
+
+    BITS_SHIFT = 0x5,
+    BITS_MASK = 0x1F,
+  };
+
+  voronoi_cell_type* cell_;
+  voronoi_vertex_type* vertex_;
+  voronoi_edge_type* twin_;
+  voronoi_edge_type* next_;
+  voronoi_edge_type* prev_;
   mutable color_type color_;
 };
 
 template <typename T>
 struct voronoi_diagram_traits {
   typedef T coordinate_type;
-  typedef struct {
-    template <typename CT>
-    coordinate_type operator()(const CT& value) {
-      return static_cast<coordinate_type>(value);
-    }
-  } ctype_converter_type;
-  typedef detail::point_2d<coordinate_type> point_type;
   typedef voronoi_cell<coordinate_type> cell_type;
   typedef voronoi_vertex<coordinate_type> vertex_type;
   typedef voronoi_edge<coordinate_type> edge_type;
   typedef class {
   public:
     enum { ULPS = 128 };
-    bool operator()(const point_type &v1, const point_type &v2) const {
+    bool operator()(const vertex_type& v1, const vertex_type& v2) const {
       return (ulp_cmp(v1.x(), v2.x(), ULPS) ==
               detail::ulp_comparison<T>::EQUAL) &&
              (ulp_cmp(v1.y(), v2.y(), ULPS) ==
@@ -265,7 +293,6 @@
 class voronoi_diagram {
 public:
   typedef typename TRAITS::coordinate_type coordinate_type;
-  typedef typename TRAITS::point_type point_type;
   typedef typename TRAITS::cell_type cell_type;
   typedef typename TRAITS::vertex_type vertex_type;
   typedef typename TRAITS::edge_type edge_type;
@@ -286,7 +313,7 @@
   // construct the Voronoi diagram.
   class voronoi_diagram_builder {
   public:
-    void vd(voronoi_diagram *vd) {
+    void vd(voronoi_diagram* vd) {
       vd_ = vd;
     }
 
@@ -299,20 +326,20 @@
     }
 
     template <typename SEvent>
-    void process_single_site(const SEvent &site) {
+    void process_single_site(const SEvent& site) {
       vd_->process_single_site(site);
     }
 
     template <typename SEvent>
     std::pair<void*, void*> insert_new_edge(
-        const SEvent &site1, const SEvent &site2) {
+        const SEvent& site1, const SEvent& site2) {
       return vd_->insert_new_edge(site1, site2);
     }
 
     template <typename SEvent, typename CEvent>
-    std::pair<void *, void *> insert_new_edge(
-        const SEvent &site1, const SEvent &site3, const CEvent &circle,
-        void *data12, void *data23) {
+    std::pair<void*, void*> insert_new_edge(
+        const SEvent& site1, const SEvent& site3, const CEvent& circle,
+        void* data12, void* data23) {
       return vd_->insert_new_edge(site1, site3, circle, data12, data23);
     }
 
@@ -322,7 +349,7 @@
     }
 
   private:
-    voronoi_diagram *vd_;
+    voronoi_diagram* vd_;
   };
 
   voronoi_diagram() {
@@ -336,31 +363,35 @@
     builder_.vd(&(*this));
   }
 
-  const cell_container_type &cells() const {
+  const cell_container_type& cells() const {
     return cells_;
   }
 
-  const vertex_container_type &vertices() const {
+  const vertex_container_type& vertices() const {
     return vertices_;
   }
 
-  const edge_container_type &edges() const {
+  const edge_container_type& edges() const {
     return edges_;
   }
 
-  unsigned int num_cells() const {
+  std::size_t num_cells() const {
     return cells_.size();
   }
 
-  unsigned int num_edges() const {
+  std::size_t num_edges() const {
     return edges_.size() >> 1;
   }
 
-  unsigned int num_vertices() const {
+  std::size_t num_half_edges() const {
+    return edges_.size();
+  }
+
+  std::size_t num_vertices() const {
     return vertices_.size();
   }
 
-  voronoi_diagram_builder *builder() {
+  voronoi_diagram_builder* builder() {
     if (builder_.done()) {
       return NULL;
     } else {
@@ -369,7 +400,6 @@
   }
 
 private:
-  typedef typename TRAITS::ctype_converter_type ctype_converter_type;
   typedef typename TRAITS::vertex_equality_predicate_type
     vertex_equality_predicate_type;
 
@@ -381,12 +411,33 @@
     edges_.reserve((num_sites << 2) + (num_sites << 1));
   }
 
-  // Update the Voronoi output in case of a single point input.
   template <typename SEvent>
-  void process_single_site(const SEvent &site) {
-    // Update cell records.
-    point_type p = prepare_point(site.point0());
-    cells_.push_back(cell_type(p, NULL));
+  bool is_primary_edge(const SEvent& site1, const SEvent& site2) const {
+    bool flag1 = site1.is_segment();
+    bool flag2 = site2.is_segment();
+    if (flag1 && !flag2) {
+      return (site1.point0() != site2.point0()) &&
+             (site1.point1() != site2.point0());
+    }
+    if (!flag1 && flag2) {
+      return (site2.point0() != site1.point0()) &&
+             (site2.point1() != site1.point0());
+    }
+    return true;
+  }
+
+  template <typename SEvent>
+  bool is_linear_edge(const SEvent& site1, const SEvent& site2) const {
+    if (!is_primary_edge(site1, site2)) {
+      return true;
+    }
+    return !(site1.is_segment() ^ site2.is_segment());
+  }
+
+  template <typename SEvent>
+  void process_single_site(const SEvent& site) {
+    cells_.push_back(cell_type(
+         site.initial_index(), site.source_category(), NULL));
   }
 
   // Insert a new half-edge into the output data structure.
@@ -394,29 +445,32 @@
   // Returns a pair of pointers to a new half-edges.
   template <typename SEvent>
   std::pair<void*, void*> insert_new_edge(
-      const SEvent &site1, const SEvent &site2) {
+      const SEvent& site1, const SEvent& site2) {
     // Get sites' indexes.
     int site_index1 = site1.sorted_index();
     int site_index2 = site2.sorted_index();
 
+    bool is_linear = is_linear_edge(site1, site2);
+    bool is_primary = is_primary_edge(site1, site2);
+
     // Create a new half-edge that belongs to the first site.
-    edges_.push_back(edge_type());
-    edge_type &edge1 = edges_.back();
+    edges_.push_back(edge_type(is_linear, is_primary));
+    edge_type& edge1 = edges_.back();
 
     // Create a new half-edge that belongs to the second site.
-    edges_.push_back(edge_type());
-    edge_type &edge2 = edges_.back();
+    edges_.push_back(edge_type(is_linear, is_primary));
+    edge_type& edge2 = edges_.back();
 
     // Add the initial cell during the first edge insertion.
     if (cells_.empty()) {
-      process_single_site(site1);
+      cells_.push_back(cell_type(
+          site1.initial_index(), site1.source_category(), NULL));
     }
 
     // The second site represents a new site during site event
     // processing. Add a new cell to the cell records.
-    cells_.push_back(cell_type(prepare_point(site2.point0()),
-                               prepare_point(site2.point1()),
-                               NULL));
+    cells_.push_back(cell_type(
+      site2.initial_index(), site2.source_category(), NULL));
 
     // Set up pointers to cells.
     edge1.cell(&cells_[site_index1]);
@@ -437,28 +491,31 @@
   // pointers to those half-edges. Half-edges' direction goes out of the
   // new Voronoi vertex point. Returns a pair of pointers to a new half-edges.
   template <typename SEvent, typename CEvent>
-  std::pair<void *, void *> insert_new_edge(
-      const SEvent &site1, const SEvent &site3, const CEvent &circle,
-      void *data12, void *data23) {
-    edge_type *edge12 = static_cast<edge_type*>(data12);
-    edge_type *edge23 = static_cast<edge_type*>(data23);
+  std::pair<void*, void*> insert_new_edge(
+      const SEvent& site1, const SEvent& site3, const CEvent& circle,
+      void* data12, void* data23) {
+    edge_type* edge12 = static_cast<edge_type*>(data12);
+    edge_type* edge23 = static_cast<edge_type*>(data23);
 
     // Add a new Voronoi vertex.
-    vertices_.push_back(vertex_type(prepare_point(circle), NULL));
-    vertex_type &new_vertex = vertices_.back();
+    vertices_.push_back(vertex_type(circle.x(), circle.y(), NULL));
+    vertex_type& new_vertex = vertices_.back();
 
     // Update vertex pointers of the old edges.
     edge12->vertex0(&new_vertex);
     edge23->vertex0(&new_vertex);
 
+    bool is_linear = is_linear_edge(site1, site3);
+    bool is_primary = is_primary_edge(site1, site3);
+
     // Add a new half-edge.
-    edges_.push_back(edge_type());
-    edge_type &new_edge1 = edges_.back();
+    edges_.push_back(edge_type(is_linear, is_primary));
+    edge_type& new_edge1 = edges_.back();
     new_edge1.cell(&cells_[site1.sorted_index()]);
 
     // Add a new half-edge.
-    edges_.push_back(edge_type());
-    edge_type &new_edge2 = edges_.back();
+    edges_.push_back(edge_type(is_linear, is_primary));
+    edge_type& new_edge2 = edges_.back();
     new_edge2.cell(&cells_[site3.sorted_index()]);
 
     // Update twin pointers.
@@ -484,14 +541,14 @@
     // Remove degenerate edges.
     edge_iterator last_edge = edges_.begin();
     for (edge_iterator it = edges_.begin(); it != edges_.end(); it += 2) {
-      const vertex_type *v1 = it->vertex0();
-      const vertex_type *v2 = it->vertex1();
-      if (v1 && v2 && vertex_equality_predicate_(v1->vertex(), v2->vertex())) {
+      const vertex_type* v1 = it->vertex0();
+      const vertex_type* v2 = it->vertex1();
+      if (v1 && v2 && vertex_equality_predicate_(*v1, *v2)) {
         remove_edge(&(*it));
       } else {
         if (it != last_edge) {
-          edge_type *e1 = &(*last_edge = *it);
-          edge_type *e2 = &(*(last_edge + 1) = *(it + 1));
+          edge_type* e1 = &(*last_edge = *it);
+          edge_type* e2 = &(*(last_edge + 1) = *(it + 1));
           e1->twin(e2);
           e2->twin(e1);
           if (e1->prev()) {
@@ -522,8 +579,8 @@
       if (it->incident_edge()) {
         if (it != last_vertex) {
           *last_vertex = *it;
-          vertex_type *v = &(*last_vertex);
-          edge_type *e = v->incident_edge();
+          vertex_type* v = &(*last_vertex);
+          edge_type* e = v->incident_edge();
           do {
             e->vertex0(v);
             e = e->rot_next();
@@ -539,7 +596,7 @@
       if (!edges_.empty()) {
         // Update prev/next pointers for the line edges.
         edge_iterator edge_it = edges_.begin();
-        edge_type *edge1 = &(*edge_it);
+        edge_type* edge1 = &(*edge_it);
         edge1->next(edge1);
         edge1->prev(edge1);
         ++edge_it;
@@ -547,7 +604,7 @@
         ++edge_it;
 
         while (edge_it != edges_.end()) {
-          edge_type *edge2 = &(*edge_it);
+          edge_type* edge2 = &(*edge_it);
           ++edge_it;
 
           edge1->next(edge2);
@@ -570,7 +627,7 @@
           continue;
         // Move to the previous edge while
         // it is possible in the CW direction.
-        edge_type *left_edge = cell_it->incident_edge();
+        edge_type* left_edge = cell_it->incident_edge();
         while (left_edge->prev() != NULL) {
           left_edge = left_edge->prev();
           // Terminate if this is not a boundary cell.
@@ -581,7 +638,7 @@
         if (left_edge->prev() != NULL)
           continue;
 
-        edge_type *right_edge = cell_it->incident_edge();
+        edge_type* right_edge = cell_it->incident_edge();
         while (right_edge->next() != NULL)
           right_edge = right_edge->next();
         left_edge->prev(right_edge);
@@ -590,31 +647,24 @@
     }
   }
 
-  template <typename P>
-  point_type prepare_point(const P& p) {
-    coordinate_type nx = convert_(p.x());
-    coordinate_type ny = convert_(p.y());
-    return point_type(nx, ny);
-  }
-
   // Remove degenerate edge.
-  void remove_edge(edge_type *edge) {
+  void remove_edge(edge_type* edge) {
     // Update the endpoints of the incident edges to the second vertex.
-    vertex_type *vertex = edge->vertex0();
-    edge_type *updated_edge = edge->twin()->rot_next();
+    vertex_type* vertex = edge->vertex0();
+    edge_type* updated_edge = edge->twin()->rot_next();
     while (updated_edge != edge->twin()) {
       updated_edge->vertex0(vertex);
       updated_edge = updated_edge->rot_next();
     }
 
-    edge_type *edge1 = edge;
-    edge_type *edge2 = edge->twin();
+    edge_type* edge1 = edge;
+    edge_type* edge2 = edge->twin();
 
-    edge_type *edge1_rot_prev = edge1->rot_prev();
-    edge_type *edge1_rot_next = edge1->rot_next();
+    edge_type* edge1_rot_prev = edge1->rot_prev();
+    edge_type* edge1_rot_next = edge1->rot_next();
 
-    edge_type *edge2_rot_prev = edge2->rot_prev();
-    edge_type *edge2_rot_next = edge2->rot_next();
+    edge_type* edge2_rot_prev = edge2->rot_prev();
+    edge_type* edge2_rot_next = edge2->rot_next();
 
     // Update prev/next pointers for the incident edges.
     edge1_rot_next->twin()->next(edge2_rot_prev);
@@ -628,7 +678,6 @@
   edge_container_type edges_;
 
   voronoi_diagram_builder builder_;
-  ctype_converter_type convert_;
   vertex_equality_predicate_type vertex_equality_predicate_;
 
   // Disallow copy constructor and operator=

Deleted: trunk/boost/polygon/voronoi_utils.hpp
==============================================================================
--- trunk/boost/polygon/voronoi_utils.hpp	2012-08-27 14:49:52 EDT (Mon, 27 Aug 2012)
+++ (empty file)
@@ -1,479 +0,0 @@
-// Boost.Polygon library voronoi_utils.hpp header file
-
-//          Copyright Andrii Sydorchuk 2010-2012.
-// Distributed under the Boost Software License, Version 1.0.
-//    (See accompanying file LICENSE_1_0.txt or copy at
-//          http://www.boost.org/LICENSE_1_0.txt)
-
-// See http://www.boost.org for updates, documentation, and revision history.
-
-#ifndef BOOST_POLYGON_VORONOI_UTILS
-#define BOOST_POLYGON_VORONOI_UTILS
-
-#include <cmath>
-#include <stack>
-#include <vector>
-
-#include "voronoi_diagram.hpp"
-
-namespace boost {
-namespace polygon {
-
-// Bounding rectangle data structure. Contains coordinates
-// of the bottom left and upper right corners of rectangle.
-template <typename T>
-class bounding_rectangle {
-public:
-  typedef T coordinate_type;
-
-  bounding_rectangle() : x_min_(1), x_max_(0) {}
-
-  bounding_rectangle(coordinate_type x1, coordinate_type y1,
-                     coordinate_type x2, coordinate_type y2) {
-    x_min_ = (std::min)(x1, x2);
-    y_min_ = (std::min)(y1, y2);
-    x_max_ = (std::max)(x1, x2);
-    y_max_ = (std::max)(y1, y2);
-  }
-
-  // Update bounding rectangle.
-  void update(coordinate_type x, coordinate_type y) {
-    if (x_min_ > x_max_) {
-      x_min_ = x_max_ = x;
-      y_min_ = y_max_ = y;
-    } else {
-      x_min_ = (std::min)(x_min_, x);
-      y_min_ = (std::min)(y_min_, y);
-      x_max_ = (std::max)(x_max_, x);
-      y_max_ = (std::max)(y_max_, y);
-    }
-  }
-
-  bool is_empty() const {
-    return x_min_ > x_max_;
-  }
-
-  void clear() {
-    x_min_ = 1;
-    x_max_ = 0;
-  }
-
-  bool contains(coordinate_type x, coordinate_type y) const {
-    return x > x_min_ && x < x_max_ && y > y_min_ && y < y_max_;
-  }
-
-  // Return the x-coordinate of the bottom left point of the rectangle.
-  coordinate_type x_min() const {
-    return x_min_;
-  }
-
-  // Return the y-coordinate of the bottom left point of the rectangle.
-  coordinate_type y_min() const {
-    return y_min_;
-  }
-
-  // Return the x-coordinate of the upper right point of the rectangle.
-  coordinate_type x_max() const {
-    return x_max_;
-  }
-
-  // Return the y-coordinate of the upper right point of the rectangle.
-  coordinate_type y_max() const {
-    return y_max_;
-  }
-
-private:
-  coordinate_type x_min_;
-  coordinate_type y_min_;
-  coordinate_type x_max_;
-  coordinate_type y_max_;
-};
-
-template <typename fpt>
-struct voronoi_utils_traits {
-  typedef fpt coordinate_type;
-  typedef detail::point_2d<coordinate_type> point_type;
-  typedef std::vector<point_type> point_set_type;
-  typedef bounding_rectangle<coordinate_type> brect_type;
-  typedef struct {
-    template <typename CT>
-      coordinate_type operator()(const CT& value) const {
-        return static_cast<coordinate_type>(value);
-    }
-  } ctype_converter_type;
-};
-
-// Voronoi output post-processing tools.
-template <typename T, typename TRAITS = voronoi_utils_traits<T> >
-class voronoi_utils {
-public:
-  typedef typename TRAITS::coordinate_type coordinate_type;
-  typedef typename TRAITS::point_type point_type;
-  typedef typename TRAITS::point_set_type point_set_type;
-  typedef typename TRAITS::brect_type brect_type;
-
-  // Get scaled by a factor bounding rectangle.
-  template <typename CT>
-  static brect_type scale(const bounding_rectangle<CT> &brect,
-      coordinate_type factor) {
-    coordinate_type x_len = to_fpt(brect.x_max()) - to_fpt(brect.x_min());
-    coordinate_type y_len = to_fpt(brect.y_max()) - to_fpt(brect.y_min());
-    coordinate_type x_mid = to_fpt(brect.x_max()) + to_fpt(brect.x_min());
-    coordinate_type y_mid = to_fpt(brect.y_max()) + to_fpt(brect.y_min());
-    x_mid *= to_fpt(0.5);
-    y_mid *= to_fpt(0.5);
-    coordinate_type offset = (std::max)(x_len, y_len) * factor * to_fpt(0.5);
-    brect_type new_brect(x_mid - offset, y_mid - offset,
-                         x_mid + offset, y_mid + offset);
-    return new_brect;
-  }
-
-  // Discretizes finite Voronoi edge.
-  // Discretization absolute error is defined by max_error value.
-  template <typename CT>
-  static void discretize(const voronoi_edge<CT> &edge,
-      coordinate_type max_error, point_set_type &discretization) {
-    // Don't process infinite edges.
-    if (!edge.is_finite())
-      return;
-
-    // Retrieve the cell to the left of the current edge.
-    const typename voronoi_edge<CT>::voronoi_cell_type *cell1 = edge.cell();
-    // Retrieve the cell to the right of the current edge.
-    const typename voronoi_edge<CT>::voronoi_cell_type *cell2 =
-        edge.twin()->cell();
-
-    discretization.push_back(get_point(edge.vertex0()->vertex()));
-    discretization.push_back(get_point(edge.vertex1()->vertex()));
-    if (edge.is_curved()) {
-      bool flag = cell1->contains_segment();
-      // point1 - site point;
-      point_type point1 = flag ?
-          get_point(cell2->point0()) : get_point(cell1->point0());
-      // point2 - start-point of the segment site;
-      point_type point2 = flag ?
-          get_point(cell1->point0()) : get_point(cell2->point0());
-      // point3 - endpoint of the segment site;
-      point_type point3 = flag ?
-          get_point(cell1->point1()) : get_point(cell2->point1());
-      fill_intermediate_points(
-          point1, point2, point3, max_error, discretization);
-    }
-  }
-
-  // Clip a linear Voronoi edge with a given rectangle.
-  template <typename CT>
-  static void clip(const voronoi_edge<CT> &edge,
-      const brect_type &brect, point_set_type &clipped_edge) {
-    // Don't process curved edges.
-    if (edge.is_curved())
-      return;
-
-    if (edge.is_finite()) {
-      clip(edge.vertex0()->vertex(), edge.vertex1()->vertex(),
-          brect, clipped_edge);
-    } else {
-      const typename voronoi_edge<CT>::voronoi_cell_type *cell1 = edge.cell();
-      const typename voronoi_edge<CT>::voronoi_cell_type *cell2 =
-          edge.twin()->cell();
-      point_type point1 = get_point(cell1->point0());
-      point_type point2 = get_point(cell2->point0());
-      if (point1 == point2) {
-        if (cell1->contains_segment())
-          point1 = get_point(cell1->point1());
-        else
-          point2 = get_point(cell2->point1());
-      }
-
-      if (edge.vertex0()) {
-        // Ray edge.
-        point_type origin = get_point(edge.vertex0()->vertex());
-        point_type direction(point1.y() - point2.y(), point2.x() - point1.x());
-        if (brect.contains(origin.x(), origin.y()))
-          clipped_edge.push_back(origin);
-        intersect(origin, direction, RAY, brect, clipped_edge);
-      } else if (edge.vertex1()) {
-        // Ray edge.
-        point_type origin = get_point(edge.vertex1()->vertex());
-        point_type direction(point2.y() - point1.y(), point1.x() - point2.x());
-        if (brect.contains(origin.x(), origin.y()))
-          clipped_edge.push_back(origin);
-        intersect(origin, direction, RAY, brect, clipped_edge);
-        // Keep correct ordering.
-        (std::reverse)(clipped_edge.begin(), clipped_edge.end());
-      } else if (cell1->contains_point() && cell2->contains_point()) {
-        // Primary line edge.
-        point_type origin((point1.x() + point2.x()) * 0.5, (point1.y() + point2.y()) * 0.5);
-        point_type direction(point1.y() - point2.y(), point2.x() - point1.x());
-        intersect(origin, direction, LINE, brect, clipped_edge);
-      } else {
-        point_type origin = cell1->contains_point() ? point1 : point2;
-        point_type direction(point1.y() - point2.y(), point2.x() - point1.x());
-        intersect(origin, direction, LINE, brect, clipped_edge);
-      }
-    }
-  }
-
-  // Clip an edge with the given rectangle.
-  template <typename PointType>
-  static void clip(const PointType &p1, const PointType &p2,
-      const brect_type &brect, point_set_type &clipped_edge) {
-    point_type start = get_point(p1);
-    point_type end = get_point(p2);
-    point_type direction(end.x() - start.x(), end.y() - start.y());
-    if (brect.contains(start.x(), start.y()))
-      clipped_edge.push_back(start);
-    if (p1 != p2)
-      intersect(start, direction, SEGMENT, brect, clipped_edge);
-    if (brect.contains(end.x(), end.y()))
-      clipped_edge.push_back(end);
-  }
-
-private:
-  typedef typename TRAITS::ctype_converter_type ctype_converter_type;
-
-  voronoi_utils();
-
-  // There are three different types of linear primitive:
-  //   1) SEGMENT - has two finite endpoints;
-  //   2) RAY - has one finite and one infinite endpoint;
-  //   3) LINE - has two infinite endpoints.
-  enum EdgeType {
-    SEGMENT = 0,
-    RAY = 1,
-    LINE = 2
-  };
-
-  template <typename P>
-  static point_type get_point(const P &point) {
-    coordinate_type x = to_fpt(point.x());
-    coordinate_type y = to_fpt(point.y());
-    return point_type(x, y);
-  }
-
-  static void intersect(
-      const point_type &origin, const point_type &direction,
-      EdgeType edge_type, const brect_type &brect,
-      point_set_type &clipped_edge) {
-    // Find the center of the rectangle.
-    coordinate_type center_x = (brect.x_min() + brect.x_max()) * to_fpt(0.5);
-    coordinate_type center_y = (brect.y_min() + brect.y_max()) * to_fpt(0.5);
-
-    // Find the half-diagonal vector of the rectangle.
-    coordinate_type len_x = brect.x_max() - center_x;
-    coordinate_type len_y = brect.y_max() - center_y;
-
-    // Find the vector between the origin and the center of the
-    // rectangle.
-    coordinate_type diff_x = origin.x() - center_x;
-    coordinate_type diff_y = origin.y() - center_y;
-
-    // Find the vector perpendicular to the direction vector.
-    coordinate_type perp_x = direction.y();
-    coordinate_type perp_y = -direction.x();
-
-    // Projection of the vector between the origin and the center of
-    // the rectangle on the line perpendicular to the direction vector.
-    coordinate_type lexpr = magnitude(perp_x * diff_x + perp_y * diff_y);
-
-    // Maximum projection among any of the half-diagonals of the
-    // rectangle on the line perpendicular to the direction vector.
-    coordinate_type rexpr =
-        magnitude(perp_x * len_x) + magnitude(perp_y * len_y);
-
-    // Intersection check. Compare projections.
-    if (lexpr > rexpr)
-      return;
-
-    // Intersection parameters. If fT[i]_used is true than:
-    // origin + fT[i] * direction = intersection point, i = 0 .. 1.
-    // For different edge types next fT values are acceptable:
-    // Segment: [0; 1];
-    // Ray: [0; infinity];
-    // Line: [-infinity; infinity].
-    bool fT0_used = false;
-    bool fT1_used = false;
-    coordinate_type fT0 = 0;
-    coordinate_type fT1 = 0;
-
-    // Check for the intersections with the lines
-    // going through the sides of the bounding rectangle.
-    clip_line(+direction.x(), -diff_x - len_x, fT0_used, fT1_used, fT0, fT1);
-    clip_line(-direction.x(), +diff_x - len_x, fT0_used, fT1_used, fT0, fT1);
-    clip_line(+direction.y(), -diff_y - len_y, fT0_used, fT1_used, fT0, fT1);
-    clip_line(-direction.y(), +diff_y - len_y, fT0_used, fT1_used, fT0, fT1);
-
-    // Update intersections vector.
-    if (fT0_used && check_extent(fT0, edge_type))
-      clipped_edge.push_back(point_type(
-          origin.x() + fT0 * direction.x(), origin.y() + fT0 * direction.y()));
-    if (fT1_used && fT0 != fT1 && check_extent(fT1, edge_type))
-      clipped_edge.push_back(point_type(
-          origin.x() + fT1 * direction.x(), origin.y() + fT1 * direction.y()));
-  }
-
-  // Find intermediate points of the parabola.
-  // Parabola is a locus of points equidistant from the point and segment
-  // sites. intermediate_points should contain two initial endpoints
-  // of the edge (Voronoi vertices). Intermediate points are inserted
-  // between the given two endpoints.
-  // Max_dist is the maximum distance allowed between parabola and line
-  // segments that discretize it.
-  static void fill_intermediate_points(
-      point_type point_site, point_type segment_site_start,
-      point_type segment_site_end, coordinate_type max_dist,
-      point_set_type &intermediate_points) {
-    // Apply the linear transformation to move start point of the
-    // segment to the point with coordinates (0, 0) and the direction
-    // of the segment to coincide the positive direction of the x-axis.
-    coordinate_type segm_vec_x =
-        segment_site_end.x() - segment_site_start.x();
-    coordinate_type segm_vec_y =
-        segment_site_end.y() - segment_site_start.y();
-    coordinate_type sqr_segment_length =
-        segm_vec_x * segm_vec_x + segm_vec_y * segm_vec_y;
-
-    // Compute x-coordinates of the endpoints of the edge
-    // in the transformed space.
-    coordinate_type projection_start = sqr_segment_length *
-        get_point_projection(
-            intermediate_points[0], segment_site_start, segment_site_end);
-    coordinate_type projection_end = sqr_segment_length *
-        get_point_projection(
-            intermediate_points[1], segment_site_start, segment_site_end);
-
-    // Compute parabola parameters in the transformed space.
-    // Parabola has next representation:
-    // f(x) = ((x-rot_x)^2 + rot_y^2) / (2.0*rot_y).
-    coordinate_type point_vec_x = point_site.x() - segment_site_start.x();
-    coordinate_type point_vec_y = point_site.y() - segment_site_start.y();
-    coordinate_type rot_x =
-        segm_vec_x * point_vec_x + segm_vec_y * point_vec_y;
-    coordinate_type rot_y =
-        segm_vec_x * point_vec_y - segm_vec_y * point_vec_x;
-
-    // Save the last point.
-    point_type last_point = intermediate_points[1];
-    intermediate_points.pop_back();
-
-    // Use stack to avoid recursion.
-    std::stack<coordinate_type> point_stack;
-    point_stack.push(projection_end);
-    coordinate_type cur_x = projection_start;
-    coordinate_type cur_y = parabola_y(cur_x, rot_x, rot_y);
-
-    // Adjust max_dist parameter in the transformed space.
-    max_dist *= max_dist * sqr_segment_length;
-
-    while (!point_stack.empty()) {
-      coordinate_type new_x = point_stack.top();
-      coordinate_type new_y = parabola_y(new_x, rot_x, rot_y);
-
-      // Compute coordinates of the point of the parabola that is
-      // furthest from the current line segment.
-      coordinate_type mid_x =
-          (new_y - cur_y) / (new_x - cur_x) * rot_y + rot_x;
-      coordinate_type mid_y = parabola_y(mid_x, rot_x, rot_y);
-
-      // Compute maximum distance between the given parabolic arc
-      // and line segment that discretize it.
-      coordinate_type dist = (new_y - cur_y) * (mid_x - cur_x) -
-          (new_x - cur_x) * (mid_y - cur_y);
-      dist = dist * dist / ((new_y - cur_y) * (new_y - cur_y) +
-          (new_x - cur_x) * (new_x - cur_x));
-      if (dist <= max_dist) {
-        // Distance between parabola and line segment is
-        // not greater than max_dist.
-        point_stack.pop();
-        coordinate_type inter_x = (segm_vec_x * new_x - segm_vec_y * new_y) /
-            sqr_segment_length + segment_site_start.x();
-        coordinate_type inter_y = (segm_vec_x * new_y + segm_vec_y * new_x) /
-            sqr_segment_length + segment_site_start.y();
-        intermediate_points.push_back(point_type(inter_x, inter_y));
-        cur_x = new_x;
-        cur_y = new_y;
-      } else {
-        point_stack.push(mid_x);
-      }
-    }
-
-    // Update last point.
-    intermediate_points.back() = last_point;
-  }
-
-  // Compute y(x) = ((x - a) * (x - a) + b * b) / (2 * b).
-  static coordinate_type parabola_y(
-      coordinate_type x, coordinate_type a, coordinate_type b) {
-    return ((x - a) * (x - a) + b * b) / (to_fpt(2.0) * b);
-  }
-
-  // Check whether extent is compatible with the edge type.
-  static bool check_extent(coordinate_type extent, EdgeType etype) {
-    switch (etype) {
-      case SEGMENT: return extent >= to_fpt(0.0) && extent <= to_fpt(1.0);
-      case RAY: return extent >= to_fpt(0.0);
-      case LINE: return true;
-    }
-    return true;
-  }
-
-  // Compute the absolute value.
-  static inline coordinate_type magnitude(coordinate_type value) {
-    return (value >= to_fpt(0.0)) ? value : -value;
-  }
-
-  // Find fT min and max values: fT = numer / denom.
-  static void clip_line(
-      coordinate_type denom, coordinate_type numer,
-      bool &fT0_used, bool &fT1_used,
-      coordinate_type &fT0, coordinate_type &fT1) {
-    if (denom > to_fpt(0.0)) {
-      if (fT1_used && numer > denom * fT1)
-        return;
-      if (!fT0_used || numer > denom * fT0) {
-        fT0_used = true;
-        fT0 = numer / denom;
-      }
-    } else if (denom < to_fpt(0.0)) {
-      if (fT0_used && numer > denom * fT0)
-        return;
-      if (!fT1_used || numer > denom * fT1) {
-        fT1_used = true;
-        fT1 = numer / denom;
-      }
-    }
-  }
-
-  // Get normalized length of the distance between:
-  //   1) point projection onto the segment;
-  //   2) start point of the segment.
-  // Return this length divided by the segment length.
-  // This is made to avoid sqrt computation during transformation from
-  // the initial space to the transformed one and vice versa.
-  // Assumption is made that projection of the point lies
-  // between the start-point and endpoint of the segment.
-  static coordinate_type get_point_projection(
-      const point_type &point,
-      const point_type &segment_start,
-      const point_type &segment_end) {
-    coordinate_type segment_vec_x = segment_end.x() - segment_start.x();
-    coordinate_type segment_vec_y = segment_end.y() - segment_start.y();
-    coordinate_type point_vec_x = point.x() - segment_start.x();
-    coordinate_type point_vec_y = point.y() - segment_start.y();
-    coordinate_type sqr_segment_length =
-        segment_vec_x * segment_vec_x + segment_vec_y * segment_vec_y;
-    coordinate_type vec_dot =
-        segment_vec_x * point_vec_x + segment_vec_y * point_vec_y;
-    return vec_dot / sqr_segment_length;
-  }
-
-  template <typename CT>
-  static coordinate_type to_fpt(const CT& value) {
-    static ctype_converter_type converter;
-    return converter(value);
-  }
-};
-}  // polygon
-}  // boost
-
-#endif  // BOOST_POLYGON_VORONOI_UTILS

Added: trunk/libs/polygon/example/voronoi_visual_utils.hpp
==============================================================================
--- (empty file)
+++ trunk/libs/polygon/example/voronoi_visual_utils.hpp	2012-08-27 14:49:52 EDT (Mon, 27 Aug 2012)
@@ -0,0 +1,187 @@
+// Boost.Polygon library voronoi_graphic_utils.hpp header file
+
+//          Copyright Andrii Sydorchuk 2010-2012.
+// Distributed under the Boost Software License, Version 1.0.
+//    (See accompanying file LICENSE_1_0.txt or copy at
+//          http://www.boost.org/LICENSE_1_0.txt)
+
+// See http://www.boost.org for updates, documentation, and revision history.
+
+#ifndef BOOST_POLYGON_VORONOI_VISUAL_UTILS
+#define BOOST_POLYGON_VORONOI_VISUAL_UTILS
+
+#include <stack>
+#include <vector>
+
+#include <boost/polygon/isotropy.hpp>
+#include <boost/polygon/point_concept.hpp>
+#include <boost/polygon/segment_concept.hpp>
+#include <boost/polygon/rectangle_concept.hpp>
+
+namespace boost{
+namespace polygon{
+
+// Utilities class, that contains set of routines handful for visualization.
+template <typename CT>
+class voronoi_visual_utils {
+public:
+  // Discretize parabolic Voronoi edge.
+  // Parabolic Voronoi edges are always formed by one point and one segment
+  // from the initial input set.
+  //
+  // Args:
+  //   point: input point.
+  //   segment: input segment.
+  //   max_dist: maximum discretization distance.
+  //   discretization: point discretization of the given Voronoi edge.
+  //
+  // Template arguments:
+  //   InCT: coordinate type of the input geometries (usually integer).
+  //   Point: point type, should model point concept.
+  //   Segment: segment type, should model segment concept.
+  // 
+  // Important:
+  //   discretization should contain both edge endpoints initially.
+  template <class InCT1, class InCT2,
+            template<class> class Point,
+            template<class> class Segment>
+  static
+  typename enable_if<
+    typename gtl_and<
+      typename gtl_if<
+        typename is_point_concept<
+          typename geometry_concept< Point<InCT1> >::type
+        >::type
+      >::type,
+      typename gtl_if<
+        typename is_segment_concept<
+          typename geometry_concept< Segment<InCT2> >::type
+        >::type
+      >::type
+    >::type,
+    void
+  >::type discretize(
+      const Point<InCT1>& point,
+      const Segment<InCT2>& segment,
+      const CT max_dist,
+      std::vector< Point<CT> >* discretization) {
+    // Apply the linear transformation to move start point of the segment to
+    // the point with coordinates (0, 0) and the direction of the segment to
+    // coincide the positive direction of the x-axis.
+    CT segm_vec_x = cast(x(high(segment))) - cast(x(low(segment)));
+    CT segm_vec_y = cast(y(high(segment))) - cast(y(low(segment)));
+    CT sqr_segment_length = segm_vec_x * segm_vec_x + segm_vec_y * segm_vec_y;
+
+    // Compute x-coordinates of the endpoints of the edge
+    // in the transformed space.
+    CT projection_start = sqr_segment_length *
+        get_point_projection((*discretization)[0], segment);
+    CT projection_end = sqr_segment_length *
+        get_point_projection((*discretization)[1], segment);
+
+    // Compute parabola parameters in the transformed space.
+    // Parabola has next representation:
+    // f(x) = ((x-rot_x)^2 + rot_y^2) / (2.0*rot_y).
+    CT point_vec_x = cast(x(point)) - cast(x(low(segment)));
+    CT point_vec_y = cast(y(point)) - cast(y(low(segment)));
+    CT rot_x = segm_vec_x * point_vec_x + segm_vec_y * point_vec_y;
+    CT rot_y = segm_vec_x * point_vec_y - segm_vec_y * point_vec_x;
+
+    // Save the last point.
+    Point<CT> last_point = (*discretization)[1];
+    discretization->pop_back();
+
+    // Use stack to avoid recursion.
+    std::stack<CT> point_stack;
+    point_stack.push(projection_end);
+    CT cur_x = projection_start;
+    CT cur_y = parabola_y(cur_x, rot_x, rot_y);
+
+    // Adjust max_dist parameter in the transformed space.
+    const CT max_dist_transformed = max_dist * max_dist * sqr_segment_length;
+    while (!point_stack.empty()) {
+      CT new_x = point_stack.top();
+      CT new_y = parabola_y(new_x, rot_x, rot_y);
+
+      // Compute coordinates of the point of the parabola that is
+      // furthest from the current line segment.
+      CT mid_x = (new_y - cur_y) / (new_x - cur_x) * rot_y + rot_x;
+      CT mid_y = parabola_y(mid_x, rot_x, rot_y);
+
+      // Compute maximum distance between the given parabolic arc
+      // and line segment that discretize it.
+      CT dist = (new_y - cur_y) * (mid_x - cur_x) -
+          (new_x - cur_x) * (mid_y - cur_y);
+      dist = dist * dist / ((new_y - cur_y) * (new_y - cur_y) +
+          (new_x - cur_x) * (new_x - cur_x));
+      if (dist <= max_dist_transformed) {
+        // Distance between parabola and line segment is less than max_dist.
+        point_stack.pop();
+        CT inter_x = (segm_vec_x * new_x - segm_vec_y * new_y) /
+            sqr_segment_length + cast(x(low(segment)));
+        CT inter_y = (segm_vec_x * new_y + segm_vec_y * new_x) /
+            sqr_segment_length + cast(y(low(segment)));
+        discretization->push_back(Point<CT>(inter_x, inter_y));
+        cur_x = new_x;
+        cur_y = new_y;
+      } else {
+        point_stack.push(mid_x);
+      }
+    }
+
+    // Update last point.
+    discretization->back() = last_point;
+  }
+
+private:
+  // Compute y(x) = ((x - a) * (x - a) + b * b) / (2 * b).
+  static CT parabola_y(CT x, CT a, CT b) {
+    return ((x - a) * (x - a) + b * b) / (b + b);
+  }
+
+  // Get normalized length of the distance between:
+  //   1) point projection onto the segment
+  //   2) start point of the segment
+  // Return this length divided by the segment length. This is made to avoid
+  // sqrt computation during transformation from the initial space to the
+  // transformed one and vice versa. The assumption is made that projection of
+  // the point lies between the start-point and endpoint of the segment.
+  template <class InCT,
+            template<class> class Point,
+            template<class> class Segment>
+  static
+  typename enable_if<
+    typename gtl_and<
+      typename gtl_if<
+        typename is_point_concept<
+          typename geometry_concept< Point<int> >::type
+        >::type
+      >::type,
+      typename gtl_if<
+        typename is_segment_concept<
+          typename geometry_concept< Segment<long> >::type
+        >::type
+      >::type
+    >::type,
+    CT
+  >::type get_point_projection(
+      const Point<CT>& point, const Segment<InCT>& segment) {
+    CT segment_vec_x = cast(x(high(segment))) - cast(x(low(segment)));
+    CT segment_vec_y = cast(y(high(segment))) - cast(y(low(segment)));
+    CT point_vec_x = x(point) - cast(x(low(segment)));
+    CT point_vec_y = y(point) - cast(y(low(segment)));
+    CT sqr_segment_length =
+        segment_vec_x * segment_vec_x + segment_vec_y * segment_vec_y;
+    CT vec_dot = segment_vec_x * point_vec_x + segment_vec_y * point_vec_y;
+    return vec_dot / sqr_segment_length;
+  }
+
+  template <typename InCT>
+  static CT cast(const InCT& value) {
+    return static_cast<CT>(value);
+  }
+};
+}
+}
+
+#endif  // BOOST_POLYGON_VORONOI_VISUAL_UTILS

Modified: trunk/libs/polygon/example/voronoi_visualizer.cpp
==============================================================================
--- trunk/libs/polygon/example/voronoi_visualizer.cpp	(original)
+++ trunk/libs/polygon/example/voronoi_visualizer.cpp	2012-08-27 14:49:52 EDT (Mon, 27 Aug 2012)
@@ -7,20 +7,24 @@
 
 // See http://www.boost.org for updates, documentation, and revision history.
 
-#include <vector>
+#include <cassert>
+#include <iostream>
 
 #include <QtOpenGL/QGLWidget>
 #include <QtGui/QtGui>
 
-#include <boost/polygon/voronoi_builder.hpp>
-#include <boost/polygon/voronoi_diagram.hpp>
-#include <boost/polygon/voronoi_utils.hpp>
+#include <boost/polygon/point_data.hpp>
+#include <boost/polygon/rectangle_data.hpp>
+#include <boost/polygon/segment_data.hpp>
+#include <boost/polygon/voronoi.hpp>
 using namespace boost::polygon;
 
+#include "voronoi_visual_utils.hpp"
+
 class GLWidget : public QGLWidget {
   Q_OBJECT
 public:
-  GLWidget(QMainWindow *parent = NULL) :
+  GLWidget(QMainWindow* parent = NULL) :
       QGLWidget(QGLFormat(QGL::SampleBuffers), parent),
       primary_edges_only_(false),
       internal_edges_only_(false) {
@@ -31,52 +35,36 @@
     return QSize(600, 600);
   }
 
-  void build(QString file_path) {
-    brect_.clear();
-    vb_.clear();
-    vd_.clear();
+  void build(const QString& file_path) {
+    // Clear all containers.
+    clear();
 
-    // Open file.
-    QFile data(file_path);
-    if (!data.open(QFile::ReadOnly)) {
-      QMessageBox::warning(this, tr("Voronoi Visualizer"),
-                           tr("Disable to open file ") + file_path);
-    }
+    // Read data.
+    read_data(file_path);
 
-    // Read points from the file.
-    QTextStream in_stream(&data);
-    int num_point_sites = 0;
-    int num_edge_sites = 0;
-    int x1, y1, x2, y2;
-    in_stream >> num_point_sites;
-    for (int i = 0; i < num_point_sites; ++i) {
-      in_stream >> x1 >> y1;
-      vb_.insert_point(x1, y1);
-      brect_.update(x1, y1);
-    }
-    in_stream >> num_edge_sites;
-    for (int i = 0; i < num_edge_sites; ++i) {
-      in_stream >> x1 >> y1 >> x2 >> y2;
-      vb_.insert_segment(x1, y1, x2, y2);
-      brect_.update(x1, y1);
-      brect_.update(x2, y2);
+    // No data, don't proceed.
+    if (!brect_initialized_) {
+      return;
     }
-    in_stream.flush();
 
-    // Build voronoi diagram.
-    vb_.construct(&vd_);
-    brect_ = voronoi_utils<coordinate_type>::scale(brect_, 1.4);
-    shift_ = point_type((brect_.x_min() + brect_.x_max()) * 0.5,
-                        (brect_.y_min() + brect_.y_max()) * 0.5);
+    // Construct bounding rectangle.
+    construct_brect();
 
+    // Construct voronoi diagram.
+    construct_voronoi(
+        point_data_.begin(), point_data_.end(),
+        segment_data_.begin(), segment_data_.end(),
+        &vd_);
+
+    // Color exterior edges.
     for (const_edge_iterator it = vd_.edges().begin();
-        it != vd_.edges().end(); ++it) {
+         it != vd_.edges().end(); ++it) {
       if (!it->is_finite()) {
-        remove_exterior(&(*it));
+        color_exterior(&(*it));
       }
     }
 
-    // Update view.
+    // Update view port.
     update_view_port();
   }
 
@@ -100,74 +88,10 @@
     qglClearColor(QColor::fromRgb(255, 255, 255));
     glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
 
-    // Draw voronoi sites.
-    {
-      glColor3f(0.0f, 0.5f, 1.0f);
-      glPointSize(9);
-      glBegin(GL_POINTS);
-      for (const_cell_iterator it = vd_.cells().begin();
-          it != vd_.cells().end(); ++it) {
-        if (!it->contains_segment()) {
-          glVertex2f(it->point0().x() - shift_.x(),
-                     it->point0().y() - shift_.y());
-        }
-      }
-      glEnd();
-      glPointSize(6);
-      glLineWidth(2.7f);
-      glBegin(GL_LINES);
-      for (const_cell_iterator it = vd_.cells().begin();
-          it != vd_.cells().end(); ++it) {
-        if (it->contains_segment()) {
-          glVertex2f(it->point0().x() - shift_.x(),
-                     it->point0().y() - shift_.y());
-          glVertex2f(it->point1().x() - shift_.x(),
-                     it->point1().y() - shift_.y());
-        }
-      }
-      glEnd();
-      glLineWidth(1.0);
-    }
-
-    // Draw voronoi vertices.
-    /*{
-      glColor3f(0.0f, 0.0f, 0.0f);
-      glBegin(GL_POINTS);
-      for (const_vertex_iterator it = vd_.vertices().begin();
-          it != vd_.vertices().end(); ++it) {
-        glVertex2f(it->vertex().x() - shift_.x(),
-                   it->vertex().y() - shift_.y());
-      }
-      glEnd();
-    }*/
-
-    // Draw voronoi edges.
-    {
-      glColor3f(0.0f, 0.0f, 0.0f);
-      glLineWidth(1.7f);
-      glBegin(GL_LINES);
-      for (const_edge_iterator it = vd_.edges().begin();
-          it != vd_.edges().end(); ++it) {
-        if (primary_edges_only_ && !it->is_primary()) {
-          continue;
-        }
-        if (internal_edges_only_ && it->color()) {
-          continue;
-        }
-        std::vector<point_type> vec;
-        if (!it->is_finite())
-          voronoi_utils<coordinate_type>::clip(*it, brect_, vec);
-        else {
-          coordinate_type max_error = 1E-3 * (brect_.x_max() - brect_.x_min());
-          voronoi_utils<coordinate_type>::discretize(*it, max_error, vec);
-        }
-        for (int i = 0; i < static_cast<int>(vec.size()) - 1; ++i) {
-          glVertex2f(vec[i].x() - shift_.x(), vec[i].y() - shift_.y());
-          glVertex2f(vec[i+1].x() - shift_.x(), vec[i+1].y() - shift_.y());
-        }
-      }
-      glEnd();
-    }
+    draw_points();
+    draw_segments();
+    draw_vertices();
+    draw_edges();
   }
 
   void resizeGL(int width, int height) {
@@ -175,14 +99,20 @@
     glViewport((width - side) / 2, (height - side) / 2, side, side);
   }
 
-  void timerEvent(QTimerEvent *) {
+  void timerEvent(QTimerEvent*) {
     update();
   }
 
 private:
   typedef double coordinate_type;
-  typedef detail::point_2d<double> point_type;
-  typedef voronoi_diagram<double> VD;
+  typedef point_data<double> point_type;
+  typedef segment_data<double> segment_type;
+  typedef rectangle_data<double> rect_type;
+  typedef voronoi_builder<int> VB;
+  typedef voronoi_diagram<coordinate_type> VD;
+  typedef VD::cell_type cell_type;
+  typedef VD::cell_type::source_index_type source_index_type;
+  typedef VD::cell_type::source_category_type source_category_type;
   typedef VD::edge_type edge_type;
   typedef VD::cell_container_type cell_container_type;
   typedef VD::cell_container_type vertex_container_type;
@@ -191,19 +121,74 @@
   typedef VD::const_vertex_iterator const_vertex_iterator;
   typedef VD::const_edge_iterator const_edge_iterator;
 
-  static const std::size_t VISITED = 1;
+  static const std::size_t EXTERNAL_COLOR = 1;
+
+  void clear() {
+    brect_initialized_ = false;
+    point_data_.clear();
+    segment_data_.clear();
+    vd_.clear();
+  }
+
+  void read_data(const QString& file_path) {
+    QFile data(file_path);
+    if (!data.open(QFile::ReadOnly)) {
+      QMessageBox::warning(
+          this, tr("Voronoi Visualizer"),
+          tr("Disable to open file ") + file_path);
+    }
+    QTextStream in_stream(&data);
+    std::size_t num_points, num_segments;
+    int x1, y1, x2, y2;
+    in_stream >> num_points;
+    for (std::size_t i = 0; i < num_points; ++i) {
+      in_stream >> x1 >> y1;
+      point_type p(x1, y1);
+      update_brect(p);
+      point_data_.push_back(p);
+    }
+    in_stream >> num_segments;
+    for (std::size_t i = 0; i < num_segments; ++i) {
+      in_stream >> x1 >> y1 >> x2 >> y2;
+      point_type lp(x1, y1);
+      point_type hp(x2, y2);
+      update_brect(lp);
+      update_brect(hp);
+      segment_data_.push_back(segment_type(lp, hp));
+    }
+    in_stream.flush();
+  }
+
+  void update_brect(const point_type& point) {
+    if (brect_initialized_) {
+      encompass(brect_, point);
+    } else {
+      set_points(brect_, point, point);
+      brect_initialized_ = true;
+    }
+  }
+
+  void construct_brect() {
+    double side = (std::max)(xh(brect_) - xl(brect_), yh(brect_) - yl(brect_));
+    center(shift_, brect_);
+    set_points(brect_, shift_, shift_);
+    bloat(brect_, side * 1.2);
+  }
 
-  void remove_exterior(const VD::edge_type* edge) {
-    if (edge->color())
+  void color_exterior(const VD::edge_type* edge) {
+    if (edge->color() == EXTERNAL_COLOR) {
       return;
-    edge->color(VISITED);
-    edge->twin()->color(VISITED);
-    const voronoi_diagram<double>::vertex_type* v = edge->vertex1();
-    if (v == NULL || !edge->is_primary())
+    }
+    edge->color(EXTERNAL_COLOR);
+    edge->twin()->color(EXTERNAL_COLOR);
+    const VD::vertex_type* v = edge->vertex1();
+    if (v == NULL || !edge->is_primary()) {
       return;
-    const voronoi_diagram<double>::edge_type* e = v->incident_edge();
+    }
+    v->color(EXTERNAL_COLOR);
+    const VD::edge_type* e = v->incident_edge();
     do {
-      remove_exterior(e);
+      color_exterior(e);
       e = e->rot_next();
     } while (e != v->incident_edge());
   }
@@ -211,16 +196,182 @@
   void update_view_port() {
     glMatrixMode(GL_PROJECTION);
     glLoadIdentity();
-    glOrtho(brect_.x_min() - shift_.x(), brect_.x_max() - shift_.x(),
-            brect_.y_min() - shift_.y(), brect_.y_max() - shift_.y(),
-            -1.0, 1.0);
+    rect_type view_rect = brect_;
+    deconvolve(view_rect, shift_);
+    glOrtho(xl(view_rect), xh(view_rect), yl(view_rect), yh(view_rect), -1.0, 1.0);
     glMatrixMode(GL_MODELVIEW);
   }
 
-  bounding_rectangle<double> brect_;
+  void draw_points() {
+    // Draw input points and endpoints of the input segments.
+    glColor3f(0.0f, 0.5f, 1.0f);
+    glPointSize(9);
+    glBegin(GL_POINTS);
+    for (std::size_t i = 0; i < point_data_.size(); ++i) {
+      point_type point = point_data_[i];
+      deconvolve(point, shift_);
+      glVertex2f(point.x(), point.y());
+    }
+    for (std::size_t i = 0; i < segment_data_.size(); ++i) {
+      point_type lp = deconvolve(low(segment_data_[i]), shift_);
+      point_type hp = deconvolve(high(segment_data_[i]), shift_);
+      glVertex2f(lp.x(), lp.y());
+      glVertex2f(hp.x(), hp.y());
+    }
+    glEnd();
+  }
+
+  void draw_segments() {
+    // Draw input segments.
+    glColor3f(0.0f, 0.5f, 1.0f);
+    glLineWidth(2.7f);
+    glBegin(GL_LINES);
+    for (std::size_t i = 0; i < segment_data_.size(); ++i) {
+      point_type lp = deconvolve(low(segment_data_[i]), shift_);
+      point_type hp = deconvolve(high(segment_data_[i]), shift_);
+      glVertex2f(lp.x(), lp.y());
+      glVertex2f(hp.x(), hp.y());
+    }
+    glEnd();
+  }
+
+  void draw_vertices() {
+    // Draw voronoi vertices.
+    glColor3f(0.0f, 0.0f, 0.0f);
+    glPointSize(6);
+    glBegin(GL_POINTS);
+    for (const_vertex_iterator it = vd_.vertices().begin();
+         it != vd_.vertices().end(); ++it) {
+      if (internal_edges_only_ && (it->color() == EXTERNAL_COLOR)) {
+        continue;
+      }
+      point_type vertex(it->x(), it->y());
+      vertex = deconvolve(vertex, shift_);
+      glVertex2f(vertex.x(), vertex.y());
+    }
+    glEnd();
+  }
+  void draw_edges() {
+    // Draw voronoi edges.
+    glColor3f(0.0f, 0.0f, 0.0f);
+    glLineWidth(1.7f);
+    for (const_edge_iterator it = vd_.edges().begin();
+         it != vd_.edges().end(); ++it) {
+      if (primary_edges_only_ && !it->is_primary()) {
+        continue;
+      }
+      if (internal_edges_only_ && (it->color() == EXTERNAL_COLOR)) {
+        continue;
+      }
+      std::vector<point_type> samples;
+      if (!it->is_finite()) {
+        clip_infinite_edge(*it, &samples);
+      } else {
+        point_type vertex0(it->vertex0()->x(), it->vertex0()->y());
+        samples.push_back(vertex0);
+        point_type vertex1(it->vertex1()->x(), it->vertex1()->y());
+        samples.push_back(vertex1);
+        if (it->is_curved()) {
+          sample_curved_edge(*it, &samples);
+        }
+      }
+      glBegin(GL_LINE_STRIP);
+      for (std::size_t i = 0; i < samples.size(); ++i) {
+        point_type vertex = deconvolve(samples[i], shift_);
+        glVertex2f(vertex.x(), vertex.y());
+      }
+      glEnd();
+    }
+  }
+
+  void clip_infinite_edge(const edge_type& edge, std::vector<point_type>* clipped_edge) {
+    const cell_type& cell1 = *edge.cell();
+    const cell_type& cell2 = *edge.twin()->cell();
+    point_type origin, direction;
+    // Infinite edges could not be created by two segment sites.
+    assert(cell1.contains_point() || cell2.contains_point());
+    if (cell1.contains_point() && cell2.contains_point()) {
+      point_type p1 = retrieve_point(cell1);
+      point_type p2 = retrieve_point(cell2);
+      origin.x((p1.x() + p2.x()) * 0.5);
+      origin.y((p1.y() + p2.y()) * 0.5);
+      direction.x(p1.y() - p2.y());
+      direction.y(p2.x() - p1.x());
+    } else {
+      origin = cell1.contains_segment() ?
+          retrieve_point(cell2) :
+          retrieve_point(cell1);
+      segment_type segment = cell1.contains_segment() ?
+          retrieve_segment(cell1) :
+          retrieve_segment(cell2);
+      coordinate_type dx = high(segment).x() - low(segment).x();
+      coordinate_type dy = high(segment).y() - low(segment).y();
+      if ((low(segment) == origin) ^ cell1.contains_point()) {
+        direction.x(dy);
+        direction.y(-dx);
+      } else {
+        direction.x(-dy);
+        direction.y(dx);
+      }
+    }
+    coordinate_type side = xh(brect_) - xl(brect_);
+    coordinate_type koef = side / (std::max)(fabs(direction.x()), fabs(direction.y()));
+    if (edge.vertex0() == NULL) {
+      clipped_edge->push_back(point_type(
+          origin.x() - direction.x() * koef,
+          origin.y() - direction.y() * koef));
+    } else {
+      clipped_edge->push_back(point_type(edge.vertex0()->x(), edge.vertex0()->y()));
+    }
+    if (edge.vertex1() == NULL) {
+      clipped_edge->push_back(point_type(
+          origin.x() + direction.x() * koef,
+          origin.y() + direction.y() * koef));
+    } else {
+      clipped_edge->push_back(point_type(edge.vertex1()->x(), edge.vertex1()->y()));
+    }
+  }
+
+  void sample_curved_edge(
+      const edge_type& edge,
+      std::vector<point_type>* sampled_edge) {
+    coordinate_type max_dist = 1E-3 * (xh(brect_) - xl(brect_));
+    point_type point = edge.cell()->contains_point() ?
+        retrieve_point(*edge.cell()) :
+        retrieve_point(*edge.twin()->cell());
+    segment_type segment = edge.cell()->contains_point() ?
+        retrieve_segment(*edge.twin()->cell()) :
+        retrieve_segment(*edge.cell());
+    voronoi_visual_utils<coordinate_type>::discretize(
+        point, segment, max_dist, sampled_edge);
+  }
+
+  point_type retrieve_point(const cell_type& cell) {
+    source_index_type index = cell.source_index();
+    source_category_type category = cell.source_category();
+    if (category == detail::SOURCE_CATEGORY_SINGLE_POINT) {
+      return point_data_[index];
+    }
+    index -= point_data_.size();
+    if (category == detail::SOURCE_CATEGORY_SEGMENT_START_POINT) {
+      return low(segment_data_[index]);
+    } else {
+      return high(segment_data_[index]);
+    }
+  }
+
+  segment_type retrieve_segment(const cell_type& cell) {
+    source_index_type index = cell.source_index() - point_data_.size();
+    return segment_data_[index];
+  }
+
   point_type shift_;
-  default_voronoi_builder vb_;
-  voronoi_diagram<coordinate_type> vd_;
+  std::vector<point_type> point_data_;
+  std::vector<segment_type> segment_data_;
+  rect_type brect_;
+  VB vb_;
+  VD vd_;
+  bool brect_initialized_;
   bool primary_edges_only_;
   bool internal_edges_only_;
 };
@@ -233,7 +384,7 @@
     file_dir_ = QDir(QDir::currentPath(), tr("*.txt"));
     file_name_ = tr("");
 
-    QHBoxLayout *centralLayout = new QHBoxLayout;
+    QHBoxLayout* centralLayout = new QHBoxLayout;
     centralLayout->addWidget(glWidget_);
     centralLayout->addLayout(create_file_layout());
     setLayout(centralLayout);
@@ -281,8 +432,8 @@
   }
 
 private:
-  QGridLayout *create_file_layout() {
-    QGridLayout *file_layout = new QGridLayout;
+  QGridLayout* create_file_layout() {
+    QGridLayout* file_layout = new QGridLayout;
 
     message_label_ = new QLabel("Double click item to build voronoi diagram:");
 
@@ -292,20 +443,20 @@
                         this,
                         SLOT(build()));
 
-    QCheckBox *primary_checkbox = new QCheckBox("Show primary edges only.");
+    QCheckBox* primary_checkbox = new QCheckBox("Show primary edges only.");
     connect(primary_checkbox, SIGNAL(clicked()),
         this, SLOT(primary_edges_only()));
 
-    QCheckBox *internal_checkbox = new QCheckBox("Show internal edges only.");
+    QCheckBox* internal_checkbox = new QCheckBox("Show internal edges only.");
     connect(internal_checkbox, SIGNAL(clicked()),
         this, SLOT(internal_edges_only()));
 
-    QPushButton *browse_button =
+    QPushButton* browse_button =
         new QPushButton(tr("Browse Input Directory"));
     connect(browse_button, SIGNAL(clicked()), this, SLOT(browse()));
     browse_button->setMinimumHeight(50);
 
-    QPushButton *print_scr_button = new QPushButton(tr("Make Screenshot"));
+    QPushButton* print_scr_button = new QPushButton(tr("Make Screenshot"));
     connect(print_scr_button, SIGNAL(clicked()), this, SLOT(print_scr()));
     print_scr_button->setMinimumHeight(50);
 
@@ -334,12 +485,12 @@
 
   QDir file_dir_;
   QString file_name_;
-  GLWidget *glWidget_;
-  QListWidget *file_list_;
-  QLabel *message_label_;
+  GLWidget* glWidget_;
+  QListWidget* file_list_;
+  QLabel* message_label_;
 };
 
-int main(int argc, char *argv[]) {
+int main(int argc, char* argv[]) {
   QApplication app(argc, argv);
   MainWindow window;
   window.show();

Modified: trunk/libs/polygon/test/voronoi_builder_test.cpp
==============================================================================
--- trunk/libs/polygon/test/voronoi_builder_test.cpp	(original)
+++ trunk/libs/polygon/test/voronoi_builder_test.cpp	2012-08-27 14:49:52 EDT (Mon, 27 Aug 2012)
@@ -26,10 +26,6 @@
 typedef vd_type::const_cell_iterator const_cell_iterator;
 typedef vd_type::const_vertex_iterator const_vertex_iterator;
 
-#define CHECK_EQUAL_POINTS(p1, p2) \
-    BOOST_CHECK(p1.x() == static_cast<T>(p2.x())); \
-    BOOST_CHECK(p1.y() == static_cast<T>(p2.y()))
-
 #define CHECK_OUTPUT_SIZE(output, cells, vertices, edges) \
     BOOST_CHECK(output.num_cells() == cells); \
     BOOST_CHECK(output.num_vertices() == vertices); \
@@ -37,8 +33,6 @@
 
 #define VERIFY_OUTPUT(output) \
     BOOST_CHECK(voronoi_test_helper::verify_output(output, \
-        voronoi_test_helper::HALF_EDGE_ORIENTATION)); \
-    BOOST_CHECK(voronoi_test_helper::verify_output(output, \
         voronoi_test_helper::CELL_CONVEXITY)); \
     BOOST_CHECK(voronoi_test_helper::verify_output(output, \
         voronoi_test_helper::INCIDENT_EDGES_CCW_ORDER)); \
@@ -77,21 +71,21 @@
   const_cell_iterator cell_it = test_output.cells().begin();
   cell_it++;
 
-  const voronoi_edge_type *edge1_1 = cell_it->incident_edge();
-  const voronoi_edge_type *edge1_2 = edge1_1->twin();
+  const voronoi_edge_type* edge1_1 = cell_it->incident_edge();
+  const voronoi_edge_type* edge1_2 = edge1_1->twin();
 
-  BOOST_CHECK_EQUAL(edge1_1->twin() == edge1_2, true);
-  BOOST_CHECK_EQUAL(edge1_2->twin() == edge1_1, true);
+  BOOST_CHECK(edge1_1->twin() == edge1_2);
+  BOOST_CHECK(edge1_2->twin() == edge1_1);
 
-  BOOST_CHECK_EQUAL(edge1_1->next() == edge1_1, true);
-  BOOST_CHECK_EQUAL(edge1_1->prev() == edge1_1, true);
-  BOOST_CHECK_EQUAL(edge1_1->rot_next() == NULL, true);
-  BOOST_CHECK_EQUAL(edge1_1->rot_prev() == NULL, true);
-
-  BOOST_CHECK_EQUAL(edge1_2->next() == edge1_2, true);
-  BOOST_CHECK_EQUAL(edge1_2->prev() == edge1_2, true);
-  BOOST_CHECK_EQUAL(edge1_2->rot_next() == NULL, true);
-  BOOST_CHECK_EQUAL(edge1_2->rot_prev() == NULL, true);
+  BOOST_CHECK(edge1_1->next() == edge1_1);
+  BOOST_CHECK(edge1_1->prev() == edge1_1);
+  BOOST_CHECK(edge1_1->rot_next() == NULL);
+  BOOST_CHECK(edge1_1->rot_prev() == NULL);
+
+  BOOST_CHECK(edge1_2->next() == edge1_2);
+  BOOST_CHECK(edge1_2->prev() == edge1_2);
+  BOOST_CHECK(edge1_2->rot_next() == NULL);
+  BOOST_CHECK(edge1_2->rot_prev() == NULL);
 }
 
 // Sites: (0, 0), (1, 1), (2, 2).
@@ -106,25 +100,25 @@
   CHECK_OUTPUT_SIZE(test_output, 3, 0, 2);
 
   const_cell_iterator cell_it = test_output.cells().begin();
-  const voronoi_edge_type *edge1_1 = cell_it->incident_edge();
-  const voronoi_edge_type *edge1_2 = edge1_1->twin();
+  const voronoi_edge_type* edge1_1 = cell_it->incident_edge();
+  const voronoi_edge_type* edge1_2 = edge1_1->twin();
   cell_it++;
   cell_it++;
-  const voronoi_edge_type *edge2_2 = cell_it->incident_edge();
-  const voronoi_edge_type *edge2_1 = edge2_2->twin();
+  const voronoi_edge_type* edge2_2 = cell_it->incident_edge();
+  const voronoi_edge_type* edge2_1 = edge2_2->twin();
 
-  BOOST_CHECK_EQUAL(edge1_1->twin() == edge1_2 && edge1_2->twin() == edge1_1, true);
-  BOOST_CHECK_EQUAL(edge2_1->twin() == edge2_2 && edge2_2->twin() == edge2_1, true);
+  BOOST_CHECK(edge1_1->twin() == edge1_2 && edge1_2->twin() == edge1_1);
+  BOOST_CHECK(edge2_1->twin() == edge2_2 && edge2_2->twin() == edge2_1);
 
-  BOOST_CHECK_EQUAL(edge1_1->next() == edge1_1 && edge1_1->prev() == edge1_1, true);
-  BOOST_CHECK_EQUAL(edge1_1->rot_next() == NULL && edge1_1->rot_prev() == NULL, true);
-  BOOST_CHECK_EQUAL(edge1_2->rot_next() == NULL && edge1_2->rot_prev() == NULL, true);
-  BOOST_CHECK_EQUAL(edge2_1->rot_next() == NULL && edge2_1->rot_prev() == NULL, true);
-  BOOST_CHECK_EQUAL(edge2_2->next() == edge2_2 && edge2_2->prev() == edge2_2, true);
-  BOOST_CHECK_EQUAL(edge2_2->rot_next() == NULL && edge2_2->rot_prev() == NULL, true);
+  BOOST_CHECK(edge1_1->next() == edge1_1 && edge1_1->prev() == edge1_1);
+  BOOST_CHECK(edge1_1->rot_next() == NULL && edge1_1->rot_prev() == NULL);
+  BOOST_CHECK(edge1_2->rot_next() == NULL && edge1_2->rot_prev() == NULL);
+  BOOST_CHECK(edge2_1->rot_next() == NULL && edge2_1->rot_prev() == NULL);
+  BOOST_CHECK(edge2_2->next() == edge2_2 && edge2_2->prev() == edge2_2);
+  BOOST_CHECK(edge2_2->rot_next() == NULL && edge2_2->rot_prev() == NULL);
 
-  BOOST_CHECK_EQUAL(edge1_2->next() == edge2_1 && edge1_2->prev() == edge2_1, true);
-  BOOST_CHECK_EQUAL(edge2_1->next() == edge1_2 && edge2_1->prev() == edge1_2, true);
+  BOOST_CHECK(edge1_2->next() == edge2_1 && edge1_2->prev() == edge2_1);
+  BOOST_CHECK(edge2_1->next() == edge1_2 && edge2_1->prev() == edge1_2);
 }
 
 // Sites: (0, 0), (0, 4), (2, 1).
@@ -142,39 +136,39 @@
   CHECK_OUTPUT_SIZE(test_output, 3, 1, 3);
 
   const_vertex_iterator it = test_output.vertices().begin();
-  BOOST_CHECK_EQUAL(it->vertex().x() == static_cast<coordinate_type>(0.25) &&
-                    it->vertex().y() == static_cast<coordinate_type>(2.0), true);
+  BOOST_CHECK_EQUAL(it->x(), 0.25);
+  BOOST_CHECK_EQUAL(it->y(), 2.0);
 
-  const voronoi_edge_type *edge1_1 = it->incident_edge();
-  const voronoi_edge_type *edge1_2 = edge1_1->twin();
-  CHECK_EQUAL_POINTS(edge1_1->cell()->point0(), point2);
-  CHECK_EQUAL_POINTS(edge1_2->cell()->point0(), point3);
-
-  const voronoi_edge_type *edge2_1 = edge1_1->rot_prev();
-  const voronoi_edge_type *edge2_2 = edge2_1->twin();
-  CHECK_EQUAL_POINTS(edge2_1->cell()->point0(), point3);
-  CHECK_EQUAL_POINTS(edge2_2->cell()->point0(), point1);
-
-  const voronoi_edge_type *edge3_1 = edge2_1->rot_prev();
-  const voronoi_edge_type *edge3_2 = edge3_1->twin();
-  CHECK_EQUAL_POINTS(edge3_1->cell()->point0(), point1);
-  CHECK_EQUAL_POINTS(edge3_2->cell()->point0(), point2);
-
-  BOOST_CHECK_EQUAL(edge1_2->twin() == edge1_1, true);
-  BOOST_CHECK_EQUAL(edge2_2->twin() == edge2_1, true);
-  BOOST_CHECK_EQUAL(edge3_2->twin() == edge3_1, true);
-
-  BOOST_CHECK_EQUAL(edge1_1->prev() == edge3_2 && edge1_1->next() == edge3_2, true);
-  BOOST_CHECK_EQUAL(edge2_1->prev() == edge1_2 && edge2_1->next() == edge1_2, true);
-  BOOST_CHECK_EQUAL(edge3_1->prev() == edge2_2 && edge3_1->next() == edge2_2, true);
-
-  BOOST_CHECK_EQUAL(edge1_2->next() == edge2_1 && edge1_2->prev() == edge2_1, true);
-  BOOST_CHECK_EQUAL(edge2_2->next() == edge3_1 && edge2_2->prev() == edge3_1, true);
-  BOOST_CHECK_EQUAL(edge3_2->next() == edge1_1 && edge3_2->prev() == edge1_1, true);
-
-  BOOST_CHECK_EQUAL(edge1_1->rot_next() == edge3_1, true);
-  BOOST_CHECK_EQUAL(edge3_1->rot_next() == edge2_1, true);
-  BOOST_CHECK_EQUAL(edge2_1->rot_next() == edge1_1, true);
+  const voronoi_edge_type* edge1_1 = it->incident_edge();
+  const voronoi_edge_type* edge1_2 = edge1_1->twin();
+  BOOST_CHECK(edge1_1->cell()->source_index() == 1);
+  BOOST_CHECK(edge1_2->cell()->source_index() == 2);
+
+  const voronoi_edge_type* edge2_1 = edge1_1->rot_prev();
+  const voronoi_edge_type* edge2_2 = edge2_1->twin();
+  BOOST_CHECK(edge2_1->cell()->source_index() == 2);
+  BOOST_CHECK(edge2_2->cell()->source_index() == 0);
+
+  const voronoi_edge_type* edge3_1 = edge2_1->rot_prev();
+  const voronoi_edge_type* edge3_2 = edge3_1->twin();
+  BOOST_CHECK(edge3_1->cell()->source_index() == 0);
+  BOOST_CHECK(edge3_2->cell()->source_index() == 1);
+
+  BOOST_CHECK(edge1_2->twin() == edge1_1);
+  BOOST_CHECK(edge2_2->twin() == edge2_1);
+  BOOST_CHECK(edge3_2->twin() == edge3_1);
+
+  BOOST_CHECK(edge1_1->prev() == edge3_2 && edge1_1->next() == edge3_2);
+  BOOST_CHECK(edge2_1->prev() == edge1_2 && edge2_1->next() == edge1_2);
+  BOOST_CHECK(edge3_1->prev() == edge2_2 && edge3_1->next() == edge2_2);
+
+  BOOST_CHECK(edge1_2->next() == edge2_1 && edge1_2->prev() == edge2_1);
+  BOOST_CHECK(edge2_2->next() == edge3_1 && edge2_2->prev() == edge3_1);
+  BOOST_CHECK(edge3_2->next() == edge1_1 && edge3_2->prev() == edge1_1);
+
+  BOOST_CHECK(edge1_1->rot_next() == edge3_1);
+  BOOST_CHECK(edge3_1->rot_next() == edge2_1);
+  BOOST_CHECK(edge2_1->rot_next() == edge1_1);
 }
 
 // Sites: (0, 1), (2, 0), (2, 4).
@@ -192,39 +186,39 @@
   CHECK_OUTPUT_SIZE(test_output, 3, 1, 3);
 
   const_vertex_iterator it = test_output.vertices().begin();
-  BOOST_CHECK_EQUAL(it->vertex().x() == static_cast<coordinate_type>(1.75) &&
-                    it->vertex().y() == static_cast<coordinate_type>(2.0), true);
+  BOOST_CHECK_EQUAL(it->x(), 1.75);
+  BOOST_CHECK_EQUAL(it->y(), 2.0);
 
-  const voronoi_edge_type *edge1_1 = it->incident_edge();
-  const voronoi_edge_type *edge1_2 = edge1_1->twin();
-  CHECK_EQUAL_POINTS(edge1_1->cell()->point0(), point3);
-  CHECK_EQUAL_POINTS(edge1_2->cell()->point0(), point2);
-
-  const voronoi_edge_type *edge2_1 = edge1_1->rot_prev();
-  const voronoi_edge_type *edge2_2 = edge2_1->twin();
-  CHECK_EQUAL_POINTS(edge2_1->cell()->point0(), point2);
-  CHECK_EQUAL_POINTS(edge2_2->cell()->point0(), point1);
-
-  const voronoi_edge_type *edge3_1 = edge2_1->rot_prev();
-  const voronoi_edge_type *edge3_2 = edge3_1->twin();
-  CHECK_EQUAL_POINTS(edge3_1->cell()->point0(), point1);
-  CHECK_EQUAL_POINTS(edge3_2->cell()->point0(), point3);
-
-  BOOST_CHECK_EQUAL(edge1_2->twin() == edge1_1, true);
-  BOOST_CHECK_EQUAL(edge2_2->twin() == edge2_1, true);
-  BOOST_CHECK_EQUAL(edge3_2->twin() == edge3_1, true);
-
-  BOOST_CHECK_EQUAL(edge1_1->prev() == edge3_2 && edge1_1->next() == edge3_2, true);
-  BOOST_CHECK_EQUAL(edge2_1->prev() == edge1_2 && edge2_1->next() == edge1_2, true);
-  BOOST_CHECK_EQUAL(edge3_1->prev() == edge2_2 && edge3_1->next() == edge2_2, true);
-
-  BOOST_CHECK_EQUAL(edge1_2->next() == edge2_1 && edge1_2->prev() == edge2_1, true);
-  BOOST_CHECK_EQUAL(edge2_2->next() == edge3_1 && edge2_2->prev() == edge3_1, true);
-  BOOST_CHECK_EQUAL(edge3_2->next() == edge1_1 && edge3_2->prev() == edge1_1, true);
-
-  BOOST_CHECK_EQUAL(edge1_1->rot_next() == edge3_1, true);
-  BOOST_CHECK_EQUAL(edge3_1->rot_next() == edge2_1, true);
-  BOOST_CHECK_EQUAL(edge2_1->rot_next() == edge1_1, true);
+  const voronoi_edge_type* edge1_1 = it->incident_edge();
+  const voronoi_edge_type* edge1_2 = edge1_1->twin();
+  BOOST_CHECK(edge1_1->cell()->source_index() == 2);
+  BOOST_CHECK(edge1_2->cell()->source_index() == 1);
+
+  const voronoi_edge_type* edge2_1 = edge1_1->rot_prev();
+  const voronoi_edge_type* edge2_2 = edge2_1->twin();
+  BOOST_CHECK(edge2_1->cell()->source_index() == 1);
+  BOOST_CHECK(edge2_2->cell()->source_index() == 0);
+
+  const voronoi_edge_type* edge3_1 = edge2_1->rot_prev();
+  const voronoi_edge_type* edge3_2 = edge3_1->twin();
+  BOOST_CHECK(edge3_1->cell()->source_index() == 0);
+  BOOST_CHECK(edge3_2->cell()->source_index() == 2);
+
+  BOOST_CHECK(edge1_2->twin() == edge1_1);
+  BOOST_CHECK(edge2_2->twin() == edge2_1);
+  BOOST_CHECK(edge3_2->twin() == edge3_1);
+
+  BOOST_CHECK(edge1_1->prev() == edge3_2 && edge1_1->next() == edge3_2);
+  BOOST_CHECK(edge2_1->prev() == edge1_2 && edge2_1->next() == edge1_2);
+  BOOST_CHECK(edge3_1->prev() == edge2_2 && edge3_1->next() == edge2_2);
+
+  BOOST_CHECK(edge1_2->next() == edge2_1 && edge1_2->prev() == edge2_1);
+  BOOST_CHECK(edge2_2->next() == edge3_1 && edge2_2->prev() == edge3_1);
+  BOOST_CHECK(edge3_2->next() == edge1_1 && edge3_2->prev() == edge1_1);
+
+  BOOST_CHECK(edge1_1->rot_next() == edge3_1);
+  BOOST_CHECK(edge3_1->rot_next() == edge2_1);
+  BOOST_CHECK(edge2_1->rot_next() == edge1_1);
 }
 
 // Sites: (0, 0), (0, 1), (1, 0), (1, 1).
@@ -245,49 +239,49 @@
 
   // Check voronoi vertex.
   const_vertex_iterator it = test_output.vertices().begin();
-  BOOST_CHECK_EQUAL(it->vertex().x() == static_cast<coordinate_type>(0.5) &&
-                    it->vertex().y() == static_cast<coordinate_type>(0.5), true);
+  BOOST_CHECK_EQUAL(it->x(), 0.5);
+  BOOST_CHECK_EQUAL(it->y(), 0.5);
 
   // Check voronoi edges.
-  const voronoi_edge_type *edge1_1 = it->incident_edge();
-  const voronoi_edge_type *edge1_2 = edge1_1->twin();
-  CHECK_EQUAL_POINTS(edge1_1->cell()->point0(), points[3]);
-  CHECK_EQUAL_POINTS(edge1_2->cell()->point0(), points[2]);
-
-  const voronoi_edge_type *edge2_1 = edge1_1->rot_prev();
-  const voronoi_edge_type *edge2_2 = edge2_1->twin();
-  CHECK_EQUAL_POINTS(edge2_1->cell()->point0(), points[2]);
-  CHECK_EQUAL_POINTS(edge2_2->cell()->point0(), points[0]);
-
-  const voronoi_edge_type *edge3_1 = edge2_1->rot_prev();
-  const voronoi_edge_type *edge3_2 = edge3_1->twin();
-  CHECK_EQUAL_POINTS(edge3_1->cell()->point0(), points[0]);
-  CHECK_EQUAL_POINTS(edge3_2->cell()->point0(), points[1]);
-
-  const voronoi_edge_type *edge4_1 = edge3_1->rot_prev();
-  const voronoi_edge_type *edge4_2 = edge4_1->twin();
-  CHECK_EQUAL_POINTS(edge4_1->cell()->point0(), points[1]);
-  CHECK_EQUAL_POINTS(edge4_2->cell()->point0(), points[3]);
-
-  BOOST_CHECK_EQUAL(edge1_2->twin() == edge1_1, true);
-  BOOST_CHECK_EQUAL(edge2_2->twin() == edge2_1, true);
-  BOOST_CHECK_EQUAL(edge3_2->twin() == edge3_1, true);
-  BOOST_CHECK_EQUAL(edge4_2->twin() == edge4_1, true);
-
-  BOOST_CHECK_EQUAL(edge1_1->prev() == edge4_2 && edge1_1->next() == edge4_2, true);
-  BOOST_CHECK_EQUAL(edge2_1->prev() == edge1_2 && edge2_1->next() == edge1_2, true);
-  BOOST_CHECK_EQUAL(edge3_1->prev() == edge2_2 && edge3_1->next() == edge2_2, true);
-  BOOST_CHECK_EQUAL(edge4_1->prev() == edge3_2 && edge4_1->next() == edge3_2, true);
-
-  BOOST_CHECK_EQUAL(edge1_2->next() == edge2_1 && edge1_2->prev() == edge2_1, true);
-  BOOST_CHECK_EQUAL(edge2_2->next() == edge3_1 && edge2_2->prev() == edge3_1, true);
-  BOOST_CHECK_EQUAL(edge3_2->next() == edge4_1 && edge3_2->prev() == edge4_1, true);
-  BOOST_CHECK_EQUAL(edge4_2->next() == edge1_1 && edge4_2->prev() == edge1_1, true);
-
-  BOOST_CHECK_EQUAL(edge1_1->rot_next() == edge4_1, true);
-  BOOST_CHECK_EQUAL(edge4_1->rot_next() == edge3_1, true);
-  BOOST_CHECK_EQUAL(edge3_1->rot_next() == edge2_1, true);
-  BOOST_CHECK_EQUAL(edge2_1->rot_next() == edge1_1, true);
+  const voronoi_edge_type* edge1_1 = it->incident_edge();
+  const voronoi_edge_type* edge1_2 = edge1_1->twin();
+  BOOST_CHECK(edge1_1->cell()->source_index() == 3);
+  BOOST_CHECK(edge1_2->cell()->source_index() == 2);
+
+  const voronoi_edge_type* edge2_1 = edge1_1->rot_prev();
+  const voronoi_edge_type* edge2_2 = edge2_1->twin();
+  BOOST_CHECK(edge2_1->cell()->source_index() == 2);
+  BOOST_CHECK(edge2_2->cell()->source_index() == 0);
+
+  const voronoi_edge_type* edge3_1 = edge2_1->rot_prev();
+  const voronoi_edge_type* edge3_2 = edge3_1->twin();
+  BOOST_CHECK(edge3_1->cell()->source_index() == 0);
+  BOOST_CHECK(edge3_2->cell()->source_index() == 1);
+
+  const voronoi_edge_type* edge4_1 = edge3_1->rot_prev();
+  const voronoi_edge_type* edge4_2 = edge4_1->twin();
+  BOOST_CHECK(edge4_1->cell()->source_index() == 1);
+  BOOST_CHECK(edge4_2->cell()->source_index() == 3);
+
+  BOOST_CHECK(edge1_2->twin() == edge1_1);
+  BOOST_CHECK(edge2_2->twin() == edge2_1);
+  BOOST_CHECK(edge3_2->twin() == edge3_1);
+  BOOST_CHECK(edge4_2->twin() == edge4_1);
+
+  BOOST_CHECK(edge1_1->prev() == edge4_2 && edge1_1->next() == edge4_2);
+  BOOST_CHECK(edge2_1->prev() == edge1_2 && edge2_1->next() == edge1_2);
+  BOOST_CHECK(edge3_1->prev() == edge2_2 && edge3_1->next() == edge2_2);
+  BOOST_CHECK(edge4_1->prev() == edge3_2 && edge4_1->next() == edge3_2);
+
+  BOOST_CHECK(edge1_2->next() == edge2_1 && edge1_2->prev() == edge2_1);
+  BOOST_CHECK(edge2_2->next() == edge3_1 && edge2_2->prev() == edge3_1);
+  BOOST_CHECK(edge3_2->next() == edge4_1 && edge3_2->prev() == edge4_1);
+  BOOST_CHECK(edge4_2->next() == edge1_1 && edge4_2->prev() == edge1_1);
+
+  BOOST_CHECK(edge1_1->rot_next() == edge4_1);
+  BOOST_CHECK(edge4_1->rot_next() == edge3_1);
+  BOOST_CHECK(edge3_1->rot_next() == edge2_1);
+  BOOST_CHECK(edge2_1->rot_next() == edge1_1);
 }
 
 #ifdef NDEBUG
@@ -303,11 +297,12 @@
     point_vec_small.clear();
     point_vec_large.clear();
     int koef = std::numeric_limits<int>::max() / max_value[k];
-    for (int i = 0; i < grid_size[k]; i++)
+    for (int i = 0; i < grid_size[k]; i++) {
       for (int j = 0; j < grid_size[k]; j++) {
         point_vec_small.push_back(point_data<T>(i, j));
         point_vec_large.push_back(point_data<T>(koef * i, koef * j));
       }
+    }
     construct_voronoi(point_vec_small.begin(), point_vec_small.end(), &test_output_small);
     construct_voronoi(point_vec_large.begin(), point_vec_large.end(), &test_output_large);
     VERIFY_OUTPUT(test_output_small);
@@ -478,6 +473,19 @@
   VERIFY_NO_HALF_EDGE_INTERSECTIONS(test_output);
 }
 
+BOOST_AUTO_TEST_CASE_TEMPLATE(segment_site_test9, T, test_types) {
+  vd_type test_output;
+  std::vector< segment_data<T> > segments;
+  point_data<T> point1(0, 0);
+  point_data<T> point2(2, 0);
+  point_data<T> point3(4, 0);
+  segments.push_back(segment_data<T>(point1, point2));
+  segments.push_back(segment_data<T>(point2, point3));
+  construct_voronoi(segments.begin(), segments.end(), &test_output);
+  CHECK_OUTPUT_SIZE(test_output, 5, 0, 4);
+  VERIFY_NO_HALF_EDGE_INTERSECTIONS(test_output);
+}
+
 #ifdef NDEBUG
 BOOST_AUTO_TEST_CASE_TEMPLATE(segment_grid_test, T, test_types) {
   vd_type test_output_small, test_output_large;

Modified: trunk/libs/polygon/test/voronoi_structures_test.cpp
==============================================================================
--- trunk/libs/polygon/test/voronoi_structures_test.cpp	(original)
+++ trunk/libs/polygon/test/voronoi_structures_test.cpp	2012-08-27 14:49:52 EDT (Mon, 27 Aug 2012)
@@ -30,35 +30,56 @@
   BOOST_CHECK_EQUAL(p.y(), 4);
 }
 
-BOOST_AUTO_TEST_CASE(site_event_test2) {
+BOOST_AUTO_TEST_CASE(source_category_test1) {
+  BOOST_CHECK(belongs(SOURCE_CATEGORY_SINGLE_POINT, GEOMETRY_CATEGORY_POINT));
+  BOOST_CHECK(belongs(SOURCE_CATEGORY_SEGMENT_START_POINT, GEOMETRY_CATEGORY_POINT));
+  BOOST_CHECK(belongs(SOURCE_CATEGORY_SEGMENT_END_POINT, GEOMETRY_CATEGORY_POINT));
+  BOOST_CHECK(!belongs(SOURCE_CATEGORY_INITIAL_SEGMENT, GEOMETRY_CATEGORY_POINT));
+  BOOST_CHECK(!belongs(SOURCE_CATEGORY_REVERSE_SEGMENT, GEOMETRY_CATEGORY_POINT));
+
+  BOOST_CHECK(!belongs(SOURCE_CATEGORY_SINGLE_POINT, GEOMETRY_CATEGORY_SEGMENT));
+  BOOST_CHECK(!belongs(SOURCE_CATEGORY_SEGMENT_START_POINT, GEOMETRY_CATEGORY_SEGMENT));
+  BOOST_CHECK(!belongs(SOURCE_CATEGORY_SEGMENT_END_POINT, GEOMETRY_CATEGORY_SEGMENT));
+  BOOST_CHECK(belongs(SOURCE_CATEGORY_INITIAL_SEGMENT, GEOMETRY_CATEGORY_SEGMENT));
+  BOOST_CHECK(belongs(SOURCE_CATEGORY_REVERSE_SEGMENT, GEOMETRY_CATEGORY_SEGMENT));
+}
+
+BOOST_AUTO_TEST_CASE(site_event_test1) {
   site_type s(1, 2);
-  BOOST_CHECK(s.x0() == s.x1() && s.x1() == 1);
-  BOOST_CHECK(s.y0() == s.y1() && s.y1() == 2);
-  BOOST_CHECK(s.sorted_index() == 0);
-  BOOST_CHECK(!s.is_segment());
-  BOOST_CHECK(!s.is_inverse());
   s.sorted_index(1);
   s.initial_index(2);
+  s.source_category(SOURCE_CATEGORY_SEGMENT_START_POINT);
+  BOOST_CHECK(s.x0() == 1 && s.x1() == 1);
+  BOOST_CHECK(s.y0() == 2 && s.y1() == 2);
+  BOOST_CHECK(s.is_point());
+  BOOST_CHECK(!s.is_segment());
+  BOOST_CHECK(!s.is_inverse());
   BOOST_CHECK(s.sorted_index() == 1);
   BOOST_CHECK(s.initial_index() == 2);
-  BOOST_CHECK(!s.is_inverse());
+  BOOST_CHECK(s.source_category() == SOURCE_CATEGORY_SEGMENT_START_POINT);
 }
 
-BOOST_AUTO_TEST_CASE(site_event_test3) {
+BOOST_AUTO_TEST_CASE(site_event_test2) {
   site_type s(1, 2, 3, 4);
+  s.sorted_index(1);
+  s.initial_index(2);
+  s.source_category(SOURCE_CATEGORY_INITIAL_SEGMENT);
   BOOST_CHECK(s.x0(true) == 1 && s.x0() == 1);
   BOOST_CHECK(s.y0(true) == 2 && s.y0() == 2);
   BOOST_CHECK(s.x1(true) == 3 && s.x1() == 3);
   BOOST_CHECK(s.y1(true) == 4 && s.y1() == 4);
-  BOOST_CHECK(s.sorted_index() == 0);
+  BOOST_CHECK(!s.is_point());
   BOOST_CHECK(s.is_segment());
   BOOST_CHECK(!s.is_inverse());
+  BOOST_CHECK(s.source_category() == SOURCE_CATEGORY_INITIAL_SEGMENT);
+  
   s.inverse();
   BOOST_CHECK(s.x1(true) == 1 && s.x0() == 1);
   BOOST_CHECK(s.y1(true) == 2 && s.y0() == 2);
   BOOST_CHECK(s.x0(true) == 3 && s.x1() == 3);
   BOOST_CHECK(s.y0(true) == 4 && s.y1() == 4);
   BOOST_CHECK(s.is_inverse());
+  BOOST_CHECK(s.source_category() == SOURCE_CATEGORY_INITIAL_SEGMENT);
 }
 
 BOOST_AUTO_TEST_CASE(circle_event_test) {
@@ -107,13 +128,13 @@
 
 BOOST_AUTO_TEST_CASE(beach_line_node_data_test) {
   node_data_type node_data(NULL);
-  BOOST_CHECK_EQUAL(node_data.edge() == NULL, true);
-  BOOST_CHECK_EQUAL(node_data.circle_event() == NULL, true);
+  BOOST_CHECK(node_data.edge() == NULL);
+  BOOST_CHECK(node_data.circle_event() == NULL);
   int data = 4;
   node_data.circle_event(&data);
-  BOOST_CHECK_EQUAL(node_data.edge() == NULL, true);
-  BOOST_CHECK_EQUAL(node_data.circle_event() == &data, true);
+  BOOST_CHECK(node_data.edge() == NULL);
+  BOOST_CHECK(node_data.circle_event() == &data);
   node_data.edge(&data);
-  BOOST_CHECK_EQUAL(node_data.edge() == &data, true);
-  BOOST_CHECK_EQUAL(node_data.circle_event() == &data, true);
+  BOOST_CHECK(node_data.edge() == &data);
+  BOOST_CHECK(node_data.circle_event() == &data);
 }

Modified: trunk/libs/polygon/test/voronoi_test_helper.hpp
==============================================================================
--- trunk/libs/polygon/test/voronoi_test_helper.hpp	(original)
+++ trunk/libs/polygon/test/voronoi_test_helper.hpp	2012-08-27 14:49:52 EDT (Mon, 27 Aug 2012)
@@ -27,56 +27,39 @@
   LEFT = 1
 };
 
-template <typename Point2D>
-kOrientation get_orientation(const Point2D &point1,
-                             const Point2D &point2,
-                             const Point2D &point3) {
-  typename Point2D::coordinate_type a = (point2.x() - point1.x()) * (point3.y() - point2.y());
-  typename Point2D::coordinate_type b = (point2.y() - point1.y()) * (point3.x() - point2.x());
-  if (a == b)
+template <typename VERTEX>
+kOrientation get_orientation(const VERTEX& v1, const VERTEX& v2, const VERTEX& v3) {
+  typename VERTEX::coordinate_type lhs = (v2.x() - v1.x()) * (v3.y() - v2.y());
+  typename VERTEX::coordinate_type rhs = (v2.y() - v1.y()) * (v3.x() - v2.x());
+  if (lhs == rhs) {
     return COLLINEAR;
-  return (a < b) ? RIGHT : LEFT;
-}
-
-template <typename Output>
-bool verify_half_edge_orientation(const Output &output) {
-  typedef typename Output::point_type point_type;
-  typename Output::const_edge_iterator edge_it;
-  for (edge_it = output.edges().begin(); 
-       edge_it != output.edges().end(); edge_it++) {
-    if (edge_it->is_finite()) {
-      const point_type &site_point = edge_it->cell()->point0();
-      const point_type &start_point = edge_it->vertex0()->vertex();
-      const point_type &end_point = edge_it->vertex1()->vertex();
-      if (get_orientation(start_point, end_point, site_point) != LEFT)
-        return false;
-    }
   }
-  return true;
+  return (lhs < rhs) ? RIGHT : LEFT;
 }
 
-template <typename Output>
-bool verify_cell_convexity(const Output &output) {
-  typedef typename Output::point_type point_type;
-  typename Output::const_cell_iterator cell_it;
+template <typename OUTPUT>
+bool verify_cell_convexity(const OUTPUT& output) {
+  typename OUTPUT::const_cell_iterator cell_it;
   for (cell_it = output.cells().begin();
        cell_it != output.cells().end(); cell_it++) {
-    const typename Output::edge_type *edge = cell_it->incident_edge();
+    const typename OUTPUT::edge_type* edge = cell_it->incident_edge();
     if (edge)
       do {
-        if (edge->next()->prev() != edge)
+        if (edge->next()->prev() != edge) {
           return false;
-        if (edge->cell() != &(*cell_it))
+        }
+        if (edge->cell() != &(*cell_it)) {
+          return false;
+        }
+        if (edge->vertex1() != edge->next()->vertex0()) {
           return false;
+        }
         if (edge->vertex0() != NULL &&
             edge->vertex1() != NULL &&
-            edge->vertex1() == edge->next()->vertex0() &&
             edge->next()->vertex1() != NULL) {
-          const point_type &vertex1 = edge->vertex0()->vertex();
-          const point_type &vertex2 = edge->vertex1()->vertex();
-          const point_type &vertex3 = edge->next()->vertex1()->vertex();
-          if (get_orientation(vertex1, vertex2, vertex3) != LEFT)
+          if (get_orientation(*edge->vertex0(), *edge->vertex1(), *edge->next()->vertex1()) != LEFT) {
             return false;
+          }
         }
         edge = edge->next();
       } while(edge != cell_it->incident_edge());
@@ -84,38 +67,36 @@
   return true;
 }
 
-template <typename Output>
-bool verify_incident_edges_ccw_order(const Output &output) {
-  typedef typename Output::point_type point_type;
-  typedef typename Output::edge_type voronoi_edge_type;
-  typename Output::const_vertex_iterator vertex_it;
+template <typename OUTPUT>
+bool verify_incident_edges_ccw_order(const OUTPUT& output) {
+  typedef typename OUTPUT::edge_type voronoi_edge_type;
+  typename OUTPUT::const_vertex_iterator vertex_it;
   for (vertex_it = output.vertices().begin();
-     vertex_it != output.vertices().end(); vertex_it++) {
+       vertex_it != output.vertices().end(); vertex_it++) {
     if (vertex_it->is_degenerate())
       continue;
-    const voronoi_edge_type *edge = vertex_it->incident_edge();
+    const voronoi_edge_type* edge = vertex_it->incident_edge();
     do {
-      const voronoi_edge_type *edge_next1 = edge->rot_next();
-      const voronoi_edge_type *edge_next2 = edge_next1->rot_next();
-      const point_type &site1 = edge->cell()->point0();
-      const point_type &site2 = edge_next1->cell()->point0();
-      const point_type &site3 = edge_next2->cell()->point0();
-
-      if (get_orientation(site1, site2, site3) != LEFT)
+      const voronoi_edge_type* next_edge = edge->rot_next();
+      if (edge->vertex0() != next_edge->vertex0()) {
         return false;
-
+      }
+      if (edge->vertex1() != NULL && next_edge->vertex1() != NULL &&
+          get_orientation(*edge->vertex1(), *edge->vertex0(), *next_edge->vertex1()) == LEFT) {
+        return false;
+      }
       edge = edge->rot_next();
     } while (edge != vertex_it->incident_edge());
   }
   return true;
 }
 
-template <typename P>
+template <typename VERTEX>
 struct cmp {
-  bool operator()(const P& p1, const P& p2) const {
-    if (p1.x() != p2.x()) 
-      return p1.x() < p2.x();
-    return p1.y() < p2.y();
+  bool operator()(const VERTEX& v1, const VERTEX& v2) const {
+    if (v1.x() != v2.x()) 
+      return v1.x() < v2.x();
+    return v1.y() < v2.y();
   }
 };
 
@@ -124,17 +105,16 @@
   // Create map from edges with first point less than the second one.
   // Key is the first point of the edge, value is a vector of second points
   // with the same first point.
-  typedef typename Output::point_type point_type;
-  cmp<point_type> comparator;
-  std::map< point_type, std::vector<point_type>, cmp<point_type> > edge_map;
+  typedef typename Output::vertex_type vertex_type;
+  cmp<vertex_type> comparator;
+  std::map< vertex_type, std::vector<vertex_type>, cmp<vertex_type> > edge_map;
   typename Output::const_edge_iterator edge_it;
   for (edge_it = output.edges().begin();
        edge_it != output.edges().end(); edge_it++) {
     if (edge_it->is_finite()) {
-        const point_type &vertex0 = edge_it->vertex0()->vertex();
-        const point_type &vertex1 = edge_it->vertex1()->vertex();
-      if (comparator(vertex0, vertex1))
-        edge_map[vertex0].push_back(vertex1);
+      if (comparator(*edge_it->vertex0(), *edge_it->vertex1())) {
+        edge_map[*edge_it->vertex0()].push_back(*edge_it->vertex1());
+      }
     }
   }
   return !intersection_check(edge_map);
@@ -192,19 +172,16 @@
 }
 
 enum kVerification {
-  HALF_EDGE_ORIENTATION = 1,
-  CELL_CONVEXITY = 2,
-  INCIDENT_EDGES_CCW_ORDER = 4,
-  NO_HALF_EDGE_INTERSECTIONS = 8,
-  FAST_VERIFICATION = 7,
-  COMPLETE_VERIFICATION = 15
+  CELL_CONVEXITY = 1,
+  INCIDENT_EDGES_CCW_ORDER = 2,
+  NO_HALF_EDGE_INTERSECTIONS = 4,
+  FAST_VERIFICATION = 3,
+  COMPLETE_VERIFICATION = 7
 };
 
 template <typename Output>
 bool verify_output(const Output &output, kVerification mask) {
   bool result = true;
-  if (mask & HALF_EDGE_ORIENTATION)
-    result &= verify_half_edge_orientation(output);
   if (mask & CELL_CONVEXITY)
     result &= verify_cell_convexity(output);
   if (mask & INCIDENT_EDGES_CCW_ORDER)
@@ -215,7 +192,7 @@
 }
 
 template <typename PointIterator>
-void save_points(PointIterator first, PointIterator last, const char *file_name) {
+void save_points(PointIterator first, PointIterator last, const char* file_name) {
   std::ofstream ofs(file_name);
   ofs << std::distance(first, last) << std::endl;
   for (PointIterator it = first; it != last; ++it) {
@@ -225,7 +202,7 @@
 }
 
 template <typename SegmentIterator>
-void save_segments(SegmentIterator first, SegmentIterator last, const char *file_name) {
+void save_segments(SegmentIterator first, SegmentIterator last, const char* file_name) {
   std::ofstream ofs(file_name);
   ofs << std::distance(first, last) << std::endl;
   for (SegmentIterator it = first; it != last; ++it) {
@@ -236,7 +213,7 @@
 }
 
 template <typename T>
-void clean_segment_set(std::vector< segment_data<T> > &data) {
+void clean_segment_set(std::vector< segment_data<T> >& data) {
   typedef T Unit;
   typedef typename scanline_base<Unit>::Point Point;
   typedef typename scanline_base<Unit>::half_edge half_edge;